Add methurator tool

[edogiuili]

In this PR, I add methurator as an additional QC tool for both the Bismark and bwa-meth subworkflows.

Methurator cannot currently be used with the TAPS/bwa-mem subworkflow because it relies on MethylDackel to extract CpG methylation information from BAM files.

Methurator is a Python package designed to estimate CpG sequencing saturation in DNA methylation sequencing data. It complements Preseq by estimating the number of sequencing reads required to achieve coverage of at least N CpGs with a minimum read depth of r. Because the estimator is based on a frequency-of-frequencies matrix derived from methylation count files, it can also be used to estimate library complexity for RRBS datasets.

Notes

• Methurator currently does not generate version YAML files. Instead, it outputs a topic channel, following the latest nf-core guidelines.
• Methurator v2.1.1 (the version used in this PR) does not yet have MultiQC support. A corresponding MultiQC PR has already been opened: Update methurator module MultiQC/MultiQC#3479.
• The fasta_index_methylseq lint test currently fails because I had to introduce a new methurator boolean parameter to trigger generation of the FASTA.FAI index. This index is required by Methurator due to its dependency on MethylDackel.

PR checklist

• This comment contains a description of changes (with reason).
• If you've fixed a bug or added code that should be tested, add tests!
• If you've added a new tool - have you followed the pipeline conventions in the contribution docs
• If necessary, also make a PR on the nf-core/methylseq branch on the nf-core/test-datasets repository.
• Make sure your code lints (nf-core pipelines lint).
• Ensure the test suite passes (nextflow run . -profile test,docker --outdir <OUTDIR>).
• Check for unexpected warnings in debug mode (nextflow run . -profile debug,test,docker --outdir <OUTDIR>).
• Usage Documentation in docs/usage.md is updated.
• Output Documentation in docs/output.md is updated.
• CHANGELOG.md is updated.
• README.md is updated (including new tool citations and authors/contributors).

[nf-core-bot]

Warning

Newer version of the nf-core template is available.

Your pipeline is using an old version of the nf-core template: 3.5.1.
Please update your pipeline to the latest version.

For more documentation on how to update your pipeline, please see the nf-core documentation and Synchronisation documentation.

[edogiuili]

Tests 6,7,8,9 and 10 (and 3 only for conda profile) fail but I haven't changed them or anything related to their workflows. Might be that these tests were already failing?

[FelixKrueger]

Thanks for this, @edogiuili — nice addition, and the module structure follows the current conventions well (topic-channel versions, on-demand .fai generation, both module and pipeline-level tests). I checked it out and ran the bismark_run_methurator pipeline scenario end-to-end (Nextflow 25.04.0, Docker); it completes and produces the .yml summaries and interactive .html plots as intended. A few things to address before merge, one of which I think is important.

1. compute_ci = true default + errorStrategy = 'ignore' → silent, non-reproducible QC gaps ⚠️

In my run, one of the four samples (SRR389222_sub1) failed inside methurator's confidence-interval bootstrap and produced no output, but because METHURATOR_GTESTIMATOR has errorStrategy = 'ignore', the pipeline still reported success. The snapshot assertion then failed because the committed snapshot does contain SRR389222_sub1.yml — i.e. this sample succeeded on your machine but not on mine.

The crash:

run_estimator.py:274 → preseqR_rSAC_bootstrap → estimator(1)
ValueError: setting an array element with a sequence. The requested array has an
inhomogeneous shape after 1 dimensions. The detected shape was (30,) + inhomogeneous part
(preceded by RuntimeWarning: invalid value encountered in divide, rational_function.py:40)

I confirmed by re-running gt-estimator on that sample's BAM with the flag toggled:

• with --compute_ci → crashes (as above)
• without --compute_ci → completes and writes the summary

So compute_ci (defaulted on in this PR) is the trigger, and it fails on sparse / low-CpG samples — which are exactly the samples a saturation analysis is aimed at. The errorStrategy = 'ignore' then hides it, so users get a silently incomplete methurator report, and the pipeline-level snapshot becomes environment-dependent (same pinned container, different numpy/platform arithmetic near a divide-by-zero → different pass/fail). That will flake in CI.

Suggestions:

• Default compute_ci = false (make CI opt-in), or gate it on a minimum coverage / CpG count.
• Please report the crash upstream to methurator — it should degrade gracefully (skip CI + warn) rather than throw on a degenerate bootstrap. (methurator:2.1.1, gt_utils/run_estimator.py:220.)
• At minimum, emit a logged warning when a sample is dropped, so the silent gap is visible.
• Reconsider snapshotting the full output file-set of an errorStrategy='ignore' process — assert on stable per-sample outputs, or exclude CI-dependent files.

2. Housekeeping before merge

• CHANGELOG.md is not updated (required).
• This branch is ~19 commits behind dev and carries a stale template artifact — nextflow.config references conf/aws/batch/nextflow.config (the aws_batch profile) but that file isn't in the branch, which is a hard config-parse error under Nextflow ≥25.10's strict parser. A rebase on current dev should clear it.
• Unrelated nextflow_schema.json churn: the diff deletes the empty bwameth_options/bwamem_options groups and reorders generic_options↔institutional_config_options. Unrelated to methurator — better split out to keep the diff clean.

3. Minor

• METHURATOR_PLOT has no errorStrategy — if plotting fails after a successful gt-estimator, the whole run fails. Consider 'ignore' for consistency with the estimator.
• error("ERROR: --run_methurator can't be running using TAPS workflow.") is triggered by params.taps || aligner == 'bwamem' but only mentions TAPS, and reads awkwardly. Suggest: "--run_methurator is not supported with the TAPS / bwa-mem workflow (methurator relies on MethylDackel)."
• Param names minimum_coverage, t_max, compute_ci are generic for global scope; a methurator_-style prefix would age better. run_methurator is also missing an fa_icon in the schema.
• docs/output.md: the summaries bullet trails off ("…extrapolation results for each"), the plots bullet calls a single .html a "folder", and the section is an ## among ### tool siblings.

[edogiuili]

Hi Felix,

Thanks for the review and the suggestions! Below you can find the points I addressed.

1. compute_ci = true

• Changed default to compute_ci = false
• Updated Methurator to v2.2.0 to report crashes upstream — now emits a warning when a sample is dropped instead of failing silently

2. Housekeeping before merge

• Branch rebased/merged and up to date with target
• CHANGELOG updated

3. Minor changes

• METHURATOR_PLOT now has an errorStrategy defined
• Fixed error message wording when activating methurator using TAPS/bwamem workflow
• Added methurator_ prefix to Methurator params for consistency
• Fixed docs — kept ## (not ###) for the Methurator section to match the heading level used by other tools (Preseq, BamQC)