Extracting unique reads from a .ma or .bam file?

[JohnK]

Hi all,

I wanted to know the most efficient way to extract uniquely mapped reads from BioScope's resequencing mapping pipeline output in either the .ma or .bam format (those are the files I have to work with)? By unique, I mean reads that mapped to one location alone. Thanks!

J

[drio]

.ma way:

They keep all the coordinates for all the mapping positions available. You can use that
to drop reads that map to multiple locations.
Actually, they should have another file (.ma.unique.*) that does that for you already.

.bam way:

Look to see if they have any tag (NH for example) to indicate how many alignments where done for that read.

You can post some records here we can take a look.

[JohnK]

.ma way:

They keep all the coordinates for all the mapping positions available. You can use that
to drop reads that map to multiple locations.
Actually, they should have another file (.ma.unique.*) that does that for you already.

.bam way:

Look to see if they have any tag (NH for example) to indicate how many alignments where done for that read.

You can post some records here we can take a look.

I was looking for the .unique.ma, but I didn't find it or any options with BioScope after mapping my resequencing sample. I understand that Corona produced that automatically. From what I understand, unique reads in bioscope are a matter of score > 5.



>2_29_53_F3,9_-110926286.224.2.15):q7
T33333203331202321313113130023213300030321303003300
>2_29_134_F3
T23331313312033202331333300333110330110312303003300
>2_29_155_F3
T11322113002233330002013020033212300020300303003302
>2_29_160_F3
T11331323300033211333223330333013320210301303003300
>2_29_200_F3,17_-90327533.241.5.0):q21,8_-43471518.229.3.0):q0
T01030013331133333333303333033303333330333303303300
>2_29_217_F3
T13330023311233323333300302133123220020302303002300



******

Is there a way to pick out a uniquely mapped read from a .bam file with SAM/Bam tools?

[drio]

Is there a way to pick out a uniquely mapped read from a .bam file with SAM/Bam tools?

No within the mandatory fields. That's why I was asking if you can post here
some SAM entries from Bioscope so we can take a look.

[lh3]

http://sourceforge.net/apps/mediawik...from_SAM.2FBAM.

[drio]

http://sourceforge.net/apps/mediawik...from_SAM.2FBAM.

He wants to get the unique mapped reads.
But yes I agree with this (from the FAQ):

We prefer to say an alignment is `reliable' rather than `unique' as `uniqueness' is not well defined in general cases. You can get reliable alignments by setting a threshold on mapping quality

[JohnK]

He wants to get the unique mapped reads.
But yes I agree with this (from the FAQ):

I had some trouble converting between bam back to sam. It appears to be incorrect. Here's my usage:

samtools view in.bam > name.sam

Probably not right usage here. Also, I'm re-formatting the matobam.ini file so it calcs 'unique' (or more reliable) reads based on the QV > 5. If anybody knows if this threshold is optimal, your help would be greatly appreciated? I'll then use bed tools to convert to bed and get the coverage off the unique. Hopefully, this will work and I can get a rough estimate of concordance between the unique reads from Corona versus the unique reads from BioScope.

[drio]

Your samtools cmd is correct. What error/problem are you getting?

[JohnK]

Your samtools cmd is correct. What error/problem are you getting?

17_1177_738 16 chr1 2999999 7 25H25M * 0 0 NNGAATTCTTTTCTATGATTTAGTT """IIIIIIIIIIIIIIIIIIIII! RG:Z:20100521202751742 NH:i:1 CS:Z:T30123003213322000220302023010232222031211331023323 CQ:Z:BB?@;A@??B@=4>B<6<<A<?=?A=96=?>289>;7/=?04<0B3+(=> CM:i:2
774_1259_2007 16 chr1 2999999 7 25H25M * 0 0 NNGAATTCTTTTCTATGATTTAGTT """IIIIIIIIIIIIIIIIIIIII! RG:Z:20100521202751742 NH:i:1 CS:Z:T30123003213322000220302023010232222031211331023123 CQ:Z:B:AABA7BBBA?@?B=;A>B;<B:A?7=.@?2<;A:=<<?;<?5>:15?@ CM:i:2

[drio]

So what is the problem? What were you expecting as a output?

[JohnK]

So what is the problem? What were you expecting as a output?

Nah- I realized it's correct now. It's just my first time really messing with and having to look at the SAM format meticulously. So this is what I'm expecting as output, but I just want to grab the unique reads. I'll keep reading over the manuals and what not and see if I can figure this out.

[JohnK]

So what is the problem? What were you expecting as a output?

Just a follow up, if anyone cares to know. You can extract unique (aka more 'reliable' reads) via the matobam.output.parameter='unique' setting. Then with BedTools convert from bam to bed and so forth... I appreciate the help, all!

[bair]

My bam file includes 2x35bp mate-pair and 2x101 pair-end, any efficient way to split it into MP.bam and PE.bam files?

[drio]

My bam file includes 2x35bp mate-pair and 2x101 pair-end, any efficient way to split it into MP.bam and PE.bam files?

If you are confortable with C or java, use samtools or picard to iterate over your bam and dump the data to PE or MP bam based on the read length.
That's the most efficient way.

There are other ports of samtools for lisp and python (maybe more) if you
want a more higher level language.

Another option would be to stream the bam data to sam with samtools and pipe it to your split script. After that you can convert back to BAM. It would be slower but it'd work just fine.

Something like:

#!/bin/bash
#
samtools view ./i.bam  | ruby -a -ne 'puts $_ if $F[9].size == 35' > 35.sam
samtools view ./i.bam  | ruby -a -ne 'puts $_ if $F[9].size == 120' > 120.sam 
samtools view -bS ./35.sam > 35.bam
samtools view -bS ./120.sam > 120.bam

[bair]

Thank you Drio,

Your suggestions and script code are very helpful, like to use the simple code.

Thanks again.