SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bwa outputs bam with less reads than in input fastq seru Bioinformatics 2 11-04-2016 05:54 AM
How many reads mapped to the genome?(about TopHat output) wisense RNA Sequencing 17 12-30-2013 01:08 AM
BWA output bitwise flag for mapped/unmapped reads wenhuang Bioinformatics 1 08-29-2011 03:54 PM
Bowtie colorspace output has less NT than input reads SongLi Bioinformatics 5 01-20-2011 07:39 AM
BWA Uniquely Mapped Reads NF_seq Bioinformatics 0 09-06-2010 03:32 AM

Reply
 
Thread Tools
Old 08-12-2013, 12:03 AM   #1
Kennels
Senior Member
 
Location: Sydney

Join Date: Feb 2011
Posts: 149
Default BWA-MEM: output mapped reads larger than input reads

I have some bam files generated from BWA-MEM using these command:
Code:
bwa mem -t 8 ref1 R1.fa.gz R2.fa.gz | samtools view -S -h -F 4 -b - | samtools sort - bwaaln-sorted
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.
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
Kennels is offline   Reply With Quote
Old 08-12-2013, 12:53 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Regarding the read number mismatch, are you counting both mates of a pair in one case and not the other? That and multimapping are the most likely causes (I haven't a clue on the XA/XT tag issue).
dpryan is offline   Reply With Quote
Old 08-12-2013, 11:59 AM   #3
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

Caused by chimeric alignments.
lh3 is offline   Reply With Quote
Old 08-12-2013, 03:56 PM   #4
Kennels
Senior Member
 
Location: Sydney

Join Date: Feb 2011
Posts: 149
Default

Thanks.

But why am I not getting any XT,XA,X0 and X1 tags?
I do have some SA tags, maybe 1% of reads.
Any recommendations to extract unique and multi mapping reads? Should I just use mapQ=0?
Kennels is offline   Reply With Quote
Old 08-12-2013, 04:04 PM   #5
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

Just check mapQ.

For long reads and local alignment, X0/X1 do not make much sense any more. If you want multiple mappings, use "-a".
lh3 is offline   Reply With Quote
Old 08-14-2013, 05:41 PM   #6
Kennels
Senior Member
 
Location: Sydney

Join Date: Feb 2011
Posts: 149
Default

Quote:
Originally Posted by lh3 View Post
Just check mapQ.

For long reads and local alignment, X0/X1 do not make much sense any more. If you want multiple mappings, use "-a".
I tried getting all alignments with mapQ>0 and I got 214,694,901 reads (alignments), which should be unique mappings (?) since bwa only outputs the primary hit by default for mapQ > 0 (default options used here).

Just to check, i coupled the read name with the flag and counted how many uniq ones there were;
samtools view -q1 sorted.bam | awk ' { print $1"-"$2 } ' | sort | uniq | wc -l
214694478

and so there are around 420 reads which appear duplicated even with mapQ>0.
Is it possible for a read with mapQ>0 to get reported twice for some reason?

Thanks.
Kennels is offline   Reply With Quote
Old 08-14-2013, 08:33 PM   #7
Kennels
Senior Member
 
Location: Sydney

Join Date: Feb 2011
Posts: 149
Default

Quote:
Originally Posted by Kennels View Post
I tried getting all alignments with mapQ>0 and I got 214,694,901 reads (alignments), which should be unique mappings (?) since bwa only outputs the primary hit by default for mapQ > 0 (default options used here).

Just to check, i coupled the read name with the flag and counted how many uniq ones there were;
samtools view -q1 sorted.bam | awk ' { print $1"-"$2 } ' | sort | uniq | wc -l
214694478

and so there are around 420 reads which appear duplicated even with mapQ>0.
Is it possible for a read with mapQ>0 to get reported twice for some reason?

Thanks.
I had a closer look at the output, and it would appear that these duplicated reads are due to chimeric reads as you mentioned Heng. I understand now what you meant. For example, the reads below were output by '-q1' option from samtools view.

Code:
HWI-ST226:220:D0AU7ACXX:5:1101:11456:146339	2193	Oa_Locus_2615_Transcript_18	2194	44	42M46H	=	2182	-54	AATTGAGCTACCAAAAACCCTAACCCAAAAATTTGTAGCGTC	*	NM:i:2	AS:i:35	XS:i:0	SA:Z:Oa_Locus_2615_Transcript_14,1923,+,21S43M24S,60,0;Oa_Locus_2615_Transcript_14,1976,-,50S38M,55,0;
HWI-ST226:220:D0AU7ACXX:5:1101:11456:146339	2193	Oa_Locus_2615_Transcript_14	1976	55	50H38M	Oa_Locus_2615_Transcript_18	2182	0	GGGGGTTTTAGTAGCTCTCTCTCTGTGGCCAGAAATGG	*	NM:i:0	AS:i:38	XS:i:0	SA:Z:Oa_Locus_2615_Transcript_14,1923,+,21S43M24S,60,0;Oa_Locus_2615_Transcript_18,2194,-,42M46S,44,2;
This the second read in a pair, but gets reported twice even though their mapQ is not 0.
I suppose using the '-M' option will flag both of these reads as secondary.

So to summarize, to get unique mappings:

1. use '-q1' option to get mapQ>0 alignments. This will include reads which are chimeric, and so will include a small level of 'multi-mapping' reads
2. from step 1, filter out reads which have "SA:" tag. This should give unique mappings in the context of bwa-mem

Is this correct?

Last edited by Kennels; 08-14-2013 at 10:30 PM.
Kennels is offline   Reply With Quote
Old 08-15-2013, 12:50 PM   #8
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

That is about right. Just beware that sometimes chimeric reads are informative to break points/SVs.
lh3 is offline   Reply With Quote
Old 09-05-2013, 05:18 PM   #9
pengchy
Senior Member
 
Location: China

Join Date: Feb 2009
Posts: 116
Default

Hi,

still question here:
for this pair reads, i have sort the bam by the query name, they were uniquely and properply mapped. However, with 0 mapQ. According to the above mentioned method, they will be filtered for unique reads calculation.

bwa-0.7.5a mem was used here.

Code:
HISEQ700708:142:C2624ACXX:1:1101:10002:122014   99      scaffold2651    547291  0       100M    =       547537  346     ATCAAGGCCGTCATTACGGCCAGACGTGATTGTTCTGGGTACTGATTTCTCTGGATCTATACACCTAAACTGCGTGAAAATGTAGTGTAGCAATTTTAAT  CCCFFFFFHHHHHJJIHIJJJJJJJJHIIJJJIIJIIJJ?DHIGGGIJJJJJJIIJJJIJIHGHHHFFFFFFEDD@DDDDDDDFEDDDEDDDDDDDDEDD  NM:i:0  AS:i:100        XS:i:100
HISEQ700708:142:C2624ACXX:1:1101:10002:122014   147     scaffold2651    547537  0       100M    =       547291  -346    ATACAGCTAGTTTCTCCTCAGTTCGAATCGGCTTTGACGCGTTGCATGCCTCACTCTCAGTTTGATTTTTCATGAGACCATGCACTTCTTGGAATTGATC  D@:DCDDDCDDBDDCDDDDDDDEDDDDDDFFFFHHGIIJJIJJJJJIIFGJHFIIJIIJJJJJJJJJJJJJIIJIHHJJJIIIHJJJGGHHHFFFFFCCC  NM:i:0  AS:i:100        XS:i:100
pengchy is offline   Reply With Quote
Old 09-05-2013, 05:27 PM   #10
pengchy
Senior Member
 
Location: China

Join Date: Feb 2009
Posts: 116
Default

I understand this, there only output one record for multiple mapped reads.

Quote:
Originally Posted by pengchy View Post
Hi,

still question here:
for this pair reads, i have sort the bam by the query name, they were uniquely and properply mapped. However, with 0 mapQ. According to the above mentioned method, they will be filtered for unique reads calculation.

bwa-0.7.5a mem was used here.

Code:
HISEQ700708:142:C2624ACXX:1:1101:10002:122014   99      scaffold2651    547291  0       100M    =       547537  346     ATCAAGGCCGTCATTACGGCCAGACGTGATTGTTCTGGGTACTGATTTCTCTGGATCTATACACCTAAACTGCGTGAAAATGTAGTGTAGCAATTTTAAT  CCCFFFFFHHHHHJJIHIJJJJJJJJHIIJJJIIJIIJJ?DHIGGGIJJJJJJIIJJJIJIHGHHHFFFFFFEDD@DDDDDDDFEDDDEDDDDDDDDEDD  NM:i:0  AS:i:100        XS:i:100
HISEQ700708:142:C2624ACXX:1:1101:10002:122014   147     scaffold2651    547537  0       100M    =       547291  -346    ATACAGCTAGTTTCTCCTCAGTTCGAATCGGCTTTGACGCGTTGCATGCCTCACTCTCAGTTTGATTTTTCATGAGACCATGCACTTCTTGGAATTGATC  D@:DCDDDCDDBDDCDDDDDDDEDDDDDDFFFFHHGIIJJIJJJJJIIFGJHFIIJIIJJJJJJJJJJJJJIIJIHHJJJIIIHJJJGGHHHFFFFFCCC  NM:i:0  AS:i:100        XS:i:100
pengchy is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 05:07 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO