Question: Mapping bisulfite-treated PE reads - why are both reads always reported to align to the same strand?
0
gravatar for adam.nunn
16 months ago by
adam.nunn0
adam.nunn0 wrote:

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 • 331 views
ADD COMMENTlink written 16 months ago by adam.nunn0

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.

ADD REPLYlink written 16 months ago by Devon Ryan92k

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.

ADD REPLYlink written 16 months ago by adam.nunn0

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.

ADD REPLYlink written 16 months ago by Devon Ryan92k

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?

ADD REPLYlink written 16 months ago by adam.nunn0

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

ADD REPLYlink written 16 months ago by Devon Ryan92k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1287 users visited in the last hour