radmeth is too slow

[guilhermesena1]

It takes several hours in the human genome to test for all 28M CpGs on a two-case test with 10 replicates each. According to perf, half of the time in the code is spent on __ieee754_log_fma, which I think has to do with log-likelihoods within GSL. We may be able to optimize these.

[moqri]

@andrewdavidsmith, I'm running this on 3 case vs 3 control whole genome on a 96 x CPU took more than an hour to finish chr1.

dnmtools radmeth -factor case design_matrix.txt proportion_table.txt >output.bed

ChatGPT suggests this parallelization approach - any thoughts?

List chromosomes, adjust as needed

chroms=("chr1" "chr2" "chr3" "chrX" "chrY")

Loop through chromosomes in parallel

parallel -j 4 '
grep {1} proportion_table.txt > pt_{1}.txt
dnmtools radmeth -factor case design_matrix.txt pt_{1}.txt > output_{1}.bed
' ::: "${chroms[@]}"

[andrewdavidsmith]

Hi @moqri This shouldn't take that long. radmeth doesn't have a multi-threaded more. You can definitely break it by chromosome, but if it's taking that long for just chr1 then there is probably something else happening. What's the total size of the data? I'm happy to look directly at your data -- just contact me by email. Overall for a 3x3 and the human genome, it shouldn't be taking as long as you are seeing. This is likely easily fixable.

[moqri]

Thank you Andrew,

The size of proportion-table.txt is 1.6 GB. The data is from 3x3 PacBio samples from All of Us and I'm processing on their platform and not able to share but will need to scale to a few thousands samples so might become an issue down the road.

I'll keep trying a few things and will post here if I find a solution.

[moqri]

Update on splitting the proportion_table.txt:

I broke the proportion_table.txt into 90 pieces by rows and 88 of them was processed in around 20 minutes but a few of the pieces are still running after an hour. In particular the piece starting at position chr1:122316104 is the slowest one.

I hope this help.

[moqri]

Additional observation:

The reason the split starting at chr1:122316104 (containing ch1 centromere) was much slower is I had lots of reads mapped to this piece (likely mapping issues) so radmeth speed seem to directly relate to the number of reads

[andrewdavidsmith]

Thanks @moqri for this issue. I improved the performance of radmeth by tightening up the log-likelihood / gradient code in the parameter estimation, and by adding support for threads. I also added arguments so users can set max iterations in the estimation and adjust numerical tolerance which determines convergence.

There are two known performance issues you might encounter:

• The threads operate on "chunks" of consecutive sites, and the threads join periodically. This means if there are places in the genome where things slow down, one thread might cause the others to wait. This relates to the next point.
• Sequencing depth. The beta-binomial regression model inherently uses the numbers of reads, which is a strength. However, at high read counts, e.g. in the hundreds, not only will it slow down, but it also is likely unnecessary. "Jackpot" situations, where some parts of the genome over-amplify and get lots of reads, might result in one thread seeing lots of extreme high-read sites.

If these become a problem, there are two apparent simple solutions to try:

• Have threads operate on interleaved sites (taking a very small hit to locality), and
• Scale down extreme high read counts to, say, a maximum of 100 reads. This should not impact accuracy.

This issue is closed with PR #279. Please open another issue if you find that performance needs further improvements.