diff --git a/CHANGELOG.md b/CHANGELOG.md index fb2c964c..85e82f50 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Added` +- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - 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, checked against a pinned `--aa_data_repo_md5`, and can be supplied with `--aa_data_repo` (the reference directory or its parent); CHM13 has no published repository, so without `--aa_data_repo` the classifier is skipped with a warning (@robert-a-forsyth). +- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - 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). +- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Added stub nf-tests for the `ECDNA` subworkflow and `ASCAT_TO_CORAL_BED` (tag `small`), covering the seeded path, the empty-seed branch, the opt-in cycle re-extraction with plotting, the `--skip_ampliconclassifier` path and a run with no data repository, plus a stub test for `AADATAREPO_DOWNLOAD` (@robert-a-forsyth). - [#197](https://github.com/IntGenomicsLab/lrsomatic/pull/197) - Added CHM13 support for ClairS-TO's Verdict module, which tags tumour-only calls as germline, somatic or subclonal somatic; its resources were GRCh38-only, so on CHM13 germline variants leaked into `somatic.vcf.gz`. With `--genome CHM13 --skip_ascat` the pipeline builds a CHM13 resource set from the ASCAT files it already downloads and passes it as `--cna_resource_dir`; a prepared directory can be given with `--clairsto_cna_resources` (validated at launch). Without `--skip_ascat` tagging comes from ASCAT's own tables instead (next entry) (@ljwharbers). - [#197](https://github.com/IntGenomicsLab/lrsomatic/pull/197) - Added `CLAIRSTO_VERDICT_TAG`: when ASCAT is in the run, Verdict's germline tagging is computed from ASCAT's purity, ploidy and segments instead of Verdict's own estimate, so `CLAIRSTO` runs with `--disable_verdict` and ASCAT runs before small variant calling. Output names are unchanged. The tables the tags were computed from are published as `_Tumor_Purity_Ploidy.txt` and `_Tumor_CNA.txt`, also on `--skip_ascat` runs (@ljwharbers). - [#197](https://github.com/IntGenomicsLab/lrsomatic/pull/197) - Added a stub nf-test for `TUMORONLY_SMALLVAR` covering both germline tagging paths (tag `small`) (@ljwharbers). @@ -20,6 +23,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Changed` +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - In tumour-only mode with both `deepvariant` and `deepsomatic` selected, DeepVariant's germline calls are now kept only where DeepSomatic's verdict is `GERMLINE` or `PON` (`INFO/DS_VERDICT`); with `deepvariant` alone they are unfiltered and may include somatic variants (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - The phased somatic VCF is now split from the merged germline+somatic VCF by an `INFO/SOMATIC` provenance flag (new `VCFTAG` module) instead of by position, so germline records at somatic coordinates no longer leak into it; ClairS-TO germline records keep their original `FILTER` in `INFO/ORIG_FILTER` (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - Added `--smallvar_filter_pass` (default `true`), which restricts each small variant caller's output to `PASS` records before the caller consensus, phasing, VEP and the report; the published per-caller VCFs are unchanged (@robert-a-forsyth). - [#201](https://github.com/IntGenomicsLab/lrsomatic/pull/201) - `CLAIRSTO` and `CLAIRSTO_VERDICT_TAG` now pull the fork image from Docker Hub: `oras://docker.io/ljwharbers/clairs-to-sif:0.5.1-verdict-chm13-c0687e8-flat` under Singularity/Apptainer and `docker.io/ljwharbers/clairs-to:0.5.1-verdict-chm13-c0687e8-flat` otherwise, instead of `ghcr.io/ljwharbers/clairs-to`. The `-cpu` SIF on ghcr failed with `PROTOCOL_ERROR` on slow links: ghcr redirects every blob download to an Azure URL that expires at the next 5-minute mark and resets a stream still open then, and Apptainer resumes neither an `oras://` nor a `docker://` download. Docker Hub's download URLs are valid for 50 minutes and only checked when the request starts. `-flat` is the same software copied into an empty image in a few layers (3.3 GB instead of 7 GB); the software and its outputs are unchanged. `docs/usage.md` describes `pullTimeout`, Docker Hub's anonymous pull limit, pre-pulling, and how to recover the remaining `oras://ghcr.io` SIFs resumably (@ljwharbers). - [#199](https://github.com/IntGenomicsLab/lrsomatic/pull/199) - `CLAIRSTO` and `CLAIRSTO_VERDICT_TAG` now run the `-cpu` rebuild of the fork image (`0.5.1-verdict-chm13-c0687e8-cpu`), which swaps PyTorch's CUDA build for the CPU build of the same version. The software is otherwise unchanged, but the Apptainer SIF drops from 6.53 GB to 3.46 GB. The old image could not be pulled on a normal VSC link: Apptainer fetches an `oras://` SIF as a single unresumable stream, and the signed blob URL ghcr redirects to expires on a 15-minute wall-clock boundary, so 6.53 GB needed 7.3 MB/s sustained and was otherwise cut mid-transfer with `PROTOCOL_ERROR` (@ljwharbers). - [#197](https://github.com/IntGenomicsLab/lrsomatic/pull/197) - `CLAIRSTO` now runs `ghcr.io/ljwharbers/clairs-to:0.5.1-verdict-chm13-c0687e8` (ClairS-TO 0.5.1) instead of `docker.io/hkubal/clairs-to:v0.4.2`: a fork that lets Verdict read its CNA resources from `--cna_resource_dir`, fixes four places where Verdict's Python port of ASCAT departed from R, and disables Verdict with a warning when its resources cannot be read. **GRCh38 results move as well as CHM13 ones.** Revert to the upstream image once HKU-BAL/ClairS-TO carries these changes. The module also selects the SIF under `-profile apptainer` and sets explicit output prefixes (@ljwharbers). @@ -44,6 +50,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - [#206](https://github.com/IntGenomicsLab/lrsomatic/pull/206) - A remote (http, https or ftp) ClinVar is now downloaded once per run by the new `VEPPLUGIN_CLINVAR` step in `PREPARE_VEP_PLUGINS`, instead of being staged by Nextflow as a foreign file; local and cloud-storage paths are staged as before. `GERMLINE_VEP` and `SOMATIC_VEP` re-checked the foreign file on its host for every sample, and NCBI answered the burst from a multi-sample GRCh38 run with HTTP 503, so `SOMATIC_VEP` failed with `Can't stage file ...clinvar_20260829.vcf.gz`; `-resume` could not recover, since the failed check changed the staging cache key. The download is checked against the new `--vep_clinvar_md5` and `--vep_clinvar_tbi_md5`, set by default to the checksums NCBI (GRCh38, VCF only) and Ensembl (CHM13, VCF and index) publish, so the pinned release cannot change silently. Resuming a run that already finished re-runs `GERMLINE_VEP` and `SOMATIC_VEP` once, since ClinVar now comes from a task rather than the stage cache. The ClinVar sizes in `docs/usage.md` are also corrected, and `docs/output.md` now documents `vep_plugins/` (@AmberVerhasselt). - [#203](https://github.com/IntGenomicsLab/lrsomatic/pull/203) - `CLAIRS` no longer runs with `--haplotagged_tumor_bam_provided_so_skip_intermediate_phasing_and_haplotagging`. Since somatic calling was moved ahead of `PHASING_HAPLOTYPING` (v1.1.0), ClairS has received the untagged minimap2 BAM, so the flag made it skip its own phasing and haplotagging and call every paired sample without haplotype information: the full-alignment model saw no `HP` tags and the haplotype filtering step had nothing to filter on, the same as `--disable_phasing`. ClairS now runs Clair3 on the normal and tumour BAMs and phases the tumour itself. **Paired somatic calls change** (fewer false positives expected), and `CLAIRS` takes longer and uses more work-directory space (@ljwharbers). - [#203](https://github.com/IntGenomicsLab/lrsomatic/pull/203) - `docs/output.md` now lists the ClairS SNV output as `snvs.vcf.gz`, the name the pipeline publishes, instead of `snv.vcf.gz` (@ljwharbers). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - `*_var_combine = 'all'` now keeps both callers' private calls; the prioritized caller's own private calls were previously dropped. Invalid `combine_method`/`prioritize_caller` values now raise an error (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - `*_var_combine = 'consensus'` now splits multi-allelic records before intersecting, so they match across callers, and rejoins them before phasing. `'all'` now keeps the prioritised caller's record wherever both callers call a position, so no record mixes two callers' alleles (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - Somatic phasing no longer includes germline records at the position of an alt-genotype somatic call, which gave the somatic record the germline genotype and phase set (new `BCFTOOLS_EXCLUDE_SITES` module); `0/0` and `./.` somatic records are published unphased (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - `LRSOMATICREPORT` now renders each sample as soon as its own inputs are ready instead of waiting for the whole batch, including samples without ASCAT's optional raw segments output (@robert-a-forsyth). +- [#200](https://github.com/IntGenomicsLab/lrsomatic/pull/200) - Bumped `WAKHAN` from 0.4.3 to 0.4.4, fixing a deterministic `StatisticsError` crash in BAF binning that aborted the run (@robert-a-forsyth). - [#196](https://github.com/IntGenomicsLab/lrsomatic/pull/196) - `LRSOMATICREPORT` now points `XDG_CACHE_HOME` at the task directory alongside `HOME` and `TMPDIR`. Singularity/Apptainer inherit the host environment, so on sites that set it outside the bind-mounted work tree the render died with `Read-only file system (os error 30): mkdir '<...>/.cache/quarto'` (@AmberVerhasselt, @ljwharbers). - [#193](https://github.com/IntGenomicsLab/lrsomatic/pull/193) - `--vep_eve https://evemodel.org/api/proteins/bulk/download/` was rejected at launch because the "needs preparing" check keyed on a `.zip` suffix; it now checks whether the value is already a prepared bgzipped file (@AmberVerhasselt). - [#188](https://github.com/IntGenomicsLab/lrsomatic/pull/188) - `MODKIT_PILEUP` now runs a patched modkit 0.6.4 ([ljwharbers/modkit@pacbio-conflict-fix](https://github.com/ljwharbers/modkit/tree/pacbio-conflict-fix)): `ghcr.io/ljwharbers/modkit:0.6.4-pacbiofix-6e0afa2` under Docker and `oras://ghcr.io/ljwharbers/modkit-sif:0.6.4-pacbiofix-6e0afa2` under Singularity/Apptainer. Stock modkit 0.4.3-0.6.4 dropped 32-65 % of reads from recent PacBio HiFi BAMs and returned empty `--cpg` pileups ([nanoporetech/modkit#612](https://github.com/nanoporetech/modkit/issues/612), fix proposed in [nanoporetech/modkit#720](https://github.com/nanoporetech/modkit/pull/720)), and ignored `--phased`/`--modified-bases` for PacBio BAMs with 6mA calls. The image is `linux/amd64` only and Conda is not supported (use `--skip_modkit` there); return to the biocontainer once a release includes the fix (@ljwharbers). diff --git a/CITATIONS.md b/CITATIONS.md index a1403be7..cbb41940 100644 --- a/CITATIONS.md +++ b/CITATIONS.md @@ -18,6 +18,10 @@ > Cheng J, Novati G, Pan J, Bycroft C, Žemgulytė A, Applebaum T, Pritzel A, Wong LH, Zielinski M, Sargeant T, Schneider RG, Senior AW, Jumper J, Hassabis D, Kohli P, Avsec Ž. Accurate proteome-wide missense variant effect prediction with AlphaMissense. Science. 2023 Sep 22;381(6664):eadg7492. doi: 10.1126/science.adg7492. +- [AmpliconClassifier](https://doi.org/10.1038/s41586-023-05937-5) + + > Luebeck, J., Ng, A.W.T., Galipeau, P.C. et al. Extrachromosomal DNA in the cancerous transformation of Barrett's oesophagus. Nature 616, 798-805 (2023). https://doi.org/10.1038/s41586-023-05937-5 + - [ASCAT](https://pubmed.ncbi.nlm.nih.gov/20837533/) > Van Loo P, Nordgard SH, Lingjærde OC, Russnes HG, Rye IH, Sun W, Weigman VJ, Marynen P, Zetterberg A, Naume B, Perou CM, Børresen-Dale AL, Kristensen VN. Allele-specific copy number analysis of tumors. Proc Natl Acad Sci U S A. 2010 Sep 28;107(39):16910-5. doi: 10.1073/pnas.1009843107. Epub 2010 Sep 13. PubMed PMID: 20837533; PubMed Central PMCID: PMC2947907. @@ -46,6 +50,10 @@ > Landrum MJ, Lee JM, Benson M, Brown GR, Chao C, Chitipiralla S, Gu B, Hart J, Hoffman D, Jang W, Karapetyan K, Katz K, Liu C, Maddipatla Z, Malheiro A, McDaniel K, Ovetsky M, Riley G, Zhou G, Holmes JB, Kattman BL, Maglott DR. ClinVar: improving access to variant interpretations and supporting evidence. Nucleic Acids Res. 2018 Jan 4;46(D1):D1062-D1067. doi: 10.1093/nar/gkx1153. +- [CoRAL](https://doi.org/10.1101/gr.279131.124) + + > Zhu, K., Jones, M.G., Luebeck, J. et al. CoRAL accurately resolves extrachromosomal DNA genome structures with long-read sequencing. Genome Res. 34, 1344-1354 (2024). https://doi.org/10.1101/gr.279131.124 + - [cramino](https://github.com/wdecoster/cramino) > De Coster W. cramino: A fast and simple tool for quality control of long read sequencing data [Software]. GitHub. https://github.com/wdecoster/cramino diff --git a/bin/ascat_to_coral_bed.py b/bin/ascat_to_coral_bed.py new file mode 100755 index 00000000..2e7e7d21 --- /dev/null +++ b/bin/ascat_to_coral_bed.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +"""Convert ASCAT's cnvs.txt into the headerless BED CoRAL's --cn-seg expects. + +CoRAL reads total copy number from the last column and builds chromosome sizes from +the BAM header, so contigs are respelled to match the reference (chr1 vs 1). +""" +import argparse +import sys + +REQUIRED = ('chr', 'startpos', 'endpos', 'nMajor', 'nMinor') + + +def fai_contigs(path): + if not path: + return [] + with open(path) as fp: + return [line.split('\t', 1)[0] for line in fp if line.strip()] + + +def spell_like_reference(chrom, contigs): + """ASCAT's '1' becomes 'chr1' when that is what the reference calls it, and the reverse.""" + if not contigs or chrom in contigs: + return chrom + if chrom.startswith('chr') and chrom[3:] in contigs: + return chrom[3:] + if 'chr' + chrom in contigs: + return 'chr' + chrom + return chrom + + +def sort_key(chrom): + """Natural contig order: 1-22, then X, Y, then anything else alphabetically.""" + bare = chrom[3:] if chrom.startswith('chr') else chrom + if bare.isdigit(): + return (0, int(bare), '') + if bare in ('X', 'Y'): + return (1, 'XY'.index(bare), '') + return (2, 0, bare) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--cnvs', required=True, help="ASCAT .cnvs.txt") + parser.add_argument('--fai', help="Reference .fai, to respell contigs to match the BAM") + parser.add_argument('--output', required=True, help="BED4 written for CoRAL --cn-seg") + args = parser.parse_args() + + contigs = fai_contigs(args.fai) + + with open(args.cnvs) as fp: + header = fp.readline().rstrip('\n').split('\t') + missing = [c for c in REQUIRED if c not in header] + if missing: + sys.exit(f"ERROR: {args.cnvs} is missing required columns: {', '.join(missing)}") + idx = {c: header.index(c) for c in REQUIRED} + + rows, dropped = [], 0 + for line in fp: + if not line.strip(): + continue + fields = line.rstrip('\n').split('\t') + start, end = int(fields[idx['startpos']]), int(fields[idx['endpos']]) + # CoRAL's segment parser rejects these, so drop them here with a count + if start >= end: + dropped += 1 + continue + total_cn = round(float(fields[idx['nMajor']])) + round(float(fields[idx['nMinor']])) + rows.append((spell_like_reference(fields[idx['chr']], contigs), start, end, total_cn)) + + if dropped: + print(f"WARNING: dropped {dropped} segments with start >= end", file=sys.stderr) + + rows.sort(key=lambda r: (sort_key(r[0]), r[1])) + with open(args.output, 'w') as out: + for chrom, start, end, total_cn in rows: + out.write(f"{chrom}\t{start}\t{end}\t{total_cn}\n") + + if not rows: + print(f"WARNING: {args.output} is empty; CoRAL will find no seeds", file=sys.stderr) + + +if __name__ == '__main__': + main() diff --git a/conf/igenomes.config b/conf/igenomes.config index 430e65a2..32ab6ddc 100644 --- a/conf/igenomes.config +++ b/conf/igenomes.config @@ -24,6 +24,13 @@ params.genomes = [ vep_species : "homo_sapiens", savana_contigs : "https://raw.githubusercontent.com/cortes-ciriano-lab/savana/main/example/contigs.chr.hg38.txt", savana_g1000_vcf : "1000g_hg38", + coral_ref : "hg38", + ac_ref : "GRCh38", + // Plain build, not GRCh38_indexed: the extra BWA index is for AmpliconArchitect's + // alignment step, which AmpliconClassifier never runs. + aa_data_repo_url : "https://refs.ampliconrepository.org/data/module_support_files/AmpliconArchitect/GRCh38.tar.gz", + // Published by the host as GRCh38_md5sum.txt; the tarball is updated in place + aa_data_repo_md5 : "2bcd1fdaed027466ed296bf78b2907eb", vep_alphamissense : "https://storage.googleapis.com/dm_alphamissense/AlphaMissense_hg38.tsv.gz", vep_alphamissense_tbi : "https://g-608c0c.273595.03c0.data.globus.org/VEP_plugins/AlphaMissense_hg38.tsv.gz.tbi", // A dated release rather than the rolling vcf_GRCh38/clinvar.vcf.gz, whose VCF and @@ -57,6 +64,10 @@ params.genomes = [ vep_species : "homo_sapiens_gca009914755v4", savana_contigs : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/savana/contigs.chr.chm13.txt", savana_g1000_vcf : "1000g_t2t", + coral_ref : "t2t", + ac_ref : "CHM13", + // No published AA data repo for CHM13; supply one with --aa_data_repo + aa_data_repo_url : null, vep_alphamissense_aa : "https://g-608c0c.273595.03c0.data.globus.org/VEP_plugins/alphamissense_protein_v2023_uniprot-2026_03.tsv.gz", vep_alphamissense_aa_tbi : "https://g-608c0c.273595.03c0.data.globus.org/VEP_plugins/alphamissense_protein_v2023_uniprot-2026_03.tsv.gz.tbi", // Pinned to a release rather than current_variation/, which moves at every Ensembl release diff --git a/conf/modules.config b/conf/modules.config index 54494857..198687f0 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -121,14 +121,28 @@ process { withName: '.*:BCFTOOLS_NORM' { ext.prefix = { "${meta.id}.${meta.caller}_norm" } - ext.args = { - "-Oz" - } + // 'consensus' splits multi-allelics so BCFTOOLS_ISEC can match them per ALT; BCFTOOLS_NORM_REJOIN rejoins them. + ext.args = { meta.split ? "-m -any -Oz" : "-Oz" } publishDir = [ enabled: false ] } + // Rejoin split multi-allelics ('consensus' only); LongPhase keys variants by position and cannot take two records at one POS. + withName: '.*:GERMLINE_CONSENSUS:BCFTOOLS_NORM_REJOIN' { + ext.prefix = { "${meta.id}_germline_sorted" } + ext.args = { '-m +any --output-type z --write-index=tbi' } + publishDir = [ enabled: false ] + } + withName: '.*:SOMATIC_CONSENSUS:BCFTOOLS_NORM_REJOIN' { + ext.prefix = { "${meta.id}_somatic_sorted" } + ext.args = { '-m +any --output-type z --write-index=tbi' } + publishDir = [ enabled: false ] + } + withName: '.*_CONSENSUS:BCFTOOLS_EXCLUDE_SITES' { + ext.prefix = { "${meta.id}_nonpriority_only" } + publishDir = [ enabled: false ] + } withName: '.*:BCFTOOLS_ISEC' { ext.prefix = { "${meta.id}_isec" } ext.args ={ @@ -138,22 +152,41 @@ process { enabled: false ] } - withName: '.*STANDARDIZE_AF' { - ext.prefix = { "${meta.id}.${meta.caller}_standardized" } - ext.args = { - meta.rename_to == 'VAF' - ? "--rename-annots <(printf 'FORMAT/AF\\tFORMAT/VAF\\n') -Oz -W=tbi" - : "--rename-annots <(printf 'FORMAT/VAF\\tFORMAT/AF\\n') -Oz -W=tbi" + withName: '.*:BCFTOOLS_ANNOTATE' { + ext.prefix = { "${meta.id}.${meta.caller}" } + // Stamp INFO/CALLER; in 'all' mode (meta.rename_to set) also unify the FORMAT allele frequency key. + ext.args = { + def rename = meta.rename_to == 'VAF' + ? "--rename-annots <(printf 'FORMAT/AF\\tFORMAT/VAF\\n') " + : meta.rename_to == 'AF' + ? "--rename-annots <(printf 'FORMAT/VAF\\tFORMAT/AF\\n') " + : "" + rename + '''-h <(echo '##INFO=') \ + -c CHROM,POS,REF,ALT,INFO/CALLER \ + -Oz \ + -W=tbi''' } publishDir = [ enabled: false ] } - withName: '.*:BCFTOOLS_ANNOTATE' { - ext.prefix = { "${meta.id}.${meta.caller}" } + // GERMLINE VERDICT TRANSFER (tumor-only): DeepSomatic's verdict filters DeepVariant's germline calls. + withName: '.*:DS_VERDICT_QUERY' { + ext.prefix = { "${meta.id}.ds_verdict" } + // Exclude PASS/RefCall rather than include GERMLINE/PON: bcftools fails on an undeclared FILTER. ext.args = { - '''-h <(echo '##INFO=') \ - -c CHROM,POS,REF,ALT,INFO/CALLER \ + "-e 'FILTER=\"PASS\" || FILTER=\"RefCall\"' -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\n'" + } + publishDir = [ + enabled: false + ] + } + + withName: '.*:DS_VERDICT_ANNOTATE' { + ext.prefix = { "${meta.id}.deepvariant_verdict" } + ext.args = { + '''-h <(echo '##INFO=') \ + -c CHROM,POS,REF,ALT,INFO/DS_VERDICT \ -Oz \ -W=tbi''' } @@ -161,6 +194,18 @@ process { enabled: false ] } + + withName: '.*:DS_GERMLINE_SELECT' { + ext.prefix = { "${meta.id}.deepvariant_germline" } + // Keep only DeepSomatic-adjudicated germline sites; all others lack DS_VERDICT or fail this test. + ext.args = { + "-i 'INFO/DS_VERDICT=\"GERMLINE\" || INFO/DS_VERDICT=\"PON\"' --output-type z --write-index=tbi" + } + publishDir = [ + enabled: false + ] + } + withName: '.*:BCFTOOLS_QUERY' { ext.args = { "-f '%CHROM\t%POS\t%REF\t%ALT\t${meta.caller}\n'" @@ -401,8 +446,47 @@ process { enabled: false ] } + // Provenance flags stamped before the phasing merge, so the somatic arm is selected by origin. + withName: '.*:TAG_SOMATIC' { + ext.prefix = { "${meta.id}_somatic_tagged" } + publishDir = [ + enabled: false + ] + } + withName: '.*:TAG_GERMLINE' { + ext.prefix = { "${meta.id}_germline_tagged" } + publishDir = [ + enabled: false + ] + } + // Only alt-genotype somatic records are phased; GERMLINE_ANCHORS drops germline records at their positions. + withName: '.*:PHASING_HAPLOTYPING:SOMATIC_ALT' { + ext.prefix = { "${meta.id}_somatic_alt" } + ext.args = { "-i 'GT=\"alt\"'" } + publishDir = [ enabled: false ] + } + withName: '.*:PHASING_HAPLOTYPING:SOMATIC_NONALT' { + ext.prefix = { "${meta.id}_somatic_nonalt" } + ext.args = { "-e 'GT=\"alt\"'" } + publishDir = [ enabled: false ] + } + withName: '.*:PHASING_HAPLOTYPING:GERMLINE_ANCHORS' { + ext.prefix = { "${meta.id}_germline_anchors" } + publishDir = [ enabled: false ] + } withName: '.*:PHASING_HAPLOTYPING:BCFTOOLS_VIEW' { + ext.prefix = { "${meta.id}_somatic_phased_alt" } + ext.args = { "-i 'INFO/SOMATIC=1'" } + publishDir = [ enabled: false ] + } + withName: '.*:PHASING_HAPLOTYPING:CONCAT_SOMATIC_UNPHASED' { + ext.prefix = { "${meta.id}_somatic_combined" } + ext.args = { '-Oz -a -W=tbi' } + publishDir = [ enabled: false ] + } + withName: '.*:PHASING_HAPLOTYPING:SORT_SOMATIC_PHASED' { ext.prefix = { "somatic_smallvariants" } + ext.args = { '-Oz -W=tbi' } publishDir = [ path: { "${params.outdir}/${meta.id}/variants/phased" }, mode: params.publish_dir_mode, @@ -510,12 +594,12 @@ process { ] } withName: '.*:GERMLINE_CONSENSUS:BCFTOOLS_SORT_CONSENSUS' { - ext.prefix = { "${meta.id}_germline_sorted" } + ext.prefix = { "${meta.id}_germline_split_sorted" } ext.args = { '-Oz -W=tbi' } publishDir = [ enabled: false ] } withName: '.*:SOMATIC_CONSENSUS:BCFTOOLS_SORT_CONSENSUS' { - ext.prefix = { "${meta.id}_somatic_sorted" } + ext.prefix = { "${meta.id}_somatic_split_sorted" } ext.args = { '-Oz -W=tbi' } publishDir = [ enabled: false ] } @@ -576,7 +660,7 @@ process { } withName: '.*:CLAIR3' { - ext.args = { "--sample_name=${meta.id}" } + // --sample_name is passed by the module itself, from ext.prefix (which defaults to meta.id). publishDir = [ path: { "${params.outdir}/${meta.id}/variants/clair3" }, mode: params.publish_dir_mode, @@ -609,6 +693,90 @@ process { ] } + withName: '.*:ASCAT_TO_CORAL_BED' { + publishDir = [ + path: { "${params.outdir}/${meta.id}/ecdna" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + withName: '.*:CORAL_.*' { + // Gurobi licences are user-supplied and mounted; SCIP, the default, needs nothing. + // Apptainer/Singularity mount with --bind, Docker/Podman with --volume. + containerOptions = { + if (params.coral_solver != 'gurobi_direct' || !params.gurobi_license) { return null } + def mount = workflow.containerEngine in ['singularity', 'apptainer'] ? '--bind' : '--volume' + "${mount} ${params.gurobi_license}:/opt/gurobi/gurobi.lic:ro --env GRB_LICENSE_FILE=/opt/gurobi/gurobi.lic" + } + publishDir = [ + path: { "${params.outdir}/${meta.id}/ecdna/coral" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + withName: '.*:CORAL_SEED' { + ext.args = { + [ + "--gain ${params.coral_gain}", + "--min-seed-size ${params.coral_min_seed_size}", + "--max-seg-gap ${params.coral_max_seg_gap}", + ].join(' ') + } + } + + withName: '.*:CORAL_RECONSTRUCT' { + // Gurobi WLS licences cap concurrent solver sessions, so serialise under Gurobi. + // 1000 stands in for "no extra limit"; the executor's queueSize already bounds this. + maxForks = params.coral_solver == 'gurobi_direct' ? 1 : 1000 + ext.args = { + [ + "--solver ${params.coral_solver}", + "--solver-threads ${task.cpus}", + "--solver-time-limit ${params.coral_solver_time_limit}", + "--global-time-limit ${params.coral_global_time_limit}", + "--min-bp-support ${params.coral_min_bp_support}", + ].join(' ') + } + } + + withName: '.*:CORAL_CYCLE' { + // Gurobi WLS licences cap concurrent solver sessions, so serialise under Gurobi. + // 1000 stands in for "no extra limit"; the executor's queueSize already bounds this. + maxForks = params.coral_solver == 'gurobi_direct' ? 1 : 1000 + ext.args = { + [ + "--solver ${params.coral_solver}", + "--threads ${task.cpus}", + "--solver-time-limit ${params.coral_solver_time_limit}", + "--global-time-limit ${params.coral_global_time_limit}", + "--alpha ${params.coral_cycle_decomp_alpha}", + ].join(' ') + } + publishDir = [ + path: { "${params.outdir}/${meta.id}/ecdna/coral/cycles" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + withName: '.*:CORAL_PLOT' { + publishDir = [ + path: { "${params.outdir}/${meta.id}/ecdna/coral/plots" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + withName: '.*:AMPLICONCLASSIFIER' { + publishDir = [ + path: { "${params.outdir}/${meta.id}/ecdna/amplicon_classifier" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + withName : '.*:WAKHAN' { ext.args = { [ @@ -656,7 +824,7 @@ process { ] } - withName : '.*:UNTAR' { + withName : '.*:UNTAR(_AA_DATA_REPO)?' { publishDir = [ enabled: false ] @@ -676,10 +844,17 @@ process { } // The wget container carries no CA bundle, as for WGET above; the pinned MD5 still checks the download - withName : '.*:VEPPLUGIN_CLINVAR' { + withName : '.*:(VEPPLUGIN_CLINVAR|AADATAREPO_DOWNLOAD)' { ext.args = { "--no-check-certificate" } } + // Unpacked into the work dir only: pass --aa_data_repo to reuse an extracted copy + withName : '.*:AADATAREPO_DOWNLOAD' { + publishDir = [ + enabled: false + ] + } + // Published so a later run can skip both the download and the reshaping by pointing // --vep_revel / --vep_eve / --vep_clinvar (and their _tbi) at these files withName : '.*:VEPPLUGIN_(REVEL|EVE|CLINVAR)' { @@ -710,6 +885,20 @@ process { ] } + // PASS-only per-caller copies for downstream steps; ClairS-TO is covered by VCFSPLIT. + // --write-index=tbi: the module's index output is optional and the join needs it. + withName: '.*:(CLAIR3|CLAIRS|DEEPVARIANT|DEEPSOMATIC)_PASS_FILTER' { + ext.args = '--apply-filters PASS --output-type z --write-index=tbi' + publishDir = [ + enabled: false + ] + } + + withName: '.*:CLAIR3_PASS_FILTER' { ext.prefix = { "${meta.id}_clair3_pass" } } + withName: '.*:CLAIRS_PASS_FILTER' { ext.prefix = { "${meta.id}_clairs_pass" } } + withName: '.*:DEEPVARIANT_PASS_FILTER' { ext.prefix = { "${meta.id}_deepvariant_pass" } } + withName: '.*:DEEPSOMATIC_PASS_FILTER' { ext.prefix = { "${meta.id}_deepsomatic_pass" } } + withName : '.*:SIGPROFILER_MATRIXGENERATOR' { ext.args = { params.sigprofiler_matrix_args ?: '' } publishDir = [ diff --git a/conf/test.config b/conf/test.config index 1829a02e..3de55ac9 100644 --- a/conf/test.config +++ b/conf/test.config @@ -69,6 +69,8 @@ params { vep_species = "caenorhabditis_elegans" skip_wakhan = true skip_ascat = true + // CoRAL seeds from ASCAT, which is skipped here, and the chr19 slice has no amplification + skip_coral = true skip_modkit = true savana_chromosomes = "19" // SAVANA's het-SNP coverage/mapq floors (--allele_min_reads default 10, --allele_mapq diff --git a/containers/ampliconclassifier.Dockerfile b/containers/ampliconclassifier.Dockerfile new file mode 100644 index 00000000..7fa4545d --- /dev/null +++ b/containers/ampliconclassifier.Dockerfile @@ -0,0 +1,48 @@ +# AmpliconClassifier v2.0.0 with T2T-CHM13 support, from a pinned fork commit. +# +# v2.0.0 is the first release with official CoRAL support. Installed with +# `pip install .` rather than a PATH symlink so the pinned BFBArchitect +# dependency, the console entry points and the bundled CHM13 lncRNA GFF3 +# resource all land correctly. +FROM docker.io/library/python:3.11-slim + +ARG AC_REPO=https://github.com/robert-a-forsyth/AmpliconClassifier.git +ARG AC_REF=cdeaa63 +ARG AC_VERSION=2.0.0 + +LABEL org.opencontainers.image.title="AmpliconClassifier" \ + org.opencontainers.image.description="AmpliconClassifier 2.0.0, T2T-CHM13 fork" \ + org.opencontainers.image.source="${AC_REPO}" \ + org.opencontainers.image.revision="${AC_REF}" \ + org.opencontainers.image.version="${AC_VERSION}" \ + org.opencontainers.image.licenses="BSD-2-Clause" + +ENV LANG=C.UTF-8 \ + PYTHONDONTWRITEBYTECODE=1 + +RUN apt-get update && apt-get install -y --no-install-recommends \ + gcc g++ git procps ca-certificates \ + zlib1g-dev libbz2-dev liblzma-dev libcurl4-openssl-dev libssl-dev \ + coinor-cbc \ + && rm -rf /var/lib/apt/lists/* + +# BFBArchitect==1.0.1 comes in as a pinned dependency and pulls PuLP, CNVkit, +# pysam and matplotlib. gurobipy arrives too but needs no licence: BFBArchitect +# falls back Gurobi -> MOSEK -> CBC, so the image runs licence-free. +RUN pip install --no-cache-dir --upgrade pip \ + && git clone "${AC_REPO}" /opt/AmpliconClassifier \ + && git -C /opt/AmpliconClassifier checkout "${AC_REF}" \ + && pip install --no-cache-dir /opt/AmpliconClassifier \ + && rm -rf /opt/AmpliconClassifier/.git + +# Mount point only. The AA data repo is ~1.1 GB of third-party-derived +# annotation with no stated licence, so it is staged at runtime. +RUN mkdir -p /opt/data_repo +ENV AA_DATA_REPO=/opt/data_repo + +RUN amplicon_classifier.py --version \ + && python -c "import ampclasslib.ac_util as u; \ +p = u.get_ncrna_file_loc('CHM13'); \ +import os; assert os.path.exists(p), p; print('CHM13 lncRNA resource OK')" + +CMD ["amplicon_classifier.py", "--help"] diff --git a/containers/coral.Dockerfile b/containers/coral.Dockerfile new file mode 100644 index 00000000..612e6ad5 --- /dev/null +++ b/containers/coral.Dockerfile @@ -0,0 +1,55 @@ +# CoRAL with T2T-CHM13 support, built from a pinned commit of the fork. +# +# micromamba rather than python:slim because SCIP -- the open-source global MINLP +# solver CoRAL needs for its non-convex MIQCP -- is only packaged on conda-forge. +# Everything is installed into the base env and PATH is set explicitly, so the +# image works without shell activation (Nextflow/Singularity run no entrypoint). +FROM docker.io/mambaorg/micromamba:2.0.5-debian12-slim + +ARG CORAL_REPO=https://github.com/robert-a-forsyth/CoRAL.git +ARG CORAL_REF=847f3d4 +ARG CORAL_VERSION=3.0.0 + +LABEL org.opencontainers.image.title="CoRAL" \ + org.opencontainers.image.description="CoRAL amplicon reconstruction, T2T-CHM13 fork" \ + org.opencontainers.image.source="${CORAL_REPO}" \ + org.opencontainers.image.revision="${CORAL_REF}" \ + org.opencontainers.image.version="${CORAL_VERSION}" \ + org.opencontainers.image.licenses="BSD-3-Clause" + +USER root +ENV PATH=/opt/conda/bin:$PATH \ + LANG=C.UTF-8 \ + PYTHONDONTWRITEBYTECODE=1 + +# Build toolchain and the headers pysam and cvxopt fail without (upstream README) +RUN apt-get update && apt-get install -y --no-install-recommends \ + gcc g++ git make pkg-config curl procps ca-certificates \ + libhdf5-dev libbz2-dev liblzma-dev zlib1g-dev \ + libcurl4-openssl-dev libssl-dev libffi-dev \ + libsuitesparse-dev \ + && rm -rf /var/lib/apt/lists/* + +# scip provides bin/scip, whose built-in AMPL reader is what Pyomo's SCIPAMPL +# plugin shells out to. pyscipopt is deliberately not installed -- it is unused. +# htslib is not needed: pysam installs from a manylinux wheel. +RUN micromamba install -y -n base -c conda-forge python=3.12 scip=10.1.0 \ + && micromamba clean --all --yes + +# CPU-only torch first, or pomegranate/cnvkit pull the multi-GB CUDA wheel +RUN pip install --no-cache-dir --upgrade pip \ + && pip install --no-cache-dir torch --index-url https://download.pytorch.org/whl/cpu + +RUN git clone "${CORAL_REPO}" /opt/CoRAL \ + && git -C /opt/CoRAL checkout "${CORAL_REF}" \ + && pip install --no-cache-dir /opt/CoRAL \ + && rm -rf /opt/CoRAL/.git + +# Mount point only. Gurobi licences are user-supplied and never baked in. +RUN mkdir -p /opt/gurobi +ENV GRB_LICENSE_FILE=/opt/gurobi/gurobi.lic + +# No `| head`: SCIP dies on SIGPIPE when the reader closes early (exit 141) +RUN coral --help > /dev/null && scip -v > /dev/null + +CMD ["coral", "--help"] diff --git a/docs/output.md b/docs/output.md index 54fb647c..f0bd0515 100644 --- a/docs/output.md +++ b/docs/output.md @@ -358,6 +358,55 @@ The germline/somatic split comes from a panel of normals and from ClairS-TO's Ve | `read_qual.txt` | file containing quality statistics about identified segements | | `severus.log` | log file | +#### `ecdna` + +Extrachromosomal DNA and focal amplification. CoRAL reconstructs amplicon structures from the tumour +BAM, seeded by ASCAT's copy-number segments, and AmpliconClassifier labels each reconstructed +amplicon (ecDNA, BFB, linear, and so on). Runs on every tumour sample with ASCAT calls; disable with +`--skip_coral`, or keep reconstruction and drop classification with `--skip_ampliconclassifier`. + +A sample with no segment above `--coral_gain` produces an empty seed BED and is skipped with a log +message rather than failing. CoRAL defaults to the open-source SCIP solver; `--coral_solver +gurobi_direct` is faster but needs `--gurobi_license`. AmpliconClassifier needs an AmpliconArchitect +data repository, downloaded automatically for GRCh38 and supplied with `--aa_data_repo` for CHM13; +a CHM13 run without one skips classification with a warning. + +``` +├── ecdna +│ ├── sample_coral_cn.bed +│ ├── coral +│ │ ├── sample_CNV_SEEDS.bed +│ │ ├── reconstruct +│ │ │ ├── sample_amplicon1_graph.txt +│ │ │ ├── sample_amplicon1_cycles.txt +│ │ │ ├── sample_summary.txt +│ │ │ └── sample_reconstruct.log +│ │ └── plots +│ │ ├── sample_amplicon1_graph.png +│ │ └── sample_amplicon1_cycles.png +│ └── amplicon_classifier +│ ├── sample_amplicon_classification_profiles.tsv +│ ├── sample_gene_list.tsv +│ ├── sample_ecDNA_counts.tsv +│ ├── sample_result_table.tsv +│ └── sample_classification_bed_files/ +``` + +| File | Description | +| --------------------------------------------- | ------------------------------------------------------------------ | +| `sample_coral_cn.bed` | ASCAT's copy number as the BED CoRAL seeds from | +| `sample_CNV_SEEDS.bed` | Amplified intervals above `--coral_gain`; empty means no amplicons | +| `sample_amplicon_graph.txt` | Breakpoint graph per amplicon, in AmpliconArchitect format | +| `sample_amplicon_cycles.txt` | Decomposed cycles and paths per amplicon | +| `sample_summary.txt` | Per-run amplicon summary; written even when no amplicon is found | +| `sample_reconstruct.log` | CoRAL reconstruction log, including solver output | +| `sample_amplicon_{graph,cycles}.png` | Per-amplicon copy-number and cycle plots | +| `sample_amplicon_classification_profiles.tsv` | The headline call per amplicon: ecDNA+, BFB+, decomposition class | +| `sample_gene_list.tsv` | Genes intersecting each classified amplicon | +| `sample_ecDNA_counts.tsv` | Number of distinct ecDNA species detected | +| `sample_result_table.tsv` | Combined per-sample table, the format AmpliconRepository ingests | +| `sample_classification_bed_files/` | Per-feature BED intervals for each classified amplicon | + #### `savana` SAVANA structural variant and copy-number calling. Runs alongside Severus/ASCAT rather than replacing @@ -453,12 +502,12 @@ Phased variant calls produced by Longphase. Present in all samples. │ ├── somatic_smallvariants.vcf.gz.tbi ``` -| File | Description | -| ----------------------------------- | ---------------------------------------------------------------- | -| `germline_smallvariants.vcf.gz` | Longphase-phased germline SNV/indel VCF with haplotype (PS) tags | -| `germline_smallvariants.vcf.gz.tbi` | Index for the phased germline VCF | -| `somatic_smallvariants.vcf.gz` | Longphase-phased somatic SNV/indel VCF with haplotype (PS) tags | -| `somatic_smallvariants.vcf.gz.tbi` | Index for the phased somatic VCF | +| File | Description | +| ----------------------------------- | ---------------------------------------------------------------------------------------------------------- | +| `germline_smallvariants.vcf.gz` | Longphase-phased germline SNV/indel VCF with haplotype (PS) tags | +| `germline_smallvariants.vcf.gz.tbi` | Index for the phased germline VCF | +| `somatic_smallvariants.vcf.gz` | Longphase-phased somatic SNV/indel VCF with haplotype (PS) tags; `0/0` and `./.` records are left unphased | +| `somatic_smallvariants.vcf.gz.tbi` | Index for the phased somatic VCF | diff --git a/docs/usage.md b/docs/usage.md index 7d822ecb..5089434b 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -135,6 +135,16 @@ For structural variants, the CHM13 panel of normals is a merged panel combining For tumour-only small variants, ClairS-TO separates germline from somatic calls with a panel of normals and with its Verdict module, which tags each call as germline, somatic or subclonal somatic from tumour purity and allele-specific copy number. `--genome CHM13` supplies five CHM13 PON VCFs (gnomAD, dbSNP, 1000 Genomes, CoLoRSdb and ASAP), which **replace** the GRCh38 databases inside the container. Unless `--skip_ascat` is set, purity and copy number come from the pipeline's own ASCAT run (`CLAIRSTO_VERDICT_TAG`); only with `--skip_ascat` does ClairS-TO estimate them itself, from assembly-specific loci, allele and GC content files. A GRCh38 resource set on a CHM13 run leaves germline variants untagged. +When `--germline_var_keep` includes `deepvariant`, the tumour-only germline arm +runs DeepVariant on the **tumour** BAM, so on its own those calls mix germline and +clonal somatic variants. When `deepsomatic` is also in `--somatic_var_keep`, the +pipeline records DeepSomatic's verdict at each site in `INFO/DS_VERDICT` and keeps +only sites it calls `GERMLINE` or `PON`; `RefCall` and sites DeepSomatic never +evaluated are dropped. Without `deepsomatic`, the DeepVariant germline calls are +used without a verdict filter and may include somatic variants. Either way the +tumour-only germline arm is a tumour-derived proxy, not a call set from normal +tissue. + With `--genome CHM13 --skip_ascat` the pipeline builds a CHM13 resource set from the ASCAT files it already downloads, so no extra setup is needed. LogR correction is GC-only, as ClairS-TO recommends for CHM13: no replication timing file is published for the assembly. Without `--skip_ascat` nothing is built, because the tagging comes from ASCAT's own tables. To use a resource set of your own — another assembly, or a CHM13 set carrying an `RT_.txt` for replication timing correction — pass `--clairsto_cna_resources` together with `--skip_ascat`. Without `--skip_ascat` it is ignored, with a warning. Expected layout: @@ -384,14 +394,74 @@ Both tools run from `ghcr.io/ljwharbers/sigprofiler`, which adds CHM13 support n These options control how variants from multiple callers are filtered and merged. -| Parameter | Description | -| ------------------------------ | --------------------------------------------------------------------------------------------------- | -| `--germline_var_keep` | Expression or threshold for retaining germline variants after calling. Default = `null` | -| `--somatic_var_keep` | Expression or threshold for retaining somatic variants after calling. Default = `null` | -| `--germline_var_combine` | Strategy for combining germline variant caller outputs (e.g. union, intersection). Default = `null` | -| `--somatic_var_combine` | Strategy for combining somatic variant caller outputs (e.g. union, intersection). Default = `null` | -| `--prioritize_caller_germline` | Comma-separated caller priority order used when combining germline calls. Default = `null` | -| `--prioritize_caller_somatic` | Comma-separated caller priority order used when combining somatic calls. Default = `null` | +| Parameter | Description | +| ------------------------------ | ----------------------------------------------------------------------------------------------------------- | +| `--germline_var_keep` | Comma-separated germline callers to run: `deepvariant`, `clair`. Default = `clair` | +| `--somatic_var_keep` | Comma-separated somatic callers to run: `deepsomatic`, `clair`. Default = `clair` | +| `--germline_var_combine` | How to combine germline caller outputs: `consensus` (shared calls only) or `all` (union). Default = `all` | +| `--somatic_var_combine` | How to combine somatic caller outputs: `consensus` (shared calls only) or `all` (union). Default = `all` | +| `--prioritize_caller_germline` | Whose record to use where both germline callers call a variant: `deepvariant` or `clair`. Default = `clair` | +| `--prioritize_caller_somatic` | Whose record to use where both somatic callers call a variant: `deepsomatic` or `clair`. Default = `clair` | +| `--smallvar_filter_pass` | Keep only PASS records from each small variant caller downstream. Default = `true` | + +DeepVariant and DeepSomatic emit a record for every site they evaluate, not only +for the variants they call, so most of their records are `RefCall`, `GERMLINE` or +`PON`; Clair3 and ClairS also keep their `LowQual` and `NonSomatic` records. Without +a `PASS` filter, `*_var_combine = 'all'` would carry all of these into the phased +VCFs and the mutation burden. + +`--smallvar_filter_pass` (`true` by default) restricts the copy of each caller's +VCF that is handed to the caller consensus, phasing, VEP and the report. In +tumor-only mode ClairS-TO is unaffected by the setting: `VCFSPLIT` already +restricts its somatic split to `PASS`, and its germline split is `PASS`-rewritten +rather than `PASS`-filtered. The per-caller VCFs published under +`//variants/` are never filtered, so no calls are lost +from the results directory. + +Set it to `false` to restore the previous unfiltered behaviour: each caller's +records are passed on with their original `FILTER`. Only the ClairS-TO germline +split is normalised to `PASS`, with its original value kept in +`INFO/ORIG_FILTER`. + +`consensus` keeps only alleles called by both callers, using the prioritised +caller's record. Multi-allelic records are split so each allele can be matched +across callers, and rejoined before phasing. A multi-allelic call of which only one +allele is shared (e.g. DeepVariant `1/2`) is therefore kept as that allele alone, +and its `PL` values for the dropped allele are lost. + +`all` keeps the union by position: every record of the prioritised caller, plus +the other caller's records at positions where the prioritised caller has none. +Where both callers call a position, even with different alleles, only the +prioritised caller's record is kept, so each output record is one caller's call. +Records are not split in this mode. + +#### Germline and somatic provenance + +Germline and somatic small variants are merged into one VCF for somatic phasing, +because Longphase needs all variant sites in a single file to produce consistent +phase blocks. The somatic arm is then recovered from the phased result by an +`INFO/SOMATIC` flag stamped on each arm before the merge, not by position, since a +germline record at the same coordinate as a somatic call would otherwise be kept. +Tagging leaves `FILTER` unchanged. + +Longphase phases by position, so a germline record at the position of a somatic +call would lend the somatic record its genotype and phase set. Germline records +at the position of any somatic call with an alternate genotype are therefore +left out of somatic phasing. Somatic records without an alternate genotype +(`0/0` or `./.`, present only with `--smallvar_filter_pass false`) are not +phased: they are added back to `somatic_smallvariants.vcf.gz` unchanged. + +Three INFO fields carry this provenance: + +| Field | Meaning | +| ------------- | --------------------------------------------------------------------------------------------------------------------- | +| `SOMATIC` | Record came from the somatic call set | +| `GERMLINE` | Record came from the germline call set | +| `ORIG_FILTER` | Original `FILTER` of ClairS-TO germline records, before normalisation to `PASS`. Multiple filters are joined with `,` | + +Germline calls dropped from `variants/phased/somatic_smallvariants.vcf.gz` are not +lost: they remain in `variants/phased/germline_smallvariants.vcf.gz`, in +`vep/germline/`, and in the unfiltered per-caller VCFs under `variants//`. #### PON Options @@ -619,6 +689,67 @@ Two of these predictors get there anyway, because they score _proteins_ rather t > for `--custom` files, but a name it cannot map annotates nothing rather than raising an error, so > confirm that `CLNSIG` values appear in an annotated VCF before trusting them. +## ecDNA and focal amplification + +CoRAL reconstructs amplicon structures and AmpliconClassifier labels them. Both run by default on +every tumour sample that has ASCAT calls. CoRAL seeds from ASCAT's copy number, so `--skip_ascat` +together with CoRAL is rejected at launch: use `--skip_coral` as well, or keep ASCAT. + +### Solver + +CoRAL's cycle decomposition is a non-convex mixed-integer quadratically-constrained problem. The +pipeline defaults to `--coral_solver scip`, which is open-source, shipped in the image and needs no +licence. Gurobi is available with: + +```bash +--coral_solver gurobi_direct --gurobi_license /path/to/gurobi.lic +``` + +The licence is mounted into the CoRAL tasks (`--bind` under Singularity/Apptainer, `--volume` under +Docker/Podman), so give an absolute path; it is never baked into the image. Under Gurobi the +reconstruction steps are serialised (`maxForks = 1`), because a Web License Service licence caps +concurrent solver sessions, and the solver uses the task's CPUs. SCIP runs single-threaded. + +Gurobi is faster, but on the models this pipeline produces that has not mattered: across 1102 solver +logs from the earlier standalone cohort, the largest model was 620 rows by 438 columns and the +slowest single solve took 0.49 s, with no model reaching the time limit. Prefer the default unless +you have measured a reason not to. + +### AmpliconClassifier reference data + +AmpliconClassifier reads an AmpliconArchitect data repository at runtime. For `--genome GRCh38` it is +downloaded automatically (about 1.1 GB) and checked against the MD5 the host publishes +(`--aa_data_repo_md5`), so a re-published repository fails the run rather than changing results +silently. The download and its ~4 GB unpack repeat on every run; to avoid both, extract +`GRCh38.tar.gz` once and pass it: + +```bash +--aa_data_repo /path/to/GRCh38 +``` + +For `--genome CHM13` there is no published repository. Without `--aa_data_repo` the classifier is +skipped with a warning and reconstruction still runs: + +```bash +--genome CHM13 --aa_data_repo /path/to/AA_DATA_REPO +``` + +`--aa_data_repo` may be the reference directory itself (`GRCh38/`, `CHM13/`) or the directory that +contains it. `--skip_ampliconclassifier` drops classification explicitly and needs no reference data. + +### Tuning + +`--coral_gain` (default 6.0) sets the total copy number a segment must reach to seed an amplicon. A +sample with nothing above it produces an empty seed file and is skipped with a log message rather +than failing. + +`--coral_min_bp_support` defaults to 1.75, the value validated on this lab's cohort. CoRAL's own +documentation recommends a considerably higher value for WGS, around 10.0; raise it if you see +spurious breakpoints. + +`--coral_run_cycle` re-extracts cycles with `coral cycle_all` after reconstruction and classifies +those instead of the originals. It is off by default. + ## Core Nextflow arguments > [!NOTE] diff --git a/modules/local/aadatarepo/download/environment.yml b/modules/local/aadatarepo/download/environment.yml new file mode 100644 index 00000000..daa07275 --- /dev/null +++ b/modules/local/aadatarepo/download/environment.yml @@ -0,0 +1,8 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + # renovate: datasource=conda depName=conda-forge/wget + - conda-forge::wget=1.21.4 diff --git a/modules/local/aadatarepo/download/main.nf b/modules/local/aadatarepo/download/main.nf new file mode 100644 index 00000000..f9056e34 --- /dev/null +++ b/modules/local/aadatarepo/download/main.nf @@ -0,0 +1,52 @@ +process AADATAREPO_DOWNLOAD { + tag "${url.toString().tokenize('/').last()}" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container + ? 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/3b/3b54fa9135194c72a18d00db6b399c03248103f87e43ca75e4b50d61179994b3/data' + : 'community.wave.seqera.io/library/wget:1.21.4--8b0fcde81c17be5e'}" + + input: + tuple val(meta), val(url), val(md5) + + output: + tuple val(meta), path("${archive_name}"), emit: archive + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + archive_name = url.toString().tokenize('/').last() + def check = md5 ? "echo '${md5} ${archive_name}' | md5sum -c -" : '' + """ + wget \\ + --no-verbose \\ + --tries=10 \\ + --waitretry=10 \\ + ${args} \\ + -O ${archive_name} \\ + ${url} + + # The tarball is updated in place, so a pinned MD5 fails the task on a re-published repo + ${check} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + wget: \$(wget --version | head -1 | cut -d ' ' -f 3) + END_VERSIONS + """ + + stub: + archive_name = url.toString().tokenize('/').last() + """ + echo "" | gzip > ${archive_name} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + wget: \$(wget --version | head -1 | cut -d ' ' -f 3) + END_VERSIONS + """ +} diff --git a/modules/local/aadatarepo/download/meta.yml b/modules/local/aadatarepo/download/meta.yml new file mode 100644 index 00000000..71fff961 --- /dev/null +++ b/modules/local/aadatarepo/download/meta.yml @@ -0,0 +1,48 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "aadatarepo_download" +description: Download the AmpliconArchitect data repository tarball that AmpliconClassifier reads, checked against a pinned MD5 +keywords: + - ecdna + - ampliconclassifier + - download +tools: + - "wget": + description: "GNU Wget is a free software package for retrieving files using HTTP, HTTPS, FTP and FTPS" + homepage: "https://www.gnu.org/software/wget/" + documentation: "https://www.gnu.org/software/wget/manual/wget.html" + licence: ["GPL-3.0-or-later"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing an id for the data repository + - url: + type: string + description: URL of the data repository tarball; its basename becomes the output name + - md5: + type: string + description: | + Expected MD5 of the tarball. When set, a download that does not match fails + the task. May be null, in which case the repository is not verified. + +output: + archive: + - - meta: + type: map + description: Groovy Map containing an id for the data repository + - "${archive_name}": + type: file + description: The data repository tarball + pattern: "*.tar.gz" + versions: + - versions.yml: + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/aadatarepo/download/tests/main.nf.test b/modules/local/aadatarepo/download/tests/main.nf.test new file mode 100644 index 00000000..579696e0 --- /dev/null +++ b/modules/local/aadatarepo/download/tests/main.nf.test @@ -0,0 +1,38 @@ +nextflow_process { + + name "Test Process AADATAREPO_DOWNLOAD" + script "../main.nf" + process "AADATAREPO_DOWNLOAD" + + tag "small" + tag "modules" + tag "modules_local" + tag "aadatarepo_download" + + // Stub only: a real run downloads the 1.1 GB repository + test("GRCh38 data repo - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id: 'aa_data_repo' ], + 'https://refs.ampliconrepository.org/data/module_support_files/AmpliconArchitect/GRCh38.tar.gz', + '2bcd1fdaed027466ed296bf78b2907eb' + ] + """ + } + } + + then { + assert process.success + assertAll( + { assert file(process.out.archive[0][1]).name == 'GRCh38.tar.gz' }, + { assert snapshot(process.out).match() } + ) + } + + } +} diff --git a/modules/local/aadatarepo/download/tests/main.nf.test.snap b/modules/local/aadatarepo/download/tests/main.nf.test.snap new file mode 100644 index 00000000..083aa38f --- /dev/null +++ b/modules/local/aadatarepo/download/tests/main.nf.test.snap @@ -0,0 +1,35 @@ +{ + "GRCh38 data repo - stub": { + "content": [ + { + "0": [ + [ + { + "id": "aa_data_repo" + }, + "GRCh38.tar.gz:md5,68b329da9893e34099c7d8ad5cb9c940" + ] + ], + "1": [ + "versions.yml:md5,3d07074650f444dad718bd48add94911" + ], + "archive": [ + [ + { + "id": "aa_data_repo" + }, + "GRCh38.tar.gz:md5,68b329da9893e34099c7d8ad5cb9c940" + ] + ], + "versions": [ + "versions.yml:md5,3d07074650f444dad718bd48add94911" + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-30T16:21:36.209439147" + } +} \ No newline at end of file diff --git a/modules/local/ampliconclassifier/main.nf b/modules/local/ampliconclassifier/main.nf new file mode 100644 index 00000000..820386c1 --- /dev/null +++ b/modules/local/ampliconclassifier/main.nf @@ -0,0 +1,62 @@ +process AMPLICONCLASSIFIER { + tag "$meta.id" + label 'process_medium' + + // No conda: bioconda's `ampliconclassifier` recipe is stuck at 0.4.14 (2023) and + // predates CoRAL support; this image carries the CHM13 fork of v2.0.0. See meta.yml + container "docker.io/robertaforsyth/ampliconclassifier:2.0.0-chm13-cdeaa63" + + input: + tuple val(meta), path(reconstruction) + tuple val(meta2), path(data_repo, stageAs: 'aa_data_repo') + val(ac_ref) + + output: + tuple val(meta), path("*_amplicon_classification_profiles.tsv"), emit: classification, optional: true + tuple val(meta), path("*_gene_list.tsv"), emit: gene_list, optional: true + tuple val(meta), path("*_ecDNA_counts.tsv"), emit: ecdna_counts, optional: true + tuple val(meta), path("*_result_table.tsv"), emit: result_table, optional: true + tuple val(meta), path("*_classification_bed_files", type: 'dir'), emit: bed_files, optional: true + tuple val(meta), path("*_SV_summaries", type: 'dir'), emit: sv_summaries, optional: true + tuple val(meta), path("*_annotated_cycles_files", type: 'dir'), emit: annotated_cycles, optional: true + tuple val(meta), path("*.log"), emit: log, optional: true + tuple val("${task.process}"), val('ampliconclassifier'), eval("amplicon_classifier.py --version"), topic: versions, emit: versions_ampliconclassifier + + when: + task.ext.when == null || task.ext.when + + script: + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "AMPLICONCLASSIFIER does not support Conda. Please use Docker / Singularity / Apptainer instead." + } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + """ + # A downloaded repo arrives as the dir itself (UNTAR strips it); AC wants its parent + if [ -e aa_data_repo/file_list.txt ]; then + mkdir repo + ln -s "\$(readlink -f aa_data_repo)" repo/${ac_ref} + export AA_DATA_REPO=\$PWD/repo + else + export AA_DATA_REPO=\$(readlink -f aa_data_repo) + fi + + amplicon_classifier.py \\ + --ref ${ac_ref} \\ + --AA_results ${reconstruction} \\ + -o ${prefix} \\ + ${args} \\ + 2>&1 | tee ${prefix}_classifier.log + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}_amplicon_classification_profiles.tsv + touch ${prefix}_gene_list.tsv + touch ${prefix}_ecDNA_counts.tsv + touch ${prefix}_result_table.tsv + touch ${prefix}_classifier.log + mkdir -p ${prefix}_classification_bed_files ${prefix}_SV_summaries ${prefix}_annotated_cycles_files + """ +} diff --git a/modules/local/ampliconclassifier/meta.yml b/modules/local/ampliconclassifier/meta.yml new file mode 100644 index 00000000..7b4f7b66 --- /dev/null +++ b/modules/local/ampliconclassifier/meta.yml @@ -0,0 +1,133 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "ampliconclassifier" +description: Classify reconstructed amplicons as ecDNA, BFB, linear or no amplification. +keywords: + - ecdna + - amplicon + - classification + - bfb +tools: + - "ampliconclassifier": + description: "Classifies focal amplifications from AmpliconArchitect or CoRAL output." + homepage: "https://github.com/AmpliconSuite/AmpliconClassifier" + documentation: "https://github.com/AmpliconSuite/AmpliconClassifier" + doi: "10.1038/s41586-023-05937-5" + licence: ["BSD-2-Clause"] + identifier: "" + +## No environment.yml: bioconda's `ampliconclassifier` recipe is stuck at 0.4.14 (2023) and +## predates CoRAL support, which arrived in v2.0.0. Built from a fork of v2.0.0 adding +## T2T-CHM13 support (https://github.com/robert-a-forsyth/AmpliconClassifier, branch +## chm13-support). See containers/ampliconclassifier.Dockerfile. + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - reconstruction: + type: directory + description: | + A CoRAL output directory holding matching *_graph.txt, *_cycles.txt and + *_summary.txt files. AmpliconClassifier pairs them by shared prefix. + - - meta2: + type: map + description: Groovy Map for the data repository + - data_repo: + type: directory + description: | + AmpliconArchitect data repository, exported as $AA_DATA_REPO. Staged under a + fixed name so the tool never writes into the shared reference directory. + - - ac_ref: + type: string + description: Reference genome name, one of hg19, GRCh37, GRCh38, mm10, CHM13 + +output: + classification: + - meta: + type: map + description: Groovy Map containing sample information + - "*_amplicon_classification_profiles.tsv": + type: file + description: Per-amplicon class, ecDNA+ and BFB+ calls + pattern: "*_amplicon_classification_profiles.tsv" + gene_list: + - meta: + type: map + description: Groovy Map containing sample information + - "*_gene_list.tsv": + type: file + description: Genes intersecting each classified amplicon + pattern: "*_gene_list.tsv" + ecdna_counts: + - meta: + type: map + description: Groovy Map containing sample information + - "*_ecDNA_counts.tsv": + type: file + description: Number of distinct ecDNA species detected + pattern: "*_ecDNA_counts.tsv" + result_table: + - meta: + type: map + description: Groovy Map containing sample information + - "*_result_table.tsv": + type: file + description: Combined per-sample table, the AmpliconRepository ingest format + pattern: "*_result_table.tsv" + bed_files: + - meta: + type: map + description: Groovy Map containing sample information + - "*_classification_bed_files": + type: directory + description: Per-feature BED intervals for each classified amplicon + sv_summaries: + - meta: + type: map + description: Groovy Map containing sample information + - "*_SV_summaries": + type: directory + description: Per-amplicon structural variant summaries + annotated_cycles: + - meta: + type: map + description: Groovy Map containing sample information + - "*_annotated_cycles_files": + type: directory + description: Cycles files annotated with classification + log: + - meta: + type: map + description: Groovy Map containing sample information + - "*.log": + type: file + description: Classifier log + pattern: "*.log" + versions_ampliconclassifier: + - - ${task.process}: + type: string + description: Process name + - ampliconclassifier: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - ampliconclassifier: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/ampliconclassifier/tests/main.nf.test b/modules/local/ampliconclassifier/tests/main.nf.test new file mode 100644 index 00000000..c6e26e88 --- /dev/null +++ b/modules/local/ampliconclassifier/tests/main.nf.test @@ -0,0 +1,35 @@ +nextflow_process { + + name "Test Process AMPLICONCLASSIFIER" + script "../main.nf" + process "AMPLICONCLASSIFIER" + + tag "modules" + tag "modules_local" + tag "ampliconclassifier" + tag "small" + + // Stub only: the real tool needs a BAM, a solver and, for the classifier, a ~1 GB data repo. + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ [ id:'sample1' ], file('reconstruct') ] + input[1] = [ [:], file('aa_data_repo') ] + input[2] = 'CHM13' + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert process.out.classification[0][1].toString().endsWith('_amplicon_classification_profiles.tsv') }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/ampliconclassifier/tests/main.nf.test.snap b/modules/local/ampliconclassifier/tests/main.nf.test.snap new file mode 100644 index 00000000..39b9dce2 --- /dev/null +++ b/modules/local/ampliconclassifier/tests/main.nf.test.snap @@ -0,0 +1,167 @@ +{ + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon_classification_profiles.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "sample1" + }, + "sample1_gene_list.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + { + "id": "sample1" + }, + "sample1_ecDNA_counts.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + [ + { + "id": "sample1" + }, + "sample1_result_table.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "4": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "5": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "6": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "7": [ + [ + { + "id": "sample1" + }, + "sample1_classifier.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "8": [ + [ + "AMPLICONCLASSIFIER", + "ampliconclassifier", + "2.0.0" + ] + ], + "annotated_cycles": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "bed_files": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "classification": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon_classification_profiles.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "ecdna_counts": [ + [ + { + "id": "sample1" + }, + "sample1_ecDNA_counts.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "gene_list": [ + [ + { + "id": "sample1" + }, + "sample1_gene_list.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "log": [ + [ + { + "id": "sample1" + }, + "sample1_classifier.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "result_table": [ + [ + { + "id": "sample1" + }, + "sample1_result_table.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "sv_summaries": [ + [ + { + "id": "sample1" + }, + [ + + ] + ] + ], + "versions_ampliconclassifier": [ + [ + "AMPLICONCLASSIFIER", + "ampliconclassifier", + "2.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:50:03.780875635" + } +} \ No newline at end of file diff --git a/modules/local/ascattocoralbed/environment.yml b/modules/local/ascattocoralbed/environment.yml new file mode 100644 index 00000000..c9b822b2 --- /dev/null +++ b/modules/local/ascattocoralbed/environment.yml @@ -0,0 +1,7 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - "conda-forge::python=3.12" diff --git a/modules/local/ascattocoralbed/main.nf b/modules/local/ascattocoralbed/main.nf new file mode 100644 index 00000000..ee75800a --- /dev/null +++ b/modules/local/ascattocoralbed/main.nf @@ -0,0 +1,38 @@ +process ASCAT_TO_CORAL_BED { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/python:3.12': + 'biocontainers/python:3.12' }" + + input: + tuple val(meta), path(cnvs) + tuple val(meta2), path(fai) + + output: + tuple val(meta), path("*_coral_cn.bed"), emit: bed + tuple val("${task.process}"), val('python'), eval("python --version | sed 's/Python //'"), topic: versions, emit: versions_python + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def fai_arg = fai ? "--fai ${fai}" : '' + """ + ascat_to_coral_bed.py \\ + --cnvs ${cnvs} \\ + ${fai_arg} \\ + --output ${prefix}_coral_cn.bed \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}_coral_cn.bed + """ +} diff --git a/modules/local/ascattocoralbed/meta.yml b/modules/local/ascattocoralbed/meta.yml new file mode 100644 index 00000000..e81343f8 --- /dev/null +++ b/modules/local/ascattocoralbed/meta.yml @@ -0,0 +1,70 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "ascat_to_coral_bed" +description: Convert ASCAT copy-number calls into the BED CoRAL seeds from. +keywords: + - copy number + - ascat + - ecdna +tools: + - "python": + description: "Python, running bin/ascat_to_coral_bed.py from this pipeline." + homepage: "https://www.python.org" + documentation: "https://docs.python.org" + licence: ["PSF-2.0"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - cnvs: + type: file + description: ASCAT cnvs.txt with chr/startpos/endpos/nMajor/nMinor columns + pattern: "*.cnvs.txt" + - - meta2: + type: map + description: Groovy Map for the reference + - fai: + type: file + description: | + Reference index, used to respell contigs. ASCAT writes "1" where the BAM may + say "chr1", and CoRAL builds chromosome sizes from the BAM header. Optional. + pattern: "*.fai" + +output: + bed: + - meta: + type: map + description: Groovy Map containing sample information + - "*_coral_cn.bed": + type: file + description: Headerless BED4 with total copy number in the last column + pattern: "*_coral_cn.bed" + versions_python: + - - ${task.process}: + type: string + description: Process name + - python: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - python: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/ascattocoralbed/tests/fixtures/chr.fai b/modules/local/ascattocoralbed/tests/fixtures/chr.fai new file mode 100644 index 00000000..a862f050 --- /dev/null +++ b/modules/local/ascattocoralbed/tests/fixtures/chr.fai @@ -0,0 +1,3 @@ +chr1 248956422 0 60 61 +chr2 248956422 0 60 61 +chrX 248956422 0 60 61 diff --git a/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt b/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt new file mode 100644 index 00000000..40394b02 --- /dev/null +++ b/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt @@ -0,0 +1,6 @@ +chr startpos endpos nMajor nMinor +1 100000 200000 1 1 +1 250000 400000 5 3 +1 500000 500000 2 1 +2 100000 300000 2 0 +X 100000 200000 1 0 diff --git a/modules/local/ascattocoralbed/tests/main.nf.test b/modules/local/ascattocoralbed/tests/main.nf.test new file mode 100644 index 00000000..560fac5c --- /dev/null +++ b/modules/local/ascattocoralbed/tests/main.nf.test @@ -0,0 +1,94 @@ +nextflow_process { + + name "Test Process ASCAT_TO_CORAL_BED" + script "../main.nf" + process "ASCAT_TO_CORAL_BED" + + tag "modules" + tag "modules_local" + tag "ascat_to_coral_bed" + tag "small" + + // Fixture: an ASCAT cnvs.txt calling its contigs "1"/"2"/"X", with one zero-length + // segment that CoRAL's own parser would reject, and a .fai that spells them "chr1" etc. + test("converts ASCAT cnvs.txt to CoRAL BED and respells contigs") { + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) + ] + input[1] = [ + [:], + file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) + ] + """ + } + } + + then { + def bed = path(process.out.bed[0][1]).readLines() + + assertAll( + { assert process.success }, + // The zero-length segment (chr1:500000-500000) is dropped + { assert bed.size() == 4 }, + // Contigs respelled to match the reference + { assert bed.every { it.startsWith('chr') } }, + // Total CN is nMajor + nMinor, in the last column + { assert bed[0].split('\t')[3] == '2' }, + { assert bed[1].split('\t')[3] == '8' }, + // Natural contig order: chr1, chr2, then chrX + { assert bed.collect { it.split('\t')[0] } == ['chr1', 'chr1', 'chr2', 'chrX'] }, + { assert snapshot(process.out.bed).match() } + ) + } + } + + test("no reference index - contigs are left as ASCAT spelled them") { + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) + ] + input[1] = [ [:], [] ] + """ + } + } + + then { + def bed = path(process.out.bed[0][1]).readLines() + + assertAll( + { assert process.success }, + { assert bed.collect { it.split('\t')[0] } == ['1', '1', '2', 'X'] } + ) + } + } + + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ [ id:'sample1' ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ] + input[1] = [ [:], [] ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/ascattocoralbed/tests/main.nf.test.snap b/modules/local/ascattocoralbed/tests/main.nf.test.snap new file mode 100644 index 00000000..70007849 --- /dev/null +++ b/modules/local/ascattocoralbed/tests/main.nf.test.snap @@ -0,0 +1,60 @@ +{ + "converts ASCAT cnvs.txt to CoRAL BED and respells contigs": { + "content": [ + [ + [ + { + "id": "sample1" + }, + "sample1_coral_cn.bed:md5,771563ffb2c54c0fddbe3aafb288e984" + ] + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:36:07.311149784" + }, + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + "sample1_coral_cn.bed:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + "ASCAT_TO_CORAL_BED", + "python", + "3.12.2" + ] + ], + "bed": [ + [ + { + "id": "sample1" + }, + "sample1_coral_cn.bed:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions_python": [ + [ + "ASCAT_TO_CORAL_BED", + "python", + "3.12.2" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:36:16.961617365" + } +} \ No newline at end of file diff --git a/modules/local/bcftools/excludesites/environment.yml b/modules/local/bcftools/excludesites/environment.yml new file mode 100644 index 00000000..cb55500b --- /dev/null +++ b/modules/local/bcftools/excludesites/environment.yml @@ -0,0 +1,9 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + # renovate: datasource=conda depName=bioconda/htslib + - bioconda::bcftools=1.22 + - bioconda::htslib=1.22.1 diff --git a/modules/local/bcftools/excludesites/main.nf b/modules/local/bcftools/excludesites/main.nf new file mode 100644 index 00000000..e0d74e41 --- /dev/null +++ b/modules/local/bcftools/excludesites/main.nf @@ -0,0 +1,41 @@ +process BCFTOOLS_EXCLUDE_SITES { + tag "${meta.id}" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container + ? 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/47/474a5ea8dc03366b04df884d89aeacc4f8e6d1ad92266888e7a8e7958d07cde8/data' + : 'community.wave.seqera.io/library/bcftools_htslib:0a3fa2654b52006f'}" + + input: + tuple val(meta), path(vcf), path(tbi), path(mask), path(mask_tbi) + + output: + tuple val(meta), path("${prefix}.vcf.gz"), emit: vcf + tuple val(meta), path("${prefix}.vcf.gz.tbi"), emit: tbi + tuple val("${task.process}"), val('bcftools'), eval("bcftools --version | sed '1!d; s/^.*bcftools //'"), topic: versions, emit: versions_bcftools + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: "${meta.id}_excluded" + """ + # Drop every record at a CHROM:POS present in the mask, whatever its alleles. + bcftools view \\ + -T ^${mask} \\ + -Oz \\ + -W=tbi \\ + ${args} \\ + -o ${prefix}.vcf.gz \\ + ${vcf} + """ + + stub: + prefix = task.ext.prefix ?: "${meta.id}_excluded" + """ + echo '' | gzip > ${prefix}.vcf.gz + touch ${prefix}.vcf.gz.tbi + """ +} diff --git a/modules/local/bcftools/excludesites/meta.yml b/modules/local/bcftools/excludesites/meta.yml new file mode 100644 index 00000000..5a7d3d92 --- /dev/null +++ b/modules/local/bcftools/excludesites/meta.yml @@ -0,0 +1,77 @@ +name: bcftools_exclude_sites +description: Remove every record of a VCF at a position (CHROM:POS) that a mask VCF also has, regardless of alleles +keywords: + - filtering + - VCF + - positions +tools: + - view: + description: VCF/BCF conversion, view, subset and filter VCF/BCF files. + homepage: http://samtools.github.io/bcftools/bcftools.html + documentation: http://www.htslib.org/doc/bcftools.html + tool_dev_url: https://github.com/samtools/bcftools + doi: "10.1093/bioinformatics/btp352" + licence: ["MIT"] + identifier: biotools:bcftools +input: + - - meta: + type: map + description: Groovy Map containing sample information e.g. [ id:'test' ] + - vcf: + type: file + description: VCF to remove records from + pattern: "*.vcf.gz" + - tbi: + type: file + description: Tabix index of the VCF + pattern: "*.tbi" + - mask: + type: file + description: VCF whose record positions are removed from the input + pattern: "*.vcf.gz" + - mask_tbi: + type: file + description: Tabix index of the mask VCF + pattern: "*.tbi" +output: + vcf: + - - meta: + type: map + description: Groovy Map containing sample information + - ${prefix}.vcf.gz: + type: file + description: Input VCF without the records at mask positions + pattern: "*.vcf.gz" + tbi: + - - meta: + type: map + description: Groovy Map containing sample information + - ${prefix}.vcf.gz.tbi: + type: file + description: Tabix index of the output VCF + pattern: "*.tbi" + versions_bcftools: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version | sed '1!d; s/^.*bcftools //'": + type: string + description: The command used to generate the version of the tool +topics: + versions: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version | sed '1!d; s/^.*bcftools //'": + type: string + description: The command used to generate the version of the tool +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/bcftools/excludesites/tests/main.nf.test b/modules/local/bcftools/excludesites/tests/main.nf.test new file mode 100644 index 00000000..be7933cb --- /dev/null +++ b/modules/local/bcftools/excludesites/tests/main.nf.test @@ -0,0 +1,55 @@ +nextflow_process { + + name "Test Process BCFTOOLS_EXCLUDE_SITES" + script "../main.nf" + process "BCFTOOLS_EXCLUDE_SITES" + + tag "modules" + tag "modules_local" + tag "bcftools_exclude_sites" + + test("drops records at mask positions by position only") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file("\${projectDir}/tests/fixtures/excludesites_input.vcf.gz", checkIfExists: true), + file("\${projectDir}/tests/fixtures/excludesites_input.vcf.gz.tbi", checkIfExists: true), + file("\${projectDir}/tests/fixtures/excludesites_mask.vcf.gz", checkIfExists: true), + file("\${projectDir}/tests/fixtures/excludesites_mask.vcf.gz.tbi", checkIfExists: true) + ] + """ + } + } + + then { + assert process.success + + def positions = path(process.out.vcf[0][1]).linesGzip + .findAll { line -> !line.startsWith('#') } + .collect { line -> line.split('\t')[1] } + + // 100 and 200 go despite different alleles; 102 stays although it lies inside the mask's deletion + assert positions == [ '102', '300' ] + } + } + + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ [ id:'test' ], [], [], [], [] ] + """ + } + } + + then { + assert process.success + } + } +} diff --git a/modules/local/bcftools/view/main.nf b/modules/local/bcftools/view/main.nf index 652da9ac..63905dd2 100644 --- a/modules/local/bcftools/view/main.nf +++ b/modules/local/bcftools/view/main.nf @@ -8,7 +8,7 @@ process BCFTOOLS_VIEW { : 'community.wave.seqera.io/library/bcftools_htslib:0a3fa2654b52006f'}" input: - tuple val(meta), path(vcf), path(tbi), path(targets), path(targets_tbi) + tuple val(meta), path(vcf), path(tbi) output: tuple val(meta), path("*.vcf.gz"), emit: vcf @@ -23,7 +23,6 @@ process BCFTOOLS_VIEW { def prefix = task.ext.prefix ?: "${meta.id}" """ bcftools view \\ - -T ${targets} \\ -Oz \\ -W=tbi \\ ${args} \\ diff --git a/modules/local/bcftools/view/meta.yml b/modules/local/bcftools/view/meta.yml index 28c0fc48..29e3646b 100644 --- a/modules/local/bcftools/view/meta.yml +++ b/modules/local/bcftools/view/meta.yml @@ -1,5 +1,5 @@ name: bcftools_view -description: Filter VCF to positions defined by a targets file using bcftools view -T +description: Filter a VCF with bcftools view, with the filter expression supplied via ext.args; outputs a bgzipped VCF and tbi index keywords: - filtering - VCF @@ -25,14 +25,6 @@ input: type: file description: Tabix index of the input VCF pattern: "*.tbi" - - targets: - type: file - description: VCF file used as position filter (-T) - pattern: "*.{vcf.gz,vcf,bcf}" - - targets_tbi: - type: file - description: Tabix index of the targets VCF - pattern: "*.tbi" output: vcf: - - meta: @@ -50,7 +42,28 @@ output: type: file description: Tabix index of filtered VCF pattern: "*.tbi" + versions_bcftools: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version | sed '1!d; s/^.*bcftools //'": + type: string + description: The command used to generate the version of the tool +topics: + versions: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version | sed '1!d; s/^.*bcftools //'": + type: string + description: The command used to generate the version of the tool authors: - - "@rforsyth" + - "@robert-a-forsyth" maintainers: - - "@rforsyth" + - "@robert-a-forsyth" diff --git a/modules/local/coral/cycle/main.nf b/modules/local/coral/cycle/main.nf new file mode 100644 index 00000000..207e9c31 --- /dev/null +++ b/modules/local/coral/cycle/main.nf @@ -0,0 +1,47 @@ +process CORAL_CYCLE { + tag "$meta.id" + label 'process_low' + + // As for CORAL_RECONSTRUCT: the subworkflow warns on the missing output + errorStrategy { task.exitStatus in 130..145 ? 'retry' : 'ignore' } + + container "docker.io/robertaforsyth/coral:3.0.0-chm13-847f3d4" + + input: + tuple val(meta), path(reconstruction) + + output: + tuple val(meta), path("cycles"), emit: reconstruction + tuple val(meta), path("cycles/*_amplicon*_cycles.txt"), emit: cycles, optional: true + tuple val("${task.process}"), val('coral'), eval("python -c 'import importlib.metadata as m; print(m.version(\"CoRAL\"))'"), topic: versions, emit: versions_coral + + when: + task.ext.when == null || task.ext.when + + script: + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "CORAL_CYCLE does not support Conda. Please use Docker / Singularity / Apptainer instead." + } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + // Graphs are copied alongside the re-extracted cycles so AmpliconClassifier + // still finds a matching graph/cycles/summary set in one directory. + """ + mkdir -p cycles + cp ${reconstruction}/*_graph.txt ${reconstruction}/*_summary.txt cycles/ 2>/dev/null || true + + coral cycle_all \\ + --bp-dir ${reconstruction} \\ + --output-prefix cycles/${prefix} \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + mkdir -p cycles + touch cycles/${prefix}_amplicon1_cycles.txt + touch cycles/${prefix}_amplicon1_graph.txt + touch cycles/${prefix}_summary.txt + """ +} diff --git a/modules/local/coral/cycle/meta.yml b/modules/local/coral/cycle/meta.yml new file mode 100644 index 00000000..d45cca19 --- /dev/null +++ b/modules/local/coral/cycle/meta.yml @@ -0,0 +1,74 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "coral_cycle" +description: Re-extract cycle decompositions from existing breakpoint graphs. +keywords: + - ecdna + - amplicon + - cycle decomposition +tools: + - "coral": + description: "Reconstructs extrachromosomal DNA structures from long-read sequencing." + homepage: "https://github.com/AmpliconSuite/CoRAL" + documentation: "https://github.com/AmpliconSuite/CoRAL" + doi: "10.1101/gr.279131.124" + licence: ["BSD-3-Clause"] + identifier: "" + +## No environment.yml: CoRAL is not on bioconda -- the `coral` recipe there is an unrelated +## RNA-seq read-bridging tool. Built from a fork adding T2T-CHM13 support and two SCIP fixes +## (https://github.com/robert-a-forsyth/CoRAL, branch chm13-support). + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - reconstruction: + type: directory + description: Output directory from `coral reconstruct` + +output: + reconstruction: + - meta: + type: map + description: Groovy Map containing sample information + - "cycles": + type: directory + description: | + Re-extracted cycles alongside the graph and summary files copied from the + reconstruction, so AmpliconClassifier still sees a complete set. + cycles: + - meta: + type: map + description: Groovy Map containing sample information + - "cycles/*_amplicon*_cycles.txt": + type: file + description: Re-extracted per-amplicon cycles + pattern: "*_amplicon*_cycles.txt" + versions_coral: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/coral/cycle/tests/main.nf.test b/modules/local/coral/cycle/tests/main.nf.test new file mode 100644 index 00000000..5402f2e7 --- /dev/null +++ b/modules/local/coral/cycle/tests/main.nf.test @@ -0,0 +1,36 @@ +nextflow_process { + + name "Test Process CORAL_CYCLE" + script "../main.nf" + process "CORAL_CYCLE" + + tag "modules" + tag "modules_local" + tag "coral_cycle" + tag "small" + + // Stub only: the real tool needs a BAM, a solver and, for the classifier, a ~1 GB data repo. + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file('reconstruct') + ] + """ + } + } + + then { + assertAll( + { assert process.success }, + + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/coral/cycle/tests/main.nf.test.snap b/modules/local/coral/cycle/tests/main.nf.test.snap new file mode 100644 index 00000000..c4ed4f81 --- /dev/null +++ b/modules/local/coral/cycle/tests/main.nf.test.snap @@ -0,0 +1,67 @@ +{ + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "1": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + "CORAL_CYCLE", + "coral", + "3.0.0" + ] + ], + "cycles": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "reconstruction": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "versions_coral": [ + [ + "CORAL_CYCLE", + "coral", + "3.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:50:09.602318674" + } +} \ No newline at end of file diff --git a/modules/local/coral/plot/main.nf b/modules/local/coral/plot/main.nf new file mode 100644 index 00000000..dc5d023b --- /dev/null +++ b/modules/local/coral/plot/main.nf @@ -0,0 +1,42 @@ +process CORAL_PLOT { + tag "$meta.id" + label 'process_medium' + + // Plotting is cosmetic: never let it fail a run that reconstructed successfully + errorStrategy { task.exitStatus in 130..145 ? 'retry' : 'ignore' } + + container "docker.io/robertaforsyth/coral:3.0.0-chm13-847f3d4" + + input: + tuple val(meta), path(reconstruction), path(bam), path(bai) + val(coral_ref) + + output: + tuple val(meta), path("*_amplicon*.{pdf,png}"), emit: plots, optional: true + tuple val("${task.process}"), val('coral'), eval("python -c 'import importlib.metadata as m; print(m.version(\"CoRAL\"))'"), topic: versions, emit: versions_coral + + when: + task.ext.when == null || task.ext.when + + script: + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "CORAL_PLOT does not support Conda. Please use Docker / Singularity / Apptainer instead." + } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + """ + coral plot_all \\ + --ref ${coral_ref} \\ + --reconstruction-dir ${reconstruction} \\ + --bam ${bam} \\ + --output-prefix ${prefix} \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}_amplicon1_graph.png + touch ${prefix}_amplicon1_cycles.png + """ +} diff --git a/modules/local/coral/plot/meta.yml b/modules/local/coral/plot/meta.yml new file mode 100644 index 00000000..754fe5b0 --- /dev/null +++ b/modules/local/coral/plot/meta.yml @@ -0,0 +1,76 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "coral_plot" +description: Render per-amplicon copy-number and cycle plots. +keywords: + - ecdna + - amplicon + - visualisation +tools: + - "coral": + description: "Reconstructs extrachromosomal DNA structures from long-read sequencing." + homepage: "https://github.com/AmpliconSuite/CoRAL" + documentation: "https://github.com/AmpliconSuite/CoRAL" + doi: "10.1101/gr.279131.124" + licence: ["BSD-3-Clause"] + identifier: "" + +## No environment.yml: CoRAL is not on bioconda -- the `coral` recipe there is an unrelated +## RNA-seq read-bridging tool. Built from a fork adding T2T-CHM13 support and two SCIP fixes +## (https://github.com/robert-a-forsyth/CoRAL, branch chm13-support). + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - reconstruction: + type: directory + description: Output directory from `coral reconstruct` + - bam: + type: file + description: Tumour BAM, for the coverage track + pattern: "*.bam" + - bai: + type: file + description: BAM index + pattern: "*.bai" + - - coral_ref: + type: string + description: Reference genome name, one of hg19, hg38, mm10, t2t + +output: + plots: + - meta: + type: map + description: Groovy Map containing sample information + - "*_amplicon*.{pdf,png}": + type: file + description: Per-amplicon graph and cycle plots + pattern: "*_amplicon*.{pdf,png}" + versions_coral: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/coral/plot/tests/main.nf.test b/modules/local/coral/plot/tests/main.nf.test new file mode 100644 index 00000000..0e0a8c94 --- /dev/null +++ b/modules/local/coral/plot/tests/main.nf.test @@ -0,0 +1,39 @@ +nextflow_process { + + name "Test Process CORAL_PLOT" + script "../main.nf" + process "CORAL_PLOT" + + tag "modules" + tag "modules_local" + tag "coral_plot" + tag "small" + + // Stub only: the real tool needs a BAM, a solver and, for the classifier, a ~1 GB data repo. + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file('reconstruct'), + file('sample1.bam'), + file('sample1.bam.bai') + ] + input[1] = 't2t' + """ + } + } + + then { + assertAll( + { assert process.success }, + + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/coral/plot/tests/main.nf.test.snap b/modules/local/coral/plot/tests/main.nf.test.snap new file mode 100644 index 00000000..608bb693 --- /dev/null +++ b/modules/local/coral/plot/tests/main.nf.test.snap @@ -0,0 +1,49 @@ +{ + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.png:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "1": [ + [ + "CORAL_PLOT", + "coral", + "3.0.0" + ] + ], + "plots": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.png:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "versions_coral": [ + [ + "CORAL_PLOT", + "coral", + "3.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:50:15.950672699" + } +} \ No newline at end of file diff --git a/modules/local/coral/reconstruct/main.nf b/modules/local/coral/reconstruct/main.nf new file mode 100644 index 00000000..76094ba8 --- /dev/null +++ b/modules/local/coral/reconstruct/main.nf @@ -0,0 +1,53 @@ +process CORAL_RECONSTRUCT { + tag "$meta.id" + label 'process_medium' + + // 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' } + + container "docker.io/robertaforsyth/coral:3.0.0-chm13-847f3d4" + + input: + tuple val(meta), path(seeds), path(cn_seg), path(bam), path(bai) + + output: + tuple val(meta), path("reconstruct"), emit: reconstruction + tuple val(meta), path("reconstruct/*_amplicon*_graph.txt"), emit: graphs, optional: true + tuple val(meta), path("reconstruct/*_amplicon*_cycles.txt"), emit: cycles, optional: true + tuple val(meta), path("reconstruct/*_summary.txt"), emit: summary, optional: true + tuple val(meta), path("reconstruct/*_reconstruct.log"), emit: log, optional: true + tuple val("${task.process}"), val('coral'), eval("python -c 'import importlib.metadata as m; print(m.version(\"CoRAL\"))'"), topic: versions, emit: versions_coral + + when: + task.ext.when == null || task.ext.when + + script: + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "CORAL_RECONSTRUCT does not support Conda. Please use Docker / Singularity / Apptainer instead." + } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + // CoRAL derives its output directory by splitting --output-prefix on '/', and + // AmpliconClassifier pairs graph/cycles/summary by prefix within one directory. + """ + mkdir -p reconstruct + + coral reconstruct \\ + --lr-bam ${bam} \\ + --cnv-seed ${seeds} \\ + --cn-seg ${cn_seg} \\ + --output-prefix reconstruct/${prefix} \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + mkdir -p reconstruct + touch reconstruct/${prefix}_amplicon1_graph.txt + touch reconstruct/${prefix}_amplicon1_cycles.txt + touch reconstruct/${prefix}_summary.txt + touch reconstruct/${prefix}_reconstruct.log + """ +} diff --git a/modules/local/coral/reconstruct/meta.yml b/modules/local/coral/reconstruct/meta.yml new file mode 100644 index 00000000..f9373390 --- /dev/null +++ b/modules/local/coral/reconstruct/meta.yml @@ -0,0 +1,112 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "coral_reconstruct" +description: Reconstruct amplicon breakpoint graphs and decompose them into cycles. +keywords: + - ecdna + - amplicon + - structural variation + - long-read +tools: + - "coral": + description: "Reconstructs extrachromosomal DNA structures from long-read sequencing." + homepage: "https://github.com/AmpliconSuite/CoRAL" + documentation: "https://github.com/AmpliconSuite/CoRAL" + doi: "10.1101/gr.279131.124" + licence: ["BSD-3-Clause"] + identifier: "" + +## No environment.yml: CoRAL is not on bioconda -- the `coral` recipe there is an unrelated +## RNA-seq read-bridging tool. Built from a fork adding T2T-CHM13 support and two SCIP fixes +## (https://github.com/robert-a-forsyth/CoRAL, branch chm13-support). + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - seeds: + type: file + description: Seed intervals from `coral seed` + pattern: "*_CNV_SEEDS.bed" + - cn_seg: + type: file + description: Copy-number segments with total CN in the last column + pattern: "*.{bed,cns}" + - bam: + type: file + description: Tumour BAM + pattern: "*.bam" + - bai: + type: file + description: BAM index + pattern: "*.bai" + +output: + reconstruction: + - meta: + type: map + description: Groovy Map containing sample information + - "reconstruct": + type: directory + description: | + Graph, cycles and summary files together. AmpliconClassifier pairs them by + shared prefix, so nothing inside may be renamed. + graphs: + - meta: + type: map + description: Groovy Map containing sample information + - "reconstruct/*_amplicon*_graph.txt": + type: file + description: Per-amplicon breakpoint graph, AmpliconArchitect format + pattern: "*_amplicon*_graph.txt" + cycles: + - meta: + type: map + description: Groovy Map containing sample information + - "reconstruct/*_amplicon*_cycles.txt": + type: file + description: Per-amplicon decomposed cycles and paths + pattern: "*_amplicon*_cycles.txt" + summary: + - meta: + type: map + description: Groovy Map containing sample information + - "reconstruct/*_summary.txt": + type: file + description: Per-run amplicon summary, written even with no amplicons + pattern: "*_summary.txt" + log: + - meta: + type: map + description: Groovy Map containing sample information + - "reconstruct/*_reconstruct.log": + type: file + description: Reconstruction log including solver output + pattern: "*_reconstruct.log" + versions_coral: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/coral/reconstruct/tests/main.nf.test b/modules/local/coral/reconstruct/tests/main.nf.test new file mode 100644 index 00000000..c746a5d5 --- /dev/null +++ b/modules/local/coral/reconstruct/tests/main.nf.test @@ -0,0 +1,42 @@ +nextflow_process { + + name "Test Process CORAL_RECONSTRUCT" + script "../main.nf" + process "CORAL_RECONSTRUCT" + + tag "modules" + tag "modules_local" + tag "coral_reconstruct" + tag "small" + + // Stub only: the real tool needs a BAM, a solver and, for the classifier, a ~1 GB data repo. + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file('sample1_CNV_SEEDS.bed'), + file('sample1_coral_cn.bed'), + file('sample1.bam'), + file('sample1.bam.bai') + ] + """ + } + } + + then { + assertAll( + { assert process.success }, + // AmpliconClassifier pairs graph/cycles/summary by prefix in one directory + { assert process.out.graphs[0][1].toString().endsWith('_amplicon1_graph.txt') }, + { assert process.out.cycles[0][1].toString().endsWith('_amplicon1_cycles.txt') }, + { assert process.out.summary[0][1].toString().endsWith('_summary.txt') }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/coral/reconstruct/tests/main.nf.test.snap b/modules/local/coral/reconstruct/tests/main.nf.test.snap new file mode 100644 index 00000000..0d60bf14 --- /dev/null +++ b/modules/local/coral/reconstruct/tests/main.nf.test.snap @@ -0,0 +1,117 @@ +{ + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_reconstruct.log:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "1": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + [ + { + "id": "sample1" + }, + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "4": [ + [ + { + "id": "sample1" + }, + "sample1_reconstruct.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "5": [ + [ + "CORAL_RECONSTRUCT", + "coral", + "3.0.0" + ] + ], + "cycles": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "graphs": [ + [ + { + "id": "sample1" + }, + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "log": [ + [ + { + "id": "sample1" + }, + "sample1_reconstruct.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "reconstruction": [ + [ + { + "id": "sample1" + }, + [ + "sample1_amplicon1_cycles.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_amplicon1_graph.txt:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_reconstruct.log:md5,d41d8cd98f00b204e9800998ecf8427e", + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "summary": [ + [ + { + "id": "sample1" + }, + "sample1_summary.txt:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions_coral": [ + [ + "CORAL_RECONSTRUCT", + "coral", + "3.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:50:23.471598843" + } +} \ No newline at end of file diff --git a/modules/local/coral/seed/main.nf b/modules/local/coral/seed/main.nf new file mode 100644 index 00000000..2ee14460 --- /dev/null +++ b/modules/local/coral/seed/main.nf @@ -0,0 +1,42 @@ +process CORAL_SEED { + tag "$meta.id" + label 'process_low' + + // No conda: CoRAL is not packaged on bioconda (the `coral` recipe there is an + // unrelated RNA-seq tool), and this image carries the CHM13 fork. See meta.yml + container "docker.io/robertaforsyth/coral:3.0.0-chm13-847f3d4" + + input: + tuple val(meta), path(cn_seg), path(bam), path(bai) + val(coral_ref) + + output: + tuple val(meta), path("*_CNV_SEEDS.bed"), emit: seeds + tuple val("${task.process}"), val('coral'), eval("python -c 'import importlib.metadata as m; print(m.version(\"CoRAL\"))'"), topic: versions, emit: versions_coral + + when: + task.ext.when == null || task.ext.when + + script: + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "CORAL_SEED does not support Conda. Please use Docker / Singularity / Apptainer instead." + } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + """ + coral seed \\ + --cn-seg ${cn_seg} \\ + --ref ${coral_ref} \\ + --lr-bam ${bam} \\ + --output-prefix ${prefix} \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + // A non-empty seed by default; meta.stub_empty_seed, set only by the ECDNA test, reaches the empty-seed branch + def seed_cmd = meta.stub_empty_seed ? "touch ${prefix}_CNV_SEEDS.bed" : "printf 'chr1\\t100000\\t400000\\t8\\n' > ${prefix}_CNV_SEEDS.bed" + """ + ${seed_cmd} + """ +} diff --git a/modules/local/coral/seed/meta.yml b/modules/local/coral/seed/meta.yml new file mode 100644 index 00000000..b75e9118 --- /dev/null +++ b/modules/local/coral/seed/meta.yml @@ -0,0 +1,81 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "coral_seed" +description: Identify amplified intervals that seed CoRAL amplicon reconstruction. +keywords: + - ecdna + - amplicon + - copy number + - long-read +tools: + - "coral": + description: "Reconstructs extrachromosomal DNA structures from long-read sequencing." + homepage: "https://github.com/AmpliconSuite/CoRAL" + documentation: "https://github.com/AmpliconSuite/CoRAL" + doi: "10.1101/gr.279131.124" + licence: ["BSD-3-Clause"] + identifier: "" + +## No environment.yml: CoRAL is not on bioconda -- the `coral` recipe there is an unrelated +## RNA-seq read-bridging tool. The container is built from a fork adding T2T-CHM13 support +## (https://github.com/robert-a-forsyth/CoRAL, branch chm13-support), which also fixes two +## upstream bugs that made `--solver scip` fail. See containers/coral.Dockerfile. + +input: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - cn_seg: + type: file + description: Copy-number segments with total CN in the last column + pattern: "*.{bed,cns}" + - bam: + type: file + description: Tumour BAM; CoRAL reads chromosome sizes from its header + pattern: "*.bam" + - bai: + type: file + description: BAM index + pattern: "*.bai" + - - coral_ref: + type: string + description: Reference genome name, one of hg19, hg38, mm10, t2t + +output: + seeds: + - meta: + type: map + description: Groovy Map containing sample information + - "*_CNV_SEEDS.bed": + type: file + description: Amplified seed intervals; empty when nothing exceeds --gain + pattern: "*_CNV_SEEDS.bed" + versions_coral: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - coral: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/coral/seed/tests/main.nf.test b/modules/local/coral/seed/tests/main.nf.test new file mode 100644 index 00000000..6f3610b2 --- /dev/null +++ b/modules/local/coral/seed/tests/main.nf.test @@ -0,0 +1,41 @@ +nextflow_process { + + name "Test Process CORAL_SEED" + script "../main.nf" + process "CORAL_SEED" + + tag "modules" + tag "modules_local" + tag "coral" + tag "coral_seed" + tag "small" + + // Stub only: seeding needs a real BAM, since CoRAL reads chromosome sizes from its header. + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'sample1' ], + file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true), + file('sample1.bam'), + file('sample1.bam.bai') + ] + input[1] = 't2t' + """ + } + } + + then { + assertAll( + { assert process.success }, + // Non-empty, so the subworkflow's empty-seed branch is not taken + { assert path(process.out.seeds[0][1]).readLines().size() == 1 }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/coral/seed/tests/main.nf.test.snap b/modules/local/coral/seed/tests/main.nf.test.snap new file mode 100644 index 00000000..bcca1d64 --- /dev/null +++ b/modules/local/coral/seed/tests/main.nf.test.snap @@ -0,0 +1,43 @@ +{ + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "sample1" + }, + "sample1_CNV_SEEDS.bed:md5,f679880d9757f33af43363b38254f3f0" + ] + ], + "1": [ + [ + "CORAL_SEED", + "coral", + "3.0.0" + ] + ], + "seeds": [ + [ + { + "id": "sample1" + }, + "sample1_CNV_SEEDS.bed:md5,f679880d9757f33af43363b38254f3f0" + ] + ], + "versions_coral": [ + [ + "CORAL_SEED", + "coral", + "3.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-24T09:50:31.018971049" + } +} \ No newline at end of file diff --git a/modules/local/vcfsplit/main.nf b/modules/local/vcfsplit/main.nf index f6156d34..55cb6300 100644 --- a/modules/local/vcfsplit/main.nf +++ b/modules/local/vcfsplit/main.nf @@ -38,10 +38,21 @@ process VCFSPLIT { bcftools concat -a -Oz -o germline_tmp.vcf.gz indels_filtered.vcf.gz snv_filtered.vcf.gz tabix -p vcf germline_tmp.vcf.gz - bcftools view germline_tmp.vcf.gz | awk 'BEGIN{FS=OFS="\t"} /^#/ {print} !/^#/ { \$7="PASS"; print }' | \ - bgzip -c > germline.vcf.gz + # Normalise FILTER to PASS, keeping the original in INFO/ORIG_FILTER (";" stored as ","). + bcftools view germline_tmp.vcf.gz | awk -v q='"' 'BEGIN{FS=OFS="\t"} + /^##/ { print; next } + /^#CHROM/ { print "##INFO="; print; next } + { of = \$7; gsub(/;/, ",", of) + \$8 = (\$8 == "." || \$8 == "") ? "ORIG_FILTER=" of : \$8 ";ORIG_FILTER=" of + \$7 = "PASS" + print } + ' | bgzip -c > germline.vcf.gz tabix -p vcf germline.vcf.gz + # Fail here, not downstream, if either header does not parse. + bcftools view -h somatic.vcf.gz > /dev/null + bcftools view -h germline.vcf.gz > /dev/null + # Cleanup intermediate files rm indels_pass.vcf.gz snv_pass.vcf.gz rm indels_pass.vcf.gz.tbi snv_pass.vcf.gz.tbi diff --git a/modules/local/vcftag/environment.yml b/modules/local/vcftag/environment.yml new file mode 100644 index 00000000..b276efd9 --- /dev/null +++ b/modules/local/vcftag/environment.yml @@ -0,0 +1,7 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - bioconda::bcftools=1.20 diff --git a/modules/local/vcftag/main.nf b/modules/local/vcftag/main.nf new file mode 100644 index 00000000..390122c6 --- /dev/null +++ b/modules/local/vcftag/main.nf @@ -0,0 +1,53 @@ +process VCFTAG { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/bcftools:1.20--h8b25389_0': + 'biocontainers/bcftools:1.20--h8b25389_0' }" + + input: + tuple val(meta), path(vcf), path(tbi) + val flag + + output: + tuple val(meta), path("${prefix}.vcf.gz") , emit: vcf + tuple val(meta), path("${prefix}.vcf.gz.tbi") , emit: tbi + + tuple val("${task.process}"), val('bcftools'), eval("bcftools --version |& sed '1!d ; s/bcftools //'"), topic: versions, emit: versions_bcftools + + when: + task.ext.when == null || task.ext.when + + script: + prefix = task.ext.prefix ?: "${meta.id}_${flag.toLowerCase()}" + """ + # Stamp a constant INFO flag marking the call set; FILTER is left as the caller emitted it. + # bcftools annotate cannot set a constant INFO field without an annotation file, hence awk. + bcftools view ${vcf} | awk -v flag="${flag}" -v q='"' 'BEGIN{FS=OFS="\t"} + /^##/ { print; next } + /^#CHROM/ { + print "##INFO=" + print + next + } + { + \$8 = (\$8 == "." || \$8 == "") ? flag : \$8 ";" flag + print + } + ' | bgzip -c > ${prefix}.vcf.gz + + # Fail here, not downstream, if the header does not parse. + bcftools view -h ${prefix}.vcf.gz > /dev/null + + tabix -p vcf ${prefix}.vcf.gz + """ + + stub: + prefix = task.ext.prefix ?: "${meta.id}_${flag.toLowerCase()}" + """ + echo "" | gzip > ${prefix}.vcf.gz + touch ${prefix}.vcf.gz.tbi + """ +} diff --git a/modules/local/vcftag/meta.yml b/modules/local/vcftag/meta.yml new file mode 100644 index 00000000..648ca689 --- /dev/null +++ b/modules/local/vcftag/meta.yml @@ -0,0 +1,75 @@ +name: vcftag +description: Stamp a constant INFO flag, named by the flag input, on every record of a VCF; FILTER is left untouched +keywords: + - vcf + - annotation + - INFO + - variant calling +tools: + - bcftools: + description: Tools for variant calling and manipulating VCFs and BCFs + homepage: http://samtools.github.io/bcftools/bcftools.html + documentation: http://www.htslib.org/doc/bcftools.html + tool_dev_url: https://github.com/samtools/bcftools + doi: "10.1093/gigascience/giab008" + licence: ["MIT"] + identifier: biotools:bcftools +input: + - - meta: + type: map + description: Groovy Map containing sample information e.g. [ id:'test' ] + - vcf: + type: file + description: Input VCF file to tag + pattern: "*.vcf.gz" + - tbi: + type: file + description: Tabix index of the input VCF + pattern: "*.tbi" + - flag: + type: string + description: | + Name of the INFO flag added to every record (e.g. `GERMLINE`); + its lowercase form is also used in the default output prefix +output: + vcf: + - - meta: + type: map + description: Groovy Map containing sample information + - ${prefix}.vcf.gz: + type: file + description: VCF with the INFO flag added to every record + pattern: "*.vcf.gz" + tbi: + - - meta: + type: map + description: Groovy Map containing sample information + - ${prefix}.vcf.gz.tbi: + type: file + description: Tabix index of the tagged VCF + pattern: "*.vcf.gz.tbi" + versions_bcftools: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version |& sed '1!d ; s/bcftools //'": + type: string + description: The command used to generate the version of the tool +topics: + versions: + - - ${task.process}: + type: string + description: The process the versions were collected from + - bcftools: + type: string + description: The tool name + - "bcftools --version |& sed '1!d ; s/bcftools //'": + type: string + description: The command used to generate the version of the tool +authors: + - "@robert-a-forsyth" +maintainers: + - "@robert-a-forsyth" diff --git a/modules/local/vcftag/tests/main.nf.test b/modules/local/vcftag/tests/main.nf.test new file mode 100644 index 00000000..c0ec94a6 --- /dev/null +++ b/modules/local/vcftag/tests/main.nf.test @@ -0,0 +1,76 @@ +nextflow_process { + + name "Test Process VCFTAG" + script "../main.nf" + process "VCFTAG" + + tag "modules" + tag "modules_local" + tag "vcftag" + + // Runs for real (no -stub): the tagging is an awk program embedded in the Nextflow script + // block, so the escaping only holds if it is actually executed. + test("stamps the flag and declares its header without touching FILTER") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file("\${projectDir}/tests/fixtures/vcftag_input.vcf.gz", checkIfExists: true), + file("\${projectDir}/tests/fixtures/vcftag_input.vcf.gz.tbi", checkIfExists: true) + ] + input[1] = 'SOMATIC' + """ + } + } + + then { + assert process.success + + // linesGzip is built into nf-test; the .vcf accessor needs an nft-vcf plugin that + // nf-test.config does not load, so this reads the records directly. + def all = path(process.out.vcf[0][1]).linesGzip + def lines = all.findAll { line -> !line.startsWith('#') } + def header = all.findAll { line -> line.startsWith('##') }.join('\n') + + def filters = lines.collectEntries { l -> def f = l.split('\t'); [ (f[1]): f[6] ] } + def infos = lines.collectEntries { l -> def f = l.split('\t'); [ (f[1]): f[7] ] } + + assertAll( + // every input record survives -- tagging must not filter + { assert lines.size() == 5 }, + // the flag is declared, so bcftools can query it downstream + { assert header.contains('ID=SOMATIC') }, + // the flag is stamped exactly once per record, including the one whose INFO was '.' + { assert infos.values().every { it.split(';').count('SOMATIC') == 1 } }, + { assert infos['300'] == 'SOMATIC' }, + // FILTER is left exactly as the caller emitted it + { assert filters == [ '100':'PASS', '200':'NonSomatic', '300':'RefCall', '400':'.', '500':'LowQual;RefCall' ] }, + // pre-existing INFO is preserved, with the flag appended + { assert infos['200'] == 'CALLER=clairs-to;EXISTING;SOMATIC' }, + { assert infos['500'] == 'CALLER=clair3;EXISTING;SOMATIC' }, + // VCFTAG no longer adds ORIG_FILTER + { assert !all.any { it.contains('ORIG_FILTER') } } + ) + } + } + + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ [ id:'test' ], [], [] ] + input[1] = 'GERMLINE' + """ + } + } + + then { + assert process.success + } + } +} diff --git a/modules/local/wakhan/environment.yml b/modules/local/wakhan/environment.yml index 6b0cdb43..33c3c873 100644 --- a/modules/local/wakhan/environment.yml +++ b/modules/local/wakhan/environment.yml @@ -4,4 +4,4 @@ channels: - conda-forge - bioconda dependencies: - - "bioconda::wakhan=0.4.3" + - "bioconda::wakhan=0.4.4" diff --git a/modules/local/wakhan/main.nf b/modules/local/wakhan/main.nf index ad7aba5e..517ae995 100644 --- a/modules/local/wakhan/main.nf +++ b/modules/local/wakhan/main.nf @@ -4,8 +4,8 @@ process WAKHAN { conda "${moduleDir}/environment.yml" container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'https://depot.galaxyproject.org/singularity/wakhan:0.4.3--pyhdfd78af_0': - 'biocontainers/wakhan:0.4.3--pyhdfd78af_0' }" + 'https://depot.galaxyproject.org/singularity/wakhan:0.4.4--pyhdfd78af_0': + 'biocontainers/wakhan:0.4.4--pyhdfd78af_0' }" input: tuple val(meta), path(tumor_input), path(tumor_index), path(normal_input), path(normal_index), path(vcf), path(breakpoints) @@ -38,9 +38,9 @@ process WAKHAN { tuple val(meta), path("solutions_ranks.tsv") , emit: solutions_ranks // Whole directories, not the plots inside: every solution's plot has the same basename, // and LRSOMATICREPORT resolves them by solution_/ path - tuple val(meta), path("solution_*", type: 'dir') , emit: solution_dirs, optional: true + tuple val(meta), path("solution_*", type: 'dir') , emit: solution_dirs // WARN: Manually update version information as tool does not provide on CLI - tuple val("${task.process}"), val('wakhan'), val("0.4.3"), topic: versions, emit: versions_wakhan + tuple val("${task.process}"), val('wakhan'), val("0.4.4"), topic: versions, emit: versions_wakhan when: task.ext.when == null || task.ext.when diff --git a/nextflow.config b/nextflow.config index 628b5ca3..437856a1 100644 --- a/nextflow.config +++ b/nextflow.config @@ -20,6 +20,7 @@ params { somatic_var_combine = 'all' prioritize_caller_germline = 'clair' prioritize_caller_somatic = 'clair' + smallvar_filter_pass = true generate_gvcf = false // Longphase options @@ -77,6 +78,21 @@ params { sigprofiler_matrix_args = '--plot' sigprofiler_assignment_args = null + // CoRAL / ecDNA options + coral_gain = 6.0 + coral_min_seed_size = 100000 + coral_max_seg_gap = 300000 + coral_solver = 'scip' + coral_solver_time_limit = 7200 + coral_global_time_limit = 21600 + coral_min_bp_support = 1.75 + coral_run_cycle = false + coral_cycle_decomp_alpha = 0.01 + coral_plot = true + gurobi_license = null + aa_data_repo = null + aa_data_repo_md5 = null + // Skip options skip_qc = false skip_cramino = false @@ -97,6 +113,8 @@ params { skip_whatshapstats = false skip_signatures = false skip_report = false + skip_coral = false + skip_ampliconclassifier = false // minimap2 options minimap2_ont_model = null diff --git a/nextflow_schema.json b/nextflow_schema.json index 47356879..d04c3942 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -91,16 +91,23 @@ }, "prioritize_caller_germline": { "type": "string", - "description": "When both germline callers are used, specifies which caller's format to use for variants called by both. Must be [deepvariant, clair].", + "description": "When both germline callers are used, specifies whose record to keep where both callers call a position. Must be [deepvariant, clair].", "default": "clair", "enum": ["deepvariant", "clair"] }, "prioritize_caller_somatic": { "type": "string", - "description": "When both somatic callers are used, specifies which caller's format to use for variants called by both. Must be [deepsomatic, clair].", + "description": "When both somatic callers are used, specifies whose record to keep where both callers call a position. Must be [deepsomatic, clair].", "default": "clair", "enum": ["deepsomatic", "clair"] }, + "smallvar_filter_pass": { + "type": "boolean", + "default": true, + "description": "Keep only PASS records from each small variant caller for downstream use.", + "help_text": "DeepVariant and DeepSomatic emit a record for every site they evaluate, so most records are RefCall (or GERMLINE/PON) rather than calls, and Clair3/ClairS keep their LowQual and NonSomatic calls. Those records otherwise flow into the caller consensus, phasing, VEP and the report. Set to false to consume each caller's unfiltered output. The per-caller VCFs published under `//variants/` are unaffected either way.", + "fa_icon": "fas fa-filter" + }, "generate_gvcf": { "type": "boolean" } @@ -597,6 +604,14 @@ "skip_report": { "type": "boolean", "description": "Skip the final per-sample HTML report" + }, + "skip_coral": { + "type": "boolean", + "description": "Skip CoRAL amplicon reconstruction and AmpliconClassifier." + }, + "skip_ampliconclassifier": { + "type": "boolean", + "description": "Run CoRAL but skip AmpliconClassifier, which avoids needing an AA data repository." } } }, @@ -749,6 +764,82 @@ "description": "Display hidden parameters in the help message (only works when --help or --help_full are provided)." } } + }, + "coral_options": { + "title": "CoRAL / ecDNA options", + "type": "object", + "description": "Amplicon reconstruction with CoRAL and classification with AmpliconClassifier.", + "default": "", + "properties": { + "coral_gain": { + "type": "number", + "default": 6.0, + "description": "Minimum total copy number for a segment to seed an amplicon (CoRAL's --gain)." + }, + "coral_min_seed_size": { + "type": "integer", + "default": 100000, + "description": "Minimum seed interval size in bp (CoRAL's --min-seed-size)." + }, + "coral_max_seg_gap": { + "type": "integer", + "default": 300000, + "description": "Largest gap in bp that CoRAL will bridge when merging seed segments (--max-seg-gap)." + }, + "coral_solver": { + "type": "string", + "default": "scip", + "enum": ["scip", "gurobi_direct"], + "description": "Optimiser for CoRAL's non-convex MIQCP cycle decomposition. SCIP is open-source and needs no licence; gurobi_direct is faster but requires --gurobi_license." + }, + "coral_solver_time_limit": { + "type": "integer", + "default": 7200, + "description": "Per-amplicon solver time limit in seconds (CoRAL's --solver-time-limit)." + }, + "coral_global_time_limit": { + "type": "integer", + "default": 21600, + "description": "Per-sample time limit in seconds across all amplicons (CoRAL's --global-time-limit)." + }, + "coral_min_bp_support": { + "type": "number", + "default": 1.75, + "description": "Minimum breakpoint read support (CoRAL's --min-bp-support). CoRAL's own README recommends a much larger value such as 10.0 for WGS; 1.75 matches the setting validated on this cohort." + }, + "coral_run_cycle": { + "type": "boolean", + "description": "Re-extract cycles with `coral cycle_all` after reconstruction, and classify those instead." + }, + "coral_cycle_decomp_alpha": { + "type": "number", + "default": 0.01, + "description": "Cycle decomposition alpha, used only with --coral_run_cycle (CoRAL's --alpha)." + }, + "coral_plot": { + "type": "boolean", + "default": true, + "description": "Render CoRAL's per-amplicon graph and cycle plots. Plot failures never fail the run." + }, + "gurobi_license": { + "type": "string", + "format": "file-path", + "description": "Path to a gurobi.lic, mounted into the CoRAL tasks. Required with --coral_solver gurobi_direct." + }, + "aa_data_repo": { + "type": "string", + "format": "directory-path", + "description": "AmpliconArchitect data repository for AmpliconClassifier: the directory holding `GRCh38/` or `CHM13/`, or that reference directory itself. Downloaded automatically for GRCh38; CHM13 has no published repo, so without one the classifier is skipped with a warning.", + "help_text": "Downloading the GRCh38 repo fetches 1.1 GB and unpacks about 4 GB on every run. Extract `GRCh38.tar.gz` once and pass the directory here to skip both." + }, + "aa_data_repo_md5": { + "type": "string", + "description": "Expected MD5 of the downloaded AmpliconArchitect data repository tarball.", + "fa_icon": "fas fa-fingerprint", + "pattern": "^[0-9a-fA-F]{32}$", + "help_text": "The GRCh38 default carries the MD5 the host publishes, and a tarball that does not match fails the run, so the repository cannot change silently. A local --aa_data_repo takes no MD5, and setting --aa_data_repo drops the default." + } + } } }, "allOf": [ @@ -797,6 +888,9 @@ { "$ref": "#/$defs/report_options" }, + { + "$ref": "#/$defs/coral_options" + }, { "$ref": "#/$defs/skip_options" }, diff --git a/nf-test.config b/nf-test.config index d1220b10..f37973a1 100644 --- a/nf-test.config +++ b/nf-test.config @@ -8,6 +8,9 @@ config { // location of an optional nextflow.config file specific for executing tests configFile "tests/nextflow.config" + // Groovy helpers shared by the pipeline tests + libDir "tests/lib" + // ignore tests coming from the nf-core/modules repo ignore = [ 'modules/nf-core/**/tests/*', diff --git a/subworkflows/local/ecdna.nf b/subworkflows/local/ecdna.nf new file mode 100644 index 00000000..2c39d5af --- /dev/null +++ b/subworkflows/local/ecdna.nf @@ -0,0 +1,154 @@ +// +// ecDNA and focal amplification: CoRAL reconstruction, AmpliconClassifier classification +// + +// IMPORT MODULES +include { ASCAT_TO_CORAL_BED } from '../../modules/local/ascattocoralbed/main' +include { CORAL_SEED } from '../../modules/local/coral/seed/main' +include { CORAL_RECONSTRUCT } from '../../modules/local/coral/reconstruct/main' +include { CORAL_CYCLE } from '../../modules/local/coral/cycle/main' +include { CORAL_PLOT } from '../../modules/local/coral/plot/main' +include { AMPLICONCLASSIFIER } from '../../modules/local/ampliconclassifier/main' + +workflow ECDNA { + + take: + tumor_bam // [meta, bam, bai] -- tumour BAMs only + ascat_cnvs // [meta, cnvs_txt] -- ASCAT.out.cnvs + fai // [[:], fai] -- value channel, reused by every sample + data_repo // [[:], dir] -- value channel; empty with --skip_ampliconclassifier or no repo + coral_ref // val 'hg38' | 't2t' + ac_ref // val 'GRCh38' | 'CHM13' + + main: + // + // MODULE: ASCAT_TO_CORAL_BED (label: process_single) + // ASCAT writes cnvs.txt as chr/startpos/endpos/nMajor/nMinor; CoRAL wants + // headerless BED with total CN last, spelled like the BAM's contigs. + // Input: [meta, cnvs_txt], [[:], fai] + // Output: .bed -- [meta, bed] + // + ASCAT_TO_CORAL_BED ( + ascat_cnvs, + fai + ) + + // + // MODULE: CORAL_SEED (label: process_low) + // Input: [meta, cn_seg_bed, bam, bai] + // Output: .seeds -- [meta, bed] -- amplified intervals above --gain + // + // multiMap: seed, reconstruct and plot each take the same cn_seg/bam in their own shape + ASCAT_TO_CORAL_BED.out.bed + .join(tumor_bam, failOnMismatch: true, failOnDuplicate: true) + .multiMap { meta, cn_seg, bam, bai -> + seed: [meta, cn_seg, bam, bai] + reconstruct: [meta, cn_seg, bam, bai] + plot: [meta, bam, bai] + } + .set { coral_inputs } + // coral_inputs.seed / .reconstruct: [meta, cn_seg_bed, bam, bai]; .plot: [meta, bam, bai] + + CORAL_SEED ( + coral_inputs.seed, + coral_ref + ) + + // + // A sample with no amplification above --gain yields an empty seed BED, and + // CoRAL reconstruct errors on one. Size, not countLines(), so the file is + // not staged just to be measured. + // + CORAL_SEED.out.seeds + .branch { _meta, bed -> + seeded: bed.size() > 0 + unseeded: true + } + .set { branched_seeds } + + branched_seeds.unseeded + .subscribe { meta, _bed -> log.info("No amplified intervals found for ${meta.id}: skipping ecDNA reconstruction.") } + + // + // MODULE: CORAL_RECONSTRUCT (label: process_medium) + // Input: [meta, seeds, cn_seg_bed, bam, bai] + // Output: .reconstruction -- [meta, dir] -- graph/cycles/summary, named as AC expects + // + // remainder: false and no failOnMismatch: an unseeded sample is dropped here by design + branched_seeds.seeded + .join(coral_inputs.reconstruct, failOnDuplicate: true) + .set { coral_reconstruct_input } + // coral_reconstruct_input: [meta, seeds, cn_seg_bed, bam, bai] + + CORAL_RECONSTRUCT ( + coral_reconstruct_input + ) + + // errorStrategy 'ignore' drops a failed sample silently, so name it here + branched_seeds.seeded + .map { meta, _seeds -> [meta] } + .join(CORAL_RECONSTRUCT.out.reconstruction, remainder: true) + .filter { _meta, dir -> dir == null } + .subscribe { meta, _dir -> log.warn("CoRAL reconstruct failed for ${meta.id}: no ecDNA results for this sample.") } + + // + // MODULE: CORAL_CYCLE (label: process_low) -- opt-in cycle re-extraction + // Input: [meta, reconstruction_dir] + // Output: .reconstruction -- [meta, dir] -- re-extracted cycles beside the copied graphs + // + if (params.coral_run_cycle) { + CORAL_CYCLE ( + CORAL_RECONSTRUCT.out.reconstruction + ) + + CORAL_RECONSTRUCT.out.reconstruction + .map { meta, _dir -> [meta] } + .join(CORAL_CYCLE.out.reconstruction, remainder: true) + .filter { _meta, dir -> dir == null } + .subscribe { meta, _dir -> log.warn("CoRAL cycle_all failed for ${meta.id}: it will not be classified.") } + ch_for_classifier = CORAL_CYCLE.out.reconstruction + } + else { + ch_for_classifier = CORAL_RECONSTRUCT.out.reconstruction + } + // ch_for_classifier: [meta, dir] + + // + // MODULE: CORAL_PLOT (label: process_medium) + // Runs beside the classifier rather than in front of it: plotting is cosmetic + // and must never gate classification. Plots the same cycles the classifier reads. + // + ch_plots = channel.empty() + if (params.coral_plot) { + ch_for_classifier + .join(coral_inputs.plot, failOnDuplicate: true) + .set { coral_plot_input } + // coral_plot_input: [meta, reconstruction_dir, bam, bai] + + CORAL_PLOT ( + coral_plot_input, + coral_ref + ) + ch_plots = CORAL_PLOT.out.plots + } + + // + // MODULE: AMPLICONCLASSIFIER (label: process_medium) + // Input: [meta, reconstruction_dir], [[:], data_repo], ac_ref + // Output: .classification -- [meta, tsv] -- amplicon_classification_profiles.tsv + // + ch_classification = channel.empty() + if (!params.skip_ampliconclassifier) { + AMPLICONCLASSIFIER ( + ch_for_classifier, + data_repo, + ac_ref + ) + ch_classification = AMPLICONCLASSIFIER.out.classification + } + + emit: + classification = ch_classification // [meta, tsv] + reconstruction = CORAL_RECONSTRUCT.out.reconstruction // [meta, dir] + plots = ch_plots // [meta, [files]] +} diff --git a/subworkflows/local/paired/paired_smallvar_germline.nf b/subworkflows/local/paired/paired_smallvar_germline.nf index 3b473006..3164e33b 100644 --- a/subworkflows/local/paired/paired_smallvar_germline.nf +++ b/subworkflows/local/paired/paired_smallvar_germline.nf @@ -1,4 +1,6 @@ // IMPORT MODULES +include { BCFTOOLS_VIEW as CLAIR3_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' +include { BCFTOOLS_VIEW as DEEPVARIANT_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' include { CLAIR3 } from '../../../modules/local/clair3/main.nf' // IMPORT SUBWORKFLOWS @@ -73,8 +75,15 @@ workflow PAIRED_SMALLVAR_GERMLINE { fai ) - CLAIR3.out.vcf - .join(CLAIR3.out.tbi) + // PASS-only copy for downstream steps; published VCFs are untouched. + def clair3_vcf = CLAIR3.out.vcf.join(CLAIR3.out.tbi) + if (params.smallvar_filter_pass) { + CLAIR3_PASS_FILTER ( clair3_vcf, [], [], [] ) + clair3_vcf = CLAIR3_PASS_FILTER.out.vcf + .join(CLAIR3_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + clair3_vcf .map { meta, vcf , tbi -> def new_meta = meta + [caller:'clair3'] return [new_meta, vcf, tbi] @@ -119,8 +128,15 @@ workflow PAIRED_SMALLVAR_GERMLINE { [[:],[]] // GFF annotation (not used) ) - DEEPVARIANT.out.vcf - .join(DEEPVARIANT.out.vcf_index) + // PASS-only copy for downstream steps; published VCFs are untouched. + def deepvariant_vcf = DEEPVARIANT.out.vcf.join(DEEPVARIANT.out.vcf_index) + if (params.smallvar_filter_pass) { + DEEPVARIANT_PASS_FILTER ( deepvariant_vcf, [], [], [] ) + deepvariant_vcf = DEEPVARIANT_PASS_FILTER.out.vcf + .join(DEEPVARIANT_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + deepvariant_vcf .map{ meta, vcf, tbi -> def new_meta = meta + [caller:'deepvariant'] return [new_meta, vcf, tbi] diff --git a/subworkflows/local/paired/paired_smallvar_somatic.nf b/subworkflows/local/paired/paired_smallvar_somatic.nf index cf5a749d..ab296b99 100644 --- a/subworkflows/local/paired/paired_smallvar_somatic.nf +++ b/subworkflows/local/paired/paired_smallvar_somatic.nf @@ -1,4 +1,6 @@ // IMPORT MODULES +include { BCFTOOLS_VIEW as CLAIRS_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' +include { BCFTOOLS_VIEW as DEEPSOMATIC_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' include { CLAIRS } from '../../../modules/local/clairs/main.nf' include { BCFTOOLS_CONCAT } from '../../../modules/nf-core/bcftools/concat' include { BCFTOOLS_SORT } from '../../../modules/nf-core/bcftools/sort' @@ -71,8 +73,15 @@ workflow PAIRED_SMALLVAR_SOMATIC { BCFTOOLS_CONCAT.out.vcf ) - BCFTOOLS_SORT.out.vcf - .join(BCFTOOLS_SORT.out.tbi) + // PASS-only copy for downstream steps; published VCFs are untouched. + def clairs_vcf = BCFTOOLS_SORT.out.vcf.join(BCFTOOLS_SORT.out.tbi) + if (params.smallvar_filter_pass) { + CLAIRS_PASS_FILTER ( clairs_vcf, [], [], [] ) + clairs_vcf = CLAIRS_PASS_FILTER.out.vcf + .join(CLAIRS_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + clairs_vcf .map { meta, vcf , tbi -> def new_meta = meta + [caller:'clairs'] return [new_meta, vcf, tbi] @@ -109,8 +118,15 @@ workflow PAIRED_SMALLVAR_SOMATIC { ds_pon_channel ) - DEEPSOMATIC.out.vcf - .join(DEEPSOMATIC.out.vcf_index) + // PASS-only copy for downstream steps; published VCFs are untouched. + def deepsomatic_vcf = DEEPSOMATIC.out.vcf.join(DEEPSOMATIC.out.vcf_index) + if (params.smallvar_filter_pass) { + DEEPSOMATIC_PASS_FILTER ( deepsomatic_vcf, [], [], [] ) + deepsomatic_vcf = DEEPSOMATIC_PASS_FILTER.out.vcf + .join(DEEPSOMATIC_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + deepsomatic_vcf .map{ meta, vcf, tbi -> def new_meta = meta + [caller:'deepsomatic'] return [new_meta, vcf, tbi] diff --git a/subworkflows/local/phasing_haplotyping.nf b/subworkflows/local/phasing_haplotyping.nf index e80561c8..f268eb8f 100644 --- a/subworkflows/local/phasing_haplotyping.nf +++ b/subworkflows/local/phasing_haplotyping.nf @@ -6,8 +6,15 @@ include { LONGPHASE_MODCALL as LONGPHASE_MODCALL_GERMLINE } from '../../module include { LONGPHASE_MODCALL as LONGPHASE_MODCALL_SOMATIC } from '../../modules/local/longphase/modcall/main.nf' include { SAMTOOLS_INDEX } from '../../modules/nf-core/samtools/index/main.nf' include { BCFTOOLS_CONCAT } from '../../modules/nf-core/bcftools/concat/main' +include { BCFTOOLS_CONCAT as CONCAT_SOMATIC_UNPHASED } from '../../modules/nf-core/bcftools/concat/main' include { BCFTOOLS_SORT } from '../../modules/nf-core/bcftools/sort/main' +include { BCFTOOLS_SORT as SORT_SOMATIC_PHASED } from '../../modules/nf-core/bcftools/sort/main' include { BCFTOOLS_VIEW } from '../../modules/local/bcftools/view/main.nf' +include { BCFTOOLS_VIEW as SOMATIC_ALT } from '../../modules/local/bcftools/view/main.nf' +include { BCFTOOLS_VIEW as SOMATIC_NONALT } from '../../modules/local/bcftools/view/main.nf' +include { BCFTOOLS_EXCLUDE_SITES as GERMLINE_ANCHORS } from '../../modules/local/bcftools/excludesites/main.nf' +include { VCFTAG as TAG_SOMATIC } from '../../modules/local/vcftag/main.nf' +include { VCFTAG as TAG_GERMLINE } from '../../modules/local/vcftag/main.nf' workflow PHASING_HAPLOTYPING { @@ -138,16 +145,55 @@ workflow PHASING_HAPLOTYPING { } + // + // MODULE: VCFTAG (label: process_single), aliased TAG_SOMATIC / TAG_GERMLINE + // Stamp each arm with an INFO provenance flag before the merge; LongPhase keeps it through phasing. + // + TAG_SOMATIC ( somatic_vcf, 'SOMATIC' ) + TAG_GERMLINE( germline_vcf, 'GERMLINE' ) + + TAG_SOMATIC.out.vcf + .join(TAG_SOMATIC.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .set{ tagged_somatic_vcf } + TAG_GERMLINE.out.vcf + .join(TAG_GERMLINE.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .set{ tagged_germline_vcf } + // tagged_*_vcf: [meta, vcf, tbi] + + // LongPhase phases by position, so a germline record at a somatic call's POS would lend it its phase. + // Only alt-genotype somatic records are phased; 0/0 and ./. records rejoin unphased after the split. + SOMATIC_ALT ( tagged_somatic_vcf ) + SOMATIC_NONALT( tagged_somatic_vcf ) + + SOMATIC_ALT.out.vcf + .join(SOMATIC_ALT.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .set{ somatic_alt_vcf } + SOMATIC_NONALT.out.vcf + .join(SOMATIC_NONALT.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .set{ somatic_nonalt_vcf } + // somatic_alt_vcf / somatic_nonalt_vcf: [meta, vcf, tbi] + + // + // MODULE: GERMLINE_ANCHORS (BCFTOOLS_EXCLUDE_SITES alias) -- germline records not at an alt somatic POS + // Input: [meta, germline_vcf, tbi, somatic_alt_vcf, tbi] + // Output: .vcf / .tbi -- [meta, vcf.gz] / [meta, tbi] + // + GERMLINE_ANCHORS ( + tagged_germline_vcf.join(somatic_alt_vcf, failOnMismatch: true, failOnDuplicate: true) + ) + // Somatic phasing needs germline and somatic sites in one VCF for consistent phase blocks - germline_vcf - .join(somatic_vcf) + GERMLINE_ANCHORS.out.vcf + .join(GERMLINE_ANCHORS.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .join(somatic_alt_vcf) .map { meta, germ_vcf, germ_tbi, som_vcf, som_tbi -> - def vcfs = [som_vcf, germ_vcf] // somatic first (higher priority in phasing) + // Order is cosmetic: BCFTOOLS_CONCAT sorts its inputs by name. + def vcfs = [som_vcf, germ_vcf] def tbis = [som_tbi, germ_tbi] return [ meta, vcfs, tbis] } .set{germline_somatic_vcfs} - // germline_somatic_vcfs (pre-concat): [meta, [somatic_vcf, germline_vcf], [somatic_tbi, germline_tbi]] + // germline_somatic_vcfs (pre-concat): [meta, [somatic_alt_vcf, germline_anchor_vcf], [tbis...]] // // MODULE: BCFTOOLS_CONCAT (label: process_medium) @@ -174,7 +220,7 @@ workflow PHASING_HAPLOTYPING { if (!params.skip_modcall) { // With modcall: include base-modification VCF as additional phasing evidence normal_bams_w_tumoronly_ch - .join(germline_vcf) + .join(tagged_germline_vcf) .join(LONGPHASE_MODCALL_GERMLINE.out.mod_vcf) .map { meta, bam, bai, vcf, _tbi, mods-> def svs = [] // SVs for phasing are not used here @@ -196,7 +242,7 @@ workflow PHASING_HAPLOTYPING { else { // Without modcall: empty lists for SVs and mods normal_bams_w_tumoronly_ch - .join(germline_vcf) + .join(tagged_germline_vcf) .map { meta, bam, bai, vcf, _tbi -> def svs = [] def mods = [] @@ -253,23 +299,29 @@ workflow PHASING_HAPLOTYPING { // phased_somatic_germline_vcf: [meta, vcf, tbi] -- Longphase-phased somatic+germline VCF (unfiltered) // - // MODULE: BCFTOOLS_VIEW (label: process_medium) -- back to somatic-only, with the original somatic VCF as -T targets; PS/HP tags survive - // Input: [meta, phased_combined_vcf, phased_combined_tbi, somatic_vcf, somatic_tbi] + // MODULE: BCFTOOLS_VIEW (label: process_medium) + // Keep the somatic arm by its INFO/SOMATIC flag, not by position; PS/HP tags survive. + // Input: [meta, phased_combined_vcf, phased_combined_tbi] // Output: .vcf -- [meta, vcf.gz] -- phased somatic-only VCF // .tbi -- [meta, tbi] // - phased_somatic_germline_vcf - .join(somatic_vcf) - .map { meta, phased_vcf, phased_tbi, som_vcf, som_tbi -> - return [ meta, phased_vcf, phased_tbi, som_vcf, som_tbi ] + BCFTOOLS_VIEW ( phased_somatic_germline_vcf ) + + // Add the unphased 0/0 and ./. somatic records back (only present with --smallvar_filter_pass false) + BCFTOOLS_VIEW.out.vcf + .join(BCFTOOLS_VIEW.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .join(somatic_nonalt_vcf, failOnMismatch: true, failOnDuplicate: true) + .map { meta, phased_vcf, phased_tbi, nonalt_vcf, nonalt_tbi -> + return [ meta, [phased_vcf, nonalt_vcf], [phased_tbi, nonalt_tbi] ] } - .set { bcftools_view_input_ch } - // bcftools_view_input_ch: [meta, phased_combined_vcf, tbi, somatic_vcf, somatic_tbi] + .set{ somatic_parts } + // somatic_parts: [meta, [phased_alt_vcf, nonalt_vcf], [tbis...]] - BCFTOOLS_VIEW ( bcftools_view_input_ch ) + CONCAT_SOMATIC_UNPHASED ( somatic_parts ) + SORT_SOMATIC_PHASED ( CONCAT_SOMATIC_UNPHASED.out.vcf ) - BCFTOOLS_VIEW.out.vcf - .join(BCFTOOLS_VIEW.out.tbi) + SORT_SOMATIC_PHASED.out.vcf + .join(SORT_SOMATIC_PHASED.out.tbi, failOnMismatch: true, failOnDuplicate: true) .set{ phased_somatic_vcf } // phased_somatic_vcf: [meta, vcf.gz, tbi] -- phased somatic-only VCF (germline removed) diff --git a/subworkflows/local/prepare_aa_data_repo.nf b/subworkflows/local/prepare_aa_data_repo.nf new file mode 100644 index 00000000..69386299 --- /dev/null +++ b/subworkflows/local/prepare_aa_data_repo.nf @@ -0,0 +1,56 @@ +// +// Stage the AmpliconArchitect data repository AmpliconClassifier reads at runtime +// + +include { AADATAREPO_DOWNLOAD } from '../../modules/local/aadatarepo/download/main' +include { UNTAR as UNTAR_AA_DATA_REPO } from '../../modules/nf-core/untar/main' + +workflow PREPARE_AA_DATA_REPO { + + take: + data_repo // path or null -- params.aa_data_repo + repo_url // URL or null -- genome attribute aa_data_repo_url + repo_md5 // MD5 or null -- --aa_data_repo_md5, else the genome attribute + + main: + ch_versions = channel.empty() + + // The MD5 is checked by the download task, so a local repo would silently skip it + if (data_repo && params.aa_data_repo_md5) { + error("--aa_data_repo_md5: only checks a repository the pipeline downloads. Drop it when --aa_data_repo is set.") + } + + if (data_repo) { + ch_data_repo = channel.value([ [ id: 'aa_data_repo' ], file(data_repo, checkIfExists: true) ]) + } + else if (repo_url) { + // + // MODULES: AADATAREPO_DOWNLOAD -> UNTAR_AA_DATA_REPO (labels: process_single) + // ~1.1 GB tarball; the plain build, not GRCh38_indexed, whose extra BWA index AC never reads + // + if (!repo_md5) { + log.warn("AmpliconClassifier: the data repository '${repo_url}' is downloaded without --aa_data_repo_md5, so the repository is not verified.") + } + AADATAREPO_DOWNLOAD ( + channel.value([ [ id: 'aa_data_repo' ], repo_url, repo_md5 ]) + ) + + UNTAR_AA_DATA_REPO ( + AADATAREPO_DOWNLOAD.out.archive + ) + + // .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() + ch_versions = ch_versions.mix(AADATAREPO_DOWNLOAD.out.versions, UNTAR_AA_DATA_REPO.out.versions) + } + else { + // An empty repo runs no classifier tasks; CoRAL reconstruction is unaffected + log.warn("No AmpliconArchitect data repository is published for ${params.genome}: skipping AmpliconClassifier. Set --aa_data_repo to classify.") + ch_data_repo = channel.empty() + } + + emit: + data_repo = ch_data_repo // [[id:'aa_data_repo'], dir], or empty when no repo is available + versions = ch_versions +} diff --git a/subworkflows/local/small_variant_consensus.nf b/subworkflows/local/small_variant_consensus.nf index 7e541aaa..aaedb58c 100644 --- a/subworkflows/local/small_variant_consensus.nf +++ b/subworkflows/local/small_variant_consensus.nf @@ -1,12 +1,13 @@ include { BCFTOOLS_NORM } from '../../modules/nf-core/bcftools/norm/main' +include { BCFTOOLS_NORM as BCFTOOLS_NORM_REJOIN } from '../../modules/nf-core/bcftools/norm/main' include { BCFTOOLS_ISEC } from '../../modules/nf-core/bcftools/isec/main' include { BCFTOOLS_QUERY } from '../../modules/nf-core/bcftools/query/main' include { BCFTOOLS_ANNOTATE } from '../../modules/nf-core/bcftools/annotate/main' -include { BCFTOOLS_ANNOTATE as STANDARDIZE_AF } from '../../modules/nf-core/bcftools/annotate/main' include { BCFTOOLS_CONCAT } from '../../modules/nf-core/bcftools/concat/main' include { BCFTOOLS_SORT } from '../../modules/nf-core/bcftools/sort/main' include { BCFTOOLS_SORT as SORT_POST_NORM } from '../../modules/nf-core/bcftools/sort/main' include { BCFTOOLS_SORT as BCFTOOLS_SORT_CONSENSUS } from '../../modules/nf-core/bcftools/sort/main' +include { BCFTOOLS_EXCLUDE_SITES } from '../../modules/local/bcftools/excludesites/main' @@ -17,16 +18,26 @@ workflow SMALL_VARIANT_CONSENSUS { fasta // [[:], fasta] _fai // [[:], fai] prioritize_caller // str: which caller's calls take priority ('deepvariant'/'deepsomatic' or 'clair') - combine_method // str: 'consensus' (intersection only) or 'all' (intersection + private calls from priority caller) + combine_method // str: 'consensus' (shared calls only) or 'all' (union of both callers' calls) main: + if (!(combine_method in ['consensus', 'all'])) { + error("combine_method must be 'consensus' or 'all', got '${combine_method}'") + } + if (!(prioritize_caller in ['deepvariant', 'deepsomatic', 'clair'])) { + error("prioritize_caller must be one of [deepvariant, deepsomatic, clair], got '${prioritize_caller}'") + } + // - // MODULE: BCFTOOLS_NORM (label: process_medium) -- left-align and normalise; sorted after, since left-alignment can reorder records - // Input: [meta, vcf, tbi] -- per-caller VCF + // MODULE: BCFTOOLS_NORM (label: process_medium) -- left-align; in 'consensus' mode also split multi-allelics for isec + // Input: [meta(+split), vcf, tbi] -- per-caller VCF // Output: .vcf -- [meta, vcf] -- left-aligned, normalised VCF (unsorted) // - BCFTOOLS_NORM(mixed_vcfs, fasta) + BCFTOOLS_NORM( + mixed_vcfs.map { meta, vcf, tbi -> [meta + [split: combine_method == 'consensus'], vcf, tbi] }, + fasta + ) // // MODULE: SORT_POST_NORM (BCFTOOLS_SORT alias, label: process_medium) -- re-sort and index after normalisation @@ -42,31 +53,10 @@ workflow SMALL_VARIANT_CONSENSUS { // normalized_vcfs: [meta(+caller), vcf.gz, tbi] -- normalised, sorted per-caller VCF // - // MODULE: STANDARDIZE_AF (BCFTOOLS_ANNOTATE alias, label: process_low) -- rename the AF FORMAT field to the priority caller's: + // ALLELE FREQUENCY KEY -- BCFTOOLS_ANNOTATE below renames the AF FORMAT field to the priority caller's: // FORMAT/AF -> FORMAT/VAF when prioritize_caller is 'deepvariant'/'deepsomatic' // FORMAT/VAF -> FORMAT/AF when prioritize_caller is 'clair' - // - if (combine_method == 'all') { - normalized_vcfs - .map { meta, vcf, tbi -> - def rename_to = prioritize_caller in ['deepvariant', 'deepsomatic'] ? 'VAF' : 'AF' - def new_meta = meta + [rename_to: rename_to] - return [new_meta, vcf, tbi, [], [], [], [], []] - } - .set { standardize_input } - - STANDARDIZE_AF(standardize_input) - - STANDARDIZE_AF.out.vcf - .join(STANDARDIZE_AF.out.tbi) - .map { meta, vcf, tbi -> - def clean_meta = meta.findAll { k, _v -> k != 'rename_to' } - return [clean_meta, vcf, tbi] - } - .set { normalized_vcfs } - // normalized_vcfs: [meta(+caller), vcf, tbi] -- normalised, AF-standardized per-caller VCF - } - // In 'consensus' mode, normalized_vcfs comes from SORT_POST_NORM (post-BCFTOOLS_NORM re-sorting) + // Only 'all' mode renames: it merges both callers, so the merged VCF needs one AF key for WAKHAN. // // MODULE: BCFTOOLS_QUERY (label: process_single) @@ -79,13 +69,17 @@ workflow SMALL_VARIANT_CONSENSUS { // Prepare BCFTOOLS_ANNOTATE input: VCF + caller-name annotation file normalized_vcfs - .join(BCFTOOLS_QUERY.out.output) - .join(BCFTOOLS_QUERY.out.index) + .join(BCFTOOLS_QUERY.out.output, failOnMismatch: true, failOnDuplicate: true) + .join(BCFTOOLS_QUERY.out.index, failOnMismatch: true, failOnDuplicate: true) .map{ meta, vcf, tbi, annotations, annotations_index -> def columns = [] // no extra column specs def header_lines = [] // no extra header lines def rename_chrs = [] // no chromosome renaming - return [ meta, vcf, tbi, annotations, annotations_index, columns, header_lines, rename_chrs ] + // 'all' mode merges both callers, so unify the AF key; 'consensus' needs no rename. + def new_meta = combine_method == 'all' + ? meta + [rename_to: (prioritize_caller in ['deepvariant', 'deepsomatic'] ? 'VAF' : 'AF')] + : meta + return [ new_meta, vcf, tbi, annotations, annotations_index, columns, header_lines, rename_chrs ] } .set{annotate_input} // annotate_input: [meta, vcf, tbi, annotations_tsv, annotations_tbi, [], [], []] @@ -100,17 +94,28 @@ workflow SMALL_VARIANT_CONSENSUS { BCFTOOLS_ANNOTATE(annotate_input) BCFTOOLS_ANNOTATE.out.vcf - .join(BCFTOOLS_ANNOTATE.out.tbi) + .join(BCFTOOLS_ANNOTATE.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .map { meta, vcf, tbi -> + def clean_meta = meta.findAll { k, _v -> !(k in ['rename_to', 'split']) } + return [clean_meta, vcf, tbi] + } .set{annotated_vcfs} // annotated_vcfs: [meta(+caller), vcf, tbi] -- VCF with CALLER INFO tag // Branch annotated VCFs by caller family for the intersection step + // `other` errors on an unrecognised meta.caller instead of silently dropping the sample. annotated_vcfs .branch { meta, _vcfs, _tbi -> deepvariant: meta.caller in [ 'deepvariant', 'deepsomatic' ] clair: meta.caller in ['clair3','clairs-to','clairs'] + other: true } .set{annotated_vcfs_branched} + + annotated_vcfs_branched.other + .map { meta, _vcfs, _tbi -> + error("SMALL_VARIANT_CONSENSUS: unrecognised meta.caller '${meta.caller}' for sample '${meta.id}'; expected one of [deepvariant, deepsomatic, clair3, clairs-to, clairs]") + } // annotated_vcfs_branched.deepvariant: [meta(caller=deepvariant/deepsomatic), vcf, tbi] // annotated_vcfs_branched.clair: [meta(caller=clair3/clairs-to/clairs), vcf, tbi] @@ -153,8 +158,9 @@ workflow SMALL_VARIANT_CONSENSUS { // deepvariant_ch: [meta (no caller), vcf, tbi] // Join DeepVariant and Clair VCFs per sample into a single tuple for BCFTOOLS_ISEC + // failOnMismatch: a sample missing one caller would otherwise be dropped silently. deepvariant_ch - .join(clair_ch) + .join(clair_ch, failOnMismatch: true, failOnDuplicate: true) .map { meta, deepvar_vcf, deepvar_tbi, clair_vcf, clair_tbi -> def vcfs = [deepvar_vcf, clair_vcf] def tbis = [deepvar_tbi, clair_tbi] @@ -163,87 +169,85 @@ workflow SMALL_VARIANT_CONSENSUS { .set{mixed_vcfs} // mixed_vcfs (re-paired): [meta, [deepvar_vcf, clair_vcf], [deepvar_tbi, clair_tbi]] - // Add empty optional fields required by BCFTOOLS_ISEC - mixed_vcfs - .map{ meta, vcfs, tbis -> - def file = [] // no regions file - def target = [] // no target sites - def regions = [] // no region string - return [meta, vcfs, tbis, file, target, regions] - } - .set{isec_input} - // isec_input: [meta, [deepvar_vcf, clair_vcf], [deepvar_tbi, clair_tbi], [], [], []] - - // - // MODULE: BCFTOOLS_ISEC (label: process_medium) -- shared and private sets of the two callers - // Input: [meta, [vcf1, vcf2], [tbi1, tbi2], [], [], []] - // Output: .deepvar_consensus_vcf / .clair_consensus_vcf -- [meta, vcf] -- shared calls, DeepVariant or Clair record - // .deepvar_private_vcf / .clair_private_vcf -- [meta, vcf] -- caller-private calls (+ .tbi for each) - // - BCFTOOLS_ISEC(isec_input) - if (combine_method == 'consensus') { - // Take only the intersection: variants called by BOTH callers - // Use the record from the prioritized caller - if (prioritize_caller in ['deepvariant', 'deepsomatic']) { - BCFTOOLS_ISEC.out.deepvar_consensus_vcf - .set{isec_consensus_vcf} - } - else if (prioritize_caller == 'clair') { - BCFTOOLS_ISEC.out.clair_consensus_vcf - .set{isec_consensus_vcf} - } + // Add empty optional fields required by BCFTOOLS_ISEC + mixed_vcfs + .map{ meta, vcfs, tbis -> + def file = [] // no regions file + def target = [] // no target sites + def regions = [] // no region string + return [meta, vcfs, tbis, file, target, regions] + } + .set{isec_input} + // isec_input: [meta, [deepvar_vcf, clair_vcf], [deepvar_tbi, clair_tbi], [], [], []] + + // + // MODULE: BCFTOOLS_ISEC (label: process_medium) -- shared and private sets of the two callers + // Input: [meta, [vcf1, vcf2], [tbi1, tbi2], [], [], []] + // Output: .deepvar_consensus_vcf / .clair_consensus_vcf -- [meta, vcf] -- shared calls, DeepVariant or Clair record + // + BCFTOOLS_ISEC(isec_input) + + // Take only the intersection: variants called by BOTH callers, from the prioritized caller's record + def isec_consensus_vcf = prioritize_caller in ['deepvariant', 'deepsomatic'] + ? BCFTOOLS_ISEC.out.deepvar_consensus_vcf + : BCFTOOLS_ISEC.out.clair_consensus_vcf // ISEC always writes 0002.vcf.gz, so germline and somatic would collide by basename in BCFTOOLS_CONCAT; // BCFTOOLS_SORT_CONSENSUS renames it per sample (conf/modules.config) BCFTOOLS_SORT_CONSENSUS(isec_consensus_vcf) - BCFTOOLS_SORT_CONSENSUS.out.vcf.set{vcf} - BCFTOOLS_SORT_CONSENSUS.out.tbi.set{tbi} - // vcf/tbi: [meta, vcf/tbi] -- consensus-only calls from the priority caller, renamed - } - else if (combine_method == 'all') { - // Intersection plus the priority caller's private calls - if (prioritize_caller in ['deepvariant', 'deepsomatic']) { - // consensus (DeepVariant record) + DeepVariant-private variants - BCFTOOLS_ISEC.out.deepvar_consensus_vcf - .join(BCFTOOLS_ISEC.out.deepvar_consensus_tbi) - .join(BCFTOOLS_ISEC.out.clair_private_vcf) - .join(BCFTOOLS_ISEC.out.clair_private_tbi) - .map{ meta, deepvar_vcf, deepvar_tbi, clair_vcf, clair_tbi -> - return[meta, [deepvar_vcf, clair_vcf], [deepvar_tbi, clair_tbi]] - } - .set{concat_input} - // concat_input: [meta, [consensus_vcf, private_vcf], [consensus_tbi, private_tbi]] - BCFTOOLS_CONCAT(concat_input) - BCFTOOLS_CONCAT.out.vcf - .set{concat_out} - } - else if (prioritize_caller == 'clair') { - // consensus (Clair record) + Clair-private variants - BCFTOOLS_ISEC.out.deepvar_private_vcf - .join(BCFTOOLS_ISEC.out.deepvar_private_tbi) - .join(BCFTOOLS_ISEC.out.clair_consensus_vcf) - .join(BCFTOOLS_ISEC.out.clair_consensus_tbi) - .map{ meta, deepvar_vcf, deepvar_tbi, clair_vcf, clair_tbi -> - return[meta, [deepvar_vcf, clair_vcf], [deepvar_tbi, clair_tbi]] - } - .set{concat_input} - // concat_input: [meta, [private_vcf, consensus_vcf], [private_tbi, consensus_tbi]] - BCFTOOLS_CONCAT(concat_input) - BCFTOOLS_CONCAT.out.vcf - .set{concat_out} - } - // concat_out: [meta, vcf] -- unsorted concatenated VCF (consensus + priority-caller-private) - BCFTOOLS_SORT(concat_out) - BCFTOOLS_SORT.out.vcf - .set{vcf} - BCFTOOLS_SORT.out.tbi - .set{tbi} - // vcf/tbi: [meta, vcf/tbi] -- sorted combined VCF + // + // MODULE: BCFTOOLS_NORM_REJOIN (BCFTOOLS_NORM alias) -- rejoin split sites (-m +any) so LongPhase and Wakhan see one record per position + // Input: [meta, vcf, tbi] -- sorted consensus VCF (one caller's records only) + // Output: .vcf -- [meta, vcf.gz] + // .tbi -- [meta, tbi] + // + BCFTOOLS_NORM_REJOIN( + BCFTOOLS_SORT_CONSENSUS.out.vcf.join(BCFTOOLS_SORT_CONSENSUS.out.tbi, failOnMismatch: true, failOnDuplicate: true), + fasta + ) + BCFTOOLS_NORM_REJOIN.out.vcf.set{ vcf } + BCFTOOLS_NORM_REJOIN.out.tbi.set{ tbi } + // vcf/tbi: [meta, vcf/tbi] -- consensus calls from the priority caller, multi-allelics rejoined + } + else { + // Union by locus: the priority caller's records, plus the other caller's at positions it has no record for. + // Records are unsplit here, so every output record is one caller's call. + mixed_vcfs + .multiMap { meta, vcfs, tbis -> + def prio = prioritize_caller in ['deepvariant', 'deepsomatic'] ? 0 : 1 + priority: [meta, vcfs[prio], tbis[prio]] + exclude: [meta, vcfs[1 - prio], tbis[1 - prio], vcfs[prio], tbis[prio]] + } + .set{ by_priority } + // by_priority.priority: [meta, prio_vcf, prio_tbi] + // by_priority.exclude: [meta, other_vcf, other_tbi, prio_vcf, prio_tbi] + + // + // MODULE: BCFTOOLS_EXCLUDE_SITES (label: process_single) -- the other caller's records at positions the priority caller lacks + // Input: [meta, other_vcf, other_tbi, prio_vcf, prio_tbi] + // Output: .vcf / .tbi -- [meta, vcf.gz] / [meta, tbi] + // + BCFTOOLS_EXCLUDE_SITES(by_priority.exclude) + + by_priority.priority + .join(BCFTOOLS_EXCLUDE_SITES.out.vcf, failOnMismatch: true, failOnDuplicate: true) + .join(BCFTOOLS_EXCLUDE_SITES.out.tbi, failOnMismatch: true, failOnDuplicate: true) + .map { meta, prio_vcf, prio_tbi, other_vcf, other_tbi -> + return [meta, [prio_vcf, other_vcf], [prio_tbi, other_tbi]] + } + .set{concat_input} + // concat_input: [meta, [prio_vcf, other_only_vcf], [tbis...]] + + BCFTOOLS_CONCAT(concat_input) + BCFTOOLS_SORT(BCFTOOLS_CONCAT.out.vcf) + BCFTOOLS_SORT.out.vcf.set{ vcf } + BCFTOOLS_SORT.out.tbi.set{ tbi } + // vcf/tbi: [meta, vcf/tbi] -- sorted union VCF, one record per position } emit: - vcf // [meta, vcf] -- final consensus/combined VCF + vcf // [meta, vcf] -- final consensus/combined VCF, one record per position tbi // [meta, tbi] } diff --git a/subworkflows/local/tests/ecdna.nf.test b/subworkflows/local/tests/ecdna.nf.test new file mode 100644 index 00000000..f0553c07 --- /dev/null +++ b/subworkflows/local/tests/ecdna.nf.test @@ -0,0 +1,195 @@ +nextflow_workflow { + + name "Test Subworkflow ECDNA" + script "../ecdna.nf" + workflow "ECDNA" + + tag "subworkflows" + tag "subworkflows_local" + tag "ecdna" + tag "small" + + // Stub-only: CoRAL and AmpliconClassifier need real BAMs, a solver and a ~1 GB data + // repo. What these prove is the wiring -- that the joins pair, that the empty-seed + // branch drops a sample before reconstruct, and that the params and repo gate the optional steps. + test("-stub - reconstructs and classifies a seeded sample") { + + options "-stub" + + when { + params { + coral_run_cycle = false + coral_plot = true + skip_ampliconclassifier = false + } + workflow { + """ + input[0] = channel.of([ [ id:'sample1' ], file('sample1.bam'), file('sample1.bam.bai') ]) + input[1] = channel.of([ [ id:'sample1' ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) ]) + input[3] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures") ]) + input[4] = 't2t' + input[5] = 'CHM13' + """ + } + } + + then { + def names = workflow.trace.tasks().collect { it.name.split(' ')[0] } + + assertAll( + { assert workflow.success }, + { assert names.contains('ECDNA:ASCAT_TO_CORAL_BED') }, + { assert names.contains('ECDNA:CORAL_SEED') }, + { assert names.contains('ECDNA:CORAL_RECONSTRUCT') }, + { assert names.contains('ECDNA:CORAL_PLOT') }, + { assert names.contains('ECDNA:AMPLICONCLASSIFIER') }, + // Cycle re-extraction is opt-in + { assert !names.contains('ECDNA:CORAL_CYCLE') }, + { assert workflow.out.classification.size() == 1 } + ) + } + } + + test("-stub - coral_run_cycle classifies and plots the re-extracted cycles") { + + options "-stub" + + when { + params { + coral_run_cycle = true + coral_plot = true + skip_ampliconclassifier = false + } + workflow { + """ + input[0] = channel.of([ [ id:'sample1' ], file('sample1.bam'), file('sample1.bam.bai') ]) + input[1] = channel.of([ [ id:'sample1' ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) ]) + input[3] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures") ]) + input[4] = 't2t' + input[5] = 'CHM13' + """ + } + } + + then { + def names = workflow.trace.tasks().collect { it.name.split(' ')[0] } + + assertAll( + { assert workflow.success }, + { assert names.contains('ECDNA:CORAL_CYCLE') }, + { assert names.contains('ECDNA:CORAL_PLOT') }, + { assert names.contains('ECDNA:AMPLICONCLASSIFIER') }, + { assert workflow.out.plots.size() == 1 } + ) + } + } + + test("-stub - skip_ampliconclassifier stops after reconstruction") { + + options "-stub" + + when { + params { + coral_run_cycle = false + coral_plot = false + skip_ampliconclassifier = true + } + workflow { + """ + input[0] = channel.of([ [ id:'sample1' ], file('sample1.bam'), file('sample1.bam.bai') ]) + input[1] = channel.of([ [ id:'sample1' ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) ]) + input[3] = channel.empty() + input[4] = 't2t' + input[5] = 'CHM13' + """ + } + } + + then { + def names = workflow.trace.tasks().collect { it.name.split(' ')[0] } + + assertAll( + { assert workflow.success }, + { assert names.contains('ECDNA:CORAL_RECONSTRUCT') }, + { assert !names.contains('ECDNA:AMPLICONCLASSIFIER') } + ) + } + } + + test("-stub - an empty seed skips reconstruction and classification") { + + options "-stub" + + when { + params { + coral_run_cycle = false + coral_plot = true + skip_ampliconclassifier = false + } + workflow { + """ + // stub_empty_seed makes the CORAL_SEED stub write an empty BED; the joins match on the whole meta + input[0] = channel.of([ [ id:'sample1', stub_empty_seed:true ], file('sample1.bam'), file('sample1.bam.bai') ]) + input[1] = channel.of([ [ id:'sample1', stub_empty_seed:true ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) ]) + input[3] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures") ]) + input[4] = 't2t' + input[5] = 'CHM13' + """ + } + } + + then { + def names = workflow.trace.tasks().collect { it.name.split(' ')[0] } + + assertAll( + { assert workflow.success }, + { assert names.contains('ECDNA:CORAL_SEED') }, + { assert !names.contains('ECDNA:CORAL_RECONSTRUCT') }, + { assert !names.contains('ECDNA:CORAL_PLOT') }, + { assert !names.contains('ECDNA:AMPLICONCLASSIFIER') }, + { assert workflow.out.reconstruction.size() == 0 }, + { assert workflow.out.classification.size() == 0 } + ) + } + } + + test("-stub - no data repo reconstructs without classifying") { + + options "-stub" + + when { + params { + coral_run_cycle = false + coral_plot = false + skip_ampliconclassifier = false + } + workflow { + """ + // PREPARE_AA_DATA_REPO emits nothing on CHM13 without --aa_data_repo + input[0] = channel.of([ [ id:'sample1' ], file('sample1.bam'), file('sample1.bam.bai') ]) + input[1] = channel.of([ [ id:'sample1' ], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/sample1.cnvs.txt", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/modules/local/ascattocoralbed/tests/fixtures/chr.fai", checkIfExists: true) ]) + input[3] = channel.empty() + input[4] = 't2t' + input[5] = 'CHM13' + """ + } + } + + then { + def names = workflow.trace.tasks().collect { it.name.split(' ')[0] } + + assertAll( + { assert workflow.success }, + { assert names.contains('ECDNA:CORAL_RECONSTRUCT') }, + { assert !names.contains('ECDNA:AMPLICONCLASSIFIER') }, + { assert workflow.out.reconstruction.size() == 1 }, + { assert workflow.out.classification.size() == 0 } + ) + } + } +} diff --git a/subworkflows/local/tumor_only/tumoronly_smallvar.nf b/subworkflows/local/tumor_only/tumoronly_smallvar.nf index b37beaf2..ab726ceb 100644 --- a/subworkflows/local/tumor_only/tumoronly_smallvar.nf +++ b/subworkflows/local/tumor_only/tumoronly_smallvar.nf @@ -1,4 +1,6 @@ // IMPORT MODULES +include { BCFTOOLS_VIEW as DEEPVARIANT_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' +include { BCFTOOLS_VIEW as DEEPSOMATIC_PASS_FILTER } from '../../../modules/nf-core/bcftools/view/main' include { CLAIRSTO } from '../../../modules/local/clairsto/main.nf' include { CLAIRSTO_VERDICT_TAG } from '../../../modules/local/clairsto/verdict_tag/main.nf' include { VCFSPLIT } from '../../../modules/local/vcfsplit/main.nf' @@ -9,6 +11,11 @@ include { DEEPSOMATIC } from '../../../subwork include { SMALL_VARIANT_CONSENSUS as GERMLINE_CONSENSUS } from '../../../subworkflows/local/small_variant_consensus.nf' include { SMALL_VARIANT_CONSENSUS as SOMATIC_CONSENSUS } from '../../../subworkflows/local/small_variant_consensus.nf' +// Germline verdict transfer: DeepSomatic adjudicates DeepVariant's tumor-derived germline calls. +include { BCFTOOLS_QUERY as DS_VERDICT_QUERY } from '../../../modules/nf-core/bcftools/query/main' +include { BCFTOOLS_ANNOTATE as DS_VERDICT_ANNOTATE } from '../../../modules/nf-core/bcftools/annotate/main' +include { BCFTOOLS_VIEW as DS_GERMLINE_SELECT } from '../../../modules/nf-core/bcftools/view/main' + workflow TUMORONLY_SMALLVAR { @@ -130,6 +137,50 @@ workflow TUMORONLY_SMALLVAR { // clairsto_somatic_ch: [meta(+caller:'clairs-to'), vcf, tbi] -- somatic variants } + // DEEPSOMATIC in tumor-only mode: normal BAM/BAI are empty lists + if(somatic_var_keep.contains('deepsomatic')) { + tumor_bams + .map { meta, tumor_bam, tumor_bai -> + def normal_bam = [] + def normal_bai = [] + return [meta,normal_bam,normal_bai,tumor_bam,tumor_bai] + } + .set{deepsomatic_input_ch} + // deepsomatic_input_ch: [meta, [], [], tumor_bam, tumor_bai] + // empty normal_bam/bai signals tumor-only mode to DEEPSOMATIC subworkflow + + // + // SUBWORKFLOW: DEEPSOMATIC (local) + // Input: [meta, [], [], tumor_bam, tumor_bai] -- tumor-only (no normal) + // [[:],[]] / fasta / fai / [[:],[]] + // Output: .vcf -- [meta, vcf] + // .vcf_index -- [meta, tbi] + // + DEEPSOMATIC ( + deepsomatic_input_ch, + [[:],[]], // intervals (empty = genome-wide) + fasta, + fai, + [[:],[]], // GZI (empty if FASTA is uncompressed) + ds_pon_channel + ) + // PASS-only copy for downstream steps; published VCFs are untouched. + def deepsomatic_vcf = DEEPSOMATIC.out.vcf.join(DEEPSOMATIC.out.vcf_index) + if (params.smallvar_filter_pass) { + DEEPSOMATIC_PASS_FILTER ( deepsomatic_vcf, [], [], [] ) + deepsomatic_vcf = DEEPSOMATIC_PASS_FILTER.out.vcf + .join(DEEPSOMATIC_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + deepsomatic_vcf + .map{ meta, vcf, tbi -> + def new_meta = meta + [caller:'deepsomatic'] + return [new_meta, vcf, tbi] + } + .set{deepsomatic_ch} + // deepsomatic_ch: [meta(+caller:'deepsomatic'), vcf, tbi] + } + // DEEPVARIANT: germline-only variant calling (no somatic mode for tumor-only) if(germline_var_keep.contains('deepvariant')) { @@ -156,14 +207,62 @@ workflow TUMORONLY_SMALLVAR { [[:],[]] // GFF annotation (not used) ) - DEEPVARIANT.out.vcf - .join(DEEPVARIANT.out.vcf_index) + // PASS-only copy for downstream steps; published VCFs are untouched. + def deepvariant_vcf = DEEPVARIANT.out.vcf.join(DEEPVARIANT.out.vcf_index) + if (params.smallvar_filter_pass) { + DEEPVARIANT_PASS_FILTER ( deepvariant_vcf, [], [], [] ) + deepvariant_vcf = DEEPVARIANT_PASS_FILTER.out.vcf + .join(DEEPVARIANT_PASS_FILTER.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + // Keep only DeepSomatic-adjudicated germline sites; skipped when deepsomatic isn't selected. + def deepvariant_germline = deepvariant_vcf + if (somatic_var_keep.contains('deepsomatic')) { + // GERMLINE VERDICT TRANSFER: DeepVariant on the tumor BAM cannot tell germline from somatic. + // + // MODULE: DS_VERDICT_QUERY (BCFTOOLS_QUERY alias, label: process_single) + // Input: [meta, deepsomatic_vcf, tbi] -- the RAW DeepSomatic VCF, before its PASS filter + // Output: .output/.index -- [meta, tsv.gz/tbi] -- CHROM POS REF ALT FILTER, non-PASS/RefCall rows only + // + DS_VERDICT_QUERY ( DEEPSOMATIC.out.vcf.join(DEEPSOMATIC.out.vcf_index), [], [], [] ) + + // + // MODULE: DS_VERDICT_ANNOTATE (BCFTOOLS_ANNOTATE alias, label: process_medium) + // Stamps INFO/DS_VERDICT on each DeepVariant record from the DeepSomatic verdict table. + // + deepvariant_vcf + .join(DS_VERDICT_QUERY.out.output, failOnMismatch: true, failOnDuplicate: true) + .join(DS_VERDICT_QUERY.out.index, failOnMismatch: true, failOnDuplicate: true) + .map { meta, vcf, tbi, annotations, annotations_index -> + def columns = [] // no extra column specs + def header_lines = [] // no extra header lines + def rename_chrs = [] // no chromosome renaming + return [ meta, vcf, tbi, annotations, annotations_index, columns, header_lines, rename_chrs ] + } + .set{ ds_verdict_annotate_input } + + DS_VERDICT_ANNOTATE ( ds_verdict_annotate_input ) + + // + // MODULE: DS_GERMLINE_SELECT (BCFTOOLS_VIEW alias, label: process_medium) + // Keeps only the positively-adjudicated germline records (see ext.args in conf/modules.config). + // + DS_GERMLINE_SELECT ( + DS_VERDICT_ANNOTATE.out.vcf.join(DS_VERDICT_ANNOTATE.out.tbi, failOnMismatch: true, failOnDuplicate: true), + [], [], [] + ) + + deepvariant_germline = DS_GERMLINE_SELECT.out.vcf + .join(DS_GERMLINE_SELECT.out.index, failOnMismatch: true, failOnDuplicate: true) + } + + deepvariant_germline .map{ meta, vcf, tbi -> def new_meta = meta + [caller:'deepvariant'] return [new_meta, vcf, tbi] } .set{deepvariant_ch} - // deepvariant_ch: [meta(+caller:'deepvariant'), vcf, tbi] + // deepvariant_ch: [meta(+caller:'deepvariant'), vcf, tbi] -- germline-adjudicated if deepsomatic ran } // COMBINE GERMLINE VARIANTS @@ -196,42 +295,6 @@ workflow TUMORONLY_SMALLVAR { .set{germline_vcf} } - // DEEPSOMATIC in tumor-only mode: normal BAM/BAI are empty lists - if(somatic_var_keep.contains('deepsomatic')) { - tumor_bams - .map { meta, tumor_bam, tumor_bai -> - def normal_bam = [] - def normal_bai = [] - return [meta,normal_bam,normal_bai,tumor_bam,tumor_bai] - } - .set{deepsomatic_input_ch} - // deepsomatic_input_ch: [meta, [], [], tumor_bam, tumor_bai] - // empty normal_bam/bai signals tumor-only mode to DEEPSOMATIC subworkflow - - // - // SUBWORKFLOW: DEEPSOMATIC (local) - // Input: [meta, [], [], tumor_bam, tumor_bai] -- tumor-only (no normal) - // [[:],[]] / fasta / fai / [[:],[]] - // Output: .vcf -- [meta, vcf] - // .vcf_index -- [meta, tbi] - // - DEEPSOMATIC ( - deepsomatic_input_ch, - [[:],[]], // intervals (empty = genome-wide) - fasta, - fai, - [[:],[]], // GZI (empty if FASTA is uncompressed) - ds_pon_channel - ) - DEEPSOMATIC.out.vcf - .join(DEEPSOMATIC.out.vcf_index) - .map{ meta, vcf, tbi -> - def new_meta = meta + [caller:'deepsomatic'] - return [new_meta, vcf, tbi] - } - .set{deepsomatic_ch} - // deepsomatic_ch: [meta(+caller:'deepsomatic'), vcf, tbi] - } // COMBINE SOMATIC VARIATION if (somatic_var_keep.size() > 1) { diff --git a/tests/.nftignore b/tests/.nftignore index dfc63dec..c835ed89 100644 --- a/tests/.nftignore +++ b/tests/.nftignore @@ -28,6 +28,7 @@ pipeline_info/*.{html,json,txt,yml} */variants/deepsomatic/*.{vcf.gz,vcf.gz.tbi} */variants/deepvariant/*.{vcf.gz,vcf.gz.tbi} */variants/savana/* +*/ecdna/* */signatures/matrices/** */signatures/assignment/** */report/*.html diff --git a/tests/chm13.nf.test b/tests/chm13.nf.test index 4e73e43e..e8020b88 100644 --- a/tests/chm13.nf.test +++ b/tests/chm13.nf.test @@ -19,6 +19,9 @@ nextflow_pipeline { then { assertAll( { assert workflow.success}, + + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3'], ['sample3']) }, { //files exist assert file("$outputDir/sample1/variants/clair3/merge_output.vcf.gz").exists() assert file("$outputDir/sample1/variants/clairs/indel.vcf.gz").exists() diff --git a/tests/chm13.nf.test.snap b/tests/chm13.nf.test.snap index a267d392..c6dfc0a5 100644 --- a/tests/chm13.nf.test.snap +++ b/tests/chm13.nf.test.snap @@ -14,6 +14,9 @@ "CLAIR3": { "clair3": "1.2.0" }, + "CLAIR3_PASS_FILTER": { + "bcftools": "1.23.1" + }, "CLAIRS": { "clairs": "0.4.4" }, @@ -23,12 +26,21 @@ "CLAIRSTO_CNA_RESOURCES": { "coreutils": 9.5 }, + "CLAIRS_PASS_FILTER": { + "bcftools": "1.23.1" + }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, "CRAMINO_PRE": { "cramino": "1.3.0" }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -91,16 +103,31 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "SORT_SOMATIC_PHASED": { + "bcftools": 1.22 + }, "SV_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -136,9 +163,9 @@ } ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:19:37.348390574" + "timestamp": "2026-09-29T14:27:45.948524511" } } \ No newline at end of file diff --git a/tests/clair_only.nf.test b/tests/clair_only.nf.test index 13f6d77c..d8f2c35d 100644 --- a/tests/clair_only.nf.test +++ b/tests/clair_only.nf.test @@ -50,17 +50,8 @@ nextflow_pipeline { } }, - // ── Phased VCFs exist and have data ────────────────────────── - { - ['sample1', 'sample2', 'sample3', 'sample4', 'sample5'].each { s -> - def germline = file("$launchDir/output/${s}/variants/phased/germline_smallvariants.vcf.gz") - def somatic = file("$launchDir/output/${s}/variants/phased/somatic_smallvariants.vcf.gz") - assert germline.exists() - assert somatic.exists() - assert germline.size() > 0 - assert somatic.size() > 0 - } - }, + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3', 'sample4', 'sample5'], ['sample3', 'sample4', 'sample5']) }, // ── BAM files exist ────────────────────────────────────────── { diff --git a/tests/clair_only.nf.test.snap b/tests/clair_only.nf.test.snap index ec1ff735..08fdf6aa 100644 --- a/tests/clair_only.nf.test.snap +++ b/tests/clair_only.nf.test.snap @@ -24,18 +24,30 @@ "CLAIR3": { "clair3": "1.2.0" }, + "CLAIR3_PASS_FILTER": { + "bcftools": "1.23.1" + }, "CLAIRS": { "clairs": "0.4.4" }, "CLAIRSTO": { "clairsto": "0.5.1" }, + "CLAIRS_PASS_FILTER": { + "bcftools": "1.23.1" + }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, "CRAMINO_PRE": { "cramino": "1.3.0" }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -104,16 +116,31 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "SORT_SOMATIC_PHASED": { + "bcftools": 1.22 + }, "SV_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -749,38 +776,38 @@ "sample5/vep/somatic/sample5_SOMATIC_VEP.vcf.gz_summary.html" ], [ - "sample1_normal.bam:md5,772e41f7cd86c03a22afbe5ec0592a6b", - "sample1_normal.bam.bai:md5,1b501f6a11efe5d2e6f47b7f1523220b", - "sample1_tumor.bam:md5,c8315c80dc92dfb5d874aef3f5dd46fb", - "sample1_tumor.bam.bai:md5,bc35f807be4b93fc795a14d701469367", + "sample1_normal.bam:md5,9a73c3f90bc4d9a140bcf8f652b0f269", + "sample1_normal.bam.bai:md5,37a866c569f24ed2b38f093f9475b9d6", + "sample1_tumor.bam:md5,dc15bf0e9ff1d401491347c15c408b6d", + "sample1_tumor.bam.bai:md5,fa57db4206692079d9b2084a866bf411", "sample1_normal.flagstat:md5,1c41ea9923945501eb7e41f83a90502d", "sample1_normal.idxstats:md5,902e503387799123ea59255e3fca172c", "sample1_normal.stats:md5,a8b3fba9c54efbc0934d6eacc1807140", "sample1_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample1_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample1_tumor.stats:md5,1c60a1d249d2e503b0678c72e851ea93", - "sample1_whatshap_stats.gtf:md5,eff050a68e36e778b06e0ec19435c569", - "sample1_whatshap_stats.log:md5,76b73731f74fe32ef2d11f6bb0a0f71a", - "sample1_whatshap_stats.tsv:md5,f566ae25b3c5a8f7e94b3d6c1b0417f8", + "sample1_whatshap_stats.gtf:md5,36bda647c08358df0eac0be321c24b20", + "sample1_whatshap_stats.log:md5,938f792fd22bb658a2c7cb1b7653035f", + "sample1_whatshap_stats.tsv:md5,cf8917be389dfef1344eeb6b7e99c7c5", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,47cb0e0bbe71abdbf4f40217dfda43f9", "read_qual.txt:md5,78247dfa2ea336eac0e128eba5e9eef4", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample2_normal.bam:md5,3157bd11ba095a884c7951aafcfcfb1c", - "sample2_normal.bam.bai:md5,edebda44c4383173caea728acde4ac43", - "sample2_tumor.bam:md5,47b2c5f86e0493ba94ff72cea77eeae3", - "sample2_tumor.bam.bai:md5,abf2c290c815f54c2b3f8179f717d9bd", + "sample2_normal.bam:md5,c18bbb1bc05b0ec830fa1ac3c6ef542f", + "sample2_normal.bam.bai:md5,de77684ef264476b3e0531b61ca23740", + "sample2_tumor.bam:md5,8b8cb5ac7668b8ac2f4097f584481d35", + "sample2_tumor.bam.bai:md5,3b46456b0e7b00688519dfb25c8e3c88", "sample2_normal.flagstat:md5,714d0cc0c213e2640e54a16f3d0e6e7e", "sample2_normal.idxstats:md5,72eb83bb11748dc863fef1a0a5497e4b", "sample2_normal.stats:md5,20c47cb94f9ac739d69c57be6daf82c5", "sample2_tumor.flagstat:md5,4344a8745efef9cc2a017024218d61c6", "sample2_tumor.idxstats:md5,69467fc02c83a30084736aeea8b785fb", "sample2_tumor.stats:md5,8635df10132c85a13f2d9878b7cf90a2", - "sample2_whatshap_stats.gtf:md5,4d8f4393e3aebe4e945c0b8236cf3b3e", - "sample2_whatshap_stats.log:md5,10bba7bae6dd99b989ece5e5dac7a8f9", - "sample2_whatshap_stats.tsv:md5,bb46226e486af9026ab76e014624e903", + "sample2_whatshap_stats.gtf:md5,35cd28699c298d99d01cee1c24c6d61b", + "sample2_whatshap_stats.log:md5,75a69a8e651979e25467d6e3c84cdf90", + "sample2_whatshap_stats.tsv:md5,92d7234e355833b1ec6f54951c38c09d", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,48baac86492026a4a7947bc708c47e6e", @@ -830,9 +857,9 @@ ] ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:27:45.728407577" + "timestamp": "2026-09-29T13:57:41.300856488" } } \ No newline at end of file diff --git a/tests/consensus.nf.test b/tests/consensus.nf.test index 9689f212..cc1fd98f 100644 --- a/tests/consensus.nf.test +++ b/tests/consensus.nf.test @@ -51,17 +51,8 @@ nextflow_pipeline { assert file("$launchDir/output/sample3/variants/clairsto/somatic.vcf.gz").exists() }, - // ── Phased consensus VCFs exist and have data ──────────────── - { - ['sample1', 'sample2', 'sample3'].each { s -> - def germline = file("$launchDir/output/${s}/variants/phased/germline_smallvariants.vcf.gz") - def somatic = file("$launchDir/output/${s}/variants/phased/somatic_smallvariants.vcf.gz") - assert germline.exists() - assert somatic.exists() - assert germline.size() > 0 - assert somatic.size() > 0 - } - }, + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3'], ['sample3']) }, // ── BAM files ──────────────────────────────────────────────── { diff --git a/tests/consensus.nf.test.snap b/tests/consensus.nf.test.snap index 6bf3a01d..7260fad2 100644 --- a/tests/consensus.nf.test.snap +++ b/tests/consensus.nf.test.snap @@ -14,6 +14,9 @@ "BCFTOOLS_NORM": { "bcftools": 1.22 }, + "BCFTOOLS_NORM_REJOIN": { + "bcftools": 1.22 + }, "BCFTOOLS_QUERY": { "bcftools": 1.22 }, @@ -29,12 +32,21 @@ "CLAIR3": { "clair3": "1.2.0" }, + "CLAIR3_PASS_FILTER": { + "bcftools": "1.23.1" + }, "CLAIRS": { "clairs": "0.4.4" }, "CLAIRSTO": { "clairsto": "0.5.1" }, + "CLAIRS_PASS_FILTER": { + "bcftools": "1.23.1" + }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, @@ -47,6 +59,9 @@ "DEEPSOMATIC_MAKEEXAMPLES": { "deepsomatic": "1.7.0" }, + "DEEPSOMATIC_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPSOMATIC_POSTPROCESSVARIANTS": { "deepsomatic": "1.7.0" }, @@ -56,9 +71,24 @@ "DEEPVARIANT_MAKEEXAMPLES": { "deepvariant": "1.9.0" }, + "DEEPVARIANT_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPVARIANT_POSTPROCESSVARIANTS": { "deepvariant": "1.9.0" }, + "DS_GERMLINE_SELECT": { + "bcftools": "1.23.1" + }, + "DS_VERDICT_ANNOTATE": { + "bcftools": 1.22 + }, + "DS_VERDICT_QUERY": { + "bcftools": 1.22 + }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -121,6 +151,12 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -129,11 +165,20 @@ "SORT_POST_NORM": { "bcftools": 1.22 }, + "SORT_SOMATIC_PHASED": { + "bcftools": 1.22 + }, "SV_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -607,38 +652,38 @@ "sample3/vep/somatic/sample3_SOMATIC_VEP.vcf.gz_summary.html" ], [ - "sample1_normal.bam:md5,93dd8ac8b67eb4eb4bf27e09c8f5f99b", - "sample1_normal.bam.bai:md5,75402ef1cc35229cc131155d9ec973e0", - "sample1_tumor.bam:md5,69cba03cad51bcc1d1ee8c48da042527", - "sample1_tumor.bam.bai:md5,5d633ed05021ad81ce24b1f18cbf38b4", + "sample1_normal.bam:md5,5d3f0615b0d8a9748e6b6bc7d76ab577", + "sample1_normal.bam.bai:md5,3ef784f2536c7506a50d0ed20c3d78fb", + "sample1_tumor.bam:md5,7317ffe61fe3c4f73961db60329adfa9", + "sample1_tumor.bam.bai:md5,6aaeaebe40c422358bb55296c508db2b", "sample1_normal.flagstat:md5,1c41ea9923945501eb7e41f83a90502d", "sample1_normal.idxstats:md5,902e503387799123ea59255e3fca172c", "sample1_normal.stats:md5,a8b3fba9c54efbc0934d6eacc1807140", "sample1_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample1_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample1_tumor.stats:md5,1c60a1d249d2e503b0678c72e851ea93", - "sample1_whatshap_stats.gtf:md5,9f09f9ad1a788384cb8e46a933f77b3b", - "sample1_whatshap_stats.log:md5,20135b4e9965a31d3f9bb0df7d2cec90", - "sample1_whatshap_stats.tsv:md5,264d2d76a9b8d34ea4933aee325ce36e", + "sample1_whatshap_stats.gtf:md5,30bde8f88b7d4e88b935e88e00997ce7", + "sample1_whatshap_stats.log:md5,a3cf683c728ce63a5c2be33edd20f9cc", + "sample1_whatshap_stats.tsv:md5,cca55ca6f99eb0857759a4007c10b346", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,47cb0e0bbe71abdbf4f40217dfda43f9", "read_qual.txt:md5,78247dfa2ea336eac0e128eba5e9eef4", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample2_normal.bam:md5,f2ba30d007c521d479c6158e1e22367a", - "sample2_normal.bam.bai:md5,c3096f52115ec1e24c46fedc41f1f3d3", - "sample2_tumor.bam:md5,1c0287d24fa5b25b86e48024f2f55031", - "sample2_tumor.bam.bai:md5,62849cea5a005e3d8dbe8f9edcefaf60", + "sample2_normal.bam:md5,ece242f245fc4e2210e13e6de93b8fdc", + "sample2_normal.bam.bai:md5,d96d0071ab25ca8dc2327acba4395515", + "sample2_tumor.bam:md5,0c538dc0e27566677313d1f9374fa2b1", + "sample2_tumor.bam.bai:md5,ab2fca59e6729e011c30310eaee2ef1c", "sample2_normal.flagstat:md5,714d0cc0c213e2640e54a16f3d0e6e7e", "sample2_normal.idxstats:md5,72eb83bb11748dc863fef1a0a5497e4b", "sample2_normal.stats:md5,20c47cb94f9ac739d69c57be6daf82c5", "sample2_tumor.flagstat:md5,4344a8745efef9cc2a017024218d61c6", "sample2_tumor.idxstats:md5,69467fc02c83a30084736aeea8b785fb", "sample2_tumor.stats:md5,8635df10132c85a13f2d9878b7cf90a2", - "sample2_whatshap_stats.gtf:md5,f15fb43f0af73d02fc73b66fdc12d5d8", - "sample2_whatshap_stats.log:md5,ca87088fc2f11665eca3fb9c80489085", - "sample2_whatshap_stats.tsv:md5,ca53f81e39bf5d46aa4f604216add1f6", + "sample2_whatshap_stats.gtf:md5,2e5ace4cac0b42bb6132513062781e47", + "sample2_whatshap_stats.log:md5,0979359e459a14037a704733f32d3b95", + "sample2_whatshap_stats.tsv:md5,8c974019af5c839c72604e7526ae8a3d", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,48baac86492026a4a7947bc708c47e6e", @@ -662,9 +707,9 @@ ] ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:40:36.642665738" + "timestamp": "2026-09-29T13:31:57.487084446" } } \ No newline at end of file diff --git a/tests/deep_only.nf.test b/tests/deep_only.nf.test index 5da64375..b5ad3347 100644 --- a/tests/deep_only.nf.test +++ b/tests/deep_only.nf.test @@ -41,17 +41,8 @@ nextflow_pipeline { assert !file("$launchDir/output/sample3/variants/clairsto/somatic.vcf.gz").exists() }, - // ── Phased VCFs exist and have data ────────────────────────── - { - ['sample1', 'sample2', 'sample3'].each { s -> - def germline = file("$launchDir/output/${s}/variants/phased/germline_smallvariants.vcf.gz") - def somatic = file("$launchDir/output/${s}/variants/phased/somatic_smallvariants.vcf.gz") - assert germline.exists() - assert somatic.exists() - assert germline.size() > 0 - assert somatic.size() > 0 - } - }, + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3'], ['sample3']) }, // ── BAM files ──────────────────────────────────────────────── { diff --git a/tests/deep_only.nf.test.snap b/tests/deep_only.nf.test.snap index 7887aaca..95f8fdb1 100644 --- a/tests/deep_only.nf.test.snap +++ b/tests/deep_only.nf.test.snap @@ -11,6 +11,9 @@ "BCFTOOLS_VIEW": { "bcftools": 1.22 }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, @@ -23,6 +26,9 @@ "DEEPSOMATIC_MAKEEXAMPLES": { "deepsomatic": "1.7.0" }, + "DEEPSOMATIC_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPSOMATIC_POSTPROCESSVARIANTS": { "deepsomatic": "1.7.0" }, @@ -32,9 +38,24 @@ "DEEPVARIANT_MAKEEXAMPLES": { "deepvariant": "1.9.0" }, + "DEEPVARIANT_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPVARIANT_POSTPROCESSVARIANTS": { "deepvariant": "1.9.0" }, + "DS_GERMLINE_SELECT": { + "bcftools": "1.23.1" + }, + "DS_VERDICT_ANNOTATE": { + "bcftools": 1.22 + }, + "DS_VERDICT_QUERY": { + "bcftools": 1.22 + }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -97,16 +118,31 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "SORT_SOMATIC_PHASED": { + "bcftools": 1.22 + }, "SV_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -563,8 +599,8 @@ "sample1_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample1_tumor.stats:md5,1c60a1d249d2e503b0678c72e851ea93", "sample1_whatshap_stats.gtf:md5,e1d0e87353a5f9aed8a9ac4bf7973427", - "sample1_whatshap_stats.log:md5,bd6b83a062e22cd3201523dc4c2c13e7", - "sample1_whatshap_stats.tsv:md5,7a1508751cb1daa841a577ae25f55586", + "sample1_whatshap_stats.log:md5,43811a1aa4726c2ef62646ff22d1fad0", + "sample1_whatshap_stats.tsv:md5,a821c5d0e3451f82327645fee89e68eb", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,47cb0e0bbe71abdbf4f40217dfda43f9", @@ -582,22 +618,22 @@ "sample2_tumor.idxstats:md5,69467fc02c83a30084736aeea8b785fb", "sample2_tumor.stats:md5,8635df10132c85a13f2d9878b7cf90a2", "sample2_whatshap_stats.gtf:md5,af33281699a1d0da83fbe7eaff198d03", - "sample2_whatshap_stats.log:md5,bbd9ab2ce07a009d9348a1d78bc6fc70", - "sample2_whatshap_stats.tsv:md5,c65436f930c23ddbfd568532d07dce70", + "sample2_whatshap_stats.log:md5,0cd536df69e244a5271e9c5440ea0f3f", + "sample2_whatshap_stats.tsv:md5,c7cc47024ef622a72e7f4f0d5fa119eb", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,48baac86492026a4a7947bc708c47e6e", "read_qual.txt:md5,8b92ff7dc4536188be159b95525511cd", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample3_tumor.bam:md5,965264ef8436cb887ad02c92d40fb50e", - "sample3_tumor.bam.bai:md5,cfc6329667a3c6c66e3c0ca0ace4c6e9", + "sample3_tumor.bam:md5,9449340acf09825767e1f5b4a3f0e5c2", + "sample3_tumor.bam.bai:md5,3be46bae7402c3865f27881457b7a466", "sample3_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample3_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample3_tumor.stats:md5,ecd5ea4fee37379dd5c5ae3e89dfddda", - "sample3_whatshap_stats.gtf:md5,f47156e18c490ff9a4e6efd04d43acc5", - "sample3_whatshap_stats.log:md5,4f7648e763004ab764143cb4f8b6499e", - "sample3_whatshap_stats.tsv:md5,4cb58bb3b663aaba23da004d69adab3e", + "sample3_whatshap_stats.gtf:md5,5a9b20b6ed25ef2ede4271156a83ce88", + "sample3_whatshap_stats.log:md5,a62c149286176f863d4293d0748a6dec", + "sample3_whatshap_stats.tsv:md5,90c20ce05eef472d9f5bb77754636d31", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,56e899f85876cee082788927d0f89c5f", @@ -607,9 +643,9 @@ ] ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:46:42.562085661" + "timestamp": "2026-09-29T15:13:24.013345908" } } \ No newline at end of file diff --git a/tests/default.nf.test b/tests/default.nf.test index 7c473eba..4e1d7d2c 100644 --- a/tests/default.nf.test +++ b/tests/default.nf.test @@ -22,6 +22,9 @@ nextflow_pipeline { def stable_path = getAllFilesFromDir(params.outdir, ignoreFile: 'tests/.nftignore') assertAll( { assert workflow.success}, + + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3'], ['sample3']) }, { //files exist assert file("$launchDir/output/sample1/variants/clair3/merge_output.vcf.gz").exists() assert file("$launchDir/output/sample1/variants/clair3/merge_output.vcf.gz.tbi").exists() diff --git a/tests/default.nf.test.snap b/tests/default.nf.test.snap index 22e72e6f..b36e5f38 100644 --- a/tests/default.nf.test.snap +++ b/tests/default.nf.test.snap @@ -14,18 +14,30 @@ "CLAIR3": { "clair3": "1.2.0" }, + "CLAIR3_PASS_FILTER": { + "bcftools": "1.23.1" + }, "CLAIRS": { "clairs": "0.4.4" }, "CLAIRSTO": { "clairsto": "0.5.1" }, + "CLAIRS_PASS_FILTER": { + "bcftools": "1.23.1" + }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, "CRAMINO_PRE": { "cramino": "1.3.0" }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -88,16 +100,31 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "SORT_SOMATIC_PHASED": { + "bcftools": 1.22 + }, "SV_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -553,38 +580,38 @@ "sample3/vep/somatic/sample3_SOMATIC_VEP.vcf.gz_summary.html" ], [ - "sample1_normal.bam:md5,772e41f7cd86c03a22afbe5ec0592a6b", - "sample1_normal.bam.bai:md5,1b501f6a11efe5d2e6f47b7f1523220b", - "sample1_tumor.bam:md5,c8315c80dc92dfb5d874aef3f5dd46fb", - "sample1_tumor.bam.bai:md5,bc35f807be4b93fc795a14d701469367", + "sample1_normal.bam:md5,9a73c3f90bc4d9a140bcf8f652b0f269", + "sample1_normal.bam.bai:md5,37a866c569f24ed2b38f093f9475b9d6", + "sample1_tumor.bam:md5,dc15bf0e9ff1d401491347c15c408b6d", + "sample1_tumor.bam.bai:md5,fa57db4206692079d9b2084a866bf411", "sample1_normal.flagstat:md5,1c41ea9923945501eb7e41f83a90502d", "sample1_normal.idxstats:md5,902e503387799123ea59255e3fca172c", "sample1_normal.stats:md5,a8b3fba9c54efbc0934d6eacc1807140", "sample1_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample1_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample1_tumor.stats:md5,1c60a1d249d2e503b0678c72e851ea93", - "sample1_whatshap_stats.gtf:md5,eff050a68e36e778b06e0ec19435c569", - "sample1_whatshap_stats.log:md5,76b73731f74fe32ef2d11f6bb0a0f71a", - "sample1_whatshap_stats.tsv:md5,f566ae25b3c5a8f7e94b3d6c1b0417f8", + "sample1_whatshap_stats.gtf:md5,36bda647c08358df0eac0be321c24b20", + "sample1_whatshap_stats.log:md5,938f792fd22bb658a2c7cb1b7653035f", + "sample1_whatshap_stats.tsv:md5,cf8917be389dfef1344eeb6b7e99c7c5", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,47cb0e0bbe71abdbf4f40217dfda43f9", "read_qual.txt:md5,78247dfa2ea336eac0e128eba5e9eef4", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample2_normal.bam:md5,3157bd11ba095a884c7951aafcfcfb1c", - "sample2_normal.bam.bai:md5,edebda44c4383173caea728acde4ac43", - "sample2_tumor.bam:md5,47b2c5f86e0493ba94ff72cea77eeae3", - "sample2_tumor.bam.bai:md5,abf2c290c815f54c2b3f8179f717d9bd", + "sample2_normal.bam:md5,c18bbb1bc05b0ec830fa1ac3c6ef542f", + "sample2_normal.bam.bai:md5,de77684ef264476b3e0531b61ca23740", + "sample2_tumor.bam:md5,8b8cb5ac7668b8ac2f4097f584481d35", + "sample2_tumor.bam.bai:md5,3b46456b0e7b00688519dfb25c8e3c88", "sample2_normal.flagstat:md5,714d0cc0c213e2640e54a16f3d0e6e7e", "sample2_normal.idxstats:md5,72eb83bb11748dc863fef1a0a5497e4b", "sample2_normal.stats:md5,20c47cb94f9ac739d69c57be6daf82c5", "sample2_tumor.flagstat:md5,4344a8745efef9cc2a017024218d61c6", "sample2_tumor.idxstats:md5,69467fc02c83a30084736aeea8b785fb", "sample2_tumor.stats:md5,8635df10132c85a13f2d9878b7cf90a2", - "sample2_whatshap_stats.gtf:md5,4d8f4393e3aebe4e945c0b8236cf3b3e", - "sample2_whatshap_stats.log:md5,10bba7bae6dd99b989ece5e5dac7a8f9", - "sample2_whatshap_stats.tsv:md5,bb46226e486af9026ab76e014624e903", + "sample2_whatshap_stats.gtf:md5,35cd28699c298d99d01cee1c24c6d61b", + "sample2_whatshap_stats.log:md5,75a69a8e651979e25467d6e3c84cdf90", + "sample2_whatshap_stats.tsv:md5,92d7234e355833b1ec6f54951c38c09d", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,48baac86492026a4a7947bc708c47e6e", @@ -608,9 +635,9 @@ ] ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:10:56.589241561" + "timestamp": "2026-09-29T13:03:39.7913195" } } \ No newline at end of file diff --git a/tests/fixtures/excludesites_input.vcf b/tests/fixtures/excludesites_input.vcf new file mode 100644 index 00000000..615bfb45 --- /dev/null +++ b/tests/fixtures/excludesites_input.vcf @@ -0,0 +1,9 @@ +##fileformat=VCFv4.2 +##contig= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT testsample +chr1 100 . A T 30 PASS . GT 0/1 +chr1 102 . G C 30 PASS . GT 0/1 +chr1 200 . T C 30 PASS . GT 1/1 +chr1 300 . C G 30 PASS . GT 0/1 +chr1 400 . G A 30 PASS . GT 0/1 diff --git a/tests/fixtures/excludesites_input.vcf.gz b/tests/fixtures/excludesites_input.vcf.gz new file mode 100644 index 00000000..a22f7cb6 Binary files /dev/null and b/tests/fixtures/excludesites_input.vcf.gz differ diff --git a/tests/fixtures/excludesites_input.vcf.gz.tbi b/tests/fixtures/excludesites_input.vcf.gz.tbi new file mode 100644 index 00000000..212a1c1e Binary files /dev/null and b/tests/fixtures/excludesites_input.vcf.gz.tbi differ diff --git a/tests/fixtures/excludesites_mask.vcf b/tests/fixtures/excludesites_mask.vcf new file mode 100644 index 00000000..bf2aaef2 --- /dev/null +++ b/tests/fixtures/excludesites_mask.vcf @@ -0,0 +1,7 @@ +##fileformat=VCFv4.2 +##contig= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT testsample +chr1 100 . ACGTA A 30 PASS . GT 0/1 +chr1 200 . T A,G 30 PASS . GT 1/2 +chr1 400 . G A 30 PASS . GT 0/1 diff --git a/tests/fixtures/excludesites_mask.vcf.gz b/tests/fixtures/excludesites_mask.vcf.gz new file mode 100644 index 00000000..ccd5fc68 Binary files /dev/null and b/tests/fixtures/excludesites_mask.vcf.gz differ diff --git a/tests/fixtures/excludesites_mask.vcf.gz.tbi b/tests/fixtures/excludesites_mask.vcf.gz.tbi new file mode 100644 index 00000000..869e7d78 Binary files /dev/null and b/tests/fixtures/excludesites_mask.vcf.gz.tbi differ diff --git a/tests/fixtures/vcftag_input.vcf b/tests/fixtures/vcftag_input.vcf new file mode 100644 index 00000000..a1945635 --- /dev/null +++ b/tests/fixtures/vcftag_input.vcf @@ -0,0 +1,13 @@ +##fileformat=VCFv4.2 +##INFO= +##INFO= +##FILTER= +##FILTER= +##FILTER= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT testsample +chr1 100 . A C 30 PASS CALLER=clairs-to GT:DP 0/1:20 +chr1 200 . G T 20 NonSomatic CALLER=clairs-to;EXISTING GT:DP 0/1:18 +chr1 300 . T A 10 RefCall . GT:DP 0/0:12 +chr1 400 . C G 40 . CALLER=deepsomatic GT:DP 1/1:25 +chr1 500 . G A 5 LowQual;RefCall CALLER=clair3;EXISTING GT:DP 0/1:9 diff --git a/tests/fixtures/vcftag_input.vcf.gz b/tests/fixtures/vcftag_input.vcf.gz new file mode 100644 index 00000000..c956c3c5 Binary files /dev/null and b/tests/fixtures/vcftag_input.vcf.gz differ diff --git a/tests/fixtures/vcftag_input.vcf.gz.tbi b/tests/fixtures/vcftag_input.vcf.gz.tbi new file mode 100644 index 00000000..75130768 Binary files /dev/null and b/tests/fixtures/vcftag_input.vcf.gz.tbi differ diff --git a/tests/lib/PhasedVcf.groovy b/tests/lib/PhasedVcf.groovy new file mode 100644 index 00000000..12161928 --- /dev/null +++ b/tests/lib/PhasedVcf.groovy @@ -0,0 +1,44 @@ +import java.util.zip.GZIPInputStream + +// Record-level checks on the published phased VCFs, shared by the pipeline tests. +class PhasedVcf { + + // Non-header records of a bgzipped VCF, each split into its tab-separated fields + static List> records(String path) { + new GZIPInputStream(new FileInputStream(path)).withReader { reader -> + reader.readLines().findAll { !it.startsWith('#') }.collect { it.split('\t') as List } + } + } + + static boolean hasFlag(List rec, String flag) { + rec[7].split(';').contains(flag) + } + + // Value of a FORMAT field in the first sample column, or null if absent + static String format(List rec, String key) { + def i = rec[8].split(':').toList().indexOf(key) + def values = rec[9].split(':') + i >= 0 && i < values.size() ? values[i] : null + } + + static boolean isAlt(List rec) { + format(rec, 'GT').split(/[\/|]/).any { it != '0' && it != '.' } + } + + // Default-mode checks: both arms all PASS, somatic records all carry SOMATIC, germline not empty + static void assertPassOnly(String outdir, List samples, List tumourOnly) { + samples.each { s -> + def germ = records("${outdir}/${s}/variants/phased/germline_smallvariants.vcf.gz") + def som = records("${outdir}/${s}/variants/phased/somatic_smallvariants.vcf.gz") + assert !germ.isEmpty() : "${s}: phased germline VCF has no records" + assert germ.every { it[6] == 'PASS' } : "${s}: non-PASS germline record with the PASS filter on" + assert som.every { it[6] == 'PASS' } : "${s}: non-PASS somatic record with the PASS filter on" + assert som.every { hasFlag(it, 'SOMATIC') } : "${s}: somatic record without INFO/SOMATIC" + } + // Paired test samples have no PASS somatic calls, so the checks above only bite on tumour-only ones + tumourOnly.each { s -> + assert !records("${outdir}/${s}/variants/phased/somatic_smallvariants.vcf.gz").isEmpty() : + "${s}: tumour-only phased somatic VCF has no records" + } + } +} diff --git a/tests/union.nf.test b/tests/union.nf.test index 2fe8925b..155ffa60 100644 --- a/tests/union.nf.test +++ b/tests/union.nf.test @@ -51,16 +51,19 @@ nextflow_pipeline { assert file("$launchDir/output/sample3/variants/clairsto/somatic.vcf.gz").exists() }, - // ── Phased union VCFs exist and have data ──────────────────── + // ── Phased VCFs: records present, all PASS, somatic arm tagged ── + { PhasedVcf.assertPassOnly("$outputDir", ['sample1', 'sample2', 'sample3'], ['sample3']) }, + // ── 'all' mode keeps one caller's record per position ──────── { ['sample1', 'sample2', 'sample3'].each { s -> - def germline = file("$launchDir/output/${s}/variants/phased/germline_smallvariants.vcf.gz") - def somatic = file("$launchDir/output/${s}/variants/phased/somatic_smallvariants.vcf.gz") - assert germline.exists() - assert somatic.exists() - assert germline.size() > 0 - assert somatic.size() > 0 + def germ = PhasedVcf.records("$outputDir/${s}/variants/phased/germline_smallvariants.vcf.gz") + def pos = germ.collect { "${it[0]}:${it[1]}" } + assert pos.size() == pos.toUnique().size() : "${s}: two germline records at one position" } + // DeepVariant calls T>A,G 1/2 here and Clair3 T>A 1/1; the prioritised Clair3 record is kept whole + def site = PhasedVcf.records("$outputDir/sample2/variants/phased/germline_smallvariants.vcf.gz") + .find { it[0] == 'chr19' && it[1] == '24897193' } + assert site && site[4] == 'A' && site[7].split(';').contains('CALLER=clair3') }, // ── BAM files ──────────────────────────────────────────────── @@ -115,4 +118,49 @@ nextflow_pipeline { ) } } + + test("-profile test, union combine mode, smallvar_filter_pass=false") { + + when { + params { + outdir = "$outputDir" + germline_var_combine = 'all' + somatic_var_combine = 'all' + germline_var_keep = 'clair, deepvariant' + somatic_var_keep = 'clair, deepsomatic' + smallvar_filter_pass = false + } + } + + then { + assertAll( + { assert workflow.success }, + + // ── No PASS filter process may run ─────────────────────────── + // With the param false no _PASS_FILTER process runs or reports a version. + { + def versions = file("$outputDir/pipeline_info/lrsomatic_software_mqc_versions.yml") + assert versions.exists() + assert !versions.text.contains('_PASS_FILTER') + }, + + // ── Phased VCFs keep non-PASS records; somatic 0/0 records are not phased ── + { + ['sample1', 'sample2', 'sample3'].each { s -> + def germ = PhasedVcf.records("$outputDir/${s}/variants/phased/germline_smallvariants.vcf.gz") + def som = PhasedVcf.records("$outputDir/${s}/variants/phased/somatic_smallvariants.vcf.gz") + assert !germ.isEmpty() : "${s}: phased germline VCF has no records" + assert !som.isEmpty() : "${s}: phased somatic VCF has no records" + assert germ.any { it[6] != 'PASS' } : "${s}: no non-PASS germline record with the filter off" + assert som.any { it[6] != 'PASS' } : "${s}: no non-PASS somatic record with the filter off" + assert som.every { PhasedVcf.hasFlag(it, 'SOMATIC') } : "${s}: somatic record without INFO/SOMATIC" + // LongPhase would copy a colliding germline record's phase onto a non-alt somatic record + assert som.findAll { !PhasedVcf.isAlt(it) }.every { PhasedVcf.format(it, 'PS') in [null, '.'] } : + "${s}: non-alt somatic record carries a phase set" + } + } + ) + } + } + } diff --git a/tests/union.nf.test.snap b/tests/union.nf.test.snap index 654a82d7..9979ecac 100644 --- a/tests/union.nf.test.snap +++ b/tests/union.nf.test.snap @@ -8,7 +8,7 @@ "BCFTOOLS_CONCAT": { "bcftools": 1.22 }, - "BCFTOOLS_ISEC": { + "BCFTOOLS_EXCLUDE_SITES": { "bcftools": 1.22 }, "BCFTOOLS_NORM": { @@ -26,12 +26,21 @@ "CLAIR3": { "clair3": "1.2.0" }, + "CLAIR3_PASS_FILTER": { + "bcftools": "1.23.1" + }, "CLAIRS": { "clairs": "0.4.4" }, "CLAIRSTO": { "clairsto": "0.5.1" }, + "CLAIRS_PASS_FILTER": { + "bcftools": "1.23.1" + }, + "CONCAT_SOMATIC_UNPHASED": { + "bcftools": 1.22 + }, "CRAMINO_POST": { "cramino": "1.3.0" }, @@ -44,6 +53,9 @@ "DEEPSOMATIC_MAKEEXAMPLES": { "deepsomatic": "1.7.0" }, + "DEEPSOMATIC_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPSOMATIC_POSTPROCESSVARIANTS": { "deepsomatic": "1.7.0" }, @@ -53,9 +65,24 @@ "DEEPVARIANT_MAKEEXAMPLES": { "deepvariant": "1.9.0" }, + "DEEPVARIANT_PASS_FILTER": { + "bcftools": "1.23.1" + }, "DEEPVARIANT_POSTPROCESSVARIANTS": { "deepvariant": "1.9.0" }, + "DS_GERMLINE_SELECT": { + "bcftools": "1.23.1" + }, + "DS_VERDICT_ANNOTATE": { + "bcftools": 1.22 + }, + "DS_VERDICT_QUERY": { + "bcftools": 1.22 + }, + "GERMLINE_ANCHORS": { + "bcftools": 1.22 + }, "GERMLINE_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -118,6 +145,12 @@ "SEVERUS": { "severus": 1.6 }, + "SOMATIC_ALT": { + "bcftools": 1.22 + }, + "SOMATIC_NONALT": { + "bcftools": 1.22 + }, "SOMATIC_VEP": { "ensemblvep": 115.2, "perl-math-cdf": 0.1, @@ -126,7 +159,7 @@ "SORT_POST_NORM": { "bcftools": 1.22 }, - "STANDARDIZE_AF": { + "SORT_SOMATIC_PHASED": { "bcftools": 1.22 }, "SV_VEP": { @@ -134,6 +167,12 @@ "perl-math-cdf": 0.1, "tabix": 1.21 }, + "TAG_GERMLINE": { + "bcftools": 1.2 + }, + "TAG_SOMATIC": { + "bcftools": 1.2 + }, "UNTAR": { "untar": 1.34 }, @@ -607,52 +646,52 @@ "sample3/vep/somatic/sample3_SOMATIC_VEP.vcf.gz_summary.html" ], [ - "sample1_normal.bam:md5,cbfc940a38c74cbe8435c18b9da5dd32", - "sample1_normal.bam.bai:md5,c1498328929d45b2898fa2265b0d617c", - "sample1_tumor.bam:md5,b78866edf991393806d37505d16f7e3d", - "sample1_tumor.bam.bai:md5,f613de14ab19fc3a85403661a4f6188c", + "sample1_normal.bam:md5,969a2009965e03c84ca3ec46c2fbbd38", + "sample1_normal.bam.bai:md5,62da15ac99896f0e90a694bf52ba19d3", + "sample1_tumor.bam:md5,6267bf3e86a69536e45968b908d1cc6f", + "sample1_tumor.bam.bai:md5,1ec6d68c92182ff9051b8707c4975445", "sample1_normal.flagstat:md5,1c41ea9923945501eb7e41f83a90502d", "sample1_normal.idxstats:md5,902e503387799123ea59255e3fca172c", "sample1_normal.stats:md5,a8b3fba9c54efbc0934d6eacc1807140", "sample1_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample1_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample1_tumor.stats:md5,1c60a1d249d2e503b0678c72e851ea93", - "sample1_whatshap_stats.gtf:md5,9ae556e13516dd47d4108acf2104bddb", - "sample1_whatshap_stats.log:md5,eaddcf6a1666d4a3c1ad3316dac24139", - "sample1_whatshap_stats.tsv:md5,c2773e011c2781160fd9a7741b10546b", + "sample1_whatshap_stats.gtf:md5,f5d331899db63b2ae21e51436bad2ddd", + "sample1_whatshap_stats.log:md5,0f455b7a59d41ee276844a1ac894971c", + "sample1_whatshap_stats.tsv:md5,41c0b3aeed7b8bc7770c7bce5d40cfcf", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,47cb0e0bbe71abdbf4f40217dfda43f9", "read_qual.txt:md5,78247dfa2ea336eac0e128eba5e9eef4", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample2_normal.bam:md5,f2ba30d007c521d479c6158e1e22367a", - "sample2_normal.bam.bai:md5,c3096f52115ec1e24c46fedc41f1f3d3", - "sample2_tumor.bam:md5,1c0287d24fa5b25b86e48024f2f55031", - "sample2_tumor.bam.bai:md5,62849cea5a005e3d8dbe8f9edcefaf60", + "sample2_normal.bam:md5,c18bbb1bc05b0ec830fa1ac3c6ef542f", + "sample2_normal.bam.bai:md5,de77684ef264476b3e0531b61ca23740", + "sample2_tumor.bam:md5,8b8cb5ac7668b8ac2f4097f584481d35", + "sample2_tumor.bam.bai:md5,3b46456b0e7b00688519dfb25c8e3c88", "sample2_normal.flagstat:md5,714d0cc0c213e2640e54a16f3d0e6e7e", "sample2_normal.idxstats:md5,72eb83bb11748dc863fef1a0a5497e4b", "sample2_normal.stats:md5,20c47cb94f9ac739d69c57be6daf82c5", "sample2_tumor.flagstat:md5,4344a8745efef9cc2a017024218d61c6", "sample2_tumor.idxstats:md5,69467fc02c83a30084736aeea8b785fb", "sample2_tumor.stats:md5,8635df10132c85a13f2d9878b7cf90a2", - "sample2_whatshap_stats.gtf:md5,f15fb43f0af73d02fc73b66fdc12d5d8", - "sample2_whatshap_stats.log:md5,a6767b3490cafdcbaf3b7114644028de", - "sample2_whatshap_stats.tsv:md5,570796e5e291229e8872733425e0b133", + "sample2_whatshap_stats.gtf:md5,35cd28699c298d99d01cee1c24c6d61b", + "sample2_whatshap_stats.log:md5,e7293b22e0b3ac0798b31d897d781336", + "sample2_whatshap_stats.tsv:md5,1e9adab54274feead2597e09f558222c", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,48baac86492026a4a7947bc708c47e6e", "read_qual.txt:md5,8b92ff7dc4536188be159b95525511cd", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "sample3_tumor.bam:md5,116b6944da4aa833a8d21c46b5f5ecfe", - "sample3_tumor.bam.bai:md5,ea9eca53bbaba26d40b791a2ea1aadf6", + "sample3_tumor.bam:md5,94e7f18f0a22d16f0868780df7975447", + "sample3_tumor.bam.bai:md5,a2cc1013d6fde4fec50608fa9681a7a7", "sample3_tumor.flagstat:md5,8ff32d733c62c4910bf185ef24bf27cf", "sample3_tumor.idxstats:md5,2de140e61f9e86c9c10af20dd565cc93", "sample3_tumor.stats:md5,ecd5ea4fee37379dd5c5ae3e89dfddda", - "sample3_whatshap_stats.gtf:md5,f47156e18c490ff9a4e6efd04d43acc5", - "sample3_whatshap_stats.log:md5,679dcfa209888a9e69a07e4c4e4b049e", - "sample3_whatshap_stats.tsv:md5,035d5aa0425ba3fc32d65268b793b424", + "sample3_whatshap_stats.gtf:md5,5a9b20b6ed25ef2ede4271156a83ce88", + "sample3_whatshap_stats.log:md5,dc9a34518331417b7b5811462b722791", + "sample3_whatshap_stats.tsv:md5,26a07a6db5f1e5bf7844aa0a8c190cb7", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", "breakpoints_double.csv:md5,56e899f85876cee082788927d0f89c5f", @@ -662,9 +701,9 @@ ] ], "meta": { - "nf-test": "0.9.0", - "nextflow": "26.04.6" + "nf-test": "0.9.3", + "nextflow": "25.10.4" }, - "timestamp": "2026-09-21T11:54:50.88115842" + "timestamp": "2026-09-29T12:18:55.043468734" } } \ No newline at end of file diff --git a/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index 37bd1c4c..b5c6b7e9 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -28,6 +28,8 @@ include { NANOPLOT as NANOPLOT_PRE } from '../modules/nf-core/nanoplot/ include { NANOPLOT as NANOPLOT_POST } from '../modules/nf-core/nanoplot/main' include { MOSDEPTH } from '../modules/nf-core/mosdepth/main' include { ASCAT } from '../modules/nf-core/ascat/main' +include { ECDNA } from '../subworkflows/local/ecdna' +include { PREPARE_AA_DATA_REPO } from '../subworkflows/local/prepare_aa_data_repo' include { SEVERUS } from '../modules/nf-core/severus/main.nf' include { METAEXTRACT } from '../modules/local/metaextract/main' include { CLAIRSTO_CNA_RESOURCES } from '../modules/local/clairsto/cna_resources/main' @@ -115,6 +117,9 @@ workflow LRSOMATIC { params.vep_species = getGenomeAttribute('vep_species') params.sigprofiler_genome = getGenomeAttribute('sigprofiler_genome') params.sigprofiler_genome_url = getGenomeAttribute('sigprofiler_genome_url') + params.coral_ref = getGenomeAttribute('coral_ref') + params.ac_ref = getGenomeAttribute('ac_ref') + params.aa_data_repo_url = getGenomeAttribute('aa_data_repo_url') // Resolved once here to avoid a HEAD request per default plugin URL, and passed straight to // the VEP tasks: conf/modules.config closures do not see a param assigned here. @@ -190,9 +195,27 @@ workflow LRSOMATIC { if (params.clairsto_cna_resources && !params.skip_ascat) { log.warn("--clairsto_cna_resources is ignored without --skip_ascat: Verdict's germline tagging then comes from ASCAT's purity and copy number.") } + // Tumour-only DeepVariant germline calls are adjudicated only by DeepSomatic's verdict; warn once if it is not run + if (params.germline_var_keep.contains('deepvariant') && !params.somatic_var_keep.contains('deepsomatic')) { + ch_samplesheet + .filter { meta, _bams -> !meta.paired_data } + .first() + .subscribe { meta, _bams -> + log.warn("Tumour-only samples (e.g. ${meta.id}) use DeepVariant germline calls without DeepSomatic's verdict, so they may include somatic variants. Add 'deepsomatic' to --somatic_var_keep to filter them.") + } + } // CHM13 has no ascat_loci_rt attribute, so the built set is GC-only by construction build_clairsto_cna = clairsto_cna_dir == null && params.genome == 'CHM13' && params.skip_ascat + // CoRAL seeds from ASCAT's copy number, so it cannot run without it + if (!params.skip_coral && params.skip_ascat) { + error("CoRAL seeds from ASCAT's copy-number segments. Remove --skip_ascat, or add --skip_coral.") + } + // Gurobi needs a licence file; SCIP, the default, needs nothing + if (!params.skip_coral && params.coral_solver == 'gurobi_direct' && !params.gurobi_license) { + error("--coral_solver gurobi_direct needs a licence: pass --gurobi_license , or use the default --coral_solver scip.") + } + // A missing set would leave the join below waiting forever, so CLAIRSTO would silently never run if (build_clairsto_cna) { def missing_ascat = ['ascat_alleles': params.ascat_allele_files, @@ -649,6 +672,8 @@ workflow LRSOMATIC { ch_ascat_files = channel.empty() ascat_tumoronly_ch = channel.empty() + ch_ascat_cnvs = channel.empty() + ch_ecdna_tumor_bam = channel.empty() if (!params.skip_ascat) { branched_minimap.tumor_only @@ -696,11 +721,23 @@ workflow LRSOMATIC { // ascat_tumoronly_ch: [meta, purityploidy, segments] // All ASCAT files per sample for the report module, which globs by suffix - ch_ascat_files = ASCAT.out.segments_raw - .mix(ASCAT.out.purityploidy, ASCAT.out.png) - .groupTuple() - .map { meta, files -> [meta, files.flatten()] } + // Joined per sample; segments_raw is optional, so it arrives as null when absent + ch_ascat_files = ASCAT.out.purityploidy + .join(ASCAT.out.png) + .join(ASCAT.out.segments_raw, remainder: true) + .map { meta, purityploidy, png, segments_raw -> [meta, [purityploidy, png, segments_raw ?: []].flatten()] } // ch_ascat_files: [meta, [file, file, ...]] + + // CoRAL seeds from ASCAT's total copy number. Built from ascat_ch so the meta + // key matches ASCAT's output exactly and the joins in ECDNA pair. + ch_ascat_cnvs = ASCAT.out.cnvs + // ch_ascat_cnvs: [meta, cnvs_txt] + + ch_ecdna_tumor_bam = ascat_ch + .map { meta, _normal_bam, _normal_bai, tumor_bam, tumor_bai -> + return [meta, tumor_bam, tumor_bai] + } + // ch_ecdna_tumor_bam: [meta, tumor_bam, tumor_bai] } // SUBWORKFLOW: TUMORONLY_SMALLVAR @@ -1255,13 +1292,51 @@ workflow LRSOMATIC { ) // The WAKHAN outputs the report renders: ranked solutions, heatmap, per-solution plots + // Joined per sample; all three outputs are required ch_wakhan_files = WAKHAN.out.solutions_ranks - .mix(WAKHAN.out.heatmap_html, WAKHAN.out.solution_dirs) - .groupTuple() - .map { meta, files -> [meta, files.flatten()] } // solution_dirs contributes a list + .join(WAKHAN.out.heatmap_html) + .join(WAKHAN.out.solution_dirs) + .map { meta, ranks, heatmap, dirs -> [meta, [ranks, heatmap, dirs].flatten()] } // dirs may be a list // ch_wakhan_files: [meta, [file_or_dir, ...]] } + // + // SUBWORKFLOW: ECDNA -- CoRAL amplicon reconstruction, then AmpliconClassifier + // Input: ch_ecdna_tumor_bam -- [meta, tumor_bam, tumor_bai] + // ch_ascat_cnvs -- [meta, cnvs_txt] + // ch_fai -- [[:], fai] + // Output: .classification -- [meta, tsv] -- amplicon_classification_profiles.tsv + // + + ch_ecdna_classification = channel.empty() + + if (!params.skip_coral) { + // The data repo is only fetched when the classifier will actually use it + ch_aa_data_repo = channel.empty() + if (!params.skip_ampliconclassifier) { + // Overriding the repo drops the default MD5, which belongs to the published tarball + PREPARE_AA_DATA_REPO ( + params.aa_data_repo, + params.aa_data_repo_url, + params.aa_data_repo_md5 ?: (params.aa_data_repo ? null : getGenomeAttribute('aa_data_repo_md5')) + ) + ch_aa_data_repo = PREPARE_AA_DATA_REPO.out.data_repo + ch_versions = ch_versions.mix(PREPARE_AA_DATA_REPO.out.versions) + } + + ECDNA ( + ch_ecdna_tumor_bam, + ch_ascat_cnvs, + ch_fai, + ch_aa_data_repo, + params.coral_ref, + params.ac_ref + ) + + ch_ecdna_classification = ECDNA.out.classification + // ch_ecdna_classification: [meta, tsv] + } + // // MODULE: LRSOMATICREPORT -- per-sample HTML report; all inputs optional, so joins use remainder: true on the tumor id // @@ -1276,13 +1351,19 @@ workflow LRSOMATIC { .set { report_id_meta } // report_id_meta: [id, meta] - ch_somatic_vep_vcf - .map { meta, vcf -> [meta.id, vcf] } - .set { report_vep_ch } + // A skipped module leaves an empty leg, and remainder: true then defers every sample to + // channel close. One [] per sample keeps each leg matched so samples report independently. + def report_empty_slot = { -> report_id_meta.map { id, _meta -> [id, []] } } - ch_sv_vep_vcf - .map { meta, vcf -> [meta.id, vcf] } - .set { report_sv_vep_ch } + def report_vep_ch = params.skip_vep + ? report_empty_slot.call() + : ch_somatic_vep_vcf.map { meta, vcf -> [meta.id, vcf] } + // report_vep_ch: [id, vcf] + + def report_sv_vep_ch = params.skip_vep + ? report_empty_slot.call() + : ch_sv_vep_vcf.map { meta, vcf -> [meta.id, vcf] } + // report_sv_vep_ch: [id, vcf] SEVERUS.out.somatic_vcf .map { meta, vcf -> [meta.id, vcf] } @@ -1292,29 +1373,63 @@ workflow LRSOMATIC { .map { meta, vcf, _tbi -> [meta.id, vcf] } .set { report_somatic_ch } - ch_ascat_files - .map { meta, files -> [meta.id, files] } - .set { report_ascat_ch } + def report_ascat_ch = params.skip_ascat + ? report_empty_slot.call() + : ch_ascat_files.map { meta, files -> [meta.id, files] } + // report_ascat_ch: [id, [files]] - ch_wakhan_files - .map { meta, files -> [meta.id, files] } - .set { report_wakhan_ch } + def report_wakhan_ch = params.skip_wakhan + ? report_empty_slot.call() + : ch_wakhan_files.map { meta, files -> [meta.id, files] } + // report_wakhan_ch: [id, [files]] + + // One emission per sample per tool, none optional: mosdepth 2, cramino 1, samtools 2. + // Adding another per-sample emission to either mix below must bump this count. + def qc_files_per_sample = params.skip_qc + ? 0 + : (params.skip_mosdepth ? 0 : 2) + (params.skip_cramino ? 0 : 1) + (params.skip_bamstats ? 0 : 2) // Tumor-side QC, keyed by the sample id (= report id) + // groupKey: emit a sample's bundle on its own files; toString() restores a plain String key ch_mosdepth_summary .mix(ch_mosdepth_global, ch_cramino_post_txt, ch_bam_stats, ch_bam_flagstat) .filter { meta, _f -> meta.type == 'tumor' } - .map { meta, f -> [meta.id, f] } + .map { meta, f -> [groupKey(meta.id, qc_files_per_sample), f] } .groupTuple() - .set { report_qc_tumor_ch } + .map { key, files -> [key.toString(), files] } + .set { report_qc_tumor_grouped } + // report_qc_tumor_grouped: [id, [qc_file, ...]] // Normal-side QC (matched mode): a pair shares meta.id, so already keyed by the report id ch_mosdepth_summary .mix(ch_mosdepth_global, ch_cramino_post_txt, ch_bam_stats, ch_bam_flagstat) .filter { meta, _f -> meta.type == 'normal' } - .map { meta, f -> [meta.id, f] } + .map { meta, f -> [groupKey(meta.id, qc_files_per_sample), f] } .groupTuple() - .set { report_qc_normal_ch } + .map { key, files -> [key.toString(), files] } + .set { report_qc_normal_grouped } + // report_qc_normal_grouped: [id, [qc_file, ...]] -- paired samples only + + def report_qc_tumor_ch = qc_files_per_sample == 0 + ? report_empty_slot.call() + : report_qc_tumor_grouped + + // Normal-side QC covers paired samples only; meta.paired_data gives the tumor-only arm + // its [] up front instead of waiting out channel close for a match that never arrives. + report_id_meta + .branch { _id, meta -> + paired: meta.paired_data + tumor_only: true + } + .set { report_roster } + + def report_qc_normal_ch = qc_files_per_sample == 0 + ? report_empty_slot.call() + : report_roster.paired + .join(report_qc_normal_grouped) + .map { id, _meta, files -> [id, files] } + .mix(report_roster.tumor_only.map { id, _meta -> [id, []] }) + // report_qc_normal_ch: [id, [files] | []] -- full roster report_id_meta .join(report_vep_ch, remainder: true)