Mapping paired-end stranded RNA-seq data

[hubin.keio]

Hello. I need to map paired-end stranded RNA-seq reads. The library was prepared using ScriptSeq. Can anyone show me how to assign transcripts with the strand information based on mapping? Filtering mapped reads with the 0x10 flag seems not working here, since I got paired-end reads. Thanks.

[swbarnes2]

I don't know why that wouldn't work, unless you are looking for reads that have a flag that is literally 16. In a paired end read, they will certainly have the 1 flag, and either the 128 or the 64. And if they are properly paired, they will have the 2 flag, and they might also have the 16 and 32 (hopefully exactly one of those, but 2 or 0 is possible), and the 4 and or 8 might be flagged too.

[alexdobin]

Most RNA-seq mappers will set the 0x10 bit if the actual sequence from .fastq files is reverse complemented with respect to the "+" genome strand. I believe this is a standard SAM convention. Illimuna PE sequencing produces sequences in such a way that one of the mates agrees with the original RNA, while the other is reverse complementary. Hence, for a correctly (concordantly) paired alignment, the mates will always have inverted 0x10 bits.

If you have stranded RNA-seq data, then you know - from the library construction - which of the mates, 1st or 2nd, agrees with the original RNA. For example, in the standard "dUTP" protocol it is the 2nd mate that is the first end of the RNA fragments. As @swbarnes2 pointed out, 0x40 bit is set for the 1st mate and 0x80 bit is set for the 2nd mate. Here is awk code:

if ( (and(FLAG,0x80)>0 && and(FLAG,0x10)==0) || ( (and(FLAG,0x80)==0 && and(FLAG,0x10)>0) )
{
strand="+";
} else {
strand="-";
};

Again, you have to check which mate has the true RNA strand in your ScriptSeq prep.

[hubin.keio]

Hello @alexdobin and @swbarnes2, thank you so much for your kind response. Your explanation are really useful. After reading the sam format documentation several times and your replies, I started to understand the problem. It looks for my bam file, only the properly-paired reads (0x2) have the right strand information for our spiking control sequence. Without filtering the reads with 0x2, 0x10 flag still gave me wrong strand information in a few cases. Again, thanks for your help.