This is a test version of Biostars. For the public version, visit https://www.biostars.org.
samtools mpileup with strand-specific paired-end reads

I have strand-specific reads from a paired-end RNAseq run. These are the SAM records for spot SRR1269551.1308177 after mapping:

SRR1269551.1308177  99  chrI    35  255 34M =   17  170 CACACCACACCACACCCACACACACACATCCTAA  ?AAAAAAADEEEEDDDGGGGGGHHHHHHIIIIII  NH:i:1  HI:i:1  AS:i:66 NM:i:0  MD:Z:34
SRR1269551.1308177  147 chrI    171 255 34M =   35  -170    CACTCCGAACCACCATCCATCCCTCTACTTACTA  HHHEHHIHHHFCGGGFGGDDDDDDDDBBB??A??  NH:i:1  HI:i:1  AS:i:66 NM:i:0  MD:Z:34

samtools mpileup gives me following result:

[...]
chrI    65      C       1       .       I
chrI    66      T       1       .       I
chrI    67      A       1       .       E
chrI    68      A       1       .$      A
chrI    171     C       1       ^~,     E
chrI    172     A       1       ,       H
chrI    173     C       1       ,       H
chrI    174     T       1       ,       E
[...]

In the pileup the two reads map to different strands. This makes sense from a mapping point of view.

Looking at it from the transcript point of view, both come from the same transcript, thus the same (biological) DNA strand - it's just the technicalities of the paired-end sequencing that make the second read map to the opposite strand.

So what I actually would like to have is a pileup that flips the strand of the second read, so that I get

[...]
chrI    65      C       1       .       I
chrI    66      T       1       .       I
chrI    67      A       1       .       E
chrI    68      A       1       .$      A
chrI    171     C       1       ^~.     E
chrI    172     A       1       .       H
chrI    173     C       1       .       H
chrI    174     T       1       .       E
[...]

I guess some hack like reverse complementing the second read before mapping would do, but to me this problem seems so basic that there must be a better way to do it.

Did I miss some obvious flags during mapping/pileup that could do this for me?

rna-seq

0 answers

No answers yet.

Log in to answer this question.