SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
[BWA] Calculate % of reads not aligned to reference qesecir Bioinformatics 18 11-04-2014 03:03 AM
Convert 1000-Genomes-proje BAM to FASTA (aligned to reference, grouped by chromosome) ce.log Bioinformatics 17 01-14-2014 12:35 AM
BWA - number of sequences aligned Cbon Bioinformatics 7 05-13-2012 04:45 PM
Sam flags for bwa-aligned paired end reads with identical + / - strand coordinates spark Bioinformatics 0 03-09-2011 05:00 AM
the number of aligned reads by maq and bwa totalnew Bioinformatics 3 12-14-2009 12:18 PM

Reply
 
Thread Tools
Old 02-29-2012, 12:32 AM   #1
phatjoe
Member
 
Location: Asia

Join Date: Aug 2011
Posts: 13
Talking Are these reads aligned to my reference (BWA)?

I am currently extracting unaligned reads with the syntax below:-

Code:
samtools view -f 4 *.bam
According to http://picard.sourceforge.net/explain-flags.html, unmapped reads are flagged as 77 or 141.

These are some of the reads that I've extracted
HTML Code:
HWI-ST715:180:D0JHKACXX:4:2307:5231:174176	93	4340.m000924	1187	23	76M	=	1187	0	CGATGGCGCCCAAGGCCGAGAAGAAGCCCGCGGAGAAGAAGCCGGCCTCCGATAAA
CCGGCGGAGGAGAAGGAGAA	<DBBB;@BDCCDBDDDBDDDDDDDDDDDCEDEEECHHGJIIIJJIIIJJJJJJJJHJJJJJJJHHHHHFFFFFCCC	XT:A:U	NM:i:1	SM:i:23	AM:i:0	X0:i:1	X1:i:1	XM:i:1	XO:i:0	XG:i:0	
MD:Z:0A75	XA:Z:4340.m000879,-671,76M,2;
HWI-ST715:180:D0JHKACXX:4:1307:20018:59034	77	4340.m000900	388	0	76M	=	388	0	GTTCTCCCCGGCGAACTCGCGAAGCACGCCGTCTCCGAGGGCACTAAGGCTGTTAC
CAAGTTCACAAGTTCTTGAT	CCCFFFFFHHHHHJJJJJJJJJJJJJJIHHFFDDEEDBDDDDBDCDDDDDDDDDDDDDDDDDDEDDCDCDEEDDDA	XT:A:R	NM:i:1	SM:i:0	AM:i:0	X0:i:2	X1:i:0	XM:i:1	XO:i:0	XG:i:0	
MD:Z:75A0	XA:Z:4340.m000899,+388,76M,1;
HWI-ST715:180:D0JHKACXX:4:1206:2590:114359	77	4340.m000900	391	0	76M	=	391	0	CTCCCCGGCGAACTCGCGAAGCACGCCGTCTCCGAGGGCACTAAGGCTGTTACCAA
GTTCACAAGTTCTTGATCCG	CCCFFFFFHHHHHJJJJJJJJIJIIIJHEFEEEDBDDDDDDDDDDDDDDDDDDDCDDCDEEDDDDCCDECDDDDDA	XT:A:R	NM:i:4	SM:i:0	AM:i:0	X0:i:2	X1:i:0	XM:i:4	XO:i:0	XG:i:0	
MD:Z:72A0T0G0A0	XA:Z:4340.m000899,+391,76M,4;
From the XT:AU and XT:AR tags, does this mean that these reads are aligned at unique or at repeated regions? Or are they really unaligned reads? I was under the impression that for unaligned reads, there will be a "*" at the reference column.


Thanks

Regards,
Joanne
phatjoe is offline   Reply With Quote
Old 02-29-2012, 10:41 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

There are lots of flags that indicate unmapped reads, not just 77 or 141. If the 4 flag is present, you can't necessarily believe anything else in the sam line. And no, many unmapped reads will NOT show a * in the reference column. SAM specs say that unmapped reads should be given the mapping coordinates of the rimapped mates, where possible.

But that's not necessarily what's happening here. The flag of 93 indicates that neither end mapped, and if that were really true, both reads would have no mapping coordinates. But that one does. The answer, I bet, is that the reads are hanging off of the sequences. It's a known issue in bwa that it concatenates all your reference sequences, and if a read hangs off of the edge, it will give the mapping coordinates based on the where the read starts, and it will also put the 4 flag in there so you know that something is off.

So they really were mapped, but that edge effect is causing the 4 flag to be thrown, and it can be hard to automatically know what exactly that means for those reads.
swbarnes2 is offline   Reply With Quote
Old 03-01-2012, 05:13 PM   #3
phatjoe
Member
 
Location: Asia

Join Date: Aug 2011
Posts: 13
Talking

Quote:
Originally Posted by swbarnes2 View Post
There are lots of flags that indicate unmapped reads, not just 77 or 141. If the 4 flag is present, you can't necessarily believe anything else in the sam line. And no, many unmapped reads will NOT show a * in the reference column. SAM specs say that unmapped reads should be given the mapping coordinates of the rimapped mates, where possible.

But that's not necessarily what's happening here. The flag of 93 indicates that neither end mapped, and if that were really true, both reads would have no mapping coordinates. But that one does. The answer, I bet, is that the reads are hanging off of the sequences. It's a known issue in bwa that it concatenates all your reference sequences, and if a read hangs off of the edge, it will give the mapping coordinates based on the where the read starts, and it will also put the 4 flag in there so you know that something is off.

So they really were mapped, but that edge effect is causing the 4 flag to be thrown, and it can be hard to automatically know what exactly that means for those reads.
Thanks swbarnes2! I think you are right, as my reference is just CDS sequences, the edge-effect would have been more apparent! The scenario was also described here (http://bio-bwa.sourceforge.net/) Is it advisable to use BWA for transcriptome mapping in that case?

Regards,
Joanne
phatjoe is offline   Reply With Quote
Reply

Tags
bwa, samtools bitwise flags, unaligned reads

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 12:59 AM.


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