This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Mapping bisulfite-treated PE reads - why are both reads always reported to align to the same strand?

I am learning about bisulfite sequencing and have been attempting to map some artificially generated reads back to the original genome using erne-bs5 (among others). For some reason this tool always reports mate 1 and mate 2 in each case as both being either on the forward strand or the reverse strand, in the same direction. The positions are usually correct when compared to other aligners but the orientation is wrong. What could be going on here?

alignment
• 1,508 views • link
written by adam.nunn • adam.nunn

It'd be helpful if you provided an example with the original fastq entries. In general, people try to follow how bismark set its flags, since there's a lot of software out there that assumes that as input.

• 0 views • link
written by Devon Ryan Devon Ryan

Sure thing! Here is an example read and resulting comparison between erne bamfile and segemehl bamfile. All the reads are orientated like this in erne.

MATE 1:

@6017897_AZNP01000001.1:2332-2536
ATGGTTTTTTTATTTTTTTATTTTATTAATTTAAAGAAATTGATGATTGTATTGATATAATTGTTTTTGTTTTTAGAGGAAATTGGTGGAAAGGTAATTTATTATATTTTA
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGGGGGGGGGGGGGFFFFFFFFEEEEEEEDDDDDCCCCCBBBAAAA@@@???>>>===<<;;:::998877655

MATE 2:

@6017897_AZNP01000001.1:2332-2536
TTAATCTTAATAATACAAAAAACCTACAATAAAAAAACAATCTATTATTAATAAATTAAAAATAATTAAATAATAAAAATTATAATAAACTATTAAAATATAATAAAATAC
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGGGGGGGGGGGGGFFFFFFFFEEEEEEEDDDDDCCCCCBBBAAAA@@@???>>>===<<;;:::998877655

SEGEMEHL

6017897_AZNP01000001.1:2332-2536    163 AZNP01000001.1  2333    255 111M    =   2426    -18 TTAATCTTAATAATACAAAAAACCTACAATAAAAAAACAATCTATTATTAATAAATTAAAAATAATTAAATAATAAAAATTATAATAAACTATTAAAATATAATAAAATAC HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGGGGGGGGGGGGGFFFFFFFFEEEEEEEDDDDDCCCCCBBBAAAA@@@???>>>===<<;;:::998877655 NM:i:0  MD:Z:111    NH:i:1  XI:i:0  XA:Z:P  XB:Z:F1/GA  XD:i:21 XF:i:0
6017897_AZNP01000001.1:2332-2536    83  AZNP01000001.1  2426    255 111M    =   2333    18  TAAAATATAATAAATTACCTTTCCACCAATTTCCTCTAAAAACAAAAACAATTATATCAATACAATCATCAATTTCTTTAAATTAATAAAATAAAAAAATAAAAAAACCAT 556778899:::;;<<===>>>???@@@AAAABBBCCCCCDDDDDEEEEEEEFFFFFFFFGGGGGGGGGGGGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH NM:i:1  MD:Z:14A96  NH:i:1  XI:i:0  XA:Z:P  XB:Z:F1/GA  XD:i:30 XF:i:0

ERNE

6017897_AZNP01000001.1:2332-2536    177 AZNP01000001.1  2333    60  111M    =   2426    -18 TTAATCTTAATAATACAAAAAACCTACAATAAAAAAACAATCTATTATTAATAAATTAAAAATAATTAAATAATAAAAATTATAATAAACTATTAAAATATAATAAAATAC HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGGGGGGGGGGGGGFFFFFFFFEEEEEEEDDDDDCCCCCBBBAAAA@@@???>>>===<<;;:::998877655 RG:Z:1530617709 PG:Z:ERNE   NM:i:0  NH:i:1  IH:i:1  HI:i:1
6017897_AZNP01000001.1:2332-2536    113 AZNP01000001.1  2426    60  111M    =   2333    -204    TAAAATATAATAAATTACCTTTCCACCAATTTCCTCTAAAAACAAAAACAATTATATCAATACAATCATCAATTTCTTTAAATTAATAAAATAAAAAAATAAAAAAACCAT 556778899:::;;<<===>>>???@@@AAAABBBCCCCCDDDDDEEEEEEEFFFFFFFFGGGGGGGGGGGGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH RG:Z:1530617709 PG:Z:ERNE   NM:i:1  NH:i:1  IH:i:1  HI:i:1

That's interesting I didn't realist that Bismark set it's sam flags differently to standard format.

• 0 views • link
written by adam.nunn • adam.nunn

I would suggest that erne5 is broken, neither it nor segemehl are producing reasonable output. I would strongly encourage you to use more standard aligners like bwa-meth or bismark.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks for the recommendation. I am also investigating the data with bwa-meth and bismark (among others). Out of interest, what is unreasonable about the segemehl output?

• 0 views • link
written by adam.nunn • adam.nunn

The TLEN columns are nonsensical, they should be 204 and -204.

• 0 views • link
written by Devon Ryan Devon Ryan

0 answers

No answers yet.

Log in to answer this question.

More posts like this