Repository navigation
Add ShatterSeek chromothripsis calling with PCAWG confidence tiers - #222
Open
robert-a-forsyth wants to merge 5 commits into
Open
robert-a-forsyth wants to merge 5 commits into
robert-a-forsyth wants to merge 5 commits into
Conversation
SAVANA writes BP_NOTATION in (_1, _2) order on both mate records (all 12,679 _2 records of B2101426 are reversed relative to their own ALT). Mates were collapsed keeping whichever record came first in the file, so a junction whose _2 record came first paired the second breakend's position with the first breakend's orientation. ReConPlot's tables for the fixtures are unchanged; the snapshot moves only with the wrapper VERSION in versions.yml. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
SHATTERSEEK_CHROMOTHRIPSIS runs ShatterSeek for every CN/SV caller pair with output for a sample (ASCAT + Severus, Wakhan + Severus, SAVANA), reading the callers with the ReConPlot wrapper's parsers, and labels each chromosome's cluster with the PCAWG HC1/HC2/HC3/LC1 tiers (BH q at 0.2). The module uses the robert-a-forsyth/ShatterSeek fork (chm13 support, statistics corrected to the PCAWG definitions, one cluster per chromosome) in docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c. New parameters: --skip_shatterseek, --shatterseek_genome (inferred from --genome) and --shatterseek_min_svlen (default 150 bp). Synthetic fixtures in tests/fixtures/shatterseek cover every caller layout; module, subworkflow and pipeline-wiring tests are tagged small. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
| // No conda: the ShatterSeek fork (chm13, corrected statistics, PCAWG tiers) ships only in this image, built from | ||
| // containers/shatterseek/Dockerfile; the wrapper is pipeline glue in assets/shatterseek/ and reads the callers with | ||
| // the ReConPlot parsers in assets/reconplot/. Fork update = new commit there, rebuild, bump this tag. | ||
| container "docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c" |
There was a problem hiding this comment.
Container convention (REVIEW.md): this is a single Docker reference with no singularity/apptainer ternary. Custom images here are meant to use the oras://docker.io/...-sif / docker.io/... pair, as padfoot, reconplot, wakhan and lrsomaticreport do. As it stands, every Singularity/Apptainer run converts the image at pull time. nf-core pipelines download --container-system singularity won't stage a SIF for it either, even though the "Offline and air-gapped systems" section of docs/usage.md says the modules use the pinned-tag pair. AmpliconClassifier is the only other module that does this, so it isn't a reason to skip it here. Please push a robertaforsyth/shatterseek-sif:1.1.1-chm13-e1cee6c twin and use the ternary.
| for (ch in unique(d$chrom)) { | ||
| x <- d[d$chrom == ch, , drop = FALSE] | ||
| # no overlaps: each segment starts after the previous one ends | ||
| for (i in seq_len(nrow(x))[-1]) if (x$start[i] <= x$end[i - 1]) x$start[i] <- x$end[i - 1] + 1L |
There was a problem hiding this comment.
Bug: overlap trimming misses nested segments. Each segment is compared only with the previous row's end, and a row that gets trimmed away (start > end) still counts as the previous row. Example: A = 1–1000, B = 500–800, C = 900–2000. B becomes 1001–800 and is dropped, but C is compared with B's end (800), so it stays 900–2000 and still overlaps A. ShatterSeek then gets overlapping CN segments, which changes the oscillation count. Comparing against the running maximum end fixes it:
Suggested change
| for (i in seq_len(nrow(x))[-1]) if (x$start[i] <= x$end[i - 1]) x$start[i] <- x$end[i - 1] + 1L | |
| for (i in seq_len(nrow(x))[-1]) { prev_end <- max(x$end[seq_len(i - 1)]); if (x$start[i] <= prev_end) x$start[i] <- prev_end + 1L } |
| @@ -0,0 +1,5 @@ | |||
| process { | |||
| withName: 'SHATTERSEEK' { | |||
| ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } | |||
There was a problem hiding this comment.
Test gap: the pipeline's default config is never run. The module tests set only ext.prefix, so every "real" test runs the wrapper with --min-svlen 0. The pipeline instead runs --min-svlen 150 by default (conf/modules.config). That makes the length filter untested, along with the claim that translocations are always kept. Please add ext.args = { params.shatterseek_min_svlen ? "--min-svlen ${params.shatterseek_min_svlen}" : '' } here, or set it in one real test. That test should include an intra-chromosomal SV shorter than 150 bp in the fixture and assert the log line that reports it was dropped.
|
Automated review: 3 findings (see the inline comments).
For a human reviewer to check by running things (I couldn't run any tests):
|
|
Replace the robert-a-forsyth ShatterSeek fork with upstream parklab/ShatterSeek (3b3f4cd, v1.1) plus a build-time patch (containers/shatterseek/patch_chm13.R) that adds T2T-CHM13 cytobands and chromosome sizes and accepts "chm13" in the three genome checks; the statistics are upstream's. The wrapper now computes the PCAWG tiers itself (HC1, HC2, HC3; issue #22) with BH q-values per sample, drops the low-confidence tier (LC1, --min-osc-lc/--max-osc-lc), and seeds R's RNG because upstream's exponential tests draw random numbers. Review fixes: the module uses the singularity/docker ternary with an oras://...-sif twin; CN overlap trimming uses the running maximum end so nested segments no longer slip through; the module tests run the pipeline's default --min-svlen 150 against a new 100 bp fixture deletion. CI: inline the wrapper digests in SHATTERSEEK_CHROMOTHRIPSIS (nextflow-lint rejected the local closure) and exclude the ShatterSeek Wakhan fixture's trailing tab from trailing-whitespace, as for the Padfoot one. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds ShatterSeek chromothripsis calling. It runs for every CN/SV caller pair that has output for a sample: ASCAT + Severus, Wakhan + Severus, and SAVANA. Each chromosome's cluster gets the PCAWG high/low-confidence tiers (HC1, HC2, HC3, LC1; Cortes-Ciriano et al. 2020, Supplementary Note 1.2, parklab/ShatterSeek#22). Results go to
<sample>/shatterseek/<pair>/.SHATTERSEEK_CHROMOTHRIPSIS: builds the same pairs asRECONPLOT_FIGURESand reads the callers with the same parsers (assets/reconplot/R/parsers), so ShatterSeek and ReConPlot see the same SVs and segments.assets/shatterseek/run_shatterseek.R: converts the caller output to ShatterSeek's input.--skip_shatterseek--shatterseek_genome(inferred from--genome: GRCh38 → hg38, CHM13 → chm13)--shatterseek_min_svlen(default 150 bp; translocations are always kept)docker.io/robertaforsyth/shatterseek:1.1.1-chm13-e1cee6c, recipe incontainers/shatterseek/. Like CoRAL, Singularity pulls the Docker image, and the module refuses-profile conda.The ShatterSeek fork
Upstream ShatterSeek computes several of its statistical criteria differently from how the PCAWG paper defines them. The module uses robert-a-forsyth/ShatterSeek@
chm13-support, which fixes them in place:pos1forpos2), and the narrowed set then replaced the translocation list for every later test. DO17373 chrXnumber_TRAgoes from 1 → 600, the published value.SVs();number_CN_segments_chris no longer 2 short; the linked-windowseq()crash (upstream add cramino to the main lr_somatic.nf workflow #19/Make sure that ASCAT allele/loci files are in correct format when downloaded #32) is guarded; and T2T-CHM13 support is added.On the bundled DO17373 example, the fork reproduces PCAWG's published counts and labels: chr2, chr21 and chrX high confidence; chr3 only linked to the chr2 event. Upstream got chr3 wrong (HC2, from pooled clusters), and reached chrX only through HC3. The fork's
tests/test_do17373.Rchecks this at every image build.Also fixed
_2breakends from the_1side before collapsing mates. SAVANA writesBP_NOTATIONin (_1, _2) order on both mate records, so a junction whose_2record came first in the file got the second breakend's position with the first breakend's orientation. ReConPlot's tables for the fixtures are unchanged; its snapshot moves only through the wrapper VERSION inversions.yml.Things reviewers should know
ext.args = '--no-hc3'drops it.Tests
tests/fixtures/shatterseek/make_fixtures.py) in each caller's own layout. All pairs must call the chr19 eventhigh (HC1,HC2)and nothing else. The SAVANA fixture puts_2mates first.tests/sv_annotation_compose.nf.testnow asserts theSHATTERSEEK_*processes are present when ShatterSeek is on and absent with--skip_shatterseek.🤖 Generated with Claude Code