Found two separate (pre-existing, unrelated to my SAVANA PR) issues in lrsomatic while testing with real patient data, flagging for whoever owns PREPARE_REFERENCE_FILES/model selection: 1. No HAC entries in the Clair3/ClairS model maps clair3_modelMap / clairs_modelMap in workflows/lrsomatic.nf only have entries for ..._sup@v4.x / v5.x basecaller models. Real data basecalled in HAC mode, e.g. dna_r10.4.1_e8.2_400bps_hac@v4.3.0, which is what our actual data uses, has no match, so auto-detection silently resolves to nothing. 2. Manual clair3_model override breaks model-download staging when auto-detection doesn't resolve In subworkflows/local/prepare_reference_files.nf52 def id_new = basecall_model_meta ? clair3_modelMap.get(basecall_model_meta) : basecall_model_meta def model = (!meta.clair3_model || ...) ? clair3_modelMap.get(basecall_model_meta) : meta.clair3_model model (used for the actual download URL) correctly falls back to the samplesheet's explicit clair3_model override. id_new (used as meta.id for UNTAR's output directory name) does not, it's only ever derived from the map lookup. So when a sample's basecall model isn't in the map (issue #1, or just missing metadata) and you supply an explicit override, the download itself succeeds but UNTAR crashes: mkdir: missing operand because the staging directory name comes out empty.

I guess these are both due hac not being in there? Initially this was purposely done by us because we do not see the benefit of basecalling with anything less than SUP. However, ONT is going to stop releasing newer sup models and the hac ones are getting closer in performance (although not close enough yet in our opinion).

You're welcome to also add a pr to add hac as possibilities, although im not sure we would recommend it especially for SNV calling etc

Makes total sense on the HAC/SUP tradeoff, that's your call to make, not something I'd push on ;). Just want to flag that issue #2 is actually independent of HAC specifically, it'll break the same way for any basecall model string that isn't (yet) in clair3_modelMap, regardless of accuracy mode. To my understanding the model var (used for the download URL) falls back to the samplesheet's explicit clair3_model override when the map lookup misses, but id_new (used as the staging directory name for UNTAR) doesn't have that same fallback, it's only ever clair3_modelMap.get(basecall_model_meta). So it'll mkdir: missing operand crash on: - a brand-new SUP model ONT ships before the map gets updated for it - any sample whose BAM header is missing/has malformed basecall_model= metadata - HAC data, same as we saw So worth a small independent fix regardless of the HAC decision, just making id_new fall back to meta.clair3_model the same way model already does, so an explicit override always works even when the map doesn't have a match. Not proposing to touch the HAC-support question at all.

Sure just make a separate PR and we'll look into it

Hey! Now that SAVANA's merged, I wanted to run something else by you: I've built a differential-methylation-at-CpG-islands module and I think it'd be a good fit as an optional step inside lrsomatic itself, alongside modkit. Curious whether you'd be open to it. What it does (standalone Nextflow entry today, -entry DMR_ONLY): 1. DMR_CALLING: takes a treatment/control bedMethyl pair, restricts to CpG islands with real read coverage in both samples (so it never compares real data against absence-of-data), then runs modkit dmr pair for the actual DMR calls and annotates each with its nearest GENCODE gene. 2. DMR_RANK_PROMOTERS: ranks promoter-proximal DMRs. 3. DMR_LOCUS_PLOTS: per-locus plots from bigWig tracks. 4. DMR_REPORT: HTML/PDF report tying it together. It's nf-test covered end to end (calling, ranking, plots, report, and the entry workflow itself) and there's a written-up method (CpG-island restriction logic, modkit dmr pair params) I can send over. Why I think it fits lrsomatic specifically: the input is exactly what modkit_pileup already produces, {tumor,normal}/modkit_pileup/sample.bed.gz, bgzip bedMethyl. DMR_CALLING/DMR_RANK_PROMOTERS could run directly off that with no format changes. The only gap is bigWig for DMR_LOCUS_PLOTS (not needed for calling/ranking, just the plots), that's a small bedMethyl→bigWig conversion step, tools for it are already in our image. One caveat worth flagging up front: this is a differential comparison between two samples, and that's only meaningful when both sides are the same tissue/cell type. In my own data the "normal" is PBMC, not tissue-matched normal, so tumor-vs-PBMC methylation differences would be dominated by cell-type identity, not tumor biology, that comparison isn't valid and shouldn't run by default. In practice, for my case, what I actually want to compare is tumor sample 6 vs tumor sample 1 (two tumor regions from the same patient), not tumor vs normal at all, so the module's "treatment vs control" framing is already more general than "tumor vs normal," and whatever we land on should keep that generality explicit rather than assuming a tumor_bam/normal_bam pairing. (Example in PDF attached). That's actually already how my DMR_ONLY entry is set up, and I think it's the right shape for lrsomatic too: a separate comparisons samplesheet (dmr_comparisons_csv on our side, columns label,treatment_bedmethyl,control_bedmethyl,treatment_bigwig,control_bigwig), decoupled from the main sample sheet, where each row names which two bedMethyl outputs to compare and under what label. That sidesteps the tumor/normal assumption entirely, a row can point at tumor6 vs tumor1, tumor vs tissue-matched normal, or whatever the user actually wants, and there's no default pairing to get wrong. Would this make sense as a new --skip_dmr-style optional module in lrsomatic, built around a comparisons samplesheet like that? Just want your read on where it belongs before I put together a PR.

Hi! I think we've been looking into these things but for between haplotypes. We'll discuss it a bit internally and come back to you! Thanks 😄

@Robert Forsyth How is it going with CoRAL + BFBArchitect? Happy to help if needed!

Hi @Yann Vanrobaeys we were thinking about how to best implement this. The way you propose with performing DMR between different tumour samples seems more like a step downstream of the pipeline to us. We want to have peolpe be able to run dozens of samples at a time and the combinatorial DMRs for all samples would be a lot of compute and output. What we would like to implement is to perform these steps between haplotypes and, if present, tumour/normal pair. Unless you have a way on how to do it between samples in a clean way.

My “clean way” was to ask the user to prepare a samplesheet.csv to setup all the DMR of interest between samples. What do you think?

Yea not sure if it will be good to add more samplesheet (columns) into an already heavy samplesheet in its current state. Think I would vote for this being more of an ad-hoc post lrsomatic analysis

Fair enough, no problem. I can start on DMR between haplotypes if you want? Have you also considered polyploidy tools (more than the two parental haplotypes)?

Yea between haplotypes would definitely be something we would like to have, so you're very welcome! If you also think that between tumour/normal (if a normal is present in samplesheet) is of interest, maybe that's also an option. Not sure what you think of that

Between haplotypes: needs zero new infrastructure. --modkit_phased already emits sample_hp1.bed.gz/sample_hp2.bed.gz per sample, that's exactly the input pair DMR_CALLING needs, so it's just modkit dmr pair -a hp1 -b hp2 with the module I already built. One comparison per sample, no combinatorics, no samplesheet changes. Happy to put together a PR. Between tumor/normal: same deal if you want it, normal_bam is already in the samplesheet, so pairing is automatic, one comparison per sample, no new columns. One caveat: only valid when normal is tissue-matched. My own "normal" is PBMC, and tumor-vs-PBMC would mostly reflect cell-type identity, not tumor biology, so I'd gate this behind an opt-in flag rather than firing whenever normal_bam is populated. Between arbitrary tumor samples (my tumor6-vs-tumor1 case): agreed, that one's combinatorial and doesn't belong in lrsomatic. Dropping that ask, it already works as an ad-hoc downstream step in our pipeline (NCP's DMR_ONLY, fed directly from lrsomatic's sample.bed.gz), no changes needed on your end. So: haplotype in lrsomatic (low-cost, yes), tumor/normal as opt-in with the tissue caveat (your call), cross-sample stays downstream. Sound right?

Maybe I would then simply say, let's do the haplotype one and leave the tumour/normal out currently?

Deal, I will start working on the PR.

Regarding polyploidy detection tool, have you guys started thinking about it at all?

Do you mean with this WGD? Or just ploidy estimation in general?

We have ploidy estimates from wakhan and/or ascat. We were testing longphase-TO to perform tumour-only optimzed haplotype phasing

I was wondering since phasing produce for example bam_1 and bam_2, could we have more than just the two parental for clonal evolution?

Do you mean like this:

That is the idea of longphase-t and longphase-to.

Thanks, I think I was confusing phasing with clone deconvolution.

There are some tools for it but none that we considered to actually run in the pipeline

did you have anything in mind?

No but it's something I am interested and will investigate and test it on my own data before proposing anything 😉

Sounds good, let us know if you have something in mind 😄

This might be of interest: github.com/atabeerk/cancer_phasing

Hey @Luuk Harbers, just a gentle nudge on lrsomatic PR #204 (haplotype DMR). I think I've addressed everything from your Sep 28 review: posted a reply on the modkit_args question (point 4), re-pointed modkit/dmr at the now-merged nf-core/modules#13021, added usage/output docs + CHANGELOG, and CI's fully green. No rush if you're busy, just flagging it's ready whenever you get a chance to take another look.

Anyone else having issues with tumor purity and ploidy assessments for ONT long read tools, I can't get the different tools to agree with each other.