Aggregation step is paused shortly after beginning

[daggoo]

Hi, I'm experiencing difficulties with the aggregation step during modification calling. I ran Megalodon twice and for both runs aggregation is paused after processing approximately a third of one chromosome (based on the resulting bed file). Megalodon doesn't crash but there is no output for more than 10 hours, bed file is not accessed and one of the the threads is using 100% of one processor. I have local ssd disks. Could you please tell what could cause this behaviour?

[marcus1487]

I have not seen this behavior in recent versions of Megalodon. I'm not sure exactly what the cause could be at this point. Most of the processing for this step is I/O dependent, so if the Megalodon output directory is on local disk I'm not currently seeing a potential source for this issue.

Have you tried to re-run aggregation after the run is stopped manually (megalodon_extras aggregate run --outputs mods)? Does the same behavior occur?

[daggoo]

Yes, I redid indexing with megalodon_extras modified_bases index_database and tried to re-run aggregation with megalodon_extras aggregate run --outputs mods but had the same problem.
I will try to run it on another machine then. Thanks for the answer anyway!

[marcus1487]

I've had a dig around and found that this may be related to a bug in python that was fixed in python 3.8 (python/cpython#10305). I'll have a go at providing a fix for python versions <3.8, but could you try to upgrade your python version to see if it fixes this issue?

[marcus1487]

I've just pushed a commit to add an option to lower the batch size in the aggregation step (and lowered the default from the previous hardcoded value). Could you test this new version and let me know if this solves the issue (install with pip install git+https://github.com/nanoporetech/megalodon). You do not have to re-run megalodon, just the re-aggregation step will suffice. If this still fails, could you try running with the setting: --aggregate-batch-size 10. Hopefully this resolves this issue.

[marcus1487]

One additional question, did you have the --write-mod-log-probs flag set? This would increase the memory accumulated in a batch and could potentially cause this issue.

It looks like the core of the issue is a memory cap on the multiprocessing module that results in an uncatchable error (as it occurs at the interface of two processes and leaves both stalled). This is fixed in python 3.8, but as a Megalodon solution the --aggregate-batch-size argument has been added for cases when this happens. Please re-open this issue if this is not fixed with either of these recommendations.

[daggoo]

Hello @marcus1487, apologies for my late response. Yes, I had --write-mod-log-probs flag set. When I upgraded to python 3.8 and re-ran aggregation step without this flag the issue persisted but the new commit has fixed it even without setting --aggregate-batch-size 10. Thank you very much for your help.

[marcus1487]

The new commit lowered the default batch size from 1000 to 200 (the recommendation of 10 was to be super cautious). Glad the new default works well for your run.

I might consider removing the --write-mod-log-probs option in a future release to avoid this issue becoming more common. Did you have a particular use case for this output? Do you think the same output from the database API would suffice for access to this data in the future?

[daggoo]

No, not really. I was working only with the coverage and modified proportions by now but was curious to see the probabilities. Do you mean megalodon_extras per_read_text modified_bases by database API?

[marcus1487]

No. I mean the megalodon.mods.ModsDB python interface. The per-read database output file (per_read_modified_base_calls.db) is an sqlite database indexed by reference position. I had thought I documented this somewhere, but it does not appear to be the case. As a quick start here is a script to interact with the database (TODO add this to the official docs):

from megalodon import mods, megalodon_helper as mh

# open DB for reading
mods_db = mods.ModsDb('mega_res/per_read_modified_base_calls.db')

# iterate over log likelihood ratios by reference position
for (chrm, strand, pos), pos_lps in mod_db.iter_pos_scores(convert_pos=True, compute_llrs=True):
     for mod_base, lps in pos_lps.items():
         print(chrm, strand, pos, mod_base, lps)

# extract statistics from specific position of interest
chrm, strand, pos = 'chr1', '+', 100000
mods_db.get_pos_stats(mods_db.get_pos_dbid(chrm, mh.str_strand_to_int(strand), pos))

[daggoo]

I understand now. Yes, I think the database API will suffice to access this data. It is convenient. Thanks for the code!