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 250d345f..f19cbc44 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` +- [#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). - [#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). 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/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/assets/shatterseek/README.md b/assets/shatterseek/README.md new file mode 100644 index 00000000..a40141f0 --- /dev/null +++ b/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 (`.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`. diff --git a/assets/shatterseek/VERSION b/assets/shatterseek/VERSION new file mode 100644 index 00000000..9084fa2f --- /dev/null +++ b/assets/shatterseek/VERSION @@ -0,0 +1 @@ +1.1.0 diff --git a/assets/shatterseek/run_shatterseek.R b/assets/shatterseek/run_shatterseek.R new file mode 100644 index 00000000..f12d5c0a --- /dev/null +++ b/assets/shatterseek/run_shatterseek.R @@ -0,0 +1,211 @@ +#!/usr/bin/env Rscript +## --------------------------------------------------------------------------- +## run_shatterseek.R -- ShatterSeek chromothripsis calls from lrsomatic caller output. +## +## --[ReConPlot wrapper parsers]--> {cn, sv} --[ShatterSeek format]--> shatterseek() +## --> 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 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({ + 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("--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("--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 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), + 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 --------------------------------------------- +## 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_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, 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 41fada0d..aaa1002e 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -1060,6 +1060,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..c330d3b9 --- /dev/null +++ b/containers/shatterseek/Dockerfile @@ -0,0 +1,37 @@ +# syntax=docker/dockerfile:1.6 +FROM docker.io/mambaorg/micromamba:2.3.2 + +# 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=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 \ + # 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 new file mode 100644 index 00000000..8b8cdb02 --- /dev/null +++ b/containers/shatterseek/README.md @@ -0,0 +1,31 @@ +# ShatterSeek image + +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: 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). 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-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" +``` + +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 new file mode 100644 index 00000000..d282f70f --- /dev/null +++ b/containers/shatterseek/environment.yml @@ -0,0 +1,21 @@ +--- +# 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 (upstream parklab/ShatterSeek, 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-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" + - "conda-forge::r-optparse=1.8.2" + - "conda-forge::procps-ng" # `ps` required by Nextflow task metrics 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 bdf69b86..70d8c0ad 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` @@ -790,6 +798,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 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. +- 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 044e9e02..1e856a14 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` | @@ -470,6 +471,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-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`, `--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: @@ -478,9 +496,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/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", diff --git a/modules/local/shatterseek/main.nf b/modules/local/shatterseek/main.nf new file mode 100644 index 00000000..10620ef4 --- /dev/null +++ b/modules/local/shatterseek/main.nf @@ -0,0 +1,73 @@ +process SHATTERSEEK { + tag "${meta.id}:${cn_source}+${sv_source}" + label 'process_low' + + // 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/. 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. + // 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: 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}") + 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_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_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..e06c16a6 --- /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 (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 + - 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/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: "" + +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) 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..c7e7cb68 --- /dev/null +++ b/modules/local/shatterseek/tests/main.nf.test @@ -0,0 +1,214 @@ +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 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") { + + 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 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' } + .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') }, + // 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() } + ) + } + } + + 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)') }, + { 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') } == [ + '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..dcd8dbf4 --- /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,243442e8d24a16e7ef70c1d6f9157ca8" + ], + "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,243442e8d24a16e7ef70c1d6f9157ca8" + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-30T10:33:37.070484047" + }, + "ascat + severus - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,c2f449257a71dfee9f946d9d61c3b838" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,982eb23095d2181864328d8ad5ff771f", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-30T10:33:52.021265099" + }, + "wakhan + severus - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,fa5e8ac3f6a5639b3b9dd44e46a13b66" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,e5fc9a40fddb5b1e0371ec6ed19356d7", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-30T10:34:01.900701475" + }, + "savana single source - real": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.shatterseek.tsv:md5,e35f12492ea69f91fa3bd2747e4bf57c" + ] + ], + [ + [ + { + "id": "test" + }, + [ + "test.shatterseek_cn_input.tsv:md5,982eb23095d2181864328d8ad5ff771f", + "test.shatterseek_sv_input.tsv:md5,b668729863186626947977e9651505ae" + ] + ] + ], + [ + "versions.yml:md5,e956f799e9d3ed36eb265ced3ab3de77" + ] + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-10-30T10:34:11.328170994" + }, + "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,243442e8d24a16e7ef70c1d6f9157ca8" + ], + "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,243442e8d24a16e7ef70c1d6f9157ca8" + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "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 new file mode 100644 index 00000000..ca057aa6 --- /dev/null +++ b/modules/local/shatterseek/tests/nextflow.config @@ -0,0 +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/nextflow.config b/nextflow.config index 836884c5..4b012cbe 100644 --- a/nextflow.config +++ b/nextflow.config @@ -108,6 +108,7 @@ params { skip_savana = false skip_padfoot = false skip_reconplot = false + skip_shatterseek = false skip_wakhan = false skip_fiber = false skip_normalfiber = false @@ -158,6 +159,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 46c7fe25..e47f6876 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -515,6 +515,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", @@ -701,6 +722,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" @@ -1015,6 +1041,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..556f04f0 --- /dev/null +++ b/subworkflows/local/shatterseek_chromothripsis.nf @@ -0,0 +1,84 @@ +// 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 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]] + + // + // 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..dfed6d69 --- /dev/null +++ b/tests/fixtures/shatterseek/make_fixtures.py @@ -0,0 +1,130 @@ +#!/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) + 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 +_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", 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])) +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) +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]) +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..e5b72960 --- /dev/null +++ b/tests/fixtures/shatterseek/savana/test_tumor.classified.somatic.vcf @@ -0,0 +1,52 @@ +##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 +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 +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 00000000..5492356e Binary files /dev/null and b/tests/fixtures/shatterseek/severus/severus_somatic.vcf.gz differ 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 e6f56c9e..70690079 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' @@ -64,6 +65,7 @@ include { PAIRED_SAVANA } from '../subworkflows/local/paired/p include { DMR } from '../subworkflows/local/dmr' include { PADFOOT_ANNOTATION } from '../subworkflows/local/padfoot_annotation' include { RECONPLOT_FIGURES } from '../subworkflows/local/reconplot_figures' +include { SHATTERSEEK_CHROMOTHRIPSIS } from '../subworkflows/local/shatterseek_chromothripsis' @@ -1290,6 +1292,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 //