Reads mapped on wrong strand

Question

[qstefano]

Hello there,

I'm working on Coli MG1655, after i have aligned the reads to reference, i checked in IGV if alignments are coherent to gene expression but i noticed that despite some genes are mapped on forward strand, the reads has been aligned to the reverse strand.

i used samtools to divide reads mapped on forward and reverse strand in this way:

samtools view -F 16 -b -o forward_strand.bam in.bam samtools view -f 16 -b -o reverse_strand.bam in.bam

I am attaching an example of groL gene, mapped to the forward strand, but reads actually mapped on reverse strand. Why? what am i missing?

In the image, from top to bottom: gff annotation, coverage reverse, coverage forward, mapped reads to reverse strand.

[Pierre Lindenbaum]

DNA is fragmented regardless of the 'genes'. "High-throughput sequencing (HTS) methods" in https://en.wikipedia.org/wiki/DNA_sequencing

[Carlo Yague]

I think the OP data is RNA-seq (based on the tags), so depending on the library preparation, he could expect strandness.

@qstefano: what library preparation was used to generate the data ? It can be unstranded, stranded, or reversely stranded. If it is unstranded - or stranded but mapped as "reversely stranded" -, then it is normal that reads map on the "wrong" strand.

[qstefano]

The library is stranded

Answer 1

[A. Domingues]

Depending on the library preparation method (dUTP based), this could be completely normal.

[qstefano]

Why? i would have expected to see these reads aligned to forward strand due to gene location, something fundamental escapes me?

[A. Domingues]

In the reverse strand library prep (also know as dUTP based) what gets sequenced is not the original RNA but the reverse-complement. Thus when it gets mapped it goes to the opposite strand. This diagram might explain better: http://onetipperday.sterding.com/2012/07/how-to-tell-which-library-type-to-use.html

[Carlo Yague]

when it gets mapped it goes to the opposite strand

Depending on the alignment program used, there can be an option to correct for that during the mapping. For instance in hisat2, we can use -rna-strandness R to specify that the single-end RNA-seq data is reversely stranded. The alignment resulting from that will be in the correct orientation.

[A. Domingues]

I was not aware of that, good to know. I normally use one of STAR/bowtie/bwa, and AFAIK, there is no option to correct on those (I am happy if someone tells me I am wrong). For pseudo-alignment, salmon will autodetect the library protocol.

[qstefano]

In this case the library has been prepared with NEBNextUltra Directional, I didn't know it was based on dUTP method, so this is the answer.

thanks a lot, now is clear.

[A. Domingues]

No problem (you can click the accept bottom if you are happy with the answer). One more note, if you are going to count reads mapping to genes with something like featureCounts or HTseq-count, those tools will have have a library flag to be set to the library protocol - to count the reads correctly.