From 30ac555b7cd033ec50b1940109ff57fcca27004d Mon Sep 17 00:00:00 2001 From: robert-a-forsyth Date: Wed, 14 Oct 2026 15:45:41 +0200 Subject: [PATCH 1/5] Express SAVANA _2 breakends from the _1 side before collapsing mates SAVANA writes BP_NOTATION in (_1, _2) order on both mate records (all 12,679 _2 records of B2101426 are reversed relative to their own ALT). Mates were collapsed keeping whichever record came first in the file, so a junction whose _2 record came first paired the second breakend's position with the first breakend's orientation. ReConPlot's tables for the fixtures are unchanged; the snapshot moves only with the wrapper VERSION in versions.yml. Co-Authored-By: Claude Opus 5.5 (1M context) --- assets/reconplot/R/parsers/savana.R | 3 +++ assets/reconplot/README.md | 2 +- assets/reconplot/VERSION | 2 +- modules/local/reconplot/tests/main.nf.test.snap | 10 +++++----- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/assets/reconplot/R/parsers/savana.R b/assets/reconplot/R/parsers/savana.R index e6e094be..145f32d3 100644 --- a/assets/reconplot/R/parsers/savana.R +++ b/assets/reconplot/R/parsers/savana.R @@ -149,6 +149,9 @@ savana_read_sv_vcf <- function(file, min_support = 0, min_svlen = 0, pass_only = support = support, stringsAsFactors = FALSE ) + ## BP_NOTATION is written in (_1, _2) order on both mates: put a _2 record's mate first so positions match it + mate2 <- grepl("_2$", vcf$id) & !is.na(mate$chr) + if (any(mate2)) sv[mate2, c("chr1", "pos1", "chr2", "pos2")] <- sv[mate2, c("chr2", "pos2", "chr1", "pos1")] sv <- dedupe_breakend_pairs(sv) # BND records come in mate pairs savana_filter_sv(sv, min_support, min_svlen) } diff --git a/assets/reconplot/README.md b/assets/reconplot/README.md index 91a294ae..e36f19d7 100644 --- a/assets/reconplot/README.md +++ b/assets/reconplot/README.md @@ -14,4 +14,4 @@ applies to SAVANA too; `--genes` symbols are resolved against ReConPlot's gene t skipped with a warning instead of losing the panel); a failing BAF-track reader skips the track instead of failing the run; `--genome` is validated up front; thousands separators are accepted in `--regions`; malformed and duplicate regions are dropped; the Wakhan BED header is found with or without the comment block; `--sv-format` defaults to `auto`, which reads -SAVANA's classified somatic VCF when no BEDPE is staged (SAVANA's PacBio classifier writes no BEDPE). +SAVANA's classified somatic VCF when no BEDPE is staged (SAVANA's PacBio classifier writes no BEDPE). SAVANA writes `BP_NOTATION` in (_1, _2) order on both mate records, so a `_2` record is re-expressed from its mate's side before mates are collapsed (positions then match the orientation). diff --git a/assets/reconplot/VERSION b/assets/reconplot/VERSION index d3d30822..b438414a 100644 --- a/assets/reconplot/VERSION +++ b/assets/reconplot/VERSION @@ -1 +1 @@ -2672ccb+lrsomatic.2026-10-07 +2672ccb+lrsomatic.2026-10-14 diff --git a/modules/local/reconplot/tests/main.nf.test.snap b/modules/local/reconplot/tests/main.nf.test.snap index 9eff3e6f..8fc3d8a5 100644 --- a/modules/local/reconplot/tests/main.nf.test.snap +++ b/modules/local/reconplot/tests/main.nf.test.snap @@ -134,7 +134,7 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,a8129e310cb03d150aa22f9a6ec0421d" ] ], "timestamp": "2026-10-07T11:27:41.013206989", @@ -157,7 +157,7 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,a8129e310cb03d150aa22f9a6ec0421d" ] ], "timestamp": "2026-10-07T11:27:53.7954896", @@ -180,7 +180,7 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,a8129e310cb03d150aa22f9a6ec0421d" ] ], "timestamp": "2026-10-07T11:28:15.015932537", @@ -203,7 +203,7 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,a8129e310cb03d150aa22f9a6ec0421d" ] ], "timestamp": "2026-10-07T11:28:35.052455573", @@ -226,7 +226,7 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,a8129e310cb03d150aa22f9a6ec0421d" ] ], "timestamp": "2026-10-07T11:28:24.926693238", From 0cecbfe6eca54585f324073b121d0ac67c07ee21 Mon Sep 17 00:00:00 2001 From: robert-a-forsyth Date: Wed, 14 Oct 2026 15:45:41 +0200 Subject: [PATCH 2/5] Add ShatterSeek chromothripsis calling with PCAWG confidence tiers SHATTERSEEK_CHROMOTHRIPSIS runs ShatterSeek for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), reading the callers with the ReConPlot wrapper's parsers, and labels each chromosome's cluster with the PCAWG HC1/HC2/HC3/LC1 tiers (BH q at 0.2). The module uses the robert-a-forsyth/ShatterSeek fork (chm13 support, statistics corrected to the PCAWG definitions, one cluster per chromosome) in docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c. New parameters: --skip_shatterseek, --shatterseek_genome (inferred from --genome) and --shatterseek_min_svlen (default 150 bp). Synthetic fixtures in tests/fixtures/shatterseek cover every caller layout; module, subworkflow and pipeline-wiring tests are tagged small. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 2 + CITATIONS.md | 4 + README.md | 14 +- assets/shatterseek/README.md | 22 ++ assets/shatterseek/VERSION | 1 + assets/shatterseek/run_shatterseek.R | 181 +++++++++++ conf/modules.config | 17 + conf/test.config | 3 +- conf/test_full.config | 3 +- containers/shatterseek/Dockerfile | 24 ++ containers/shatterseek/README.md | 29 ++ containers/shatterseek/environment.yml | 20 ++ docs/output.md | 74 ++++- docs/usage.md | 23 +- modules/local/shatterseek/main.nf | 71 +++++ modules/local/shatterseek/meta.yml | 103 ++++++ modules/local/shatterseek/tests/main.nf.test | 208 ++++++++++++ .../local/shatterseek/tests/main.nf.test.snap | 301 ++++++++++++++++++ .../local/shatterseek/tests/nextflow.config | 5 + nextflow.config | 5 + nextflow_schema.json | 29 ++ .../local/shatterseek_chromothripsis.nf | 85 +++++ .../tests/shatterseek_chromothripsis.nf.test | 82 +++++ .../utils_nfcore_lrsomatic_pipeline/main.nf | 15 + .../shatterseek/ascat/test.segments.txt | 15 + tests/fixtures/shatterseek/make_fixtures.py | 127 ++++++++ .../savana/test_tumor.classified.somatic.vcf | 50 +++ ...t_tumor_segmented_absolute_copy_number.tsv | 15 + .../severus/severus_somatic.vcf.gz | Bin 0 -> 1027 bytes .../solution_rank_1/integer_profile.bed | 22 ++ tests/sv_annotation_compose.nf.test | 38 ++- workflows/lrsomatic.nf | 17 + 32 files changed, 1583 insertions(+), 22 deletions(-) create mode 100644 assets/shatterseek/README.md create mode 100644 assets/shatterseek/VERSION create mode 100644 assets/shatterseek/run_shatterseek.R create mode 100644 containers/shatterseek/Dockerfile create mode 100644 containers/shatterseek/README.md create mode 100644 containers/shatterseek/environment.yml create mode 100644 modules/local/shatterseek/main.nf create mode 100644 modules/local/shatterseek/meta.yml create mode 100644 modules/local/shatterseek/tests/main.nf.test create mode 100644 modules/local/shatterseek/tests/main.nf.test.snap create mode 100644 modules/local/shatterseek/tests/nextflow.config create mode 100644 subworkflows/local/shatterseek_chromothripsis.nf create mode 100644 subworkflows/local/tests/shatterseek_chromothripsis.nf.test create mode 100644 tests/fixtures/shatterseek/ascat/test.segments.txt create mode 100644 tests/fixtures/shatterseek/make_fixtures.py create mode 100644 tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf create mode 100644 tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv create mode 100644 tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz create mode 100644 tests/fixtures/shatterseek/wakhan/solution_rank_1/integer_profile.bed diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ea96028..5d4c408d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Added` +- Added ShatterSeek chromothripsis calling for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), published to `/shatterseek//` with a per-chromosome summary, the exact ShatterSeek inputs and a figure per called chromosome; each cluster is labelled with the PCAWG high/low-confidence tiers (HC1, HC2, HC3, LC1; Cortes-Ciriano et al. 2020, ShatterSeek issue #22), with BH q-values at 0.2 and cut-offs overridable through `ext.args`. The pipeline runs the fork [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), which adds T2T-CHM13 support and corrects the statistical criteria to the PCAWG definitions (translocations selected on the right breakend, intra-SV fragment-joins test, all-breakpoint one-sided enrichment, exponential tests on inter-breakpoint gaps, one cluster per chromosome, the strand swap in `SVs()`), inside `docker.io/robertaforsyth/shatterseek` (recipe in `containers/shatterseek/`); the callers are read by the ReConPlot wrapper's parsers. `--skip_shatterseek`, `--shatterseek_genome` and `--shatterseek_min_svlen` (intra-chromosomal SVs shorter than 150 bp are left out by default) control it (@robert-a-forsyth). - [#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). @@ -58,6 +59,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Fixed` +- The ReConPlot wrapper's SAVANA VCF parser re-expresses a `_2` breakend from its mate's side before collapsing mates: SAVANA writes `BP_NOTATION` in (_1, _2) order on both records, so a junction whose `_2` record came first in the file paired the second breakend's position with the first breakend's orientation (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `PREPARE_AA_DATA_REPO` no longer reads a `versions` output `UNTAR` does not have, which crashed every default GRCh38 run with CoRAL; overriding `--aa_data_repo_url` drops the genome's default MD5 (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `-profile test_full` starts again: it skips CoRAL and signatures and carries the chr19 SAVANA workarounds from `test` (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `CLAIRSTO_VERDICT_TAG` no longer crashes when ASCAT finds no solution and writes an empty segments table (@robert-a-forsyth). diff --git a/CITATIONS.md b/CITATIONS.md index 068c997e..216a60f0 100644 --- a/CITATIONS.md +++ b/CITATIONS.md @@ -122,6 +122,10 @@ > Espejo Valle-Inclán J, Cortés-Ciriano I. ReConPlot: an R package for the visualization and interpretation of genomic rearrangements. Bioinformatics. 2023 Dec 1;39(12):btad719. doi: 10.1093/bioinformatics/btad719. PubMed PMID: 38058190; PubMed Central PMCID: PMC10710371. +- [ShatterSeek](https://pubmed.ncbi.nlm.nih.gov/32025003/) + + > Cortés-Ciriano I, Lee JJ, Xi R, Jain D, Jung YL, Yang L, Gordenin D, Klimczak LJ, Zhang CZ, Pellman DS; PCAWG Structural Variation Working Group; Park PJ; PCAWG Consortium. Comprehensive analysis of chromothripsis in 2,658 human cancers using whole-genome sequencing. Nat Genet. 2020 Mar;52(3):331-341. doi: 10.1038/s41588-019-0576-7. PubMed PMID: 32025003; PubMed Central PMCID: PMC7058534. + - [RepeatMasker](https://www.repeatmasker.org/) > Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0 [Software]. https://www.repeatmasker.org/ diff --git a/README.md b/README.md index 761ae695..75210b44 100644 --- a/README.md +++ b/README.md @@ -72,6 +72,10 @@ c. Somatic SV and CNA functional annotation ([`Padfoot`](https://github.com/Kolm a. Rearrangement and copy-number figures per CN/SV caller pair ([`ReConPlot`](https://github.com/cortes-ciriano-lab/ReConPlot)) +**6) Chromothripsis:** + +a. Chromothripsis calls with PCAWG confidence tiers per CN/SV caller pair ([`ShatterSeek`](https://github.com/parklab/ShatterSeek)) + @@ -114,7 +118,7 @@ IntGenomicsLab/lr_somatic was originally written by Luuk Harbers, Robert Forsyth The main output is an aligned and phased tumour BAM, per-sample QC from `cramino`, `mosdepth`, `samtools` and optionally `fibertools`, a MultiQC report, and a self-contained per-sample HTML report (`/report/_report.html`; disable it with `--skip_report`). -Variant and copy number callers (`clairS`, `clairS-TO`, `severus`, `savana`, `ascat`) write to their own folders; Padfoot annotations and ReConPlot figures go to `padfoot/` and `reconplot/`, one subfolder per caller pair. See the [output documentation](/docs/output.md). +Variant and copy number callers (`clairS`, `clairS-TO`, `severus`, `savana`, `ascat`) write to their own folders; Padfoot annotations, ReConPlot figures and ShatterSeek chromothripsis calls go to `padfoot/`, `reconplot/` and `shatterseek/`, one subfolder per caller pair. See the [output documentation](/docs/output.md). Example output directory structure: @@ -136,6 +140,10 @@ Example output directory structure: │ │ ├── savana │ │ ├── severus_ascat │ │ └── severus_wakhan +│ ├── shatterseek +│ │ ├── savana +│ │ ├── severus_ascat +│ │ └── severus_wakhan │ ├── variants │ │ ├──clairS-TO │ │ ├──severus @@ -169,6 +177,10 @@ Example output directory structure: │ │ ├── savana │ │ ├── severus_ascat │ │ └── severus_wakhan +│ ├── shatterseek +│ │ ├── savana +│ │ ├── severus_ascat +│ │ └── severus_wakhan │ ├── variants │ │ ├── clair3 │ │ ├── clairS diff --git a/assets/shatterseek/README.md b/assets/shatterseek/README.md new file mode 100644 index 00000000..23f5fb28 --- /dev/null +++ b/assets/shatterseek/README.md @@ -0,0 +1,22 @@ +# ShatterSeek wrapper + +`run_shatterseek.R` turns lrsomatic caller output into ShatterSeek input, runs ShatterSeek, labels each chromosome's +cluster with the PCAWG confidence tiers and draws a figure for every called chromosome. The `SHATTERSEEK` module +stages this directory, plus `assets/reconplot/` whose `R/parsers/` read ASCAT, Wakhan, Severus and SAVANA, so +ShatterSeek and ReConPlot see the same SVs and segments. + +ShatterSeek is the [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support) +fork inside the module's image (recipe in `containers/shatterseek/`): T2T-CHM13 support, the statistical criteria fixed +to the PCAWG definitions, one cluster per chromosome, and `classify_chromothripsis()`. + +Conversion to ShatterSeek's format: + +- SVs: insertions, single breakends and anything without two `+`/`-` strands are dropped, as are SVs with a breakend + outside chromosomes 1-22 and X. Intra-chromosomal SVs are typed by strand (`+-` DEL, `-+` DUP, `++` h2hINV, + `--` t2tINV), inter-chromosomal ones are TRA. Breakends are ordered so pos1 <= pos2, strands moving with them. +- Copy number: total copy number rounded to an integer (SAVANA's is fractional), chromosomes 1-22 and X, overlaps + trimmed, and neighbouring segments with the same copy number merged (ShatterSeek counts copy-number states row by row). + +Both tables are written next to the summary (`.shatterseek_{sv,cn}_input.tsv`). Tier cut-offs default to PCAWG's +and can be changed with `--use`, `--alpha`, `--min-sv-hc1`, `--min-osc-hc`, `--min-sv-hc2`, `--min-tra-hc2`, +`--min-sv-hc3`, `--min-osc-lc`, `--max-osc-lc` and `--no-hc3`; see `Rscript run_shatterseek.R --help`. diff --git a/assets/shatterseek/VERSION b/assets/shatterseek/VERSION new file mode 100644 index 00000000..3eefcb9d --- /dev/null +++ b/assets/shatterseek/VERSION @@ -0,0 +1 @@ +1.0.0 diff --git a/assets/shatterseek/run_shatterseek.R b/assets/shatterseek/run_shatterseek.R new file mode 100644 index 00000000..65dc5afe --- /dev/null +++ b/assets/shatterseek/run_shatterseek.R @@ -0,0 +1,181 @@ +#!/usr/bin/env Rscript +## --------------------------------------------------------------------------- +## run_shatterseek.R -- ShatterSeek chromothripsis calls from lrsomatic caller output. +## +## --[ReConPlot wrapper parsers]--> {cn, sv} --[ShatterSeek format]--> shatterseek() +## --> classify_chromothripsis() (PCAWG HC1/HC2/HC3/LC1) --> summary TSV + per-chromosome figures +## +## Callers are read by the parsers in assets/reconplot/R/parsers (--reconplot-dir), so ShatterSeek and +## ReConPlot see the same SVs and segments. ShatterSeek itself is the robert-a-forsyth fork (chm13 support, +## corrected statistics) installed in the module's container. +## --------------------------------------------------------------------------- + +suppressPackageStartupMessages({ + library(optparse) + library(data.table) +}) + +option_list <- list( + make_option("--reconplot-dir", type = "character", default = NULL, + help = "ReConPlot wrapper directory whose R/ parsers read the caller output (required)"), + make_option("--cn-source", type = "character", default = NULL, help = "Copy number parser: ascat, wakhan or savana"), + make_option("--sv-source", type = "character", default = NULL, help = "SV parser: severus or savana"), + make_option("--cn-input", type = "character", default = NULL, help = "Directory with the CN caller's files"), + make_option("--sv-input", type = "character", default = NULL, help = "Directory with the SV caller's files [default: --cn-input]"), + make_option("--sample", type = "character", default = NULL, help = "Sample name (also pins caller files)"), + make_option("--prefix", type = "character", default = NULL, help = "Output file prefix [default: --sample]"), + make_option("--outdir", type = "character", default = ".", help = "Output directory [default %default]"), + make_option("--genome", type = "character", default = "hg38", help = "hg38, chm13 or hg19 [default %default]"), + make_option("--min-support", type = "double", default = 0, help = "Drop SVs with tumour read support below this [default %default]"), + make_option("--min-svlen", type = "double", default = 0, + help = "Drop intra-chromosomal SVs shorter than this many bp [default %default]"), + make_option("--exclude-vntr", action = "store_true", default = FALSE, help = "Drop Severus SVs inside a VNTR"), + make_option("--sv-format", type = "character", default = "auto", help = "SAVANA SV input: bedpe, vcf or auto [default %default]"), + make_option("--ascat-cn-mode", type = "character", default = "segments", + help = "ASCAT table: segments, raw or raw-fractional [default %default]"), + make_option("--min-size", type = "integer", default = 1, help = "shatterseek() min.Size [default %default]"), + make_option("--use", type = "character", default = "q", help = "Tier tests on BH q (q) or raw P values (p) [default %default]"), + make_option("--alpha", type = "double", default = 0.2, help = "Significance level for the tiers [default %default]"), + make_option("--min-sv-hc1", type = "double", default = 6, help = "HC1/LC1 interleaved intra SVs [default %default]"), + make_option("--min-osc-hc", type = "double", default = 7, help = "HC1/HC2 oscillating segments [default %default]"), + make_option("--min-sv-hc2", type = "double", default = 3, help = "HC2 interleaved intra SVs [default %default]"), + make_option("--min-tra-hc2", type = "double", default = 4, help = "HC2 translocations [default %default]"), + make_option("--min-sv-hc3", type = "double", default = 40, help = "HC3 SVs incl. translocations [default %default]"), + make_option("--min-osc-lc", type = "double", default = 4, help = "LC1 minimum oscillating segments [default %default]"), + make_option("--max-osc-lc", type = "double", default = 6, help = "LC1 maximum oscillating segments [default %default]"), + make_option("--no-hc3", action = "store_true", default = FALSE, help = "Skip the >= 40 SV tier (PCAWG: sarcomas only)"), + make_option("--plot-all", action = "store_true", default = FALSE, + help = "Draw every chromosome with a cluster, not only those with a tier") +) +args <- parse_args(OptionParser(usage = "%prog --reconplot-dir DIR --cn-source SRC --sv-source SRC --cn-input DIR [options]", + option_list = option_list)) +names(args) <- gsub("-", "_", names(args)) + +if (is.null(args$reconplot_dir)) stop("--reconplot-dir is required") +for (f in c("utils_common.R", "utils_vcf.R", file.path("parsers", "registry.R"))) source(file.path(args$reconplot_dir, "R", f)) +load_parsers(file.path(args$reconplot_dir, "R", "parsers")) +suppressPackageStartupMessages(library(ShatterSeek)) + +args$genome <- match.arg(args$genome, c("hg38", "chm13", "hg19")) +args$use <- match.arg(args$use, c("q", "p")) +args$sv_input <- args$sv_input %||% args$cn_input +prefix <- args$prefix %||% args$sample %||% "sample" +dir.create(args$outdir, showWarnings = FALSE, recursive = TRUE) + +ANALYSED <- c(as.character(1:22), "X") # ShatterSeek's chromosomes (Ensembl names) +STRAND_PAIRS <- c("+-", "-+", "++", "--") +SVTYPE_OF <- c("+-" = "DEL", "-+" = "DUP", "++" = "h2hINV", "--" = "t2tINV") + +## --- 1. parse with the ReConPlot parsers ------------------------------------- +parsed <- run_parsers(args) +cn <- validate_cn(parsed$cn) +sv <- validate_sv(parsed$sv, interchrom_as_tra = FALSE) # keep each translocation's strands + +## --- 2. SVs in ShatterSeek's format -------------------------------------------- +to_shatterseek_sv <- function(sv) { + if (nrow(sv) == 0) { + return(data.frame(chrom1 = character(), pos1 = integer(), chrom2 = character(), pos2 = integer(), + SVtype = character(), strand1 = character(), strand2 = character())) + } + strands <- sv$strands + by_class <- c(DEL = "+-", DUP = "-+", h2hINV = "++", t2tINV = "--") # generic parser may give classes + strands[strands %in% names(by_class)] <- by_class[strands[strands %in% names(by_class)]] + keep <- strands %in% STRAND_PAIRS + if (any(!keep)) { + tab <- table(strands[!keep]) + log_msg(" dropped ", sum(!keep), " SVs without a two-breakend orientation (", + paste(sprintf("%s=%d", names(tab), as.integer(tab)), collapse = ", "), ")") + } + out <- data.frame(chrom1 = sub("^chr", "", sv$chr1[keep]), pos1 = as.integer(sv$pos1[keep]), + chrom2 = sub("^chr", "", sv$chr2[keep]), pos2 = as.integer(sv$pos2[keep]), + strand1 = substr(strands[keep], 1, 1), strand2 = substr(strands[keep], 2, 2), + stringsAsFactors = FALSE) + analysed <- out$chrom1 %in% ANALYSED & out$chrom2 %in% ANALYSED + if (any(!analysed)) log_msg(" dropped ", sum(!analysed), " SVs with a breakend outside chromosomes 1-22 and X") + out <- out[analysed, , drop = FALSE] + # pos1 <= pos2, strands moving with their breakends (the order SVs() would impose) + flip <- out$pos1 > out$pos2 + out[flip, c("chrom1", "pos1", "strand1", "chrom2", "pos2", "strand2")] <- + out[flip, c("chrom2", "pos2", "strand2", "chrom1", "pos1", "strand1")] + out$SVtype <- ifelse(out$chrom1 == out$chrom2, unname(SVTYPE_OF[paste0(out$strand1, out$strand2)]), "TRA") + out[order(match(out$chrom1, ANALYSED), out$pos1), c("chrom1", "pos1", "chrom2", "pos2", "SVtype", "strand1", "strand2")] +} + +## --- 3. copy number in ShatterSeek's format ----------------------------------- +to_shatterseek_cn <- function(cn) { + d <- data.frame(chrom = sub("^chr", "", cn$chr), start = as.integer(cn$start), end = as.integer(cn$end), + total_cn = round(as.numeric(cn$copyNumber)), stringsAsFactors = FALSE) + if (any(cn$copyNumber != d$total_cn, na.rm = TRUE)) { + log_msg(" rounded ", sum(cn$copyNumber != d$total_cn, na.rm = TRUE), " non-integer total copy numbers") + } + d <- d[d$chrom %in% ANALYSED & !is.na(d$total_cn), , drop = FALSE] + d <- d[order(match(d$chrom, ANALYSED), d$start, d$end), , drop = FALSE] + out <- list() + for (ch in unique(d$chrom)) { + x <- d[d$chrom == ch, , drop = FALSE] + # no overlaps: each segment starts after the previous one ends + for (i in seq_len(nrow(x))[-1]) if (x$start[i] <= x$end[i - 1]) x$start[i] <- x$end[i - 1] + 1L + x <- x[x$start <= x$end, , drop = FALSE] + # neighbours with equal copy number become one segment (ShatterSeek counts states row by row) + run <- cumsum(c(TRUE, diff(x$total_cn) != 0)) + out[[ch]] <- data.frame(chrom = ch, start = tapply(x$start, run, min), end = tapply(x$end, run, max), + total_cn = tapply(x$total_cn, run, `[`, 1), stringsAsFactors = FALSE) + } + out <- do.call(rbind, out) + rownames(out) <- NULL + if (is.null(out) || nrow(out) == 0) stop("No copy number segments on chromosomes 1-22 and X") + out +} + +sv_ss <- to_shatterseek_sv(sv) +cn_ss <- to_shatterseek_cn(cn) +log_msg(sprintf("ShatterSeek input: %d SVs (%s), %d CN segments after merging equal neighbours", + nrow(sv_ss), paste(sprintf("%s=%d", names(table(sv_ss$SVtype)), as.integer(table(sv_ss$SVtype))), collapse = " "), + nrow(cn_ss))) +fwrite(sv_ss, file.path(args$outdir, paste0(prefix, ".shatterseek_sv_input.tsv")), sep = "\t") +fwrite(cn_ss, file.path(args$outdir, paste0(prefix, ".shatterseek_cn_input.tsv")), sep = "\t") + +## --- 4. ShatterSeek + PCAWG tiers --------------------------------------------- +res <- shatterseek( + SV.sample = SVs(chrom1 = sv_ss$chrom1, pos1 = as.numeric(sv_ss$pos1), chrom2 = sv_ss$chrom2, pos2 = as.numeric(sv_ss$pos2), + SVtype = sv_ss$SVtype, strand1 = sv_ss$strand1, strand2 = sv_ss$strand2), + seg.sample = CNVsegs(chrom = cn_ss$chrom, start = as.numeric(cn_ss$start), end = as.numeric(cn_ss$end), + total_cn = as.numeric(cn_ss$total_cn)), + min.Size = args$min_size, genome = args$genome) +res@chromSummary <- classify_chromothripsis( + res@chromSummary, use = args$use, alpha = args$alpha, min_sv_hc1 = args$min_sv_hc1, min_osc_hc = args$min_osc_hc, + min_sv_hc2 = args$min_sv_hc2, min_tra_hc2 = args$min_tra_hc2, min_sv_hc3 = args$min_sv_hc3, + min_osc_lc = args$min_osc_lc, max_osc_lc = args$max_osc_lc, hc3 = !isTRUE(args$no_hc3)) + +summary <- res@chromSummary +summary$inter_other_chroms_coords_all <- trimws(gsub("[\n;]+", ";", summary$inter_other_chroms_coords_all)) +out <- cbind(data.frame(sample = args$sample %||% prefix, cn_source = args$cn_source, sv_source = args$sv_source, + stringsAsFactors = FALSE), summary) +out$chrom <- paste0("chr", out$chrom) +out$start <- as.integer(out$start) +out$end <- as.integer(out$end) +fwrite(out, file.path(args$outdir, paste0(prefix, ".shatterseek.tsv")), sep = "\t", na = "NA", quote = FALSE) +called <- summary[summary$chromothripsis_confidence != "none", , drop = FALSE] +log_msg(sprintf("tiers: %s", if (nrow(called)) paste(sprintf("chr%s %s (%s)", called$chrom, called$chromothripsis_confidence, + called$chromothripsis_criteria), collapse = "; ") else "none")) + +## --- 5. figures ----------------------------------------------------------------- +draw <- summary$chrom[if (isTRUE(args$plot_all)) !is.na(summary$start) else summary$chromothripsis_confidence != "none"] +if (length(draw)) dir.create(file.path(args$outdir, "plots"), showWarnings = FALSE) +for (ch in draw) { + f <- file.path(args$outdir, "plots", sprintf("%s_chr%s.pdf", prefix, ch)) + ok <- tryCatch({ + p <- plot_chromothripsis(ShatterSeek_output = res, chr = ch, sample_name = prefix, genome = args$genome) + pdf(f, width = 9, height = 11) + gridExtra::grid.arrange(p[[1]], p[[2]], p[[3]], p[[4]], nrow = 4, ncol = 1, heights = c(0.2, 0.4, 0.4, 0.4)) + dev.off() + TRUE + }, error = function(e) { + if (dev.cur() > 1) dev.off() + unlink(f) + log_msg("ERROR: figure for chr", ch, " not drawn: ", conditionMessage(e)) + FALSE + }) + if (ok) log_msg(" drew ", basename(f)) +} +log_msg("done") diff --git a/conf/modules.config b/conf/modules.config index d56f3ffe..c6cdda35 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -1010,6 +1010,23 @@ process { ] } + // + // Chromothripsis + // + + withName : '.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + ext.args = { params.shatterseek_min_svlen ? "--min-svlen ${params.shatterseek_min_svlen}" : '' } + // Downstream of every caller: a sample ShatterSeek cannot process must not stop the run. + // Resource kills are retried as in conf/base.config first. + errorStrategy = { task.exitStatus in ((130..145) + 104 + (175..177)) && task.attempt <= task.maxRetries ? 'retry' : 'ignore' } + publishDir = [ + path: { "${params.outdir}/${meta.id}/shatterseek" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + ///////// // VEP // ///////// diff --git a/conf/test.config b/conf/test.config index ce2cb581..f7587f44 100644 --- a/conf/test.config +++ b/conf/test.config @@ -72,9 +72,10 @@ params { // CoRAL seeds from ASCAT, which is skipped here, and the chr19 slice has no amplification skip_coral = true skip_modkit = true - // Padfoot/ReConPlot need whole-genome caller output; covered by the module- and subworkflow-level nf-tests + // Padfoot/ReConPlot/ShatterSeek need whole-genome caller output; covered by the module- and subworkflow-level nf-tests skip_padfoot = true skip_reconplot = true + skip_shatterseek = true savana_chromosomes = "19" // SAVANA's het-SNP coverage/mapq floors (--allele_min_reads default 10, --allele_mapq // default 5) aren't met on this minimal chr19 slice; lowering them gets real allele diff --git a/conf/test_full.config b/conf/test_full.config index 7fab9491..ddcd5e1f 100644 --- a/conf/test_full.config +++ b/conf/test_full.config @@ -51,9 +51,10 @@ params { savana_chromosomes = "19" savana_allele_min_reads = 2 savana_allele_mapq = 0 - // Padfoot/ReConPlot need whole-genome caller output; covered by the module- and subworkflow-level nf-tests + // Padfoot/ReConPlot/ShatterSeek need whole-genome caller output; covered by the module- and subworkflow-level nf-tests skip_padfoot = true skip_reconplot = true + skip_shatterseek = true // Both callers, union combine mode germline_var_keep = 'clair, deepvariant' diff --git a/containers/shatterseek/Dockerfile b/containers/shatterseek/Dockerfile new file mode 100644 index 00000000..1a194caa --- /dev/null +++ b/containers/shatterseek/Dockerfile @@ -0,0 +1,24 @@ +# syntax=docker/dockerfile:1.6 +FROM docker.io/mambaorg/micromamba:2.3.2 + +# R runtime for the ShatterSeek wrapper, with the robert-a-forsyth/ShatterSeek fork (chm13 support, statistical +# criteria fixed to the PCAWG definitions, PCAWG confidence tiers) installed from a pinned commit. The wrapper ships +# with the pipeline in assets/shatterseek/. +COPY --chown=$MAMBA_USER containers/shatterseek/environment.yml /tmp/environment.yml + +ARG SHATTERSEEK_COMMIT=e1cee6c8ee75ea1a1e9f4d1c57b7a8a0a199ae91 +ADD --checksum=sha256:74b4258a0e6b509aeba4976348feebf910c098b9556e85008795d1b4e6286b69 --chown=$MAMBA_USER \ + https://github.com/robert-a-forsyth/ShatterSeek/archive/${SHATTERSEEK_COMMIT}.tar.gz /tmp/shatterseek.tar.gz + +RUN micromamba create --yes --name shatterseek --file /tmp/environment.yml \ + && micromamba clean --all --yes \ + && export PATH=/opt/conda/envs/shatterseek/bin:${PATH} \ + && mkdir -p /tmp/shatterseek-src \ + && tar -xzf /tmp/shatterseek.tar.gz -C /tmp/shatterseek-src --strip-components 1 \ + && R CMD INSTALL --no-docs --no-html /tmp/shatterseek-src \ + # the fork's regression test: DO17373 against the PCAWG-published values, plus synthetic hg38 / chm13 events + && Rscript /tmp/shatterseek-src/tests/test_do17373.R \ + && rm -rf /tmp/shatterseek.tar.gz /tmp/shatterseek-src + +ENV PATH=/opt/conda/envs/shatterseek/bin:${PATH} \ + SHATTERSEEK_COMMIT=${SHATTERSEEK_COMMIT} diff --git a/containers/shatterseek/README.md b/containers/shatterseek/README.md new file mode 100644 index 00000000..12adceae --- /dev/null +++ b/containers/shatterseek/README.md @@ -0,0 +1,29 @@ +# ShatterSeek image + +R runtime for the ShatterSeek wrapper used by the `SHATTERSEEK` module, with the +[robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support) fork installed from +a pinned, checksum-verified commit (`SHATTERSEEK_COMMIT` in the Dockerfile). Relative to upstream +[parklab/ShatterSeek](https://github.com/parklab/ShatterSeek) the fork adds T2T-CHM13 (`genome = "chm13"`), corrects +the statistical criteria to the PCAWG definitions, keeps one cluster per chromosome and labels clusters with the PCAWG +confidence tiers (`classify_chromothripsis()`); its README lists every change. The build runs the fork's regression +test (`tests/test_do17373.R`: the DO17373 example against the PCAWG-published values, plus synthetic hg38 and chm13 +events) and fails if it does not pass. The wrapper is part of the pipeline in `assets/shatterseek/`. The module does +not support `-profile conda`: the image ships the package. + +Licence: the fork's `DESCRIPTION` declares GPL (>= 3), as upstream does, while the upstream README states that +ShatterSeek is free for academic use only. + +Build and publish from the pipeline root (on the VSC login nodes: podman with the cgroupfs/file flags, TMPDIR on local +scratch). After a new fork commit, update `SHATTERSEEK_COMMIT` and the `ADD --checksum` (the sha256 of +`https://github.com/robert-a-forsyth/ShatterSeek/archive/.tar.gz`) first: + +```bash +TAG=1.1.1-chm13- +podman --cgroup-manager=cgroupfs --events-backend=file build --format docker \ + -f containers/shatterseek/Dockerfile -t "docker.io/robertaforsyth/shatterseek:$TAG" . +podman --cgroup-manager=cgroupfs --events-backend=file push "docker.io/robertaforsyth/shatterseek:$TAG" +``` + +There is no SIF twin: like CoRAL, Singularity/Apptainer pull the Docker image. The module pins the tag in +`modules/local/shatterseek/main.nf`; override per site via +`process { withName: '.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = ... } }`. diff --git a/containers/shatterseek/environment.yml b/containers/shatterseek/environment.yml new file mode 100644 index 00000000..fc9fbf52 --- /dev/null +++ b/containers/shatterseek/environment.yml @@ -0,0 +1,20 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +# R runtime for the ShatterSeek wrapper (assets/shatterseek), used only to build containers/shatterseek/Dockerfile; +# ShatterSeek itself (the robert-a-forsyth fork, not on conda) is installed into the image from a pinned commit. +channels: + - conda-forge + - bioconda +dependencies: + - "conda-forge::r-base=4.4.3" + - "bioconda::bioconductor-biocgenerics=0.52.0" + - "bioconda::bioconductor-s4vectors=0.44.0" + - "bioconda::bioconductor-iranges=2.40.0" + - "bioconda::bioconductor-genomicranges=1.58.0" + - "bioconda::bioconductor-graph=1.84.0" + - "conda-forge::r-foreach=1.5.2" + - "conda-forge::r-ggplot2=4.0.3" + - "conda-forge::r-gridextra=2.3" + - "conda-forge::r-data.table=1.18.4" + - "conda-forge::r-optparse=1.8.2" + - "conda-forge::procps-ng" # `ps` required by Nextflow task metrics diff --git a/docs/output.md b/docs/output.md index 9487ba34..65c35dd9 100644 --- a/docs/output.md +++ b/docs/output.md @@ -36,6 +36,10 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ ├── severus_ascat │ │ ├── severus_wakhan │ │ └── savana +│ ├── shatterseek +│ │ ├── severus_ascat +│ │ ├── severus_wakhan +│ │ └── savana │ ├── signatures │ │ ├── assignment │ │ └── matrices @@ -89,6 +93,10 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ ├── severus_ascat │ │ ├── severus_wakhan │ │ └── savana +│ ├── shatterseek +│ │ ├── severus_ascat +│ │ ├── severus_wakhan +│ │ └── savana │ ├── signatures │ │ ├── assignment │ │ └── matrices @@ -110,7 +118,7 @@ The pipeline produces per-sample output directories. Two modes exist depending o └── multiqc ``` -The `padfoot` and `reconplot` directories are only present when the corresponding step is enabled (`--skip_padfoot`, `--skip_reconplot`); SAVANA's own output lives under `variants/savana`. Within them, each caller-pair subdirectory requires both of its callers to have produced output for that sample: `severus_wakhan` needs `--skip_wakhan false`, `severus_ascat` needs `--skip_ascat false` (ASCAT runs for paired and tumour-only samples, so this pair forms for both), and the `savana` subdirectories additionally need SAVANA copy number, which is only produced when an SNP source is available (the phased germline VCF for paired samples, or the bundled 1000G panel for tumour-only samples) and SAVANA finds an acceptable purity/ploidy fit. +The `padfoot`, `reconplot` and `shatterseek` directories are only present when the corresponding step is enabled (`--skip_padfoot`, `--skip_reconplot`, `--skip_shatterseek`); SAVANA's own output lives under `variants/savana`. Within them, each caller-pair subdirectory requires both of its callers to have produced output for that sample: `severus_wakhan` needs `--skip_wakhan false`, `severus_ascat` needs `--skip_ascat false` (ASCAT runs for paired and tumour-only samples, so this pair forms for both), and the `savana` subdirectories additionally need SAVANA copy number, which is only produced when an SNP source is available (the phased germline VCF for paired samples, or the bundled 1000G panel for tumour-only samples) and SAVANA finds an acceptable purity/ploidy fit. ### `ascat` @@ -769,6 +777,70 @@ A figure ReConPlot cannot draw is skipped and the task keeps the others: ReConPl +### `shatterseek` + +
+Output files + +``` +├── shatterseek +│ ├── severus_ascat +│ │ ├── sample.shatterseek.tsv +│ │ ├── sample.shatterseek_sv_input.tsv +│ │ ├── sample.shatterseek_cn_input.tsv +│ │ ├── plots/sample_chr.pdf +│ │ └── shatterseek.log +│ ├── severus_wakhan +│ │ └── (same layout) +│ └── savana +│ └── (same layout) +``` + +| File | Description | +| -------------------------------------- | --------------------------------------------------------------------------------------------------------------------------------- | +| `sample.shatterseek.tsv` | One row per chromosome (chr1-22, chrX): the largest cluster of interleaved SVs, the statistical criteria and the confidence tiers | +| `sample.shatterseek_{sv,cn}_input.tsv` | The SV and copy-number tables exactly as given to ShatterSeek (Ensembl chromosome names, integer copy number) | +| `plots/*.pdf` | ShatterSeek figure (ideogram, SVs, copy number, summary) for every chromosome with a confidence tier (optional) | +| `shatterseek.log` | Wrapper log: files read, SVs dropped and why, copy numbers rounded, the tiers called | + +[ShatterSeek](https://github.com/parklab/ShatterSeek) (Cortes-Ciriano et al. 2020) finds, per chromosome, the largest cluster of interleaved intra-chromosomal SVs and tests it against the hallmarks of chromothripsis. The pipeline runs a fork ([robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support)) whose statistics follow the definitions in the PCAWG Supplementary Note; they differ from upstream ShatterSeek: + +| Column | Meaning | +| ------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------ | +| `start`, `end`, `clusterSize` | Span and number of intra-chromosomal SVs of the chromosome's largest cluster (one cluster; NA/0 without one) | +| `number_DEL`/`_DUP`/`_h2hINV`/`_t2tINV` | The cluster's SVs by type | +| `number_TRA`, `clusterSize_including_TRA` | Translocations with a breakend inside the cluster, and the cluster size including them | +| `pval_fragment_joins` | Chi-square test of the cluster's four SV types against equal proportions. Chromothripsis expects random joins, so this should **not** be significant | +| `chr_breakpoint_enrichment` | One-sided binomial test: more breakpoints on this chromosome than its share of the genome predicts (hg38: non-N bases; chm13: length) | +| `pval_exp_chr`, `pval_exp_cluster` | Kolmogorov-Smirnov test of the gaps between consecutive breakpoints (chromosome / cluster) against an exponential; significant = clustered breakpoints | +| `max_number_oscillating_CN_segments_2_states` | Longest run of adjacent copy-number segments in the cluster region alternating between two states | +| `q_*` | Benjamini-Hochberg q-values of the four tests across the sample's chromosomes | +| `inter_*` | Other chromosomes linked to this one by translocations, and the fragment-joins test over the linked clusters | +| `HC1`, `HC2`, `HC3`, `LC1` | PCAWG confidence tiers met (see below) | +| `chromothripsis_criteria`, `chromothripsis_confidence` | The tiers met, e.g. `HC1,HC2`, and `high` (any HC), `low` (LC1 only) or `none` | + +The tiers use the PCAWG cut-offs (Cortes-Ciriano et al. 2020, Supplementary Note 1.2; ShatterSeek issue #22), with significance at q ≤ 0.2 and NA counting as a fail: + +| Tier | Interleaved intra SVs | Other | Oscillating CN segments (2 states) | Fragment joins | Enrichment or exponential test | +| ---- | --------------------- | ----------------------------- | ---------------------------------- | -------------- | ------------------------------ | +| HC1 | ≥ 6 | | ≥ 7 | q > 0.2 | q ≤ 0.2 | +| HC2 | ≥ 3 | ≥ 4 translocations | ≥ 7 | q > 0.2 | | +| HC3 | | ≥ 40 SVs incl. translocations | | q > 0.2 | | +| LC1 | ≥ 6 | | 4-6 | q > 0.2 | q ≤ 0.2 | + +Points to keep in mind when reading the calls: + +- PCAWG applied HC3 only to osteosarcomas, leiomyosarcomas and liposarcomas, where amplification erodes the oscillating pattern; here it applies to every tumour. Drop it with `ext.args = '--no-hc3'`. +- PCAWG corrected the P values once across every chromosome of all 2,658 tumours; within one sample the correction covers 23 chromosomes. The raw P values are kept, so a cohort-wide correction can be recomputed from the summaries (each row carries `sample`, `cn_source` and `sv_source`). +- The cut-offs were set by visual inspection of PCAWG's short-read consensus calls. Long-read callers report many small SVs, whose two breakpoints a few hundred bp apart make the exponential test significant on most chromosomes; intra-chromosomal SVs shorter than `--shatterseek_min_svlen` (default 150 bp) are therefore left out, and raising it tests larger rearrangements only. +- ShatterSeek does not mask centromeres or satellite arrays. A cluster spanning a centromere or heterochromatin block (e.g. chr16 q11.2, chr1 q12, chr9 q12) is more likely mapping artefact than chromothripsis; check the figure. +- Tumour-only Severus calls include germline SVs that the panel of normals did not remove. +- Insertions and single breakends are not used, nor are SVs touching chrY, the mitochondrion or unplaced contigs. SAVANA's absolute copy number and fractional Wakhan states are rounded to integers, and neighbouring segments with equal copy number are merged, because ShatterSeek counts copy-number states row by row. + +Tier cut-offs, the q/P choice and SV filters can be changed with `ext.args` on `.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)` (see `Rscript assets/shatterseek/run_shatterseek.R --help`). A sample ShatterSeek cannot process is skipped without stopping the run. + +
+ ### `wakhan`
diff --git a/docs/usage.md b/docs/usage.md index 3a64a4ac..8146be08 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -163,7 +163,7 @@ Keep one resource set per directory: ClairS-TO derives the per-contig prefix fro If the loci cannot belong to the reference, ClairS-TO disables Verdict with a warning rather than applying the wrong coordinates. On a `--skip_ascat` run, look for `VERDICT CNA RESOURCE DIRECTORY` in the ClairS-TO log to confirm which set was used; otherwise Verdict runs outside ClairS-TO and the tables the tags came from are published next to the VCFs. -Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. With `--genome CHM13` it is skipped with a warning unless you also pass `--padfoot_gff` and `--padfoot_rm` (see [Padfoot Options](#padfoot-options)). ReConPlot needs nothing extra on CHM13 (`T2T` preset). +Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. With `--genome CHM13` it is skipped with a warning unless you also pass `--padfoot_gff` and `--padfoot_rm` (see [Padfoot Options](#padfoot-options)). ReConPlot (`T2T` preset) and ShatterSeek (`chm13`) need nothing extra on CHM13. ### Pipeline options @@ -185,6 +185,7 @@ Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. With `-- | `--skip_savana` | A boolean to skip `savana` (SV + copy-number calling). Default = `false` | | `--skip_padfoot` | A boolean to skip `padfoot` SV/CNA annotation. Default = `false` | | `--skip_reconplot` | A boolean to skip `reconplot` SV/CNA figures. Default = `false` | +| `--skip_shatterseek` | A boolean to skip `shatterseek` chromothripsis calling. Default = `false` | | `--skip_bamstats` | A boolean to skip `bamstats`. Default = `false` | | `--skip_wakhan` | A boolean to skip `wakhan`. Default = `false` | | `--skip_vep` | A boolean to skip `vep`. Default = `false` | @@ -457,6 +458,23 @@ Each pair produces `per_chromosome/` (one figure per chromosome), `genome_wide/` | `--reconplot_baf_track` | Add a het-SNP BAF track to the focus panel (ASCAT and SAVANA only). Default = `false` | | `--reconplot_format` | Output formats: `pdf`, `png` or `pdf,png`. Default = `pdf,png` | +#### ShatterSeek Options + +[ShatterSeek](https://github.com/parklab/ShatterSeek) chromothripsis calls are made for every CN/SV caller pair available for a sample (`severus_ascat`, `severus_wakhan`, `savana`, the same pairs as ReConPlot), into `shatterseek//`, and each chromosome's cluster is labelled with the PCAWG high/low-confidence tiers (see [output](output.md#shatterseek)). The pipeline runs a fork, [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), with T2T-CHM13 support and the statistical criteria corrected to the PCAWG definitions; it ships inside the module's image (`docker.io/robertaforsyth/shatterseek`, recipe in `containers/shatterseek/`), so the module is not available under `-profile conda`. The callers are read by the same parsers as ReConPlot (`assets/reconplot/R/parsers`), through the wrapper in `assets/shatterseek/`. + +| Parameter | Description | +| ------------------------- | --------------------------------------------------------------------------------------------------------------------------- | +| `--shatterseek_genome` | ShatterSeek genome (`hg38`, `chm13`, `hg19`). Default = `null` (inferred from `--genome`; skipped with a warning otherwise) | +| `--shatterseek_min_svlen` | Drop intra-chromosomal SVs shorter than this (bp) before ShatterSeek; translocations kept. Default = `150` | + +Tier cut-offs (`--min-sv-hc1`, `--min-osc-hc`, `--min-sv-hc2`, `--min-tra-hc2`, `--min-sv-hc3`, `--min-osc-lc`, `--max-osc-lc`, `--no-hc3`), the significance level (`--alpha`, default 0.2), q-values or raw P values (`--use q|p`) and the other SV filters (`--min-support`, `--exclude-vntr`) are wrapper options, set through `ext.args` (keep `--min-svlen` in it, or use `--shatterseek_min_svlen`, since `ext.args` replaces the pipeline's value): + +```groovy +process { + withName: '.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { ext.args = '--min-svlen 150 --no-hc3' } +} +``` + ##### Offline and air-gapped systems Both modules use the same pinned-tag pair as the pipeline's other large custom images (ClairS-TO): a prebuilt SIF (`oras://docker.io/timmy9527/-sif:`) under Singularity/Apptainer and the Docker image (`docker.io/timmy9527/:`) otherwise, both on Docker Hub, so `nf-core pipelines download --container-system singularity` stages them like every other container in the pipeline. To use a different image (e.g. a local mirror) override it in a config file: @@ -465,9 +483,12 @@ Both modules use the same pinned-tag pair as the pipeline's other large custom i process { withName: '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/padfoot-repeatmasker.sif' } withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/reconplot.sif' } + withName: '.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/shatterseek.sif' } } ``` +ShatterSeek has no SIF twin: like CoRAL, its Docker image is used under every engine (Singularity/Apptainer convert it when pulling). + Neither step downloads anything at run time: Padfoot and the ReConPlot R package are inside the two images, and the ReConPlot wrapper is part of the pipeline. Updating either tool means rebuilding its image, pushing the Docker image and its SIF twin, and bumping the tags in the module (see the READMEs under `containers/`). #### Variant Filtering and Combining Options diff --git a/modules/local/shatterseek/main.nf b/modules/local/shatterseek/main.nf new file mode 100644 index 00000000..e2675a78 --- /dev/null +++ b/modules/local/shatterseek/main.nf @@ -0,0 +1,71 @@ +process SHATTERSEEK { + tag "${meta.id}:${cn_source}+${sv_source}" + label 'process_low' + + // No conda: the ShatterSeek fork (chm13, corrected statistics, PCAWG tiers) ships only in this image, built from + // containers/shatterseek/Dockerfile; the wrapper is pipeline glue in assets/shatterseek/ and reads the callers with + // the ReConPlot parsers in assets/reconplot/. Fork update = new commit there, rebuild, bump this tag. + container "docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c" + + input: + // cn_files / sv_files are the caller output files the parsers discover by name. + // When cn_source == sv_source (savana) put everything in cn_files and leave sv_files empty. + tuple val(meta), val(cn_source), path(cn_files, stageAs: 'cn_input/*'), val(sv_source), path(sv_files, stageAs: 'sv_input/*') + tuple val(meta2), path(shatterseek_src) // ShatterSeek wrapper (run_shatterseek.R + VERSION) + tuple val(meta3), path(reconplot_src) // ReConPlot wrapper, for its caller parsers (R/) + val(genome) // hg38 | chm13 | hg19 + + output: + tuple val(meta), path("${prefix}/${meta.id}.shatterseek.tsv") , emit: summary + tuple val(meta), path("${prefix}/${meta.id}.shatterseek_{sv,cn}_input.tsv") , emit: inputs + tuple val(meta), path("${prefix}/plots/*.pdf") , emit: plots, optional: true + tuple val(meta), path("${prefix}/shatterseek.log") , emit: log + tuple val(meta), path("${prefix}/staged.txt") , emit: staged, optional: true // stub only: the staged inputs, for the tests + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + // Exit if running this module with -profile conda / -profile mamba + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "SHATTERSEEK does not support Conda: the ShatterSeek fork ships only inside its container. Use Docker / Singularity / Apptainer, or --skip_shatterseek." + } + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}") + def sv_in = cn_source == sv_source ? '' : '--sv-input sv_input' + """ + mkdir -p ${prefix} + + Rscript ${shatterseek_src}/run_shatterseek.R \\ + --reconplot-dir ${reconplot_src} \\ + --cn-source ${cn_source} --sv-source ${sv_source} --cn-input cn_input ${sv_in} \\ + --sample ${meta.id} --prefix ${meta.id} --genome ${genome} --outdir ${prefix} \\ + ${args} 2>&1 | tee ${prefix}/shatterseek.log + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + shatterseek: \$(Rscript -e 'cat(as.character(packageVersion("ShatterSeek")))' 2>/dev/null) + shatterseek_fork_commit: \${SHATTERSEEK_COMMIT:-unknown} + shatterseek_wrapper: \$(cat ${shatterseek_src}/VERSION 2>/dev/null || echo unknown) + r-base: \$(Rscript -e 'cat(R.version\$major, R.version\$minor, sep=".")' 2>/dev/null) + END_VERSIONS + """ + + stub: + prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}") + """ + mkdir -p ${prefix}/plots + touch ${prefix}/${meta.id}.shatterseek.tsv ${prefix}/${meta.id}.shatterseek_sv_input.tsv ${prefix}/${meta.id}.shatterseek_cn_input.tsv + touch ${prefix}/plots/${meta.id}_chr1.pdf ${prefix}/shatterseek.log + { ls -1 cn_input/* sv_input/* 2>/dev/null || true; } | sort > ${prefix}/staged.txt # sv_input is absent in single-source mode + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + shatterseek: stub + shatterseek_fork_commit: stub + shatterseek_wrapper: stub + r-base: stub + END_VERSIONS + """ +} diff --git a/modules/local/shatterseek/meta.yml b/modules/local/shatterseek/meta.yml new file mode 100644 index 00000000..b08929ad --- /dev/null +++ b/modules/local/shatterseek/meta.yml @@ -0,0 +1,103 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "shatterseek" +description: Chromothripsis calls with ShatterSeek (robert-a-forsyth fork with chm13 support and corrected statistics), labelled with the PCAWG high/low-confidence tiers, from ASCAT/Wakhan/SAVANA copy number and Severus/SAVANA SVs. +keywords: + - chromothripsis + - structural-variation + - copy-number + - somatic + - long-read +tools: + - "shatterseek": + description: "ShatterSeek: detection of chromothripsis from SV and copy-number calls" + homepage: "https://github.com/parklab/ShatterSeek" + documentation: "https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support" + tool_dev_url: "https://github.com/robert-a-forsyth/ShatterSeek" + doi: "10.1038/s41588-019-0576-7" + licence: ["GPL-3.0-or-later"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - cn_source: + type: string + description: CN parser name (`ascat`, `wakhan`, `savana`) + - cn_files: + type: file + description: Caller output files the CN parser discovers by name (staged into `cn_input/`). ASCAT `segments.txt`; Wakhan's top-ranked `integer_profile.bed`; SAVANA's segmented CN table plus its classified somatic VCF. + - sv_source: + type: string + description: SV parser name (`severus`, `savana`) + - sv_files: + type: file + description: SV caller output files (staged into `sv_input/`); empty list when `cn_source == sv_source` + - - meta2: + type: map + description: Groovy Map for the wrapper source + - shatterseek_src: + type: directory + description: ShatterSeek wrapper directory (`assets/shatterseek`, contains `run_shatterseek.R` and `VERSION`) + - - meta3: + type: map + description: Groovy Map for the parser source + - reconplot_src: + type: directory + description: ReConPlot wrapper directory (`assets/reconplot`), whose `R/` parsers read the caller output + - - genome: + type: string + description: ShatterSeek genome (`hg38`, `chm13`, `hg19`) + +output: + - summary: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/${meta.id}.shatterseek.tsv": + type: file + description: One row per chromosome (1-22, X) with the cluster, the statistical criteria, BH q-values, the PCAWG tiers (HC1, HC2, HC3, LC1) and chromothripsis_confidence + pattern: "*.shatterseek.tsv" + - inputs: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/${meta.id}.shatterseek_{sv,cn}_input.tsv": + type: file + description: The exact SV and copy-number tables given to ShatterSeek + pattern: "*.shatterseek_{sv,cn}_input.tsv" + - plots: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/plots/*.pdf": + type: file + description: ShatterSeek figure for each chromosome with a confidence tier (optional) + pattern: "*.pdf" + - log: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/shatterseek.log": + type: file + description: Wrapper log (parsing, dropped SVs, tiers) + pattern: "shatterseek.log" + - staged: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/staged.txt": + type: file + description: Stub runs only -- the CN and SV file names the task received, asserted on by the subworkflow tests (optional) + pattern: "staged.txt" + - 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/shatterseek/tests/main.nf.test b/modules/local/shatterseek/tests/main.nf.test new file mode 100644 index 00000000..34e33745 --- /dev/null +++ b/modules/local/shatterseek/tests/main.nf.test @@ -0,0 +1,208 @@ +nextflow_process { + + name "Test Process SHATTERSEEK" + script "../main.nf" + process "SHATTERSEEK" + config "./nextflow.config" + + tag "modules" + tag "modules_local" + tag "shatterseek" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // The stub tests check the I/O contract; the "real" tests run the wrapper and the ShatterSeek fork on the synthetic + // event in tests/fixtures/shatterseek (chr19 chromothripsis linked to chr7, an isolated chr2 deletion), written in + // each caller's own layout. Every pair must call chr19 high confidence (HC1, HC2) and nothing else. + + test("ascat + severus - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'ascat', + [ file("${projectDir}/tests/fixtures/shatterseek/ascat/test.segments.txt", checkIfExists: true) ], + 'severus', + [ file("${projectDir}/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz", checkIfExists: true) ] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("savana single source - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'savana', + [ file("${projectDir}/tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv", checkIfExists: true), + file("${projectDir}/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf", checkIfExists: true) ], + 'savana', + [] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("ascat + severus - real") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'ascat', + [ file("${projectDir}/tests/fixtures/shatterseek/ascat/test.segments.txt", checkIfExists: true) ], + 'severus', + [ file("${projectDir}/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz", checkIfExists: true) ] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + def rows = path(process.out.summary.get(0).get(1)).readLines() + def head = rows[0].split('\t') as List + def called = rows.drop(1).collect { r -> r.split('\t', -1) } + .findAll { r -> r[head.indexOf('chromothripsis_confidence')] != 'none' } + .collect { r -> "${r[head.indexOf('chrom')]}:${r[head.indexOf('chromothripsis_confidence')]}:${r[head.indexOf('chromothripsis_criteria')]}" } + def log_text = path(process.out.log.get(0).get(1)).text + assertAll( + { assert process.success }, + { assert called == ['chr19:high:HC1,HC2'] }, + { assert rows.size() == 24 }, // header + chromosomes 1-22 and X + { assert log_text.contains('dropped 1 non-PASS VCF records') && log_text.contains('INS=1') }, + { assert process.out.plots.get(0).get(1).toString().endsWith('test_chr19.pdf') }, + { assert snapshot(process.out.summary, process.out.inputs, process.out.versions).match() } + ) + } + } + + test("wakhan + severus - real") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'wakhan', + [ file("${projectDir}/tests/fixtures/shatterseek/wakhan/solution_rank_1/integer_profile.bed", checkIfExists: true) ], + 'severus', + [ file("${projectDir}/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz", checkIfExists: true) ] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + def log_text = path(process.out.log.get(0).get(1)).text + assertAll( + { assert process.success }, + { assert log_text.contains('tiers: chr19 high (HC1,HC2)') }, + // the subclonal chr7 haplotype state (1 + 1.2) is rounded to an integer total + { assert log_text.contains('rounded 1 non-integer total copy numbers') }, + { assert snapshot(process.out.summary, process.out.inputs, process.out.versions).match() } + ) + } + } + + test("savana single source - real") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'savana', + [ file("${projectDir}/tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv", checkIfExists: true), + file("${projectDir}/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf", checkIfExists: true) ], + 'savana', + [] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + def sv_rows = path(process.out.inputs.get(0).get(1).find { f -> f.toString().endsWith('_sv_input.tsv') }).readLines() + assertAll( + { assert process.success }, + { assert path(process.out.log.get(0).get(1)).text.contains('tiers: chr19 high (HC1,HC2)') }, + // SAVANA writes BP_NOTATION in (_1, _2) order; the chr7 _2 records come first in the file, and the + // translocations must keep the chr19 side's orientation (the same table as from Severus) + { assert sv_rows.findAll { r -> r.contains('TRA') } == [ + '19\t20500000\t7\t50100000\tTRA\t+\t-', '19\t22500000\t7\t50300000\tTRA\t-\t+', + '19\t24500000\t7\t50500000\tTRA\t+\t+', '19\t26500000\t7\t50700000\tTRA\t-\t-', + '19\t28500000\t7\t50900000\tTRA\t+\t-' ] }, + { assert snapshot(process.out.summary, process.out.inputs, process.out.versions).match() } + ) + } + } + + test("ascat + severus - real, chm13") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'ascat', + [ file("${projectDir}/tests/fixtures/shatterseek/ascat/test.segments.txt", checkIfExists: true) ], + 'severus', + [ file("${projectDir}/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz", checkIfExists: true) ] + ] + input[1] = [ [ id:'shatterseek' ], file("${projectDir}/assets/shatterseek", checkIfExists: true) ] + input[2] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[3] = 'chm13' + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert path(process.out.log.get(0).get(1)).text.contains('tiers: chr19 high (HC1,HC2)') }, + { assert process.out.plots.get(0).get(1).toString().endsWith('test_chr19.pdf') } + ) + } + } +} diff --git a/modules/local/shatterseek/tests/main.nf.test.snap b/modules/local/shatterseek/tests/main.nf.test.snap new file mode 100644 index 00000000..a284d8e3 --- /dev/null +++ b/modules/local/shatterseek/tests/main.nf.test.snap @@ -0,0 +1,301 @@ +{ + "ascat + severus - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.shatterseek_sv_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "2": [ + [ + { + "id": "test" + }, + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + [ + { + "id": "test" + }, + "shatterseek.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "4": [ + [ + { + "id": "test" + }, + "staged.txt:md5,3f1ba489ac329fd2881e318803d3068a" + ] + ], + "5": [ + "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + ], + "inputs": [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.shatterseek_sv_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "log": [ + [ + { + "id": "test" + }, + "shatterseek.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "plots": [ + [ + { + "id": "test" + }, + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "staged": [ + [ + { + "id": "test" + }, + "staged.txt:md5,3f1ba489ac329fd2881e318803d3068a" + ] + ], + "summary": [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions": [ + "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-14T15:40:32.880945189" + }, + "ascat + severus - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,52bdbcbc33b8fceb7cce44a19ec47ae7" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,982eb23095d2181864328d8ad5ff771f", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-14T15:40:48.576974826" + }, + "wakhan + severus - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,4298d3136a0d42899397d807b4d7a516" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,e5fc9a40fddb5b1e0371ec6ed19356d7", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-14T15:40:58.358617435" + }, + "savana single source - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,9dfb723da4741ead534812fa6e212799" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,982eb23095d2181864328d8ad5ff771f", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-14T15:41:08.198924416" + }, + "savana single source - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.shatterseek_sv_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "2": [ + [ + { + "id": "test" + }, + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + [ + { + "id": "test" + }, + "shatterseek.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "4": [ + [ + { + "id": "test" + }, + "staged.txt:md5,4464eb9289c1b04cbe9f13b2fa94befc" + ] + ], + "5": [ + "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + ], + "inputs": [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.shatterseek_sv_input.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "log": [ + [ + { + "id": "test" + }, + "shatterseek.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "plots": [ + [ + { + "id": "test" + }, + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "staged": [ + [ + { + "id": "test" + }, + "staged.txt:md5,4464eb9289c1b04cbe9f13b2fa94befc" + ] + ], + "summary": [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions": [ + "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-14T15:40:38.682355173" + } +} \ No newline at end of file diff --git a/modules/local/shatterseek/tests/nextflow.config b/modules/local/shatterseek/tests/nextflow.config new file mode 100644 index 00000000..8c5ffc88 --- /dev/null +++ b/modules/local/shatterseek/tests/nextflow.config @@ -0,0 +1,5 @@ +process { + withName: 'SHATTERSEEK' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + } +} diff --git a/nextflow.config b/nextflow.config index 3a2a0f3d..a325cc62 100644 --- a/nextflow.config +++ b/nextflow.config @@ -103,6 +103,7 @@ params { skip_savana = false skip_padfoot = false skip_reconplot = false + skip_shatterseek = false skip_wakhan = false skip_fiber = false skip_normalfiber = false @@ -153,6 +154,10 @@ params { reconplot_baf_track = false // het-SNP BAF track on the focus panel (ASCAT/SAVANA) reconplot_format = 'pdf,png' + // ShatterSeek options + shatterseek_genome = null // hg38 | chm13 | hg19; null = inferred from --genome + shatterseek_min_svlen = 150 // drop intra-chromosomal SVs shorter than this (bp); 0 = keep all + // ASCAT options ascat_ploidy = null ascat_min_base_qual = 20 diff --git a/nextflow_schema.json b/nextflow_schema.json index a7234ecf..8b3e2330 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -501,6 +501,27 @@ } } }, + "shatterseek_options": { + "title": "ShatterSeek options", + "type": "object", + "description": "Options for ShatterSeek chromothripsis calling", + "default": "", + "properties": { + "shatterseek_genome": { + "type": "string", + "enum": ["hg38", "chm13", "hg19"], + "help_text": "Genome for ShatterSeek's chromosome lengths (breakpoint enrichment test) and ideograms: `hg38`, `chm13` or `hg19`. Inferred from `--genome` (GRCh38 -> hg38, CHM13 -> chm13); for any other genome ShatterSeek is skipped with a warning unless this is set.", + "description": "ShatterSeek genome. Default: inferred from --genome (GRCh38 -> hg38, CHM13 -> chm13)." + }, + "shatterseek_min_svlen": { + "type": "integer", + "default": 150, + "minimum": 0, + "help_text": "Intra-chromosomal SVs shorter than this many bp are not given to ShatterSeek; translocations are always kept. Small SVs (e.g. VNTR indels) add pairs of breakpoints a few hundred bp apart, which make the exponential-distribution test read clustering where there is none. 0 keeps every SV.", + "description": "Drop intra-chromosomal SVs shorter than this many bp before ShatterSeek (translocations are kept). 0 keeps everything." + } + } + }, "ascat_parameters": { "title": "ASCAT parameters", "type": "object", @@ -683,6 +704,11 @@ "help_text": "Skips the ReConPlot rearrangement and copy-number figures (ASCAT + Severus, Wakhan + Severus and SAVANA pairs). ReConPlot runs only inside its container; skip it under `-profile conda`.", "description": "Skips ReConPlot SV/CNA figures" }, + "skip_shatterseek": { + "type": "boolean", + "help_text": "Skips ShatterSeek chromothripsis calling (ASCAT + Severus, Wakhan + Severus and SAVANA pairs). ShatterSeek runs only inside its container; skip it under `-profile conda`.", + "description": "Skips ShatterSeek chromothripsis calling" + }, "skip_m6a": { "type": "boolean", "description": "Skip m6a calling by Fibertools" @@ -997,6 +1023,9 @@ { "$ref": "#/$defs/reconplot_options" }, + { + "$ref": "#/$defs/shatterseek_options" + }, { "$ref": "#/$defs/ascat_parameters" }, diff --git a/subworkflows/local/shatterseek_chromothripsis.nf b/subworkflows/local/shatterseek_chromothripsis.nf new file mode 100644 index 00000000..6dcd5727 --- /dev/null +++ b/subworkflows/local/shatterseek_chromothripsis.nf @@ -0,0 +1,85 @@ +// IMPORT MODULES +include { SHATTERSEEK as SHATTERSEEK_SEVERUS_ASCAT } from '../../modules/local/shatterseek/main' +include { SHATTERSEEK as SHATTERSEEK_SEVERUS_WAKHAN } from '../../modules/local/shatterseek/main' +include { SHATTERSEEK as SHATTERSEEK_SAVANA } from '../../modules/local/shatterseek/main' + +// +// ShatterSeek chromothripsis calls with PCAWG confidence tiers for every CN/SV caller pair that produced output for a +// sample: ASCAT + Severus, Wakhan + Severus, and SAVANA on its own. Pass channel.empty() for a caller that did not run. +// The pairs are built as in RECONPLOT_FIGURES, and the callers are read by the same parsers (assets/reconplot/R/parsers). +// +workflow SHATTERSEEK_CHROMOTHRIPSIS { + + take: + severus_vcf // [meta, severus_somatic.vcf.gz] + ascat_segments // [meta, segments.txt] + wakhan_solution_dirs // [meta, [solution_rank_N/, ...]] + savana_cna // [meta, segmented_absolute_copy_number.tsv] + savana_vcf // [meta, classified.somatic.vcf] + genome // ShatterSeek genome: hg38, chm13 or hg19 + + main: + ch_versions = channel.empty() + + // A digest of every wrapper file goes into the meta map, so editing either wrapper invalidates -resume + def wrapper_src = { name -> + def dir = file("${projectDir}/assets/${name}", type: 'dir', checkIfExists: true) + def digest = files("${projectDir}/assets/${name}/**", type: 'file').sort { f -> f.toString() }.collect { f -> f.text }.join('\n').md5() + return channel.value([[id: name, digest: digest], dir]) + } + shatterseek_src = wrapper_src('shatterseek') + reconplot_src = wrapper_src('reconplot') + + severus_sv_files = severus_vcf.map { meta, vcf -> [meta, [vcf]] } + // severus_sv_files: [meta, [severus_somatic.vcf.gz]] + + // + // MODULE: SHATTERSEEK_SEVERUS_ASCAT (label: process_low) + // Input: [meta, 'ascat', [segments.txt], 'severus', [vcf]] + // + // ASCAT writes an empty segments.txt when it finds no solution: nothing to test for that pair + ascat_segments + .filter { _meta, seg -> seg.size() > 0 && seg.countLines() > 1 } + .map { meta, seg -> [meta, [seg]] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'ascat', cn, 'severus', sv] } + .set { ascat_input } + + SHATTERSEEK_SEVERUS_ASCAT( ascat_input, shatterseek_src, reconplot_src, genome ) + ch_versions = ch_versions.mix(SHATTERSEEK_SEVERUS_ASCAT.out.versions) + + // + // MODULE: SHATTERSEEK_SEVERUS_WAKHAN (label: process_low) + // Input: [meta, 'wakhan', [integer_profile.bed], 'severus', [vcf]] -- the top-ranked solution (solution_rank_1/) + // + wakhan_solution_dirs + .map { meta, dirs -> + def best = [dirs].flatten().find { dir -> dir.name == 'solution_rank_1' } + def bed = best ? best.resolve('integer_profile.bed') : null + return [meta, bed?.exists() ? bed : null] + } + .filter { _meta, bed -> bed != null } + .map { meta, bed -> [meta, [bed]] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'wakhan', cn, 'severus', sv] } + .set { wakhan_input } + + SHATTERSEEK_SEVERUS_WAKHAN( wakhan_input, shatterseek_src, reconplot_src, genome ) + ch_versions = ch_versions.mix(SHATTERSEEK_SEVERUS_WAKHAN.out.versions) + + // + // MODULE: SHATTERSEEK_SAVANA (label: process_low) + // Input: [meta, 'savana', [cna.tsv, somatic.vcf], 'savana', []] -- single-source mode + // + savana_cna + .join(savana_vcf) + .map { meta, cna, sv_vcf -> [meta, 'savana', [cna, sv_vcf], 'savana', []] } + .set { savana_input } + + SHATTERSEEK_SAVANA( savana_input, shatterseek_src, reconplot_src, genome ) + ch_versions = ch_versions.mix(SHATTERSEEK_SAVANA.out.versions) + + emit: + summary = SHATTERSEEK_SEVERUS_ASCAT.out.summary.mix(SHATTERSEEK_SEVERUS_WAKHAN.out.summary, SHATTERSEEK_SAVANA.out.summary) // [meta, tsv] + versions = ch_versions // [versions.yml] +} diff --git a/subworkflows/local/tests/shatterseek_chromothripsis.nf.test b/subworkflows/local/tests/shatterseek_chromothripsis.nf.test new file mode 100644 index 00000000..d9a00fa3 --- /dev/null +++ b/subworkflows/local/tests/shatterseek_chromothripsis.nf.test @@ -0,0 +1,82 @@ +nextflow_workflow { + + name "Test Workflow SHATTERSEEK_CHROMOTHRIPSIS" + script "../shatterseek_chromothripsis.nf" + workflow "SHATTERSEEK_CHROMOTHRIPSIS" + + tag "subworkflows" + tag "subworkflows_local" + tag "shatterseek_chromothripsis" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // ShatterSeek needs whole-genome caller output, so the pipeline tests skip it. Under -stub what is tested is the + // pairing: ASCAT + Severus, the top-ranked Wakhan solution + Severus, SAVANA alone, and that an ASCAT run without + // a solution or a Wakhan run without a rank-1 profile forms no pair. + + test("each CN/SV caller pair forms per sample") { + + options "-stub" + + when { + params { + outdir = "$outputDir" + } + workflow { + """ + // Empty stand-ins carrying the real callers' file names; the stub does not read them + def mk = { String rel, String text = '' -> def f = file("\${launchDir}/fixtures/\${rel}"); f.parent.mkdirs(); f.text = text; return f } + def rank = { String rel -> mk.call(rel + '/integer_profile.bed'); return file("\${launchDir}/fixtures/\${rel}") } + def seg = 'sample\\tchr\\tstartpos\\tendpos\\tnMajor\\tnMinor\\nt\\t1\\t1\\t1000\\t1\\t1\\n' + def s1 = [ id:'sample1', paired_data:true, platform:'ont' ] + def s2 = [ id:'sample2', paired_data:true, platform:'pb' ] + def s3 = [ id:'sample3', paired_data:false, platform:'ont' ] + + input[0] = channel.of( + [ s1, mk.call('severus/sample1/severus_somatic.vcf.gz') ], + [ s2, mk.call('severus/sample2/severus_somatic.vcf.gz') ], + [ s3, mk.call('severus/sample3/severus_somatic.vcf.gz') ] + ) + // sample2's ASCAT found no solution (empty segments.txt): its ASCAT pair must not form + input[1] = channel.of([ s1, mk.call('ascat/sample1.segments.txt', seg) ], [ s2, mk.call('ascat/sample2.segments.txt') ], [ s3, mk.call('ascat/sample3.segments.txt', seg) ]) + // Wakhan >= 0.5: sample1 has two ranked solutions (rank 1 is used); sample2 has only a rank-2 profile + input[2] = channel.of( + [ s1, [ rank.call('wakhan/sample1/solution_rank_1'), rank.call('wakhan/sample1/solution_rank_2') ] ], + [ s2, [ rank.call('wakhan/sample2/solution_rank_2') ] ] + ) + input[3] = channel.of( + [ s1, mk.call('savana/sample1/sample1_segmented_absolute_copy_number.tsv') ], + [ s3, mk.call('savana/sample3/sample3_segmented_absolute_copy_number.tsv') ] + ) + input[4] = channel.of( + [ s1, mk.call('savana/sample1/sample1.classified.somatic.vcf') ], + [ s3, mk.call('savana/sample3/sample3.classified.somatic.vcf') ] + ) + input[5] = 'hg38' + """ + } + } + + then { + def tasks = workflow.trace.tasks().collect { task -> task.name } + // conf/modules.config publishes each task's / under //shatterseek/; + // the stub lists what was staged into cn_input/ and sv_input/ + def staged = { String rel -> path("$outputDir/${rel}/staged.txt").text.trim().split('\n') as List } + + assertAll( + { assert workflow.success }, + { assert tasks.count { name -> name.contains('SHATTERSEEK_SEVERUS_ASCAT') } == 2 }, + { assert tasks.any { name -> name.contains('SHATTERSEEK_SEVERUS_ASCAT (sample1:') } }, + { assert tasks.any { name -> name.contains('SHATTERSEEK_SEVERUS_ASCAT (sample3:') } }, + { assert staged('sample1/shatterseek/severus_ascat') == ['cn_input/sample1.segments.txt', 'sv_input/severus_somatic.vcf.gz'] }, + { assert tasks.count { name -> name.contains('SHATTERSEEK_SEVERUS_WAKHAN') } == 1 }, + { assert staged('sample1/shatterseek/severus_wakhan') == ['cn_input/integer_profile.bed', 'sv_input/severus_somatic.vcf.gz'] }, + { assert tasks.count { name -> name.contains('SHATTERSEEK_SAVANA') } == 2 }, + { assert staged('sample3/shatterseek/savana') == ['cn_input/sample3.classified.somatic.vcf', 'cn_input/sample3_segmented_absolute_copy_number.tsv'] }, + { assert !tasks.any { name -> name.contains('(sample2:') } }, + { assert workflow.out.summary.size() == 5 }, + { assert workflow.out.versions.size() >= 1 } + ) + } + } +} diff --git a/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf b/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf index 532f7082..8eacad3f 100644 --- a/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf +++ b/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf @@ -380,6 +380,13 @@ def padfootGenome() { return params.padfoot_genome ?: (params.genome == 'GRCh38' ? 'hg38' : params.genome == 'CHM13' ? 'chm13' : null) } +// +// ShatterSeek genome: an explicit --shatterseek_genome wins, else inferred from --genome (hg38 / chm13), else null (skipped) +// +def shatterseekGenome() { + return params.shatterseek_genome ?: (params.genome == 'GRCh38' ? 'hg38' : params.genome == 'CHM13' ? 'chm13' : null) +} + // // Whether Padfoot can annotate on that preset: it bundles annotations for hg38 and mm10 only; any other // preset needs both --padfoot_gff and --padfoot_rm, and a null preset disables Padfoot even with them. @@ -407,6 +414,14 @@ def validateSvAnnotationParams() { "--padfoot_gff and --padfoot_rm (and --padfoot_genome for a genome lrsomatic cannot infer)." } + // ShatterSeek's chromosome lengths and ideograms must be those of the alignments' assembly + if (!params.skip_shatterseek && genome_preset && params.shatterseek_genome && params.shatterseek_genome != genome_preset) { + error("--shatterseek_genome ${params.shatterseek_genome} does not match --genome ${params.genome} (expected ${genome_preset}).") + } + if (!params.skip_shatterseek && !shatterseekGenome()) { + log.warn "ShatterSeek will be skipped: genome could not be inferred from '${params.genome}'. Set --shatterseek_genome (hg38, chm13, hg19)." + } + if (!params.skip_reconplot && !params.reconplot_genome && !(params.genome in ['GRCh38', 'CHM13'])) { log.warn "ReConPlot: genome could not be inferred from '${params.genome}'; falling back to hg38 gene/chromosome annotations. " + "Set --reconplot_genome (hg38, hg19, T2T, mm10, mm39) to override." diff --git a/tests/fixtures/shatterseek/ascat/test.segments.txt b/tests/fixtures/shatterseek/ascat/test.segments.txt new file mode 100644 index 00000000..30a9b95e --- /dev/null +++ b/tests/fixtures/shatterseek/ascat/test.segments.txt @@ -0,0 +1,15 @@ +sample chr startpos endpos nMajor nMinor +test.tumour 2 1 242193529 1 1 +test.tumour 7 1 159345973 1 1 +test.tumour 19 1 20000000 1 1 +test.tumour 19 20000001 21000000 1 0 +test.tumour 19 21000001 22000000 1 1 +test.tumour 19 22000001 23000000 1 0 +test.tumour 19 23000001 24000000 1 1 +test.tumour 19 24000001 25000000 1 0 +test.tumour 19 25000001 26000000 1 1 +test.tumour 19 26000001 27000000 1 0 +test.tumour 19 27000001 28000000 1 1 +test.tumour 19 28000001 29000000 1 0 +test.tumour 19 29000001 30000000 1 1 +test.tumour 19 30000001 58617616 1 1 diff --git a/tests/fixtures/shatterseek/make_fixtures.py b/tests/fixtures/shatterseek/make_fixtures.py new file mode 100644 index 00000000..c9ff4654 --- /dev/null +++ b/tests/fixtures/shatterseek/make_fixtures.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 +"""Generate the synthetic ASCAT / Wakhan / Severus / SAVANA fixtures for the SHATTERSEEK real-data tests. + +One invented chromothripsis event in every caller's own file layout (GRCh38 coordinates, sample `test`): + chr19:20.1-29.4 Mb 12 interleaved SVs, three each of DEL (+-), DUP (-+), h2hINV (++) and t2tINV (--) + (SV k joins breakpoints k and k+12, so every SV overlaps every other without nesting) + chr19 <-> chr7 5 translocations from inside the cluster to chr7:50.1-50.9 Mb + chr19:20-30 Mb copy number alternating 1/2 over ten 1 Mb segments (2 elsewhere) + chr2:10.00-10.05 Mb one isolated deletion (a chromosome with an SV but no event) +Expected ShatterSeek result: chr19 high confidence (HC1, HC2), chr2 and chr7 none. +Also present so the conversion is exercised: an insertion and a non-PASS deletion (both dropped), a Wakhan segment +with a fractional haplotype state and fractional SAVANA copy numbers (rounded), and SAVANA translocations whose +_2 mate (on chr7) comes first in the file while BP_NOTATION is written in (_1, _2) order. +""" +import gzip, os + +F = os.path.dirname(os.path.abspath(__file__)) +LEN = {"chr2": 242193529, "chr7": 159345973, "chr19": 58617616} +ORDER = {"chr2": 2, "chr7": 7, "chr19": 19} + + +def w(name, lines, gz=False): + path = os.path.join(F, name) + os.makedirs(os.path.dirname(path), exist_ok=True) + with (gzip.open(path, "wt") if gz else open(path, "w")) as fh: + fh.write("\n".join(lines) + "\n") + print("wrote", name, len(lines), "lines") + + +def alt(strands, chrom, pos): + """VCF breakend ALT for a record whose own strand is strands[0] and its mate's strands[1].""" + own, mate = strands + bracket = "[" if mate == "-" else "]" + return f"N{bracket}{chrom}:{pos}{bracket}" if own == "+" else f"{bracket}{chrom}:{pos}{bracket}N" + + +# ---------------------------------------------------------------- the event +bps = [20_100_000 + k * 400_000 + (k % 3) * 7_000 for k in range(24)] +types = ["DEL", "DUP", "h2hINV", "t2tINV"] * 3 +STRANDS = {"DEL": "+-", "DUP": "-+", "h2hINV": "++", "t2tINV": "--"} +intra = [("chr19", bps[k], bps[k + 12], types[k]) for k in range(12)] +tra = [(p19, p7, s) for p19, p7, s in zip([20_500_000, 22_500_000, 24_500_000, 26_500_000, 28_500_000], + [50_100_000, 50_300_000, 50_500_000, 50_700_000, 50_900_000], + ["+-", "-+", "++", "--", "+-"])] + +# copy number (start, end, total): ten alternating 1 Mb segments inside chr19:20-30 Mb +cn = [("chr19", 1, 20_000_000, 2)] +cn += [("chr19", 20_000_001 + i * 1_000_000, 20_000_000 + (i + 1) * 1_000_000, 1 if i % 2 == 0 else 2) for i in range(10)] +cn += [("chr19", 30_000_001, LEN["chr19"], 2), ("chr2", 1, LEN["chr2"], 2), ("chr7", 1, LEN["chr7"], 2)] +cn.sort(key=lambda r: (ORDER[r[0]], r[1])) + +# ---------------------------------------------------------------- Severus (somatic SV VCF) +sev_hdr = ["##fileformat=VCFv4.2", "##source=Severus_v1.6", + '##FILTER=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##FORMAT=', + '##FORMAT='] + \ + [f"##contig=" for c, n in LEN.items()] + \ + ["#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ttest_tumor"] +rec = [] +def sev(chrom, pos, vid, alt_, info, flt="PASS"): + rec.append((chrom, pos, f"{chrom}\t{pos}\t{vid}\tN\t{alt_}\t60.0\t{flt}\tPRECISE;{info}\tGT:DV\t0/1:12")) +for i, (c, p1, p2, t) in enumerate(intra, 1): + s = STRANDS[t] + if t in ("DEL", "DUP"): # Severus writes deletions and duplications as symbolic alleles + sev(c, p1, f"severus_{t}{i}", f"<{t}>", f"SVTYPE={t};SVLEN={p2 - p1};END={p2};STRANDS={s}") + else: # and inversion-like junctions as mirrored breakend pairs + sev(c, p1, f"severus_BND{i}_1", alt(s, c, p2), f"SVTYPE=BND;SVLEN={p2 - p1};MATE_ID=severus_BND{i}_2;STRANDS={s}") + sev(c, p2, f"severus_BND{i}_2", alt(s[::-1], c, p1), f"SVTYPE=BND;SVLEN={p2 - p1};MATE_ID=severus_BND{i}_1;STRANDS={s[::-1]}") +for j, (p19, p7, s) in enumerate(tra, 1): + sev("chr19", p19, f"severus_BND{100 + j}_1", alt(s, "chr7", p7), f"SVTYPE=BND;SVLEN=0;MATE_ID=severus_BND{100 + j}_2;STRANDS={s}") + sev("chr7", p7, f"severus_BND{100 + j}_2", alt(s[::-1], "chr19", p19), f"SVTYPE=BND;SVLEN=0;MATE_ID=severus_BND{100 + j}_1;STRANDS={s[::-1]}") +sev("chr2", 10_000_000, "severus_DEL200", "", "SVTYPE=DEL;SVLEN=50000;END=10050000;STRANDS=+-") +sev("chr2", 20_000_000, "severus_DEL201", "", "SVTYPE=DEL;SVLEN=80000;END=20080000;STRANDS=+-", flt="FAIL_LOWSUPPORT") +sev("chr19", 25_000_100, "severus_INS1", "", "SVTYPE=INS;SVLEN=300") +rec.sort(key=lambda r: (ORDER[r[0]], r[1])) +w("severus/severus_somatic.vcf.gz", sev_hdr + [r[2] for r in rec], gz=True) + +# ---------------------------------------------------------------- ASCAT (sample test.tumour; chromosomes without chr) +w("ascat/test.segments.txt", ["sample\tchr\tstartpos\tendpos\tnMajor\tnMinor"] + + [f"test.tumour\t{c[3:]}\t{s}\t{e}\t1\t{t - 1}" for c, s, e, t in cn]) + +# ---------------------------------------------------------------- Wakhan 0.5.0 integer profile (top-ranked solution) +bed_hdr = ["#chr: chromosome number", "#start: start address for CN segment", "#end: end address for CN segment", + "#hp1_coverage / hp2_coverage: median coverage for this segment, per haplotype", + "#hp1_copynumber_state / hp2_copynumber_state: detected copy number state (integer/fraction), per haplotype", + "#hp1_confidence / hp2_confidence: confidence score, per haplotype", + "#svs_breakpoints_ids: corresponding structural variations (breakpoints) IDs from VCF file (union across haplotypes)", + "#chr\tstart\tend\thp1_coverage\thp1_copynumber_state\thp1_confidence\thp2_coverage\thp2_copynumber_state\thp2_confidence\tsvs_breakpoints_ids"] +def hp(c, t): + return (1.0, 1.2) if c == "chr7" else (1.0, float(t - 1)) # chr7: a subclonal 1.2 state (total 2.2 -> 2) +w("wakhan/solution_rank_1/integer_profile.bed", bed_hdr + + [f"{c}\t{s - 1}\t{e}\t{15 * hp(c, t)[0]:.2f}\t{hp(c, t)[0]}\t0.9\t{15 * hp(c, t)[1]:.2f}\t{hp(c, t)[1]}\t0.9\t" for c, s, e, t in cn]) + +# ---------------------------------------------------------------- SAVANA (prefix test_tumor; BND pairs, notation in (_1, _2) order) +sav_hdr = ["##fileformat=VCFv4.2", "##source=SAVANAv1.3.8", + '##FILTER=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##FORMAT='] + \ + [f"##contig=" for c, n in LEN.items()] + \ + ["#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ttest_tumor"] +rec = [] +def sav_pair(n, c1, p1, c2, p2, note, svlen): + info = f"SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN={svlen};BP_NOTATION={note}" + rec.append((c1, p1, f"{c1}\t{p1}\tID_{n}_1\tN\t{alt(note, c2, p2)}\t.\tPASS\tMATEID=ID_{n}_2;{info}\tGT\t0/1")) + rec.append((c2, p2, f"{c2}\t{p2}\tID_{n}_2\tN\t{alt(note[::-1], c1, p1)}\t.\tPASS\tMATEID=ID_{n}_1;{info}\tGT\t0/1")) +for i, (c, p1, p2, t) in enumerate(intra, 1): + sav_pair(i, c, p1, c, p2, STRANDS[t], p2 - p1) +for j, (p19, p7, s) in enumerate(tra, 1): # _1 on chr19, so the chr7 _2 record is written first + sav_pair(100 + j, "chr19", p19, "chr7", p7, s, 0) +sav_pair(200, "chr2", 10_000_000, "chr2", 10_050_000, "+-", 50000) +rec.append(("chr19", 25_000_100, "chr19\t25000100\tID_300_1\tN\t\t.\tPASS\tSVTYPE=INS;TUMOUR_READ_SUPPORT=12;SVLEN=300;BP_NOTATION=\tGT\t0/1")) +rec.sort(key=lambda r: (ORDER[r[0]], r[1])) +w("savana/test_tumor.classified.somatic.vcf", sav_hdr + [r[2] for r in rec]) +frac = {1: 1.08, 2: 1.94} # SAVANA's absolute copy number is not an integer +w("savana/test_tumor_segmented_absolute_copy_number.tsv", + ["chromosome\tstart\tend\tsegment_id\tbin_count\tsum_of_bin_lengths\tweight\tcopyNumber\tminorAlleleCopyNumber\tmeanBAF\tno_hetSNPs"] + + [f"{c}\t{s}\t{e}\t{c}_seg{i + 1}\t10\t{e - s + 1}\t10.0\t{frac[t]}\t{0.0 if t == 1 else 0.97}\t0.5\t100" for i, (c, s, e, t) in enumerate(cn)]) +print("done") diff --git a/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf b/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf new file mode 100644 index 00000000..5dbe3820 --- /dev/null +++ b/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf @@ -0,0 +1,50 @@ +##fileformat=VCFv4.2 +##source=SAVANAv1.3.8 +##FILTER= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##FORMAT= +##contig= +##contig= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT test_tumor +chr2 10000000 ID_200_1 N N[chr2:10050000[ . PASS MATEID=ID_200_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=50000;BP_NOTATION=+- GT 0/1 +chr2 10050000 ID_200_2 N ]chr2:10000000]N . PASS MATEID=ID_200_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=50000;BP_NOTATION=+- GT 0/1 +chr7 50100000 ID_101_2 N ]chr19:20500000]N . PASS MATEID=ID_101_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=+- GT 0/1 +chr7 50300000 ID_102_2 N N[chr19:22500000[ . PASS MATEID=ID_102_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=-+ GT 0/1 +chr7 50500000 ID_103_2 N N]chr19:24500000] . PASS MATEID=ID_103_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=++ GT 0/1 +chr7 50700000 ID_104_2 N [chr19:26500000[N . PASS MATEID=ID_104_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=-- GT 0/1 +chr7 50900000 ID_105_2 N ]chr19:28500000]N . PASS MATEID=ID_105_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=+- GT 0/1 +chr19 20100000 ID_1_1 N N[chr19:24900000[ . PASS MATEID=ID_1_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 20500000 ID_101_1 N N[chr7:50100000[ . PASS MATEID=ID_101_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=+- GT 0/1 +chr19 20507000 ID_2_1 N ]chr19:25307000]N . PASS MATEID=ID_2_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 20914000 ID_3_1 N N]chr19:25714000] . PASS MATEID=ID_3_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 21300000 ID_4_1 N [chr19:26100000[N . PASS MATEID=ID_4_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 +chr19 21707000 ID_5_1 N N[chr19:26507000[ . PASS MATEID=ID_5_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 22114000 ID_6_1 N ]chr19:26914000]N . PASS MATEID=ID_6_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 22500000 ID_7_1 N N]chr19:27300000] . PASS MATEID=ID_7_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 22500000 ID_102_1 N ]chr7:50300000]N . PASS MATEID=ID_102_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=-+ GT 0/1 +chr19 22907000 ID_8_1 N [chr19:27707000[N . PASS MATEID=ID_8_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 +chr19 23314000 ID_9_1 N N[chr19:28114000[ . PASS MATEID=ID_9_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 23700000 ID_10_1 N ]chr19:28500000]N . PASS MATEID=ID_10_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 24107000 ID_11_1 N N]chr19:28907000] . PASS MATEID=ID_11_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 24500000 ID_103_1 N N]chr7:50500000] . PASS MATEID=ID_103_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=++ GT 0/1 +chr19 24514000 ID_12_1 N [chr19:29314000[N . PASS MATEID=ID_12_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 +chr19 24900000 ID_1_2 N ]chr19:20100000]N . PASS MATEID=ID_1_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 25000100 ID_300_1 N . PASS SVTYPE=INS;TUMOUR_READ_SUPPORT=12;SVLEN=300;BP_NOTATION= GT 0/1 +chr19 25307000 ID_2_2 N N[chr19:20507000[ . PASS MATEID=ID_2_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 25714000 ID_3_2 N N]chr19:20914000] . PASS MATEID=ID_3_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 26100000 ID_4_2 N [chr19:21300000[N . PASS MATEID=ID_4_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 +chr19 26500000 ID_104_1 N [chr7:50700000[N . PASS MATEID=ID_104_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=-- GT 0/1 +chr19 26507000 ID_5_2 N ]chr19:21707000]N . PASS MATEID=ID_5_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 26914000 ID_6_2 N N[chr19:22114000[ . PASS MATEID=ID_6_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 27300000 ID_7_2 N N]chr19:22500000] . PASS MATEID=ID_7_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 27707000 ID_8_2 N [chr19:22907000[N . PASS MATEID=ID_8_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 +chr19 28114000 ID_9_2 N ]chr19:23314000]N . PASS MATEID=ID_9_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=+- GT 0/1 +chr19 28500000 ID_10_2 N N[chr19:23700000[ . PASS MATEID=ID_10_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-+ GT 0/1 +chr19 28500000 ID_105_1 N N[chr7:50900000[ . PASS MATEID=ID_105_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=+- GT 0/1 +chr19 28907000 ID_11_2 N N]chr19:24107000] . PASS MATEID=ID_11_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=++ GT 0/1 +chr19 29314000 ID_12_2 N [chr19:24514000[N . PASS MATEID=ID_12_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=4800000;BP_NOTATION=-- GT 0/1 diff --git a/tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv b/tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv new file mode 100644 index 00000000..7a2579db --- /dev/null +++ b/tests/fixtures/shatterseek/savana/test_tumor_segmented_absolute_copy_number.tsv @@ -0,0 +1,15 @@ +chromosome start end segment_id bin_count sum_of_bin_lengths weight copyNumber minorAlleleCopyNumber meanBAF no_hetSNPs +chr2 1 242193529 chr2_seg1 10 242193529 10.0 1.94 0.97 0.5 100 +chr7 1 159345973 chr7_seg2 10 159345973 10.0 1.94 0.97 0.5 100 +chr19 1 20000000 chr19_seg3 10 20000000 10.0 1.94 0.97 0.5 100 +chr19 20000001 21000000 chr19_seg4 10 1000000 10.0 1.08 0.0 0.5 100 +chr19 21000001 22000000 chr19_seg5 10 1000000 10.0 1.94 0.97 0.5 100 +chr19 22000001 23000000 chr19_seg6 10 1000000 10.0 1.08 0.0 0.5 100 +chr19 23000001 24000000 chr19_seg7 10 1000000 10.0 1.94 0.97 0.5 100 +chr19 24000001 25000000 chr19_seg8 10 1000000 10.0 1.08 0.0 0.5 100 +chr19 25000001 26000000 chr19_seg9 10 1000000 10.0 1.94 0.97 0.5 100 +chr19 26000001 27000000 chr19_seg10 10 1000000 10.0 1.08 0.0 0.5 100 +chr19 27000001 28000000 chr19_seg11 10 1000000 10.0 1.94 0.97 0.5 100 +chr19 28000001 29000000 chr19_seg12 10 1000000 10.0 1.08 0.0 0.5 100 +chr19 29000001 30000000 chr19_seg13 10 1000000 10.0 1.94 0.97 0.5 100 +chr19 30000001 58617616 chr19_seg14 10 28617616 10.0 1.94 0.97 0.5 100 diff --git a/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz b/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz new file mode 100644 index 0000000000000000000000000000000000000000..1d8532e0cc64935e25f1b09a44b956e70160ef5b GIT binary patch literal 1027 zcmV+e1pNCSiwFquOV4To|8r$_WpZ_MUvqD5VRUI@E_P#P0IXMAZ<|UGe%}5H)I2uH zVHTE4QYyItr$`iY4E9kKRhDaS5~~hafFt$apXC;@_S(3?eu)XYGvCbj&BYfNkE=C( zOwyNmZcXoouj)1aaB-0(+jK##kiOD%o6TO~T4Vnmx?bdrt04asA<_yb(^LAAxh}Z(n}V)Pl?C&xSe3(^dS$8!BE% z$)l9z>2{HC)A?F@ou{jLobMX+9Dm228|U}6^0K03kDJ&lu`z1+QQPftlwb)-m7&6vc%)JPUdoUfirV~gh$U#?45 zEUaQM6@tzEGkS9dh^eT=)Rp}+U^a({p@FV}#?IV*8;`!DU=$+PMq_7)23~~zoeVrw zDHs(gMpbkur&&JBw=YThp}>Veu0ANM#cw+vW{*Dl%)Wg=ntUyzLjEWiJ9lpA+*ZoA z*xTDulNKrN-!frY%vkkMOT+qfh3=ypdy3?L;0AEz5(`}d?{YYBy_q-qFPsFyXdE5( zZoGMCc8xCYIwCUg-oN@b$S}jq|IaeTH@J#I{^=hMpWUhf$hWm}tsRT-2IT6K5QriG zf=eL&LyZJiNFI6s!7U&MGC~AElnO|b617fA4IhBVe*xfxTLRzbQ_&EQz{6)GDlFawvTgVTC{=O#P-QqK94uf=>XV zbHPsnTn0Z0js-snM8OlKG}f7L6V#pHl`8axd7_smPAY{4fgl4)JAN@1*zXW(iRG;t z;4iA-2o|-4hj%U7Ghz9VYI3Pr2 zXD~5|f&6YI)}x~<=_>3c9e|gjQSs7OiU8eORK-ie8=OLn=hzWW0J7Jh0drzBD~0qH zBj9c|79m3mjq8Wyd}e=-2!p+?-Lv`?+tmeQh}+2OLDhc@L085#5R9SrQdZBYGmIhV z%C|;@D`G=`p+fx}8)7H%x+0hp96^|zX-;$!r0atA+$N~*-)Csgodg*~xH@yoquG0T xB;CuSVT#%(|Dy)Z)jm5ZG8IwRYI9WgwR3c>!zBF=009600{}-$QgHzj002bK>LdUF literal 0 HcmV?d00001 diff --git a/tests/fixtures/shatterseek/wakhan/solution_rank_1/integer_profile.bed b/tests/fixtures/shatterseek/wakhan/solution_rank_1/integer_profile.bed new file mode 100644 index 00000000..db01f2b9 --- /dev/null +++ b/tests/fixtures/shatterseek/wakhan/solution_rank_1/integer_profile.bed @@ -0,0 +1,22 @@ +#chr: chromosome number +#start: start address for CN segment +#end: end address for CN segment +#hp1_coverage / hp2_coverage: median coverage for this segment, per haplotype +#hp1_copynumber_state / hp2_copynumber_state: detected copy number state (integer/fraction), per haplotype +#hp1_confidence / hp2_confidence: confidence score, per haplotype +#svs_breakpoints_ids: corresponding structural variations (breakpoints) IDs from VCF file (union across haplotypes) +#chr start end hp1_coverage hp1_copynumber_state hp1_confidence hp2_coverage hp2_copynumber_state hp2_confidence svs_breakpoints_ids +chr2 0 242193529 15.00 1.0 0.9 15.00 1.0 0.9 +chr7 0 159345973 15.00 1.0 0.9 18.00 1.2 0.9 +chr19 0 20000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 20000000 21000000 15.00 1.0 0.9 0.00 0.0 0.9 +chr19 21000000 22000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 22000000 23000000 15.00 1.0 0.9 0.00 0.0 0.9 +chr19 23000000 24000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 24000000 25000000 15.00 1.0 0.9 0.00 0.0 0.9 +chr19 25000000 26000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 26000000 27000000 15.00 1.0 0.9 0.00 0.0 0.9 +chr19 27000000 28000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 28000000 29000000 15.00 1.0 0.9 0.00 0.0 0.9 +chr19 29000000 30000000 15.00 1.0 0.9 15.00 1.0 0.9 +chr19 30000000 58617616 15.00 1.0 0.9 15.00 1.0 0.9 diff --git a/tests/sv_annotation_compose.nf.test b/tests/sv_annotation_compose.nf.test index 3f0ee513..92a1428e 100644 --- a/tests/sv_annotation_compose.nf.test +++ b/tests/sv_annotation_compose.nf.test @@ -1,27 +1,28 @@ nextflow_pipeline { - name "Test pipeline - Padfoot / ReConPlot wiring with every caller on (-preview)" + name "Test pipeline - Padfoot / ReConPlot / ShatterSeek wiring with every caller on (-preview)" script "../main.nf" tag "pipeline" tag "small" - // The test profiles skip Wakhan, ASCAT, Padfoot and ReConPlot (no whole-genome caller output on the test data), - // so no other pipeline test evaluates how the two subworkflows are wired in. -preview builds the whole workflow + // The test profiles skip Wakhan, ASCAT, Padfoot, ReConPlot and ShatterSeek (no whole-genome caller output on the test + // data), so no other pipeline test evaluates how these subworkflows are wired in. -preview builds the whole workflow // graph without running a task, which catches launch errors (e.g. an output a module no longer emits), and // -with-dag writes that graph, whose process nodes ( v12(["NAME"]) ) show which steps the run would contain: a - // regression that silently drops Padfoot or ReConPlot (a gate wrongly false, a skip flag ignored) fails here too. + // regression that silently drops Padfoot, ReConPlot or ShatterSeek (a gate wrongly false, a skip flag ignored) fails here too. - test("-profile test with Wakhan, ASCAT, Padfoot and ReConPlot on - preview") { + test("-profile test with Wakhan, ASCAT, Padfoot, ReConPlot and ShatterSeek on - preview") { options "-preview -with-dag dag.mmd" when { params { - outdir = "$outputDir" - skip_wakhan = false - skip_ascat = false - skip_padfoot = false - skip_reconplot = false + outdir = "$outputDir" + skip_wakhan = false + skip_ascat = false + skip_padfoot = false + skip_reconplot = false + skip_shatterseek = false } } @@ -31,22 +32,24 @@ nextflow_pipeline { { assert workflow.success }, { assert processes.containsAll(['SEVERUS', 'WAKHAN', 'ASCAT']) }, { assert processes.containsAll(['PADFOOT_SEVERUS_WAKHAN', 'PADFOOT_SAVANA']) }, - { assert processes.containsAll(['RECONPLOT_SEVERUS_ASCAT', 'RECONPLOT_SEVERUS_WAKHAN', 'RECONPLOT_SAVANA']) } + { assert processes.containsAll(['RECONPLOT_SEVERUS_ASCAT', 'RECONPLOT_SEVERUS_WAKHAN', 'RECONPLOT_SAVANA']) }, + { assert processes.containsAll(['SHATTERSEEK_SEVERUS_ASCAT', 'SHATTERSEEK_SEVERUS_WAKHAN', 'SHATTERSEEK_SAVANA']) } ) } } - test("-profile test with Wakhan and ReConPlot on, Padfoot skipped - preview") { + test("-profile test with Wakhan and ReConPlot on, Padfoot and ShatterSeek skipped - preview") { options "-preview -with-dag dag.mmd" when { params { - outdir = "$outputDir" - skip_wakhan = false - skip_ascat = false - skip_padfoot = true - skip_reconplot = false + outdir = "$outputDir" + skip_wakhan = false + skip_ascat = false + skip_padfoot = true + skip_reconplot = false + skip_shatterseek = true } } @@ -55,6 +58,7 @@ nextflow_pipeline { assertAll( { assert workflow.success }, { assert !processes.any { name -> name.startsWith('PADFOOT_') } }, + { assert !processes.any { name -> name.startsWith('SHATTERSEEK_') } }, { assert processes.containsAll(['RECONPLOT_SEVERUS_ASCAT', 'RECONPLOT_SEVERUS_WAKHAN', 'RECONPLOT_SAVANA']) } ) } diff --git a/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index 4a606c6f..01f5499d 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -14,6 +14,7 @@ include { reportGenePanelTokens } from '../subworkflows/local/utils_nfcore_lrso include { reportGenePanelIsFile } from '../subworkflows/local/utils_nfcore_lrsomatic_pipeline' include { resolveVepPlugins; validateVepPluginParams } from '../subworkflows/local/utils_nfcore_lrsomatic_pipeline' include { padfootGenome; padfootAnnotationOk } from '../subworkflows/local/utils_nfcore_lrsomatic_pipeline' +include { shatterseekGenome } from '../subworkflows/local/utils_nfcore_lrsomatic_pipeline' include { validateClairstoCnaResources } from '../subworkflows/local/utils_nfcore_lrsomatic_pipeline' include { PREPARE_VEP_PLUGINS } from '../subworkflows/local/prepare_vep_plugins' @@ -63,6 +64,7 @@ include { TUMORONLY_SAVANA } from '../subworkflows/local/tumor_on include { PAIRED_SAVANA } from '../subworkflows/local/paired/paired_savana' include { PADFOOT_ANNOTATION } from '../subworkflows/local/padfoot_annotation' include { RECONPLOT_FIGURES } from '../subworkflows/local/reconplot_figures' +include { SHATTERSEEK_CHROMOTHRIPSIS } from '../subworkflows/local/shatterseek_chromothripsis' @@ -1286,6 +1288,21 @@ workflow LRSOMATIC { ch_versions = ch_versions.mix(RECONPLOT_FIGURES.out.versions) } + // + // SUBWORKFLOW: SHATTERSEEK_CHROMOTHRIPSIS -- ShatterSeek chromothripsis tiers per CN/SV caller pair (same pairs as ReConPlot) + // + if (!params.skip_shatterseek && sv_annot_container_ok && shatterseekGenome()) { + SHATTERSEEK_CHROMOTHRIPSIS ( + SEVERUS.out.somatic_vcf, + params.skip_ascat ? channel.empty() : ASCAT.out.segments, + params.skip_wakhan ? channel.empty() : WAKHAN.out.solution_dirs, + savana_cna, + savana_somatic_vcf, + shatterseekGenome() + ) + ch_versions = ch_versions.mix(SHATTERSEEK_CHROMOTHRIPSIS.out.versions) + } + // // MODULE: LRSOMATICREPORT -- per-sample HTML report; all inputs optional, so joins use remainder: true on the tumor id // From 86e74d5fb188a1a2a78fc1814f521e9d4d0b3d51 Mon Sep 17 00:00:00 2001 From: robert-a-forsyth Date: Wed, 14 Oct 2026 15:46:32 +0200 Subject: [PATCH 3/5] Link the ShatterSeek CHANGELOG entries to #222 Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d4c408d..0177b381 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Added` -- Added ShatterSeek chromothripsis calling for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), published to `/shatterseek//` with a per-chromosome summary, the exact ShatterSeek inputs and a figure per called chromosome; each cluster is labelled with the PCAWG high/low-confidence tiers (HC1, HC2, HC3, LC1; Cortes-Ciriano et al. 2020, ShatterSeek issue #22), with BH q-values at 0.2 and cut-offs overridable through `ext.args`. The pipeline runs the fork [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), which adds T2T-CHM13 support and corrects the statistical criteria to the PCAWG definitions (translocations selected on the right breakend, intra-SV fragment-joins test, all-breakpoint one-sided enrichment, exponential tests on inter-breakpoint gaps, one cluster per chromosome, the strand swap in `SVs()`), inside `docker.io/robertaforsyth/shatterseek` (recipe in `containers/shatterseek/`); the callers are read by the ReConPlot wrapper's parsers. `--skip_shatterseek`, `--shatterseek_genome` and `--shatterseek_min_svlen` (intra-chromosomal SVs shorter than 150 bp are left out by default) control it (@robert-a-forsyth). +- [#222](https://github.com/IntGenomicsLab/lrsomatic/pull/222) - Added ShatterSeek chromothripsis calling for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), published to `/shatterseek//` with a per-chromosome summary, the exact ShatterSeek inputs and a figure per called chromosome; each cluster is labelled with the PCAWG high/low-confidence tiers (HC1, HC2, HC3, LC1; Cortes-Ciriano et al. 2020, ShatterSeek issue #22), with BH q-values at 0.2 and cut-offs overridable through `ext.args`. The pipeline runs the fork [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), which adds T2T-CHM13 support and corrects the statistical criteria to the PCAWG definitions (translocations selected on the right breakend, intra-SV fragment-joins test, all-breakpoint one-sided enrichment, exponential tests on inter-breakpoint gaps, one cluster per chromosome, the strand swap in `SVs()`), inside `docker.io/robertaforsyth/shatterseek` (recipe in `containers/shatterseek/`); the callers are read by the ReConPlot wrapper's parsers. `--skip_shatterseek`, `--shatterseek_genome` and `--shatterseek_min_svlen` (intra-chromosomal SVs shorter than 150 bp are left out by default) control it (@robert-a-forsyth). - [#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). @@ -59,7 +59,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Fixed` -- The ReConPlot wrapper's SAVANA VCF parser re-expresses a `_2` breakend from its mate's side before collapsing mates: SAVANA writes `BP_NOTATION` in (_1, _2) order on both records, so a junction whose `_2` record came first in the file paired the second breakend's position with the first breakend's orientation (@robert-a-forsyth). +- [#222](https://github.com/IntGenomicsLab/lrsomatic/pull/222) - The ReConPlot wrapper's SAVANA VCF parser re-expresses a `_2` breakend from its mate's side before collapsing mates: SAVANA writes `BP_NOTATION` in (_1, _2) order on both records, so a junction whose `_2` record came first in the file paired the second breakend's position with the first breakend's orientation (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `PREPARE_AA_DATA_REPO` no longer reads a `versions` output `UNTAR` does not have, which crashed every default GRCh38 run with CoRAL; overriding `--aa_data_repo_url` drops the genome's default MD5 (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `-profile test_full` starts again: it skips CoRAL and signatures and carries the chr19 SAVANA workarounds from `test` (@robert-a-forsyth). - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `CLAIRSTO_VERDICT_TAG` no longer crashes when ASCAT finds no solution and writes an empty segments table (@robert-a-forsyth). From 23dc8b213f6080497eeeb80d5024e1192bf297ce Mon Sep 17 00:00:00 2001 From: robert-a-forsyth Date: Fri, 30 Oct 2026 10:48:04 +0100 Subject: [PATCH 4/5] Run upstream ShatterSeek with PCAWG high-confidence tiers only Replace the robert-a-forsyth ShatterSeek fork with upstream parklab/ShatterSeek (3b3f4cd, v1.1) plus a build-time patch (containers/shatterseek/patch_chm13.R) that adds T2T-CHM13 cytobands and chromosome sizes and accepts "chm13" in the three genome checks; the statistics are upstream's. The wrapper now computes the PCAWG tiers itself (HC1, HC2, HC3; issue #22) with BH q-values per sample, drops the low-confidence tier (LC1, --min-osc-lc/--max-osc-lc), and seeds R's RNG because upstream's exponential tests draw random numbers. Review fixes: the module uses the singularity/docker ternary with an oras://...-sif twin; CN overlap trimming uses the running maximum end so nested segments no longer slip through; the module tests run the pipeline's default --min-svlen 150 against a new 100 bp fixture deletion. CI: inline the wrapper digests in SHATTERSEEK_CHROMOTHRIPSIS (nextflow-lint rejected the local closure) and exclude the ShatterSeek Wakhan fixture's trailing tab from trailing-whitespace, as for the Padfoot one. Co-Authored-By: Claude Opus 5.5 (1M context) --- .pre-commit-config.yaml | 3 +- CHANGELOG.md | 2 +- assets/shatterseek/README.md | 12 ++-- assets/shatterseek/VERSION | 2 +- assets/shatterseek/run_shatterseek.R | 54 ++++++++++++---- containers/shatterseek/Dockerfile | 31 ++++++--- containers/shatterseek/README.md | 34 +++++----- containers/shatterseek/environment.yml | 3 +- containers/shatterseek/patch_chm13.R | 61 ++++++++++++++++++ docs/output.md | 48 +++++++------- docs/usage.md | 4 +- modules/local/shatterseek/main.nf | 14 ++-- modules/local/shatterseek/meta.yml | 8 +-- modules/local/shatterseek/tests/main.nf.test | 12 +++- .../local/shatterseek/tests/main.nf.test.snap | 30 ++++----- .../local/shatterseek/tests/nextflow.config | 2 + .../local/shatterseek_chromothripsis.nf | 13 ++-- tests/fixtures/shatterseek/make_fixtures.py | 3 + .../savana/test_tumor.classified.somatic.vcf | 2 + .../severus/severus_somatic.vcf.gz | Bin 1027 -> 1040 bytes 20 files changed, 231 insertions(+), 107 deletions(-) create mode 100644 containers/shatterseek/patch_chm13.R diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a83210cc..eca6b76e 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -18,7 +18,8 @@ repos: modules/(?!local/).*| subworkflows/(?!local/).*| .*\.snap$| - tests/fixtures/padfoot_reconplot/wakhan/.* + tests/fixtures/padfoot_reconplot/wakhan/.*| + tests/fixtures/shatterseek/wakhan/.* )$ - id: end-of-file-fixer exclude: | diff --git a/CHANGELOG.md b/CHANGELOG.md index 0177b381..c1ef9254 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Added` -- [#222](https://github.com/IntGenomicsLab/lrsomatic/pull/222) - Added ShatterSeek chromothripsis calling for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), published to `/shatterseek//` with a per-chromosome summary, the exact ShatterSeek inputs and a figure per called chromosome; each cluster is labelled with the PCAWG high/low-confidence tiers (HC1, HC2, HC3, LC1; Cortes-Ciriano et al. 2020, ShatterSeek issue #22), with BH q-values at 0.2 and cut-offs overridable through `ext.args`. The pipeline runs the fork [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), which adds T2T-CHM13 support and corrects the statistical criteria to the PCAWG definitions (translocations selected on the right breakend, intra-SV fragment-joins test, all-breakpoint one-sided enrichment, exponential tests on inter-breakpoint gaps, one cluster per chromosome, the strand swap in `SVs()`), inside `docker.io/robertaforsyth/shatterseek` (recipe in `containers/shatterseek/`); the callers are read by the ReConPlot wrapper's parsers. `--skip_shatterseek`, `--shatterseek_genome` and `--shatterseek_min_svlen` (intra-chromosomal SVs shorter than 150 bp are left out by default) control it (@robert-a-forsyth). +- [#222](https://github.com/IntGenomicsLab/lrsomatic/pull/222) - Added ShatterSeek chromothripsis calling for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), published to `/shatterseek//` with a per-chromosome summary, the exact ShatterSeek inputs and a figure per called chromosome; each cluster is labelled with the PCAWG high-confidence tiers (HC1, HC2, HC3; Cortes-Ciriano et al. 2020, ShatterSeek issue #22), with BH q-values at 0.2 and cut-offs overridable through `ext.args`. The pipeline runs upstream [parklab/ShatterSeek](https://github.com/parklab/ShatterSeek) with a build-time patch adding T2T-CHM13, inside `docker.io/robertaforsyth/shatterseek` (SIF twin `oras://docker.io/robertaforsyth/shatterseek-sif`, recipe in `containers/shatterseek/`); the callers are read by the ReConPlot wrapper's parsers. `--skip_shatterseek`, `--shatterseek_genome` and `--shatterseek_min_svlen` (intra-chromosomal SVs shorter than 150 bp are left out by default) control it (@robert-a-forsyth). - [#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). diff --git a/assets/shatterseek/README.md b/assets/shatterseek/README.md index 23f5fb28..a40141f0 100644 --- a/assets/shatterseek/README.md +++ b/assets/shatterseek/README.md @@ -5,9 +5,11 @@ cluster with the PCAWG confidence tiers and draws a figure for every called chro stages this directory, plus `assets/reconplot/` whose `R/parsers/` read ASCAT, Wakhan, Severus and SAVANA, so ShatterSeek and ReConPlot see the same SVs and segments. -ShatterSeek is the [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support) -fork inside the module's image (recipe in `containers/shatterseek/`): T2T-CHM13 support, the statistical criteria fixed -to the PCAWG definitions, one cluster per chromosome, and `classify_chromothripsis()`. +ShatterSeek is upstream [parklab/ShatterSeek](https://github.com/parklab/ShatterSeek) inside the module's image (recipe +in `containers/shatterseek/`), with a build-time patch adding T2T-CHM13 (`--genome chm13`). The PCAWG high-confidence +tiers (HC1, HC2, HC3; Cortes-Ciriano et al. 2020, Supplementary Note 1.2, and ShatterSeek issue #22) and their BH +q-values are computed here, from ShatterSeek's per-chromosome summary. R's generator is seeded (`--seed`) because +upstream's exponential tests draw random numbers. Conversion to ShatterSeek's format: @@ -15,8 +17,8 @@ Conversion to ShatterSeek's format: outside chromosomes 1-22 and X. Intra-chromosomal SVs are typed by strand (`+-` DEL, `-+` DUP, `++` h2hINV, `--` t2tINV), inter-chromosomal ones are TRA. Breakends are ordered so pos1 <= pos2, strands moving with them. - Copy number: total copy number rounded to an integer (SAVANA's is fractional), chromosomes 1-22 and X, overlaps - trimmed, and neighbouring segments with the same copy number merged (ShatterSeek counts copy-number states row by row). + trimmed (each segment starts after every earlier one ends), and neighbouring segments with the same copy number merged (ShatterSeek counts copy-number states row by row). Both tables are written next to the summary (`.shatterseek_{sv,cn}_input.tsv`). Tier cut-offs default to PCAWG's and can be changed with `--use`, `--alpha`, `--min-sv-hc1`, `--min-osc-hc`, `--min-sv-hc2`, `--min-tra-hc2`, -`--min-sv-hc3`, `--min-osc-lc`, `--max-osc-lc` and `--no-hc3`; see `Rscript run_shatterseek.R --help`. +`--min-sv-hc3` and `--no-hc3`; see `Rscript run_shatterseek.R --help`. diff --git a/assets/shatterseek/VERSION b/assets/shatterseek/VERSION index 3eefcb9d..9084fa2f 100644 --- a/assets/shatterseek/VERSION +++ b/assets/shatterseek/VERSION @@ -1 +1 @@ -1.0.0 +1.1.0 diff --git a/assets/shatterseek/run_shatterseek.R b/assets/shatterseek/run_shatterseek.R index 65dc5afe..f12d5c0a 100644 --- a/assets/shatterseek/run_shatterseek.R +++ b/assets/shatterseek/run_shatterseek.R @@ -3,11 +3,12 @@ ## run_shatterseek.R -- ShatterSeek chromothripsis calls from lrsomatic caller output. ## ## --[ReConPlot wrapper parsers]--> {cn, sv} --[ShatterSeek format]--> shatterseek() -## --> classify_chromothripsis() (PCAWG HC1/HC2/HC3/LC1) --> summary TSV + per-chromosome figures +## --> PCAWG high-confidence tiers (HC1/HC2/HC3) --> summary TSV + per-chromosome figures ## ## Callers are read by the parsers in assets/reconplot/R/parsers (--reconplot-dir), so ShatterSeek and -## ReConPlot see the same SVs and segments. ShatterSeek itself is the robert-a-forsyth fork (chm13 support, -## corrected statistics) installed in the module's container. +## ReConPlot see the same SVs and segments. ShatterSeek itself is upstream parklab/ShatterSeek, installed in the +## module's container with a build-time patch adding T2T-CHM13 (containers/shatterseek/). +## The tiers follow Cortes-Ciriano et al. 2020, Supplementary Note 1.2, as spelled out in parklab/ShatterSeek#22. ## --------------------------------------------------------------------------- suppressPackageStartupMessages({ @@ -36,13 +37,13 @@ option_list <- list( make_option("--min-size", type = "integer", default = 1, help = "shatterseek() min.Size [default %default]"), make_option("--use", type = "character", default = "q", help = "Tier tests on BH q (q) or raw P values (p) [default %default]"), make_option("--alpha", type = "double", default = 0.2, help = "Significance level for the tiers [default %default]"), - make_option("--min-sv-hc1", type = "double", default = 6, help = "HC1/LC1 interleaved intra SVs [default %default]"), + make_option("--seed", type = "integer", default = 1, + help = "RNG seed (ShatterSeek's exponential tests draw random numbers) [default %default]"), + make_option("--min-sv-hc1", type = "double", default = 6, help = "HC1 interleaved intra SVs [default %default]"), make_option("--min-osc-hc", type = "double", default = 7, help = "HC1/HC2 oscillating segments [default %default]"), make_option("--min-sv-hc2", type = "double", default = 3, help = "HC2 interleaved intra SVs [default %default]"), make_option("--min-tra-hc2", type = "double", default = 4, help = "HC2 translocations [default %default]"), make_option("--min-sv-hc3", type = "double", default = 40, help = "HC3 SVs incl. translocations [default %default]"), - make_option("--min-osc-lc", type = "double", default = 4, help = "LC1 minimum oscillating segments [default %default]"), - make_option("--max-osc-lc", type = "double", default = 6, help = "LC1 maximum oscillating segments [default %default]"), make_option("--no-hc3", action = "store_true", default = FALSE, help = "Skip the >= 40 SV tier (PCAWG: sarcomas only)"), make_option("--plot-all", action = "store_true", default = FALSE, help = "Draw every chromosome with a cluster, not only those with a tier") @@ -113,9 +114,12 @@ to_shatterseek_cn <- function(cn) { out <- list() for (ch in unique(d$chrom)) { x <- d[d$chrom == ch, , drop = FALSE] - # no overlaps: each segment starts after the previous one ends - for (i in seq_len(nrow(x))[-1]) if (x$start[i] <= x$end[i - 1]) x$start[i] <- x$end[i - 1] + 1L - x <- x[x$start <= x$end, , drop = FALSE] + # no overlaps: each segment starts after every earlier one ends (running maximum, so nested segments drop out) + if (nrow(x) > 1) { + prev_end <- c(-Inf, cummax(x$end)[-nrow(x)]) + x$start <- as.integer(pmax(x$start, prev_end + 1)) + x <- x[x$start <= x$end, , drop = FALSE] + } # neighbours with equal copy number become one segment (ShatterSeek counts states row by row) run <- cumsum(c(TRUE, diff(x$total_cn) != 0)) out[[ch]] <- data.frame(chrom = ch, start = tapply(x$start, run, min), end = tapply(x$end, run, max), @@ -136,16 +140,42 @@ fwrite(sv_ss, file.path(args$outdir, paste0(prefix, ".shatterseek_sv_input.tsv") fwrite(cn_ss, file.path(args$outdir, paste0(prefix, ".shatterseek_cn_input.tsv")), sep = "\t") ## --- 4. ShatterSeek + PCAWG tiers --------------------------------------------- +## PCAWG high-confidence tiers (Cortes-Ciriano et al. 2020, Supplementary Note 1.2; parklab/ShatterSeek#22). Tests +## are significant at <= alpha on BH q (across the sample's chromosomes) or raw P; NA counts as a fail. +classify_pcawg <- function(d, use, alpha, min_sv_hc1, min_osc_hc, min_sv_hc2, min_tra_hc2, min_sv_hc3, hc3) { + for (stem in c("fragment_joins", "chr_breakpoint_enrichment", "exp_chr", "exp_cluster")) { + p_col <- if (stem == "chr_breakpoint_enrichment") stem else paste0("pval_", stem) + d[[paste0("q_", stem)]] <- p.adjust(d[[p_col]], method = "BH") + } + test <- function(stem) d[[if (use == "q") paste0("q_", stem) else if (stem == "chr_breakpoint_enrichment") stem + else paste0("pval_", stem)]] + ok <- function(x) !is.na(x) & x + + n_intra <- d$number_DEL + d$number_DUP + d$number_h2hINV + d$number_t2tINV + osc <- d$max_number_oscillating_CN_segments_2_states + joins_random <- ok(test("fragment_joins") > alpha) + clustered <- ok(test("chr_breakpoint_enrichment") <= alpha) | ok(test("exp_chr") <= alpha) + + d$HC1 <- ok(n_intra >= min_sv_hc1) & ok(osc >= min_osc_hc) & joins_random & clustered + d$HC2 <- ok(n_intra >= min_sv_hc2) & ok(d$number_TRA >= min_tra_hc2) & ok(osc >= min_osc_hc) & joins_random + d$HC3 <- hc3 & ok(d$clusterSize_including_TRA >= min_sv_hc3) & joins_random + + tiers <- c("HC1", "HC2", "HC3") + d$chromothripsis_criteria <- apply(as.matrix(d[, tiers]), 1, function(r) paste(tiers[r], collapse = ",")) + d$chromothripsis_confidence <- ifelse(d$HC1 | d$HC2 | d$HC3, "high", "none") + d +} + +set.seed(args$seed) res <- shatterseek( SV.sample = SVs(chrom1 = sv_ss$chrom1, pos1 = as.numeric(sv_ss$pos1), chrom2 = sv_ss$chrom2, pos2 = as.numeric(sv_ss$pos2), SVtype = sv_ss$SVtype, strand1 = sv_ss$strand1, strand2 = sv_ss$strand2), seg.sample = CNVsegs(chrom = cn_ss$chrom, start = as.numeric(cn_ss$start), end = as.numeric(cn_ss$end), total_cn = as.numeric(cn_ss$total_cn)), min.Size = args$min_size, genome = args$genome) -res@chromSummary <- classify_chromothripsis( +res@chromSummary <- classify_pcawg( res@chromSummary, use = args$use, alpha = args$alpha, min_sv_hc1 = args$min_sv_hc1, min_osc_hc = args$min_osc_hc, - min_sv_hc2 = args$min_sv_hc2, min_tra_hc2 = args$min_tra_hc2, min_sv_hc3 = args$min_sv_hc3, - min_osc_lc = args$min_osc_lc, max_osc_lc = args$max_osc_lc, hc3 = !isTRUE(args$no_hc3)) + min_sv_hc2 = args$min_sv_hc2, min_tra_hc2 = args$min_tra_hc2, min_sv_hc3 = args$min_sv_hc3, hc3 = !isTRUE(args$no_hc3)) summary <- res@chromSummary summary$inter_other_chroms_coords_all <- trimws(gsub("[\n;]+", ";", summary$inter_other_chroms_coords_all)) diff --git a/containers/shatterseek/Dockerfile b/containers/shatterseek/Dockerfile index 1a194caa..c330d3b9 100644 --- a/containers/shatterseek/Dockerfile +++ b/containers/shatterseek/Dockerfile @@ -1,24 +1,37 @@ # syntax=docker/dockerfile:1.6 FROM docker.io/mambaorg/micromamba:2.3.2 -# R runtime for the ShatterSeek wrapper, with the robert-a-forsyth/ShatterSeek fork (chm13 support, statistical -# criteria fixed to the PCAWG definitions, PCAWG confidence tiers) installed from a pinned commit. The wrapper ships -# with the pipeline in assets/shatterseek/. +# R runtime for the ShatterSeek wrapper, with upstream parklab/ShatterSeek (not on conda) installed from a pinned commit +# plus a build-time patch adding T2T-CHM13 (patch_chm13.R). The wrapper ships with the pipeline in assets/shatterseek/. COPY --chown=$MAMBA_USER containers/shatterseek/environment.yml /tmp/environment.yml +COPY --chown=$MAMBA_USER containers/shatterseek/patch_chm13.R /tmp/patch_chm13.R -ARG SHATTERSEEK_COMMIT=e1cee6c8ee75ea1a1e9f4d1c57b7a8a0a199ae91 -ADD --checksum=sha256:74b4258a0e6b509aeba4976348feebf910c098b9556e85008795d1b4e6286b69 --chown=$MAMBA_USER \ - https://github.com/robert-a-forsyth/ShatterSeek/archive/${SHATTERSEEK_COMMIT}.tar.gz /tmp/shatterseek.tar.gz +ARG SHATTERSEEK_COMMIT=3b3f4cddb23ff06a5d14d2ab5633d184e4f8d122 +ADD --checksum=sha256:85d17017e144784f2bf8714c1746b91d2238d96bc3e5ae79195ae4c6bbeb7e4d --chown=$MAMBA_USER \ + https://github.com/parklab/ShatterSeek/archive/${SHATTERSEEK_COMMIT}.tar.gz /tmp/shatterseek.tar.gz +# T2T-CHM13 v2.0 cytobands, for the chm13 ideogram and chromosome sizes +ADD --checksum=sha256:09c7fa3ca3e222ebaabbecedc3086aa2543999824c0966f96b0c3f31c1126ea2 --chown=$MAMBA_USER \ + https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/annotation/chm13v2.0_cytobands_allchrs.bed /tmp/chm13_cytobands.bed RUN micromamba create --yes --name shatterseek --file /tmp/environment.yml \ && micromamba clean --all --yes \ && export PATH=/opt/conda/envs/shatterseek/bin:${PATH} \ && mkdir -p /tmp/shatterseek-src \ && tar -xzf /tmp/shatterseek.tar.gz -C /tmp/shatterseek-src --strip-components 1 \ + && (cd /tmp/shatterseek-src && Rscript /tmp/patch_chm13.R /tmp/chm13_cytobands.bed) \ && R CMD INSTALL --no-docs --no-html /tmp/shatterseek-src \ - # the fork's regression test: DO17373 against the PCAWG-published values, plus synthetic hg38 / chm13 events - && Rscript /tmp/shatterseek-src/tests/test_do17373.R \ - && rm -rf /tmp/shatterseek.tar.gz /tmp/shatterseek-src + # smoke test: the bundled DO17373 example runs under every supported genome + && Rscript -e 'suppressPackageStartupMessages(library(ShatterSeek)); data(DO17373); \ + for (g in c("hg19", "hg38", "chm13")) { set.seed(1); \ + r <- shatterseek(SV.sample = SVs(chrom1 = as.character(SV_DO17373$chrom1), pos1 = as.numeric(SV_DO17373$start1), \ + chrom2 = as.character(SV_DO17373$chrom2), pos2 = as.numeric(SV_DO17373$end2), \ + SVtype = as.character(SV_DO17373$svclass), strand1 = as.character(SV_DO17373$strand1), \ + strand2 = as.character(SV_DO17373$strand2)), \ + seg.sample = CNVsegs(chrom = as.character(SCNA_DO17373$chromosome), start = SCNA_DO17373$start, \ + end = SCNA_DO17373$end, total_cn = SCNA_DO17373$total_cn), genome = g); \ + stopifnot(nrow(r@chromSummary) == 23); invisible(plot_chromothripsis(r, chr = "2", genome = g)) }; \ + cat("ShatterSeek", as.character(packageVersion("ShatterSeek")), "OK\n")' \ + && rm -rf /tmp/shatterseek.tar.gz /tmp/shatterseek-src /tmp/chm13_cytobands.bed /tmp/patch_chm13.R ENV PATH=/opt/conda/envs/shatterseek/bin:${PATH} \ SHATTERSEEK_COMMIT=${SHATTERSEEK_COMMIT} diff --git a/containers/shatterseek/README.md b/containers/shatterseek/README.md index 12adceae..8b8cdb02 100644 --- a/containers/shatterseek/README.md +++ b/containers/shatterseek/README.md @@ -1,29 +1,31 @@ # ShatterSeek image -R runtime for the ShatterSeek wrapper used by the `SHATTERSEEK` module, with the -[robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support) fork installed from -a pinned, checksum-verified commit (`SHATTERSEEK_COMMIT` in the Dockerfile). Relative to upstream -[parklab/ShatterSeek](https://github.com/parklab/ShatterSeek) the fork adds T2T-CHM13 (`genome = "chm13"`), corrects -the statistical criteria to the PCAWG definitions, keeps one cluster per chromosome and labels clusters with the PCAWG -confidence tiers (`classify_chromothripsis()`); its README lists every change. The build runs the fork's regression -test (`tests/test_do17373.R`: the DO17373 example against the PCAWG-published values, plus synthetic hg38 and chm13 -events) and fails if it does not pass. The wrapper is part of the pipeline in `assets/shatterseek/`. The module does -not support `-profile conda`: the image ships the package. +R runtime for the ShatterSeek wrapper used by the `SHATTERSEEK` module, with upstream +[parklab/ShatterSeek](https://github.com/parklab/ShatterSeek) (v1.1) installed from a pinned, checksum-verified commit +(`SHATTERSEEK_COMMIT` in the Dockerfile). Upstream supports hg19 and hg38 only; `patch_chm13.R` adds T2T-CHM13 +(`genome = "chm13"`) at build time: the CHM13 v2.0 cytobands for the ideogram, the chromosome lengths for the +breakpoint-enrichment test, and `chm13` in the three genome checks. The statistics are upstream's, unchanged. The build +fails if the patch does not apply exactly, and runs the bundled DO17373 example under hg19, hg38 and chm13. The wrapper +is part of the pipeline in `assets/shatterseek/`. The module does not support `-profile conda`: the image ships the package. -Licence: the fork's `DESCRIPTION` declares GPL (>= 3), as upstream does, while the upstream README states that -ShatterSeek is free for academic use only. +Licence: ShatterSeek's `DESCRIPTION` declares GPL (>= 3), while the upstream README states that it is free for academic +use only. Build and publish from the pipeline root (on the VSC login nodes: podman with the cgroupfs/file flags, TMPDIR on local -scratch). After a new fork commit, update `SHATTERSEEK_COMMIT` and the `ADD --checksum` (the sha256 of -`https://github.com/robert-a-forsyth/ShatterSeek/archive/.tar.gz`) first: +scratch). For a new upstream commit, update `SHATTERSEEK_COMMIT` and the `ADD --checksum` (the sha256 of +`https://github.com/parklab/ShatterSeek/archive/.tar.gz`) first. The SIF twin is what Singularity/Apptainer +pull, as for the ReConPlot and Padfoot images: ```bash -TAG=1.1.1-chm13- +TAG=1.1-3b3f4cd-chm13 podman --cgroup-manager=cgroupfs --events-backend=file build --format docker \ -f containers/shatterseek/Dockerfile -t "docker.io/robertaforsyth/shatterseek:$TAG" . podman --cgroup-manager=cgroupfs --events-backend=file push "docker.io/robertaforsyth/shatterseek:$TAG" +# the SIF twin, built from the image just pushed and published with ORAS +# (`singularity remote login --username oras://docker.io` first) +singularity build "shatterseek-$TAG.sif" "docker://docker.io/robertaforsyth/shatterseek:$TAG" +singularity push "shatterseek-$TAG.sif" "oras://docker.io/robertaforsyth/shatterseek-sif:$TAG" ``` -There is no SIF twin: like CoRAL, Singularity/Apptainer pull the Docker image. The module pins the tag in -`modules/local/shatterseek/main.nf`; override per site via +The module pins both tags in `modules/local/shatterseek/main.nf`; override per site via `process { withName: '.*:SHATTERSEEK_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = ... } }`. diff --git a/containers/shatterseek/environment.yml b/containers/shatterseek/environment.yml index fc9fbf52..d282f70f 100644 --- a/containers/shatterseek/environment.yml +++ b/containers/shatterseek/environment.yml @@ -1,7 +1,7 @@ --- # yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json # R runtime for the ShatterSeek wrapper (assets/shatterseek), used only to build containers/shatterseek/Dockerfile; -# ShatterSeek itself (the robert-a-forsyth fork, not on conda) is installed into the image from a pinned commit. +# ShatterSeek itself (upstream parklab/ShatterSeek, not on conda) is installed into the image from a pinned commit. channels: - conda-forge - bioconda @@ -13,6 +13,7 @@ dependencies: - "bioconda::bioconductor-genomicranges=1.58.0" - "bioconda::bioconductor-graph=1.84.0" - "conda-forge::r-foreach=1.5.2" + - "conda-forge::r-mass=7.3_66" - "conda-forge::r-ggplot2=4.0.3" - "conda-forge::r-gridextra=2.3" - "conda-forge::r-data.table=1.18.4" diff --git a/containers/shatterseek/patch_chm13.R b/containers/shatterseek/patch_chm13.R new file mode 100644 index 00000000..e64479b3 --- /dev/null +++ b/containers/shatterseek/patch_chm13.R @@ -0,0 +1,61 @@ +#!/usr/bin/env Rscript +# Build-time patch adding T2T-CHM13 to upstream parklab/ShatterSeek. Run from the unpacked package root: +# Rscript patch_chm13.R +# +# Upstream only knows hg19 and hg38, through three genome checks and two per-genome tables in R/sysdata.rda: +# info_mappa_hg38 (chromosome sizes for the breakpoint-enrichment test) and chr_info_hg38 (ideogram bands). +# This adds info_mappa_chm13 (CHM13 has no gaps, so the sizes are the chromosome lengths) and chr_info_chm13 in the same +# layout, and lets the three checks accept "chm13". The statistics are untouched. Fails if any edit does not apply. + +args <- commandArgs(trailingOnly = TRUE) +stopifnot(length(args) == 1, file.exists(args[1])) + +cyto <- read.delim(args[1], header = FALSE, stringsAsFactors = FALSE, + col.names = c("seqnames", "start", "end", "name", "gieStain")) + +env <- new.env() +load("R/sysdata.rda", envir = env) +stopifnot(exists("info_mappa_hg38", envir = env), exists("chr_info_hg38", envir = env)) + +# chr_info_chm13: the bands, in chr_info_hg38's columns and chromosome naming +hg38 <- env$chr_info_hg38 +cyto <- cyto[sub("^chr", "", cyto$seqnames) %in% sub("^chr", "", unique(as.character(hg38$seqnames))), ] +if (!any(grepl("^chr", hg38$seqnames))) cyto$seqnames <- sub("^chr", "", cyto$seqnames) +cyto$width <- cyto$end - cyto$start + 1L # as in chr_info_hg38 +cyto$strand <- "*" +cyto$y <- 1L +chr_info_chm13 <- cyto[, intersect(names(hg38), names(cyto)), drop = FALSE] +stopifnot(identical(names(chr_info_chm13), names(hg38))) +if (is.factor(hg38$seqnames)) chr_info_chm13$seqnames <- factor(chr_info_chm13$seqnames, levels = unique(chr_info_chm13$seqnames)) + +# info_mappa_chm13: one row per chromosome of info_mappa_hg38 (1-22, X, Y, M), tot = chromosome length +mappa_hg38 <- env$info_mappa_hg38 +chm13_len <- c(tapply(cyto$end, sub("^chr", "", cyto$seqnames), max), M = 16569) # chrM (rCRS) has no bands +info_mappa_chm13 <- mappa_hg38 +info_mappa_chm13$tot <- as.numeric(chm13_len[sub("^chr", "", as.character(mappa_hg38$V1))]) +stopifnot(!anyNA(info_mappa_chm13$tot), all(info_mappa_chm13$tot > 0)) +for (col in setdiff(names(info_mappa_chm13), c("V1", "tot"))) info_mappa_chm13[[col]] <- NA + +assign("chr_info_chm13", chr_info_chm13, envir = env) +assign("info_mappa_chm13", info_mappa_chm13, envir = env) +save(list = ls(env, all.names = TRUE), envir = env, file = "R/sysdata.rda", compress = "bzip2") + +# the three genome checks: (file, exact upstream text, replacement) +edits <- list( + list("R/shatterseek.R", + 'if ( !(as.character(genome) %in% c("hg19","hg38"))){stop("Reference genome assembly is not supported (Use hg19 or hg38)")}', + 'if ( !(as.character(genome) %in% c("hg19","hg38","chm13"))){stop("Reference genome assembly is not supported (Use hg19, hg38 or chm13)")}'), + list("R/function_criteria.R", + 'if(genome == "hg38"){', + 'if(genome == "chm13"){\n prob_cand = info_mappa_chm13$tot[which(info_mappa_chm13$V1 == cand)] / sum(as.numeric(info_mappa_chm13$tot))\n }else if(genome == "hg38"){'), + list("R/plot_function.R", + '}else if(genome == "hg38"){\n\t chr_info <- chr_info_hg38', + '}else if(genome == "chm13"){\n\t chr_info <- chr_info_chm13\n\t}else if(genome == "hg38"){\n\t chr_info <- chr_info_hg38') +) +for (e in edits) { + src <- paste(readLines(e[[1]], warn = FALSE), collapse = "\n") + n <- lengths(regmatches(src, gregexpr(e[[2]], src, fixed = TRUE))) + if (n != 1) stop(sprintf("%s: expected the upstream text once, found it %d times", e[[1]], n)) + writeLines(sub(e[[2]], e[[3]], src, fixed = TRUE), e[[1]]) +} +cat("ShatterSeek patched for chm13:", nrow(chr_info_chm13), "bands,", nrow(info_mappa_chm13), "chromosome sizes\n") diff --git a/docs/output.md b/docs/output.md index 65c35dd9..4ecd9612 100644 --- a/docs/output.md +++ b/docs/output.md @@ -803,35 +803,35 @@ A figure ReConPlot cannot draw is skipped and the task keeps the others: ReConPl | `plots/*.pdf` | ShatterSeek figure (ideogram, SVs, copy number, summary) for every chromosome with a confidence tier (optional) | | `shatterseek.log` | Wrapper log: files read, SVs dropped and why, copy numbers rounded, the tiers called | -[ShatterSeek](https://github.com/parklab/ShatterSeek) (Cortes-Ciriano et al. 2020) finds, per chromosome, the largest cluster of interleaved intra-chromosomal SVs and tests it against the hallmarks of chromothripsis. The pipeline runs a fork ([robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support)) whose statistics follow the definitions in the PCAWG Supplementary Note; they differ from upstream ShatterSeek: - -| Column | Meaning | -| ------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------ | -| `start`, `end`, `clusterSize` | Span and number of intra-chromosomal SVs of the chromosome's largest cluster (one cluster; NA/0 without one) | -| `number_DEL`/`_DUP`/`_h2hINV`/`_t2tINV` | The cluster's SVs by type | -| `number_TRA`, `clusterSize_including_TRA` | Translocations with a breakend inside the cluster, and the cluster size including them | -| `pval_fragment_joins` | Chi-square test of the cluster's four SV types against equal proportions. Chromothripsis expects random joins, so this should **not** be significant | -| `chr_breakpoint_enrichment` | One-sided binomial test: more breakpoints on this chromosome than its share of the genome predicts (hg38: non-N bases; chm13: length) | -| `pval_exp_chr`, `pval_exp_cluster` | Kolmogorov-Smirnov test of the gaps between consecutive breakpoints (chromosome / cluster) against an exponential; significant = clustered breakpoints | -| `max_number_oscillating_CN_segments_2_states` | Longest run of adjacent copy-number segments in the cluster region alternating between two states | -| `q_*` | Benjamini-Hochberg q-values of the four tests across the sample's chromosomes | -| `inter_*` | Other chromosomes linked to this one by translocations, and the fragment-joins test over the linked clusters | -| `HC1`, `HC2`, `HC3`, `LC1` | PCAWG confidence tiers met (see below) | -| `chromothripsis_criteria`, `chromothripsis_confidence` | The tiers met, e.g. `HC1,HC2`, and `high` (any HC), `low` (LC1 only) or `none` | - -The tiers use the PCAWG cut-offs (Cortes-Ciriano et al. 2020, Supplementary Note 1.2; ShatterSeek issue #22), with significance at q ≤ 0.2 and NA counting as a fail: - -| Tier | Interleaved intra SVs | Other | Oscillating CN segments (2 states) | Fragment joins | Enrichment or exponential test | -| ---- | --------------------- | ----------------------------- | ---------------------------------- | -------------- | ------------------------------ | -| HC1 | ≥ 6 | | ≥ 7 | q > 0.2 | q ≤ 0.2 | -| HC2 | ≥ 3 | ≥ 4 translocations | ≥ 7 | q > 0.2 | | -| HC3 | | ≥ 40 SVs incl. translocations | | q > 0.2 | | -| LC1 | ≥ 6 | | 4-6 | q > 0.2 | q ≤ 0.2 | +[ShatterSeek](https://github.com/parklab/ShatterSeek) (Cortes-Ciriano et al. 2020) finds, per chromosome, the largest cluster of interleaved intra-chromosomal SVs and tests it against the hallmarks of chromothripsis. The pipeline runs upstream ShatterSeek (v1.1, pinned commit) unchanged, apart from a build-time patch adding T2T-CHM13 chromosome sizes and cytobands. The columns are upstream's: + +| Column | Meaning | +| ------------------------------------------------------ | -------------------------------------------------------------------------------------------------------------------------------------------------- | +| `start`, `end` | Span of the chromosome's largest cluster of interleaved intra-chromosomal SVs (NA without one) | +| `number_DEL`/`_DUP`/`_h2hINV`/`_t2tINV` | The cluster's SVs by type | +| `number_TRA`, `clusterSize_including_TRA` | Translocations with a breakend inside the cluster, and the cluster size including them | +| `pval_fragment_joins` | Chi-square test of the cluster's SV types (translocations typed by strand) against equal proportions; chromothripsis should **not** be significant | +| `chr_breakpoint_enrichment` | Binomial test of the translocations touching this chromosome against its share of the genome | +| `pval_exp_chr`, `pval_exp_cluster` | Kolmogorov-Smirnov tests of the breakpoints (chromosome / cluster) against an exponential distribution; significant = clustered breakpoints | +| `max_number_oscillating_CN_segments_2_states` | Longest run of adjacent copy-number segments in the cluster region alternating between two states | +| `q_*` | Benjamini-Hochberg q-values of the four tests across the sample's chromosomes (added by the wrapper) | +| `inter_*` | Other chromosomes linked to this one by translocations, and the fragment-joins test over the linked clusters | +| `HC1`, `HC2`, `HC3` | PCAWG high-confidence tiers met (see below) | +| `chromothripsis_criteria`, `chromothripsis_confidence` | The tiers met, e.g. `HC1,HC2`, and `high` (any tier) or `none` | + +The tiers use the PCAWG high-confidence cut-offs (Cortes-Ciriano et al. 2020, Supplementary Note 1.2; ShatterSeek issue #22), with significance at q ≤ 0.2 and NA counting as a fail. Interleaved intra SVs is the sum of the four SV-type columns. PCAWG's low-confidence tier (4-6 oscillating segments) is not reported. + +| Tier | Interleaved intra SVs | Other | Oscillating CN segments (2 states) | Fragment joins | Enrichment or `exp_chr` test | +| ---- | --------------------- | ----------------------------- | ---------------------------------- | -------------- | ---------------------------- | +| HC1 | ≥ 6 | | ≥ 7 | q > 0.2 | q ≤ 0.2 | +| HC2 | ≥ 3 | ≥ 4 translocations | ≥ 7 | q > 0.2 | | +| HC3 | | ≥ 40 SVs incl. translocations | | q > 0.2 | | Points to keep in mind when reading the calls: - PCAWG applied HC3 only to osteosarcomas, leiomyosarcomas and liposarcomas, where amplification erodes the oscillating pattern; here it applies to every tumour. Drop it with `ext.args = '--no-hc3'`. - PCAWG corrected the P values once across every chromosome of all 2,658 tumours; within one sample the correction covers 23 chromosomes. The raw P values are kept, so a cohort-wide correction can be recomputed from the summaries (each row carries `sample`, `cn_source` and `sv_source`). +- Upstream's exponential tests compare against random draws; the wrapper seeds R's generator (`--seed`, default 1), so a rerun gives the same P values. Clusters of equal size are pooled before the statistics are computed, so a chromosome whose SVs are all singletons is summarised over all of them. - The cut-offs were set by visual inspection of PCAWG's short-read consensus calls. Long-read callers report many small SVs, whose two breakpoints a few hundred bp apart make the exponential test significant on most chromosomes; intra-chromosomal SVs shorter than `--shatterseek_min_svlen` (default 150 bp) are therefore left out, and raising it tests larger rearrangements only. - ShatterSeek does not mask centromeres or satellite arrays. A cluster spanning a centromere or heterochromatin block (e.g. chr16 q11.2, chr1 q12, chr9 q12) is more likely mapping artefact than chromothripsis; check the figure. - Tumour-only Severus calls include germline SVs that the panel of normals did not remove. diff --git a/docs/usage.md b/docs/usage.md index 8146be08..76142fcc 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -460,14 +460,14 @@ Each pair produces `per_chromosome/` (one figure per chromosome), `genome_wide/` #### ShatterSeek Options -[ShatterSeek](https://github.com/parklab/ShatterSeek) chromothripsis calls are made for every CN/SV caller pair available for a sample (`severus_ascat`, `severus_wakhan`, `savana`, the same pairs as ReConPlot), into `shatterseek//`, and each chromosome's cluster is labelled with the PCAWG high/low-confidence tiers (see [output](output.md#shatterseek)). The pipeline runs a fork, [robert-a-forsyth/ShatterSeek](https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support), with T2T-CHM13 support and the statistical criteria corrected to the PCAWG definitions; it ships inside the module's image (`docker.io/robertaforsyth/shatterseek`, recipe in `containers/shatterseek/`), so the module is not available under `-profile conda`. The callers are read by the same parsers as ReConPlot (`assets/reconplot/R/parsers`), through the wrapper in `assets/shatterseek/`. +[ShatterSeek](https://github.com/parklab/ShatterSeek) chromothripsis calls are made for every CN/SV caller pair available for a sample (`severus_ascat`, `severus_wakhan`, `savana`, the same pairs as ReConPlot), into `shatterseek//`, and each chromosome's cluster is labelled with the PCAWG high-confidence tiers (see [output](output.md#shatterseek)). The pipeline runs upstream ShatterSeek with a build-time patch adding T2T-CHM13; it ships inside the module's image (`oras://docker.io/robertaforsyth/shatterseek-sif` under Singularity/Apptainer, `docker.io/robertaforsyth/shatterseek` otherwise, recipe in `containers/shatterseek/`), so the module is not available under `-profile conda`. The callers are read by the same parsers as ReConPlot (`assets/reconplot/R/parsers`), through the wrapper in `assets/shatterseek/`. | Parameter | Description | | ------------------------- | --------------------------------------------------------------------------------------------------------------------------- | | `--shatterseek_genome` | ShatterSeek genome (`hg38`, `chm13`, `hg19`). Default = `null` (inferred from `--genome`; skipped with a warning otherwise) | | `--shatterseek_min_svlen` | Drop intra-chromosomal SVs shorter than this (bp) before ShatterSeek; translocations kept. Default = `150` | -Tier cut-offs (`--min-sv-hc1`, `--min-osc-hc`, `--min-sv-hc2`, `--min-tra-hc2`, `--min-sv-hc3`, `--min-osc-lc`, `--max-osc-lc`, `--no-hc3`), the significance level (`--alpha`, default 0.2), q-values or raw P values (`--use q|p`) and the other SV filters (`--min-support`, `--exclude-vntr`) are wrapper options, set through `ext.args` (keep `--min-svlen` in it, or use `--shatterseek_min_svlen`, since `ext.args` replaces the pipeline's value): +Tier cut-offs (`--min-sv-hc1`, `--min-osc-hc`, `--min-sv-hc2`, `--min-tra-hc2`, `--min-sv-hc3`, `--no-hc3`), the significance level (`--alpha`, default 0.2), q-values or raw P values (`--use q|p`) and the other SV filters (`--min-support`, `--exclude-vntr`) are wrapper options, set through `ext.args` (keep `--min-svlen` in it, or use `--shatterseek_min_svlen`, since `ext.args` replaces the pipeline's value): ```groovy process { diff --git a/modules/local/shatterseek/main.nf b/modules/local/shatterseek/main.nf index e2675a78..10620ef4 100644 --- a/modules/local/shatterseek/main.nf +++ b/modules/local/shatterseek/main.nf @@ -2,10 +2,12 @@ process SHATTERSEEK { tag "${meta.id}:${cn_source}+${sv_source}" label 'process_low' - // No conda: the ShatterSeek fork (chm13, corrected statistics, PCAWG tiers) ships only in this image, built from + // No conda: ShatterSeek (upstream, plus a build-time chm13 patch) ships only in this image, built from // containers/shatterseek/Dockerfile; the wrapper is pipeline glue in assets/shatterseek/ and reads the callers with - // the ReConPlot parsers in assets/reconplot/. Fork update = new commit there, rebuild, bump this tag. - container "docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c" + // the ReConPlot parsers in assets/reconplot/. Package update = new commit there, rebuild, bump these two tags. + container "${(workflow.containerEngine == 'singularity' || workflow.containerEngine == 'apptainer') && !task.ext.singularity_pull_docker_container + ? 'oras://docker.io/robertaforsyth/shatterseek-sif:1.1-3b3f4cd-chm13' + : 'docker.io/robertaforsyth/shatterseek:1.1-3b3f4cd-chm13'}" input: // cn_files / sv_files are the caller output files the parsers discover by name. @@ -29,7 +31,7 @@ process SHATTERSEEK { script: // Exit if running this module with -profile conda / -profile mamba if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { - error "SHATTERSEEK does not support Conda: the ShatterSeek fork ships only inside its container. Use Docker / Singularity / Apptainer, or --skip_shatterseek." + error "SHATTERSEEK does not support Conda: ShatterSeek ships only inside its container. Use Docker / Singularity / Apptainer, or --skip_shatterseek." } def args = task.ext.args ?: '' prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}") @@ -46,7 +48,7 @@ process SHATTERSEEK { cat <<-END_VERSIONS > versions.yml "${task.process}": shatterseek: \$(Rscript -e 'cat(as.character(packageVersion("ShatterSeek")))' 2>/dev/null) - shatterseek_fork_commit: \${SHATTERSEEK_COMMIT:-unknown} + shatterseek_commit: \${SHATTERSEEK_COMMIT:-unknown} shatterseek_wrapper: \$(cat ${shatterseek_src}/VERSION 2>/dev/null || echo unknown) r-base: \$(Rscript -e 'cat(R.version\$major, R.version\$minor, sep=".")' 2>/dev/null) END_VERSIONS @@ -63,7 +65,7 @@ process SHATTERSEEK { cat <<-END_VERSIONS > versions.yml "${task.process}": shatterseek: stub - shatterseek_fork_commit: stub + shatterseek_commit: stub shatterseek_wrapper: stub r-base: stub END_VERSIONS diff --git a/modules/local/shatterseek/meta.yml b/modules/local/shatterseek/meta.yml index b08929ad..e06c16a6 100644 --- a/modules/local/shatterseek/meta.yml +++ b/modules/local/shatterseek/meta.yml @@ -1,7 +1,7 @@ --- # yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json name: "shatterseek" -description: Chromothripsis calls with ShatterSeek (robert-a-forsyth fork with chm13 support and corrected statistics), labelled with the PCAWG high/low-confidence tiers, from ASCAT/Wakhan/SAVANA copy number and Severus/SAVANA SVs. +description: Chromothripsis calls with ShatterSeek (upstream, with a build-time chm13 patch), labelled with the PCAWG high-confidence tiers, from ASCAT/Wakhan/SAVANA copy number and Severus/SAVANA SVs. keywords: - chromothripsis - structural-variation @@ -12,8 +12,8 @@ tools: - "shatterseek": description: "ShatterSeek: detection of chromothripsis from SV and copy-number calls" homepage: "https://github.com/parklab/ShatterSeek" - documentation: "https://github.com/robert-a-forsyth/ShatterSeek/tree/chm13-support" - tool_dev_url: "https://github.com/robert-a-forsyth/ShatterSeek" + documentation: "https://github.com/parklab/ShatterSeek/blob/master/tutorial.pdf" + tool_dev_url: "https://github.com/parklab/ShatterSeek" doi: "10.1038/s41588-019-0576-7" licence: ["GPL-3.0-or-later"] identifier: "" @@ -57,7 +57,7 @@ output: description: Groovy Map containing sample information - "${prefix}/${meta.id}.shatterseek.tsv": type: file - description: One row per chromosome (1-22, X) with the cluster, the statistical criteria, BH q-values, the PCAWG tiers (HC1, HC2, HC3, LC1) and chromothripsis_confidence + description: One row per chromosome (1-22, X) with the cluster, the statistical criteria, BH q-values, the PCAWG tiers (HC1, HC2, HC3) and chromothripsis_confidence pattern: "*.shatterseek.tsv" - inputs: - meta: diff --git a/modules/local/shatterseek/tests/main.nf.test b/modules/local/shatterseek/tests/main.nf.test index 34e33745..c7e7cb68 100644 --- a/modules/local/shatterseek/tests/main.nf.test +++ b/modules/local/shatterseek/tests/main.nf.test @@ -11,9 +11,10 @@ nextflow_process { // "small" is what .github/workflows/nf-test.yml selects on for pull_request tag "small" - // The stub tests check the I/O contract; the "real" tests run the wrapper and the ShatterSeek fork on the synthetic - // event in tests/fixtures/shatterseek (chr19 chromothripsis linked to chr7, an isolated chr2 deletion), written in - // each caller's own layout. Every pair must call chr19 high confidence (HC1, HC2) and nothing else. + // The stub tests check the I/O contract; the "real" tests run the wrapper and ShatterSeek on the synthetic event in + // tests/fixtures/shatterseek (chr19 chromothripsis linked to chr7, an isolated chr2 deletion, a 100 bp chr2 deletion), + // written in each caller's own layout, with the pipeline's default --min-svlen 150 (./nextflow.config). Every pair + // must call chr19 high confidence (HC1, HC2) and nothing else. test("ascat + severus - stub") { @@ -95,6 +96,7 @@ nextflow_process { then { def rows = path(process.out.summary.get(0).get(1)).readLines() + def sv_rows = path(process.out.inputs.get(0).get(1).find { f -> f.toString().endsWith('_sv_input.tsv') }).readLines() def head = rows[0].split('\t') as List def called = rows.drop(1).collect { r -> r.split('\t', -1) } .findAll { r -> r[head.indexOf('chromothripsis_confidence')] != 'none' } @@ -105,6 +107,9 @@ nextflow_process { { assert called == ['chr19:high:HC1,HC2'] }, { assert rows.size() == 24 }, // header + chromosomes 1-22 and X { assert log_text.contains('dropped 1 non-PASS VCF records') && log_text.contains('INS=1') }, + // the 100 bp deletion falls below --min-svlen 150; all five translocations are kept + { assert log_text.contains('dropped 1 junctions: intra-chromosomal span < 150 bp') }, + { assert sv_rows.count { r -> r.contains('\tTRA\t') } == 5 && !sv_rows.any { r -> r.startsWith('2\t30000000\t') } }, { assert process.out.plots.get(0).get(1).toString().endsWith('test_chr19.pdf') }, { assert snapshot(process.out.summary, process.out.inputs, process.out.versions).match() } ) @@ -167,6 +172,7 @@ nextflow_process { assertAll( { assert process.success }, { assert path(process.out.log.get(0).get(1)).text.contains('tiers: chr19 high (HC1,HC2)') }, + { assert path(process.out.log.get(0).get(1)).text.contains('dropped 1 junctions: intra-chromosomal span < 150 bp') }, // SAVANA writes BP_NOTATION in (_1, _2) order; the chr7 _2 records come first in the file, and the // translocations must keep the chr19 side's orientation (the same table as from Severus) { assert sv_rows.findAll { r -> r.contains('TRA') } == [ diff --git a/modules/local/shatterseek/tests/main.nf.test.snap b/modules/local/shatterseek/tests/main.nf.test.snap index a284d8e3..dcd8dbf4 100644 --- a/modules/local/shatterseek/tests/main.nf.test.snap +++ b/modules/local/shatterseek/tests/main.nf.test.snap @@ -46,7 +46,7 @@ ] ], "5": [ - "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + "versions.yml:md5,243442e8d24a16e7ef70c1d6f9157ca8" ], "inputs": [ [ @@ -92,7 +92,7 @@ ] ], "versions": [ - "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + "versions.yml:md5,243442e8d24a16e7ef70c1d6f9157ca8" ] } ], @@ -100,7 +100,7 @@ "nf-test": "0.9.3", "nextflow": "25.10.4" }, - "timestamp": "2026-10-14T15:40:32.880945189" + "timestamp": "2026-10-30T10:33:37.070484047" }, "ascat + severus - real": { "content": [ @@ -109,7 +109,7 @@ { "id": "test" }, - "test.shatterseek.tsv:md5,52bdbcbc33b8fceb7cce44a19ec47ae7" + "test.shatterseek.tsv:md5,c2f449257a71dfee9f946d9d61c3b838" ] ], [ @@ -124,14 +124,14 @@ ] ], [ - "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" ] ], "meta": { "nf-test": "0.9.3", "nextflow": "25.10.4" }, - "timestamp": "2026-10-14T15:40:48.576974826" + "timestamp": "2026-10-30T10:33:52.021265099" }, "wakhan + severus - real": { "content": [ @@ -140,7 +140,7 @@ { "id": "test" }, - "test.shatterseek.tsv:md5,4298d3136a0d42899397d807b4d7a516" + "test.shatterseek.tsv:md5,fa5e8ac3f6a5639b3b9dd44e46a13b66" ] ], [ @@ -155,14 +155,14 @@ ] ], [ - "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" ] ], "meta": { "nf-test": "0.9.3", "nextflow": "25.10.4" }, - "timestamp": "2026-10-14T15:40:58.358617435" + "timestamp": "2026-10-30T10:34:01.900701475" }, "savana single source - real": { "content": [ @@ -171,7 +171,7 @@ { "id": "test" }, - "test.shatterseek.tsv:md5,9dfb723da4741ead534812fa6e212799" + "test.shatterseek.tsv:md5,e35f12492ea69f91fa3bd2747e4bf57c" ] ], [ @@ -186,14 +186,14 @@ ] ], [ - "versions.yml:md5,0ddceb0c12aa3dbc8fa9648c158e1770" + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" ] ], "meta": { "nf-test": "0.9.3", "nextflow": "25.10.4" }, - "timestamp": "2026-10-14T15:41:08.198924416" + "timestamp": "2026-10-30T10:34:11.328170994" }, "savana single source - stub": { "content": [ @@ -242,7 +242,7 @@ ] ], "5": [ - "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + "versions.yml:md5,243442e8d24a16e7ef70c1d6f9157ca8" ], "inputs": [ [ @@ -288,7 +288,7 @@ ] ], "versions": [ - "versions.yml:md5,8cbcc93057dfcfad0d7cd89f1b035179" + "versions.yml:md5,243442e8d24a16e7ef70c1d6f9157ca8" ] } ], @@ -296,6 +296,6 @@ "nf-test": "0.9.3", "nextflow": "25.10.4" }, - "timestamp": "2026-10-14T15:40:38.682355173" + "timestamp": "2026-10-30T10:33:42.290507479" } } \ No newline at end of file diff --git a/modules/local/shatterseek/tests/nextflow.config b/modules/local/shatterseek/tests/nextflow.config index 8c5ffc88..ca057aa6 100644 --- a/modules/local/shatterseek/tests/nextflow.config +++ b/modules/local/shatterseek/tests/nextflow.config @@ -1,5 +1,7 @@ process { withName: 'SHATTERSEEK' { ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + // as in conf/modules.config, so the real tests run the pipeline's default length filter (150 bp) + ext.args = { params.shatterseek_min_svlen ? "--min-svlen ${params.shatterseek_min_svlen}" : '' } } } diff --git a/subworkflows/local/shatterseek_chromothripsis.nf b/subworkflows/local/shatterseek_chromothripsis.nf index 6dcd5727..556f04f0 100644 --- a/subworkflows/local/shatterseek_chromothripsis.nf +++ b/subworkflows/local/shatterseek_chromothripsis.nf @@ -22,13 +22,12 @@ workflow SHATTERSEEK_CHROMOTHRIPSIS { ch_versions = channel.empty() // A digest of every wrapper file goes into the meta map, so editing either wrapper invalidates -resume - def wrapper_src = { name -> - def dir = file("${projectDir}/assets/${name}", type: 'dir', checkIfExists: true) - def digest = files("${projectDir}/assets/${name}/**", type: 'file').sort { f -> f.toString() }.collect { f -> f.text }.join('\n').md5() - return channel.value([[id: name, digest: digest], dir]) - } - shatterseek_src = wrapper_src('shatterseek') - reconplot_src = wrapper_src('reconplot') + def shatterseek_dir = file("${projectDir}/assets/shatterseek", type: 'dir', checkIfExists: true) + def shatterseek_digest = files("${projectDir}/assets/shatterseek/**", type: 'file').sort { f -> f.toString() }.collect { f -> f.text }.join('\n').md5() + def reconplot_dir = file("${projectDir}/assets/reconplot", type: 'dir', checkIfExists: true) + def reconplot_digest = files("${projectDir}/assets/reconplot/**", type: 'file').sort { f -> f.toString() }.collect { f -> f.text }.join('\n').md5() + shatterseek_src = channel.value([[id: 'shatterseek', digest: shatterseek_digest], shatterseek_dir]) + reconplot_src = channel.value([[id: 'reconplot', digest: reconplot_digest], reconplot_dir]) severus_sv_files = severus_vcf.map { meta, vcf -> [meta, [vcf]] } // severus_sv_files: [meta, [severus_somatic.vcf.gz]] diff --git a/tests/fixtures/shatterseek/make_fixtures.py b/tests/fixtures/shatterseek/make_fixtures.py index c9ff4654..dfed6d69 100644 --- a/tests/fixtures/shatterseek/make_fixtures.py +++ b/tests/fixtures/shatterseek/make_fixtures.py @@ -7,6 +7,7 @@ chr19 <-> chr7 5 translocations from inside the cluster to chr7:50.1-50.9 Mb chr19:20-30 Mb copy number alternating 1/2 over ten 1 Mb segments (2 elsewhere) chr2:10.00-10.05 Mb one isolated deletion (a chromosome with an SV but no event) + chr2:30.0000-30.0001 Mb a 100 bp deletion, dropped by the pipeline's default --min-svlen 150 Expected ShatterSeek result: chr19 high confidence (HC1, HC2), chr2 and chr7 none. Also present so the conversion is exercised: an insertion and a non-PASS deletion (both dropped), a Wakhan segment with a fractional haplotype state and fractional SAVANA copy numbers (rounded), and SAVANA translocations whose @@ -75,6 +76,7 @@ def sev(chrom, pos, vid, alt_, info, flt="PASS"): sev("chr19", p19, f"severus_BND{100 + j}_1", alt(s, "chr7", p7), f"SVTYPE=BND;SVLEN=0;MATE_ID=severus_BND{100 + j}_2;STRANDS={s}") sev("chr7", p7, f"severus_BND{100 + j}_2", alt(s[::-1], "chr19", p19), f"SVTYPE=BND;SVLEN=0;MATE_ID=severus_BND{100 + j}_1;STRANDS={s[::-1]}") sev("chr2", 10_000_000, "severus_DEL200", "", "SVTYPE=DEL;SVLEN=50000;END=10050000;STRANDS=+-") +sev("chr2", 30_000_000, "severus_DEL202", "", "SVTYPE=DEL;SVLEN=100;END=30000100;STRANDS=+-") sev("chr2", 20_000_000, "severus_DEL201", "", "SVTYPE=DEL;SVLEN=80000;END=20080000;STRANDS=+-", flt="FAIL_LOWSUPPORT") sev("chr19", 25_000_100, "severus_INS1", "", "SVTYPE=INS;SVLEN=300") rec.sort(key=lambda r: (ORDER[r[0]], r[1])) @@ -117,6 +119,7 @@ def sav_pair(n, c1, p1, c2, p2, note, svlen): for j, (p19, p7, s) in enumerate(tra, 1): # _1 on chr19, so the chr7 _2 record is written first sav_pair(100 + j, "chr19", p19, "chr7", p7, s, 0) sav_pair(200, "chr2", 10_000_000, "chr2", 10_050_000, "+-", 50000) +sav_pair(202, "chr2", 30_000_000, "chr2", 30_000_100, "+-", 100) rec.append(("chr19", 25_000_100, "chr19\t25000100\tID_300_1\tN\t\t.\tPASS\tSVTYPE=INS;TUMOUR_READ_SUPPORT=12;SVLEN=300;BP_NOTATION=\tGT\t0/1")) rec.sort(key=lambda r: (ORDER[r[0]], r[1])) w("savana/test_tumor.classified.somatic.vcf", sav_hdr + [r[2] for r in rec]) diff --git a/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf b/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf index 5dbe3820..e5b72960 100644 --- a/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf +++ b/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf @@ -13,6 +13,8 @@ #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT test_tumor chr2 10000000 ID_200_1 N N[chr2:10050000[ . PASS MATEID=ID_200_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=50000;BP_NOTATION=+- GT 0/1 chr2 10050000 ID_200_2 N ]chr2:10000000]N . PASS MATEID=ID_200_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=50000;BP_NOTATION=+- GT 0/1 +chr2 30000000 ID_202_1 N N[chr2:30000100[ . PASS MATEID=ID_202_2;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=100;BP_NOTATION=+- GT 0/1 +chr2 30000100 ID_202_2 N ]chr2:30000000]N . PASS MATEID=ID_202_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=100;BP_NOTATION=+- GT 0/1 chr7 50100000 ID_101_2 N ]chr19:20500000]N . PASS MATEID=ID_101_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=+- GT 0/1 chr7 50300000 ID_102_2 N N[chr19:22500000[ . PASS MATEID=ID_102_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=-+ GT 0/1 chr7 50500000 ID_103_2 N N]chr19:24500000] . PASS MATEID=ID_103_1;SVTYPE=BND;TUMOUR_READ_SUPPORT=12;SVLEN=0;BP_NOTATION=++ GT 0/1 diff --git a/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz b/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz index 1d8532e0cc64935e25f1b09a44b956e70160ef5b..5492356ec18734471cf7e469d5e9fcdb2c6e1532 100644 GIT binary patch literal 1040 zcmV+r1n>JFiwFpiV&rN9|8r$_WpZ_MUvqD5VRUI@E_P#P0IXMCZ<|UGeWt$xHIGeF zW?}hAN+mbo6p3K2!LF{N%5rUzSl0myaHRhGyRcxf_S%qOzr=*ynKN_FeE8zxX|;}? zlJs?+o721Dn{tglU0h_zHeEz!5WPj|Hk-Y{wMzawbllJ$o8KJE^anw3X+_y0U2XDJ z5}Tg~>$OO1DxQv+3c~Tgvw{vW{z#+wf16|#=S7(Fco{9@D9c!L4};L2IaUXA#YxezM8c!^ zr!XDGOKNO58k3i$x%;q#rTZvOa`JNjvNi3JBzmMsPsNx>z*I;UNt~~qOJj?d6kV=M zR!pQIFl3C3!!ufa2CyM3*wEy|Ghoz*u&#opf=bTZ{Tz=T1b-9=jwOukp)hbm;h)LC z6{v!RBEhz25;VKdfE>vP<^Yi?kUAwwm69Gn6W2MA^~qMtz{__RjAAjJs>ymJ7nH`eFW6N^Bed z=urF`%fBXfZ-pmvuC?=zFn2;g)rN{Xl>2uf)Fl515Q8P3+|wB$YE&J)g{Zd+QI_{a zVY|kMg{T87iX%dl_XZP@7)T#hVl_J2lCHo((g9d0>XethQUqw%qQYMic5n(Yo?}Po z0OX)S17gHzP=)jsBcOIQV!Ctg&^TUD&S&!XgfPh4$~~!HzFmzohM0}48dUwq5VU1n z1Q)(wtTCYyCT-~7b?`wu_3k+ugRP_!4QP-ndU?*L7K*C&rO1={(Xk_ z+)9v+xvMj?JQ}^1N8G+V>IScUvOj9zT69547z4_w+ literal 1027 zcmV+e1pNCSiwFquOV4To|8r$_WpZ_MUvqD5VRUI@E_P#P0IXMAZ<|UGe%}5H)I2uH zVHTE4QYyItr$`iY4E9kKRhDaS5~~hafFt$apXC;@_S(3?eu)XYGvCbj&BYfNkE=C( zOwyNmZcXoouj)1aaB-0(+jK##kiOD%o6TO~T4Vnmx?bdrt04asA<_yb(^LAAxh}Z(n}V)Pl?C&xSe3(^dS$8!BE% z$)l9z>2{HC)A?F@ou{jLobMX+9Dm228|U}6^0K03kDJ&lu`z1+QQPftlwb)-m7&6vc%)JPUdoUfirV~gh$U#?45 zEUaQM6@tzEGkS9dh^eT=)Rp}+U^a({p@FV}#?IV*8;`!DU=$+PMq_7)23~~zoeVrw zDHs(gMpbkur&&JBw=YThp}>Veu0ANM#cw+vW{*Dl%)Wg=ntUyzLjEWiJ9lpA+*ZoA z*xTDulNKrN-!frY%vkkMOT+qfh3=ypdy3?L;0AEz5(`}d?{YYBy_q-qFPsFyXdE5( zZoGMCc8xCYIwCUg-oN@b$S}jq|IaeTH@J#I{^=hMpWUhf$hWm}tsRT-2IT6K5QriG zf=eL&LyZJiNFI6s!7U&MGC~AElnO|b617fA4IhBVe*xfxTLRzbQ_&EQz{6)GDlFawvTgVTC{=O#P-QqK94uf=>XV zbHPsnTn0Z0js-snM8OlKG}f7L6V#pHl`8axd7_smPAY{4fgl4)JAN@1*zXW(iRG;t z;4iA-2o|-4hj%U7Ghz9VYI3Pr2 zXD~5|f&6YI)}x~<=_>3c9e|gjQSs7OiU8eORK-ie8=OLn=hzWW0J7Jh0drzBD~0qH zBj9c|79m3mjq8Wyd}e=-2!p+?-Lv`?+tmeQh}+2OLDhc@L085#5R9SrQdZBYGmIhV z%C|;@D`G=`p+fx}8)7H%x+0hp96^|zX-;$!r0atA+$N~*-)Csgodg*~xH@yoquG0T xB;CuSVT#%(|Dy)Z)jm5ZG8IwRYI9WgwR3c>!zBF=009600{}-$QgHzj002bK>LdUF From b25bf2b74c08e4f98ab19e975693f758074752f2 Mon Sep 17 00:00:00 2001 From: ljwharbers Date: Tue, 3 Nov 2026 16:48:22 +0100 Subject: [PATCH 5/5] Move SHATTERSEEK_CHROMOTHRIPSIS to subworkflows/local//main.nf (follows #227) Co-Authored-By: Claude Opus 5.5 --- .../main.nf} | 6 +++--- subworkflows/local/tests/shatterseek_chromothripsis.nf.test | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) rename subworkflows/local/{shatterseek_chromothripsis.nf => shatterseek_chromothripsis/main.nf} (97%) diff --git a/subworkflows/local/shatterseek_chromothripsis.nf b/subworkflows/local/shatterseek_chromothripsis/main.nf similarity index 97% rename from subworkflows/local/shatterseek_chromothripsis.nf rename to subworkflows/local/shatterseek_chromothripsis/main.nf index 556f04f0..933aa1fb 100644 --- a/subworkflows/local/shatterseek_chromothripsis.nf +++ b/subworkflows/local/shatterseek_chromothripsis/main.nf @@ -1,7 +1,7 @@ // IMPORT MODULES -include { SHATTERSEEK as SHATTERSEEK_SEVERUS_ASCAT } from '../../modules/local/shatterseek/main' -include { SHATTERSEEK as SHATTERSEEK_SEVERUS_WAKHAN } from '../../modules/local/shatterseek/main' -include { SHATTERSEEK as SHATTERSEEK_SAVANA } from '../../modules/local/shatterseek/main' +include { SHATTERSEEK as SHATTERSEEK_SEVERUS_ASCAT } from '../../../modules/local/shatterseek/main' +include { SHATTERSEEK as SHATTERSEEK_SEVERUS_WAKHAN } from '../../../modules/local/shatterseek/main' +include { SHATTERSEEK as SHATTERSEEK_SAVANA } from '../../../modules/local/shatterseek/main' // // ShatterSeek chromothripsis calls with PCAWG confidence tiers for every CN/SV caller pair that produced output for a diff --git a/subworkflows/local/tests/shatterseek_chromothripsis.nf.test b/subworkflows/local/tests/shatterseek_chromothripsis.nf.test index d9a00fa3..79fb64ab 100644 --- a/subworkflows/local/tests/shatterseek_chromothripsis.nf.test +++ b/subworkflows/local/tests/shatterseek_chromothripsis.nf.test @@ -1,7 +1,7 @@ nextflow_workflow { name "Test Workflow SHATTERSEEK_CHROMOTHRIPSIS" - script "../shatterseek_chromothripsis.nf" + script "../shatterseek_chromothripsis/main.nf" workflow "SHATTERSEEK_CHROMOTHRIPSIS" tag "subworkflows"