Methyl Dackel results compared to Hisat3n-table

[prairie-guy]

I’m using Version: 0.6.1 (using HTSlib version 1.21)

I'm doing BS analysis and using MethylDackel extract to call methylated/unmethylated regions. I'm looking at rRNA, looking at chr 28SrRNA where there are two known positive sites, 3779 and 4444 using hg38.

When I use hisat3n-table and samtools mpileup I get the predicted and exact results. Moreover, I can see the actual reads with IGV.

When I run MethylDackel with the most permissive settings, it reports 0 C and 0 T at each position.

$ rga '3779|4444' md_test.cytosine_report.txt 
28SrRNA 3779    +       0       0       CHH     CAA
28SrRNA 4444    +       0       0       CG      CGA

MethylDackel extract  --keepSingleton -q 1 -p 1 --keepDupes --keepDiscordant --cytosine_report --CHH --CHG 

I’m wondering if this is how MD treats overlapping PE reads.

Do you think this may be the problem? Do you have any suggestions?

Thanks for your help in advance.

[dpryan79]

What do the MAPQ scores look like there? For overlapping PE reads it will alter the base call phred score up or down depending on whether the overlapping portions agree or disagree (this is similar to what samtools does).

[prairie-guy]

Thank you for your response. I'm new to MethylDackel, so this is probably due to some mistake I am making.
(For my analysis, 28SrRNA:4444 is important as it is a known methylated C and is an important positive control)

Below are the quality scores. I've also attached a small subsampled bam file so you can verify if need be.

When I use samtools mpileup , histat-3n-table and IGV, I can see the reads, yet I'm not getting a signal from MethylDackel extract. What am I missing?

Thanks in advance.

$ samtools mpileup -f ../../../reference/fasta/rRNA.fa --output-MQ  --output-BP -r 28SrRNA:4444-4444 3.sub.bam

28SrRNA 4444 C 12 .,...,,t,,T^], qaqqaqqqaqIE ]]]]]]]]]]]] 104,79,68,65,38,38,30,25,23,23,13,1

MethylDackel extract ../../../reference/fasta/rRNA.fa 3.sub.bam -o 3.sub.MD --cytosine_report --CHH --CHG
rga 4444 3.sub.MD.cytosine_report.txt

28SrRNA 4444 + 0 0 CG CGA

3.sub.bam.gz

[prairie-guy]

I am trying to use MethylDackel extract instead of hisat3n-table for my BS pipelines, due not only to an order of magnitude increase in speed, but to MD’s apparently broader usage in other published pipelines.

Before doing so, I have been back testing MethylDackel against results obtained with hisat3n-table. I’ve posted one such result above, where MD does not pick up two key positive converted control reads in human 28sRNA sequences.

As I work through this on a Cytosine by Cytosine basis, comparing exact Position and Converted and Unconverted Counts from each extractor, the overall results are generally close, but not as close as I would have expected.

I’m beginning to believe that this must be due to a fundamentally different counting algorithms. The hisat3n aligner adds tags to the bam file Yf:i: and Xf:i: to capture Converted and Unconverted Reads per read, respectively. These are then used by host3n-table to extract the counts. http://daehwankimlab.github.io/hisat2/hisat-3n/

As far as I can tell, MethylDackel uses a different algorithm and does not utilize these tags.

Am I correct in this assessment? I would assume that this would be a fairly big change to MethylDackel. If not, it would be a great addition. I haven’t programmed in C for 20+ years, but I would be happy to help out to the extent I can.

Thanks for your help and a great program!

Bryan