Coral - #205
Open
robert-a-forsyth wants to merge 21 commits into
Open
Coral#205robert-a-forsyth wants to merge 21 commits into
robert-a-forsyth wants to merge 21 commits into
Conversation
…ysis in Deepvariant by default (toggleable)
# Conflicts: # CHANGELOG.md # subworkflows/local/phasing_haplotyping.nf # subworkflows/local/small_variant_consensus.nf # tests/clair_only.nf.test.snap # tests/consensus.nf.test.snap # tests/default.nf.test.snap # tests/union.nf.test.snap
…l joins BCFTOOLS_NORM ran with -Oz only, so multi-allelic records were never split. BCFTOOLS_ISEC matches on exact CHROM/POS/REF/ALT, so a site one caller reports as A>G,GT and the other as A>G never intersected: in consensus mode the variant was dropped, in union mode it appeared twice. ClairS-TO emits no multi-allelic records while DeepVariant emits ~2.5% (128,115 of 5,058,527 on B1975944), so every DeepVariant multi-allelic site was systematically excluded from the consensus. Adding -m -any splits them first; verified on real DeepSomatic output (45,638 -> 45,677 records over chr1:1-5Mb, 39 split, none left multi-allelic). This also makes the AF declaration mismatch harmless: ClairS-TO declares FORMAT/AF as Number=1 while DeepVariant and DeepSomatic declare Number=A, and bcftools concat only warns before keeping the first file's definition. Once every record is biallelic the two declarations are equivalent. The comment on STANDARDIZE_AF is corrected accordingly: it claimed DeepVariant emits VAF, but all four callers emit FORMAT/AF and none emits VAF, so the rename is a no-op under the default prioritize_caller='clair' and is kept as a guarantee for WAKHAN rather than as a live conversion. Two silent-failure paths in SMALL_VARIANT_CONSENSUS are closed. The caller branch had no `other` arm, so an unrecognised meta.caller removed the sample from every downstream result with a successful exit; it now errors. The DeepVariant/Clair join was a plain join, so a sample present for one caller but not the other was dropped just as quietly; it now uses failOnMismatch and failOnDuplicate, matching PHASING_HAPLOTYPING. Record counts in the merged output change, so the snapshots need regenerating. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
With --germline_var_keep deepvariant the tumour-only germline arm runs DeepVariant
on the tumour BAM. DeepVariant is a germline caller with no somatic
discrimination -- its FILTER vocabulary is only PASS/RefCall/LowQual/NoCall -- so
those calls mix germline and clonal somatic variants and cannot be separated on
their own. Published unchanged, germline_smallvariants.vcf.gz and vep/germline/
carried most of the somatic call set.
DeepSomatic evaluates the same sites and does emit a verdict: FILTER=GERMLINE
("Non somatic variants"), PON, RefCall or PASS. That verdict is now transferred
onto the DeepVariant records as INFO/DS_VERDICT and used to select the germline
arm, mirroring how ClairS-TO adjudicates its own calls through NonSomatic and
VCFSPLIT. This is a caller's own adjudication rather than positional subtraction
of the somatic call set: on B1975944 only 1.14% (57,684) of DeepVariant's
5,058,527 PASS calls carry a positive somatic verdict, against ~87% positional
overlap with the union somatic set.
Only positively adjudicated germline sites are kept (GERMLINE or PON). RefCall and
sites DeepSomatic never evaluated are dropped rather than assumed germline, so the
arm is not overpopulated with unadjudicated calls; it keeps 83.1% of the input.
The verdict stays in INFO/DS_VERDICT so the decision is auditable in the published
VCF.
The transfer is three independent bcftools invocations, so it is three aliased
instances of the existing upstream modules (DS_VERDICT_QUERY, DS_VERDICT_ANNOTATE,
DS_GERMLINE_SELECT) rather than a new bespoke process. DeepSomatic FILTER is
single-valued in practice (RefCall/GERMLINE/PON/PASS over 13.7M records), so
transferring it as a string cannot inject the ";" that would break INFO parsing.
DEEPSOMATIC now runs before DEEPVARIANT in the subworkflow because the verdict is
built from its raw VCF, before the PASS filter discards the GERMLINE records. The
deep family must therefore be enabled as a pair: --germline_var_keep deepvariant
without deepsomatic in --somatic_var_keep is rejected at launch rather than
silently producing an unadjudicated germline arm.
Paired mode is untouched -- both Clair3 and DeepVariant already run on the normal
BAM there, so the germline arm needs no adjudication.
Also corrects the --smallvar_filter_pass docs: VCFTAG normalises FILTER to PASS,
so `false` does not restore the pre-filter behaviour as previously claimed.
Germline record counts change, so the snapshots need regenerating.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…bites VCFTAG rewrites FILTER to PASS on both arms before phasing, which is needed so that downstream tools filtering on PASS still see the germline records ClairS-TO marks NonSomatic. But it did so destructively, and every record reaching the published phased VCFs, VEP and the report therefore read FILTER=PASS with the caller's verdict gone. That also silently disabled the signature-stage filter. SIGNATURES_BCFTOOLS_VIEW ran --apply-filters PASS on a call set whose FILTER had already been normalised to PASS, so it could never remove anything: with --smallvar_filter_pass false, RefCall/LowQual/GERMLINE/PON records reached SigProfilerMatrixGenerator that previously could not, and the matrices were computed over mostly-reference sites. The parameter was documented as restoring the previous behaviour, which it did not. VCFTAG now records the caller's FILTER in INFO/ORIG_FILTER before overwriting it, as VCFSPLIT already did for the ClairS-TO arm, and the signature filter tests that field instead. VCFSPLIT's stamp is respected rather than duplicated: both the header line and the per-record field are added only when absent, since a duplicate INFO key makes the record unparseable. Multi-valued FILTER is joined with "," because ";" separates INFO fields, and a FILTER of "." is recorded as "." rather than skipped. bcftools accepts only one of -i/-e, so the ALT="*" exclusion is folded into the same include expression. Verified against the module fixture with bcftools 1.20, including a record that already carried ORIG_FILTER and one with a multi-valued FILTER: exactly one ORIG_FILTER header line, no double stamping, and htslib parses the result. The module test is extended to cover these cases. Note that the module test cannot currently run locally: nf-test 0.9.3 generates `VCFTAG(*input)` and the available Nextflow 26.04.1 rejects the spread operator under its strict syntax, while the pipeline requires >=25.10.4 so 25.04.6 will not run either. This is pre-existing and reproduces unchanged at 5b38dcf. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The test read its output through path(...).vcf, an accessor provided by an nft-vcf plugin that nf-test.config does not load, so it failed with MissingPropertyException on every run and had never passed. It now uses the built-in linesGzip accessor instead, which needs no plugin. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
STANDARDIZE_AF and BCFTOOLS_ANNOTATE were two back-to-back bcftools annotate calls on the same file with disjoint options -- one renaming the allele frequency FORMAT key, the next stamping INFO/CALLER. --rename-annots composes with -a/-c/-h, verified against a real VCF, so a single invocation does both and the alias is removed along with its modules.config block. The rename is selected by meta.rename_to, set only when combine_method is 'all'; in 'consensus' mode every surviving record comes from one caller and needs no rename. Ordering is safe because BCFTOOLS_QUERY reads only CHROM/POS/REF/ALT and does not care whether the rename has happened. The rendered ext.args was checked for all three meta cases (VAF, AF, key absent): the escapes reach printf intact and an absent key yields no --rename-annots. The joins around the annotate step now use failOnMismatch/failOnDuplicate for the same reason as the caller join: a missing annotation table should stop the run, not quietly drop the sample. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
CLAIR3 received --sample_name twice: the module passes it from ext.prefix, which defaults to meta.id, and conf/modules.config passed the identical value again. Harmless with argparse, but it reads as if the name were configurable. Removed the config copy, matching how CLAIRSTO is already handled on dev. The comment in PHASING_HAPLOTYPING claiming the somatic VCF is passed first "(higher priority in phasing)" was false: BCFTOOLS_CONCAT sorts its input file list alphabetically, so with the _germline_tagged/_somatic_tagged prefixes the germline file is always passed first. The ordering is cosmetic under -a, which emits in coordinate order; the comment now says so rather than describing an intent the module discards. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Picks up TAG_GERMLINE and TAG_SOMATIC, and restores the clairsto and lrsomatic_report versions that the previous snapshot run had rolled back. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
# Conflicts: # tests/clair_only.nf.test.snap # tests/consensus.nf.test.snap # tests/deep_only.nf.test.snap # tests/default.nf.test.snap # tests/union.nf.test.snap
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
# Conflicts: # CHANGELOG.md
Nextflow's strict syntax resolves `report_empty_slot()` against declared functions, not local closure variables, so it read as undefined. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Folds the standalone Downstream pipeline's amplicon analysis into lrsomatic as an ECDNA subworkflow, on by default for tumour samples with ASCAT calls. CoRAL seeds from ASCAT's cnvs.txt, which the nf-core ASCAT module already emits and nothing consumed. bin/ascat_to_coral_bed.py converts it to the headerless BED CoRAL wants, respelling contigs to match the reference: ASCAT writes "1" where the BAM may say "chr1", and CoRAL builds chromosome sizes from the BAM header. Zero-length and inverted segments are dropped, which CoRAL's own parser would otherwise reject. The two tools were previously local .sif files built from uncommitted working trees, so the images could not be rebuilt by anyone. They are now built from pinned commits of public forks: - docker.io/robertaforsyth/coral:3.0.0-chm13-847f3d4 - docker.io/robertaforsyth/ampliconclassifier:2.0.0-chm13-cdeaa63 with Dockerfiles under containers/. AmpliconClassifier moves from the v1.5.2 fork to upstream v2.0.0, the first release with official CoRAL support; the CHM13 patches applied unchanged. On B2096749 the ecDNA call is reproduced exactly, and two amplicons v1.5.2 reported as "No amp/Invalid" now resolve as Complex-non-cyclic and Linear. CoRAL defaults to the open-source SCIP solver rather than Gurobi. Across the 1102 solver logs from the previous cohort the largest model was 620 rows by 438 columns and the slowest solve took 0.49 s, so a licensed solver buys nothing at this scale. --coral_solver gurobi_direct with --gurobi_license remains available, and serialises reconstruction because WLS licences cap concurrent sessions. Also fixes three bugs carried over from the Downstream modules: publishDir paths were doubled by mkdir'ing a directory named after the publish target, the empty-seed filter staged files just to measure them, and errorStrategy 'ignore' hid a real reconstruction failure without surfacing it. Not included: an ecDNA section in LRSOMATICREPORT. That module renders through render_report.R inside ghcr.io/ljwharbers/lrsomatic-report, so displaying ecDNA results needs a change in ljwharbers/lrsomatic_report first. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
ljwharbers
requested changes
Sep 28, 2026
ljwharbers
left a comment
ljwharbers
left a comment
Collaborator
There was a problem hiding this comment.
Thanks for folding the ecDNA analysis in. The port is careful: the conversion script, the pinned image builds and the join hardening all look good. There is one blocker and two defaults that need changing before this can merge.
Scope note: this branch is stacked on #200 (18 of its 21 commits are variant_filter's), and #200 has moved on since (now at ff88896). So I reviewed only the ecDNA commits (1376557, 001926c). Once #200 merges, please merge dev into coral so the diff here is ecDNA-only.
Blocker: AmpliconClassifier cannot find the downloaded GRCh38 data repo
GRCh38.tar.gzhas a single top-level directory,GRCh38/. I streamed the full listing: 25 entries, includingGRCh38/file_list.txt.- nf-core
UNTARapplies--strip-components 1whenever there is one top-level dir, soUNTAR_AA_DATA_REPOemitsaa_data_repo/withfile_list.txtdirectly inside it. - AC (
amplicon_classifier.py:2565-2568at cdeaa63) opens$AA_DATA_REPO + "/" + args.ref + "/file_list.txt", i.e.aa_data_repo/GRCh38/file_list.txt. That path does not exist. - The resulting
FileNotFoundErroris not caught (onlyKeyErroris). AMPLICONCLASSIFIER keeps the basefinishstrategy, so every default GRCh38 run with at least one seeded sample fails. - The failure is also invisible: stdout and stderr are redirected into
${prefix}_classifier.log, so.command.erris empty. - The docs (a
CHM13/subdirectory under--aa_data_repo) describe the parent-dir layout, so I suspect only the local-repo path was exercised. No test covers the download path, because the pipeline tests run withskip_coral = trueand the module and subworkflow tests are stubs.
Default CHM13 runs now fail at launch
With the feature on by default and aa_data_repo_url: null for CHM13, PREPARE_AA_DATA_REPO calls error(). Every existing --genome CHM13 command line breaks until the user adds --aa_data_repo or --skip_ampliconclassifier, and there is no CHM13 repo to point at. I'd skip the classifier with a log.warn on CHM13 when no repo is given, so reconstruction still runs.
CoRAL failures are silently dropped
CORAL_RECONSTRUCT and CORAL_CYCLE use errorStrategy ... : 'ignore'. The comment says the subworkflow warns on the missing output, but ecdna.nf has no such warning. The commit message lists this exact problem as fixed. As it stands, any real CoRAL error (a bad --ref, a solver crash) removes the sample from the results while the run finishes green. Please either drop ignore, or join branched_seeds.seeded against CORAL_RECONSTRUCT.out.reconstruction with remainder: true and log.warn the misses.
Smaller points
- CHANGELOG:
#XXXshould be#205. The test entry claims coverage of "the empty-seed branch", but none of the threeECDNAtests reaches it, because theCORAL_SEEDstub always writes a non-empty seed. - Data repo is not pinned: the tarball is updated in place (it has a
GRCh38/last_updated.txt) and nothing checks it. An MD5 as for--vep_clinvar_md5(#206) would keep runs reproducible. Please also document that--aa_data_repoavoids re-fetching 1.1 GB (and unpacking ~4 GB) on every run. - Gurobi
containerOptions:--bind … --env …is Apptainer syntax, and Docker rejects--bind. - Resources:
CORAL_RECONSTRUCTandCORAL_CYCLEareprocess_high(12 CPUs / 72 GB) for every tumour sample, while your own numbers show sub-second solves.--solver-threads -1is also not tied totask.cpus, which matters under Gurobi. - Plots vs classification: with
--coral_run_cycle,CORAL_PLOTstill plots the original reconstruction, not the re-extracted cycles that AC classifies. - Image and fork namespaces: the images and forks live in personal namespaces (
docker.io/robertaforsyth,github.com/robert-a-forsyth). For long-term maintenance, a lab-owned namespace would be safer, as was done for ClairS-TO. - Nit:
bed.toFile().length()→bed.size(), which also works on non-local work dirs. - pre-commit: it is red only on prettier (
docs/output.md,nextflow_schema.json);pre-commit run --all-filesfixes it.
Checked and fine
- ASCAT's
cnvs.txtcolumns match the converter (segments[2:6]). coral seedalways writes the seed file, possibly empty, so the size branch is sound.- Both Docker Hub tags exist.
- The AA URL is live, and the wget image fetches it over HTTPS.
- The meta keys from
ascat_chline up with ASCAT's output, so the joins pair.
|
|
||
| // .first(): UNTAR emits a queue channel, and every sample's classifier task | ||
| // needs the same repo -- without this only the first sample would get it. | ||
| ch_data_repo = UNTAR_AA_DATA_REPO.out.untar.first() |
Collaborator
There was a problem hiding this comment.
GRCh38.tar.gz has a single top-level GRCh38/, so UNTAR strips it and this dir holds file_list.txt directly. AC reads $AA_DATA_REPO/GRCh38/file_list.txt, so every GRCh38 classifier task fails. One fix is to stage it under the ref name (e.g. set ext.prefix so the untarred dir is GRCh38, and stage it in AMPLICONCLASSIFIER as aa_data_repo/${ac_ref}). Another is a small shim in the module: if aa_data_repo/file_list.txt exists, mkdir repo && ln -s ../aa_data_repo repo/${ac_ref} and export AA_DATA_REPO=$PWD/repo. The shim keeps the documented parent-dir layout working for --aa_data_repo too.
| ch_versions = ch_versions.mix(WGET_AA_DATA_REPO.out.versions, UNTAR_AA_DATA_REPO.out.versions) | ||
| } | ||
| else { | ||
| error("AmpliconClassifier needs an AmpliconArchitect data repository, which is not published for ${params.genome}. Set --aa_data_repo <path>, or use --skip_ampliconclassifier.") |
Collaborator
There was a problem hiding this comment.
The feature is on by default and CHM13 has no URL, so this now aborts every --genome CHM13 run that doesn't add --aa_data_repo or --skip_ampliconclassifier. Could this warn and skip the classifier instead, leaving reconstruction running?
|
|
||
| // A single unsolvable amplicon should not fail a whole cohort; the subworkflow | ||
| // warns on the missing output rather than letting the report slot go silent. | ||
| errorStrategy { task.exitStatus in 130..145 ? 'retry' : 'ignore' } |
Collaborator
There was a problem hiding this comment.
The comment above says the subworkflow warns on the missing output, but ecdna.nf has no such warning, so a failed reconstruction just disappears (same in CORAL_CYCLE). Either drop ignore, or add the warning: join the seeded samples against CORAL_RECONSTRUCT.out.reconstruction with remainder: true and log.warn where the reconstruction is null.
| --AA_results ${reconstruction} \\ | ||
| -o ${prefix} \\ | ||
| ${args} \\ | ||
| > ${prefix}_classifier.log 2>&1 |
Collaborator
There was a problem hiding this comment.
This hides AC's traceback from .command.err, so a failure shows as a bare non-zero exit. 2>&1 | tee ${prefix}_classifier.log keeps the published log and still surfaces the error (Nextflow runs with pipefail).
| // Gurobi licences are user-supplied and mounted; SCIP, the default, needs nothing. | ||
| containerOptions = { | ||
| params.coral_solver == 'gurobi_direct' && params.gurobi_license | ||
| ? "--bind ${params.gurobi_license}:/opt/gurobi/gurobi.lic:ro --env GRB_LICENSE_FILE=/opt/gurobi/gurobi.lic" |
Collaborator
There was a problem hiding this comment.
--bind is Apptainer/Singularity-only; under -profile docker this needs -v ${params.gurobi_license}:/opt/gurobi/gurobi.lic:ro. Branch on workflow.containerEngine.
|
|
||
| - [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - Added ecDNA and focal amplification analysis: CoRAL reconstructs amplicon structures from each tumour BAM, seeded by ASCAT's copy-number calls, and AmpliconClassifier labels each amplicon as ecDNA, BFB or linear. Runs by default on tumour samples with ASCAT calls; turn it off with `--skip_coral`, or keep reconstruction without classification with `--skip_ampliconclassifier`. CoRAL's optimisation defaults to the open-source SCIP solver rather than Gurobi: across the 1102 solver logs from the previous standalone cohort the largest model was 620 rows by 438 columns and the slowest solve took 0.49 s, so a licensed solver buys nothing at this scale, though `--coral_solver gurobi_direct` with `--gurobi_license` remains available. AmpliconClassifier's AmpliconArchitect data repository is downloaded automatically for GRCh38 and must be supplied with `--aa_data_repo` for CHM13, which has no published repository (@robert-a-forsyth). | ||
| - [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - Added `ascat_to_coral_bed.py` and the `ASCAT_TO_CORAL_BED` module, converting ASCAT's `cnvs.txt` into the headerless BED CoRAL seeds from. Contigs are respelled to match the reference, since ASCAT writes `1` where the BAM may say `chr1` and CoRAL builds its chromosome sizes from the BAM header; zero-length and inverted segments are dropped, which CoRAL's own parser would otherwise reject (@robert-a-forsyth). | ||
| - [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - Added stub nf-tests for the `ECDNA` subworkflow and `ASCAT_TO_CORAL_BED` (tag `small`), covering the empty-seed branch, the opt-in cycle re-extraction and the `--skip_ampliconclassifier` path (@robert-a-forsyth). |
Collaborator
There was a problem hiding this comment.
#XXX → #205 (also on the two lines above). The claim about covering the empty-seed branch doesn't hold: the CORAL_SEED stub always writes a non-empty seed, and none of the three ECDNA tests forces an empty one.
8 tasks
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
PR checklist
nf-core pipelines lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).