Question: Counting and extracting in a new file the unique mapped reads from BWA 0.7.15-r1140 mapping
0
gravatar for valopes
2.3 years ago by
valopes30
valopes30 wrote:

Hello guys.

I am trying to have a bam file containing just the unique mapped reads from BWA.

I've used to check the mapping results:

samtools flagstat input.bam

and I got

151399675 + 0 in total (QC-passed reads + QC-failed reads)
1299223 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
150524060 + 0 mapped (99.42% : N/A)
150100452 + 0 paired in sequencing
75050226 + 0 read1
75050226 + 0 read2
141573428 + 0 properly paired (94.32% : N/A)
148964374 + 0 with itself and mate mapped
260463 + 0 singletons (0.17% : N/A)
6491734 + 0 with mate mapped to a different chr
5052530 + 0 with mate mapped to a different chr (mapQ>=5)

I know 150524060 means the unique and multiple mapped reads, but I want to make sure how many unique mapped I have.

So I used:

samtools view -F 0x40 -c input.bam

and I got 75664312

I am not sure if this number above is right.

Could someone please help me with this? Which command should I use to count unique mapped reads? and which to extract in a new bam file just the unique mapped (excluding the multi and unmapped ones)?

Thanks in advance.

assembly • 1.5k views
ADD COMMENTlink modified 2.3 years ago by Kevin Blighe53k • written 2.3 years ago by valopes30

I added markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

101010 Button

ADD REPLYlink written 2.3 years ago by WouterDeCoster42k

Thanks. I didn't know that.

ADD REPLYlink written 2.3 years ago by valopes30

Hi, all I know for now is that 0x40 is not the correct flag. Take a look at the SAM specification here: https://samtools.github.io/hts-specs/SAMv1.pdf (go to page 5). I don't believe that BWA records anything specific about unique mappings, but it does retain information about primary (best) and secondary alignments. However, I know that it's possible to infer uniquely mapped reads if you align with Bowtie (bowtie 1), using --best and -m 1 command-line parameters. Bowtie2 doesn't have these parameters but it still retains unique mapping information. When you run Bowtie/Bowtie2, the alignment report gives this information. I've just done some mappings right now (see the part in bold):

5157135 (100.00%) were paired; of these:

805804 (15.63%) aligned concordantly 0 times

2102530 (40.77%) aligned concordantly exactly 1 time

2248801 (43.61%) aligned concordantly >1 times

...

86.78% overall alignment rate

ADD REPLYlink written 2.3 years ago by Kevin Blighe53k

Thank you. As I was used to work with Bowtie2, I was trying to get the same info with BWA. But now I see it is different.

ADD REPLYlink written 2.3 years ago by valopes30
1
gravatar for ATpoint
2.3 years ago by
ATpoint28k
Germany
ATpoint28k wrote:

Here is a tool that can help you handle the bitwise flags. 0x40 stands for 'first in pair', and samtools view -F stands for discarding those reads flagged as such. Your settings do not really make sense here. If unique means that multimappers are to be excluded and only high-quality alignments are to be counted, I would go with:

samtools view -bh -q 30 -f 3 -F 2316 -c in.bam

It keeps (-f) mapped and properly paired reads that are not multimappers (MAPQ=0), discarding (-F) non-primary alignment, unmapped reads and mates. Play around with the flags to fit your needs.

ADD COMMENTlink modified 2.3 years ago • written 2.3 years ago by ATpoint28k

Thank you. I will try it and let you know how it work for me.

ADD REPLYlink written 2.3 years ago by valopes30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1015 users visited in the last hour