Mapping bisulfite-treated PE reads - why are both reads always reported to align to the same strand?

Question

[adam.nunn]

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?

[Devon Ryan]

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.

[adam.nunn]

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.

[Devon Ryan]

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.

[adam.nunn]

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?

[Devon Ryan]

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