Add ShatterSeek chromothripsis calling with PCAWG confidence tiers

.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: |

CHANGELOG.md

@@ -7,6 +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 `<sample>/shatterseek/<pair>/` 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).
 - [#204](https://github.com/IntGenomicsLab/lrsomatic/pull/204) - Added haplotype differential methylation region (DMR) calling: with `--modkit_phased` set, `modkit dmr pair` compares methylation between the `_hp1`/`_hp2` bedMethyl pileups over CpG islands covered in both haplotypes, and each DMR call is annotated with its nearest gene via `bedtools closest`. Runs automatically whenever `--modkit_phased` is set; turn it off with `--skip_dmr`. `--dmr_cpg_islands_bed`/`--dmr_gencode_gene_bed` have built-in defaults for `--genome GRCh38`/`CHM13` and must be supplied explicitly for any other genome, or DMR is skipped with a warning rather than failing the run. Tumour/normal (cross-sample) DMR is explicitly out of scope for this PR (@YannVRB).
 - [#204](https://github.com/IntGenomicsLab/lrsomatic/pull/204) - Added a real, public subworkflow-level nf-test for the `DMR` subworkflow (tag `small`), chaining `MODKIT_PILEUP` on the public phased, methylation-tagged `test.sorted.phased.bam` fixture through the full DMR chain against a small synthetic CpG-island/gene region (@YannVRB).
 - [#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).
@@ -61,6 +62,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Fixed`
 
+- [#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).
 - [#223](https://github.com/IntGenomicsLab/lrsomatic/pull/223) - Fixed the Severus input channel (@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).

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/

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))
+
 <!-- TODO nf-core: Include a figure that guides the user through the major workflow steps. Many nf-core
      workflows use the "tube map" design for that. See https://nf-co.re/docs/guidelines/graphic_design/workflow_diagrams#examples for examples.   -->
 
@@ -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 (`<sample>/report/<sample>_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

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)
 }

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).

assets/reconplot/VERSION

@@ -1 +1 @@
-2672ccb+lrsomatic.2026-10-07
+2672ccb+lrsomatic.2026-10-14

assets/shatterseek/README.md

@@ -0,0 +1,24 @@
+# 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 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:
+
+- 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 (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 (`<sample>.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` and `--no-hc3`; see `Rscript run_shatterseek.R --help`.

assets/shatterseek/VERSION

@@ -0,0 +1 @@
+1.1.0