I have some bam files generated from BWA-MEM using these command:
So i should get only mapped reads in my output. I'm using bwa v0.7.5a.
Using samtools flagstat, i retrieved the number of mapped reads from the sorted bam file, and in some cases, this was higher than the number of original input paired-end reads.
Now this could be due to reads aligning to multiple places, but I don't have any "XA" tags, nor "XT" tags for that matter.
I already saw this post about getting uniquely mapping reads from BWA.
http://seqanswers.com/forums/showthread.php?t=30120
I want to get stats on unique and multi-mapping reads.
Is using MapQ=0 a reliable way of flagging multi-mapping reads?
And why am I not getting XA and XT tags at all?
Here are a couple of lines of my output - I am not seeing any XT/XA tags, nor X0 or X1 tags.
thanks
Code:
bwa mem -t 8 ref1 R1.fa.gz R2.fa.gz | samtools view -S -h -F 4 -b - | samtools sort - bwaaln-sorted
Using samtools flagstat, i retrieved the number of mapped reads from the sorted bam file, and in some cases, this was higher than the number of original input paired-end reads.
Now this could be due to reads aligning to multiple places, but I don't have any "XA" tags, nor "XT" tags for that matter.
I already saw this post about getting uniquely mapping reads from BWA.
http://seqanswers.com/forums/showthread.php?t=30120
I want to get stats on unique and multi-mapping reads.
Is using MapQ=0 a reliable way of flagging multi-mapping reads?
And why am I not getting XA and XT tags at all?
Here are a couple of lines of my output - I am not seeing any XT/XA tags, nor X0 or X1 tags.
Code:
HWI-ST226:220:D0AU7ACXX:5:1306:2199:194196 83 ref1 1 58 10S78M = 25 -54 CAGACGTGTGCTCTTCCGATCTCGTAAATTTAGTACAACCAAGGGGTAAAAATCCACACTACAACAAGCCGCTCACTATTAGACCCAA * NM:i:0 AS:i:78 XS:i:67 HWI-ST226:220:D0AU7ACXX:5:2206:3059:40310 163 ref1 21 0 88M = 42 109 TAGTACAACCAAGGGGTAAAAATCCAGACTACAACAAGCCGCTCACTATTAGACCCAAATGTCCAAAAACCCTTCACAAAATCGTCAA * NM:i:1 AS:i:83 XS:i:83 HWI-ST226:220:D0AU7ACXX:5:1306:2199:194196 163 ref1 25 40 66M22S = 1 54 ACAACCAAGGGGTAAAAATCCACACTACAACAAGCCGCTCACTATTAGACCCAAATGTCCAAAAACAGATCGGAAGAGCGTCGTGTAG * NM:i:0 AS:i:66 XS:i:66 HWI-ST226:220:D0AU7ACXX:5:2206:3059:40310 83 ref1 42 0 88M = 21 -109 ATCCACACTACAACAAGCCGCTCACTATTAGACCCAAATGTCCAAAAACCCTTCACAAAATCGTCAATCTTCTCACTTTCCTGGCGTC * NM:i:0 AS:i:88 XS:i:88
thanks
Comment