diff --git a/CHANGELOG.md b/CHANGELOG.md index 250d345f..8251c268 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -61,6 +61,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### `Fixed` +- [#225](https://github.com/IntGenomicsLab/lrsomatic/issues/225) - PacBio tumour-only SAVANA calls are now filtered against the genome's Severus germline SV panel of normals before VEP, Padfoot and ReConPlot. SAVANA's PacBio classifier only rejects a call on matched-normal evidence, so without a normal its somatic VCF held nearly every germline SV (27,591 events on a simulated tumour-only C8, against 2 records with the matched normal) and `RECONPLOT_SAVANA` took 6-8 h; the filter, Severus's own PoN matching rule in `bin/savana_pon_filter.py` / `SAVANA_PON_FILTER`, removes 97% of them (761 kept), and `RECONPLOT_SAVANA` then finishes in 21 min with every figure drawn (unfiltered: not finished after 2 h). The 761 left are still germline by C8's matched normal, so PacBio tumour-only SAVANA calls remain a fallback. The unfiltered VCF is still published; `--savana_pb_pon_filter false` turns it off (@ljwharbers). +- [#225](https://github.com/IntGenomicsLab/lrsomatic/issues/225) - ReConPlot figures draw at most `--reconplot_max_svs` SVs (default 1000, best supported first; the title then says `top N of M SVs`), since ReConPlot adds one plot layer per SV: tens of thousands take hours, and the genome-wide panel already overflows R's node stack at ~2,000 (no genome-wide figure was drawn for BL1, 2,084 SAVANA junctions, or for 2,500 capped C8 junctions; 761 drew). The harmonised `*.reconplot_sv.tsv` table keeps every SV (@ljwharbers). - [#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/assets/reconplot/R/plotting.R b/assets/reconplot/R/plotting.R index 8373be41..25dc581c 100644 --- a/assets/reconplot/R/plotting.R +++ b/assets/reconplot/R/plotting.R @@ -34,6 +34,7 @@ build_title <- function(user_title, meta, regions) { paste(chrs, collapse = ", ") } bits <- c(bits, where) + if (!is.null(meta$sv_note)) bits <- c(bits, meta$sv_note) paste(bits, collapse = " | ") } diff --git a/assets/reconplot/R/utils_common.R b/assets/reconplot/R/utils_common.R index 01ffe0e0..8cfaa481 100644 --- a/assets/reconplot/R/utils_common.R +++ b/assets/reconplot/R/utils_common.R @@ -110,6 +110,26 @@ validate_cn <- function(cn, drop_na_minor = FALSE) { cn } +#' Keep the `max_svs` best-supported junctions for drawing. +#' +#' ReConPlot adds one ggplot layer per SV, so a few thousand junctions take minutes and tens of thousands +#' take hours (e.g. SAVANA's PacBio tumour-only calls without a panel of normals). Junctions are ranked by +#' tumour read support, then by span (translocations first, insertions by their length); a caller without +#' a support column is ranked on span alone. +cap_svs <- function(sv, max_svs) { + if (is.null(max_svs) || is.na(max_svs) || max_svs <= 0 || nrow(sv) <= max_svs) return(sv) + support <- if (!is.null(sv$support)) suppressWarnings(as.numeric(sv$support)) else rep(NA_real_, nrow(sv)) + svlen <- if (!is.null(sv$svlen)) suppressWarnings(abs(as.numeric(sv$svlen))) else rep(NA_real_, nrow(sv)) + span <- ifelse(sv$chr1 != sv$chr2, Inf, + ifelse(sv$strands %in% "INS" & !is.na(svlen), svlen, abs(sv$pos2 - sv$pos1))) + keep <- sort(order(-support, -span, na.last = TRUE)[seq_len(max_svs)]) + log_msg(sprintf(" WARNING: %d SV junctions exceed --max-svs %d; drawing the %d best supported", + nrow(sv), max_svs, max_svs)) + out <- sv[keep, , drop = FALSE] + rownames(out) <- NULL + out +} + validate_sv <- function(sv, interchrom_as_tra = TRUE) { req <- c("chr1", "pos1", "chr2", "pos2", "strands") missing <- setdiff(req, names(sv)) diff --git a/assets/reconplot/README.md b/assets/reconplot/README.md index 91a294ae..e2faba1c 100644 --- a/assets/reconplot/README.md +++ b/assets/reconplot/README.md @@ -14,4 +14,6 @@ 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); `--max-svs` draws +only the best-supported junctions (ReConPlot adds a plot layer per SV, so tens of thousands take hours) while the +`--write-tables` table keeps them all. diff --git a/assets/reconplot/VERSION b/assets/reconplot/VERSION index d3d30822..4fee3b5d 100644 --- a/assets/reconplot/VERSION +++ b/assets/reconplot/VERSION @@ -1 +1 @@ -2672ccb+lrsomatic.2026-10-07 +2672ccb+lrsomatic.2026-10-29 diff --git a/assets/reconplot/run_reconplot.R b/assets/reconplot/run_reconplot.R index ad8a94bf..26a63679 100755 --- a/assets/reconplot/run_reconplot.R +++ b/assets/reconplot/run_reconplot.R @@ -67,6 +67,9 @@ option_list <- list( make_option("--min-svlen", type = "double", default = 0, help = paste("Drop intra-chromosomal SVs shorter than this many bp;", "translocations are never dropped [default %default]")), + make_option("--max-svs", type = "integer", default = 0, + help = paste("Draw at most this many SV junctions, the best supported first; the", + "--write-tables table keeps them all. 0 = no limit [default %default]")), make_option("--exclude-vntr", action = "store_true", default = FALSE, help = "Drop SVs flagged as inside a VNTR (Severus)"), make_option("--cluster-id", type = "character", default = NULL, @@ -170,6 +173,11 @@ if (isTRUE(args$write_tables)) { log_msg(" wrote harmonised tables to ", args$outdir) } +## after the tables, which keep every junction +n_sv <- nrow(sv) +sv <- cap_svs(sv, args$max_svs) +if (nrow(sv) < n_sv) meta$sv_note <- sprintf("top %d of %d SVs", nrow(sv), n_sv) + ## --- 3. regions ------------------------------------------------------------- regions <- if (!is.null(args$regions_file)) { read_regions_file(args$regions_file) diff --git a/bin/savana_pon_filter.py b/bin/savana_pon_filter.py new file mode 100755 index 00000000..9ae354ea --- /dev/null +++ b/bin/savana_pon_filter.py @@ -0,0 +1,163 @@ +#!/usr/bin/env python3 +"""Drop SAVANA SVs that match a germline SV panel of normals (PoN). + +SAVANA's PacBio classifier (`savana to --pb`) has no tumour-only mode: it rejects a call +only on normal-read evidence, so without a normal every germline SV with enough tumour +reads is written to the somatic VCF (~20-30k per genome). This filters that VCF against +the Severus PoN the pipeline already uses, with Severus's own matching rule (extract_pon +and add_pon in severus/breakpoint_finder.py), so both callers call germline the same way. + +PoN rows: chr1,pos1,chr2,pos2,ci1,ci2,svtype,fraction -- for INS, pos2 is the length. +A header row (the CHM13 PoN has one) is skipped because its pos1 is not a number. + +An event is matched on its lower breakpoint, and both BND mates are kept or dropped together. +""" +import argparse +import bisect +import gzip +import re +import sys +from collections import defaultdict + +# Severus add_pon constants +BUFF = 2000 +CLUST_LEN = 250 +MAX_LEN_DIFF = 50 +MAX_SUM_DIFF = 1000 + +MATE_RE = re.compile(r'[\[\]]([^:\[\]]+):(\d+)[\[\]]') + + +def open_text(path): + return gzip.open(path, 'rt') if path.endswith('.gz') else open(path) + + +def read_pon(path): + """Key (chrom, or (chrom1, chrom2) between chromosomes) -> sorted columns pos1, ci1, end, ci2, partner chrom.""" + rows = defaultdict(list) + with open_text(path) as fp: + for line in fp: + fields = line.strip().split(',') + if len(fields) < 7: + continue + c1, p1, c2, p2, ci1, ci2, svtype = fields[:7] + try: + p1, p2, ci1, ci2 = int(p1), int(p2), int(ci1), int(ci2) + except ValueError: + continue + if svtype == 'INS': + rows[c1].append((p1, ci1, p1 + p2, ci2, c1)) + elif c1 == c2: + rows[c1].append((p1, ci1, p2, ci2, c2)) + else: + rows[(c1, c2)].append((p1, ci1, p2, ci2, c2)) + rows[(c2, c1)].append((p2, ci2, p1, ci1, c1)) + pon = {} + for key, entries in rows.items(): + entries.sort() + pon[key] = list(zip(*entries)) + return pon + + +def info_dict(info): + out = {} + for kv in info.split(';'): + key, _, value = kv.partition('=') + out[key] = value + return out + + +def breakpoints(fields): + """chrom, pos, partner chrom, partner pos, length, is_ins for one VCF record.""" + chrom, pos, alt = fields[0], int(fields[1]), fields[4] + info = info_dict(fields[7]) + if info.get('BP_NOTATION') == '' or info.get('SVTYPE') == 'INS' or alt == '': + try: + length = abs(int(float(info.get('SVLEN', '0')))) + except ValueError: + length = 0 + return chrom, pos, chrom, pos, length, True + mate = MATE_RE.search(alt) + if mate: + chrom2, pos2 = mate.group(1), int(mate.group(2)) + else: + chrom2, pos2 = chrom, int(info.get('END', pos)) + return chrom, pos, chrom2, pos2, abs(pos2 - pos), False + + +def in_pon(pon, chrom, pos, chrom2, pos2, length): + """Severus add_pon, with no breakpoint confidence interval for the query (SAVANA reports none).""" + entries = pon.get(chrom if chrom == chrom2 else (chrom, chrom2)) + if not entries: + return False + starts, cis1, ends, cis2, partners = entries + lo = bisect.bisect_left(starts, pos - BUFF) + hi = bisect.bisect_right(starts, pos + BUFF) + for i in range(lo, hi): + if partners[i] != chrom2: + continue + if chrom == chrom2: + len_diff = abs(ends[i] - starts[i] - length) + else: + len_diff = abs(ends[i] - pos2) + sum_diff = abs(starts[i] - pos) + len_diff + if (len_diff <= MAX_LEN_DIFF and sum_diff <= MAX_SUM_DIFF) or sum_diff <= max(cis1[i], cis2[i]) + CLUST_LEN: + return True + return False + + +def event_id(record_id): + return re.sub(r'_[12]$', '', record_id) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--vcf', required=True, help='SAVANA classified somatic VCF') + parser.add_argument('--pon', required=True, help='Severus SV panel of normals (.tsv or .tsv.gz)') + parser.add_argument('--output', required=True, help='Filtered VCF') + parser.add_argument('--summary', required=True, help='TSV with the event counts before and after') + parser.add_argument('--sample', default='', help='Sample name for the summary') + args = parser.parse_args() + + pon = read_pon(args.pon) + + header, records = [], [] + with open_text(args.vcf) as fp: + for line in fp: + if line.startswith('#'): + header.append(line) + elif line.strip(): + records.append(line) + + # one decision per event, made on its lower breakpoint (Severus's bp_1) + lowest = {} + for line in records: + fields = line.rstrip('\n').split('\t') + bp = breakpoints(fields) + eid = event_id(fields[2]) + if eid not in lowest or (bp[0], bp[1]) < (lowest[eid][0], lowest[eid][1]): + lowest[eid] = bp + germline = {eid for eid, (c, p, c2, p2, length, _ins) in lowest.items() if in_pon(pon, c, p, c2, p2, length)} + + kept = 0 + with open(args.output, 'w') as out: + for line in header: + if line.startswith('#CHROM'): + out.write('##savana_pon_filter=\n') + out.write(line) + for line in records: + if event_id(line.split('\t', 3)[2]) not in germline: + out.write(line) + kept += 1 + + n_events = len(lowest) + with open(args.summary, 'w') as fp: + fp.write('sample\tevents_in\tevents_pon_matched\tevents_out\trecords_in\trecords_out\n') + fp.write(f'{args.sample}\t{n_events}\t{len(germline)}\t{n_events - len(germline)}\t{len(records)}\t{kept}\n') + print(f'{args.sample}: {len(germline)} of {n_events} SAVANA events match the PoN and were removed; ' + f'{n_events - len(germline)} kept ({kept} of {len(records)} records)', file=sys.stderr) + + +if __name__ == '__main__': + main() diff --git a/conf/modules.config b/conf/modules.config index 41fada0d..09c5a2c9 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -1040,6 +1040,7 @@ process { "--max-cn ${params.reconplot_max_cn}", "--format '${params.reconplot_format}'", params.reconplot_min_svlen ? "--min-svlen ${params.reconplot_min_svlen}" : '', + params.reconplot_max_svs ? "--max-svs ${params.reconplot_max_svs}" : '', (params.reconplot_exclude_vntr && sv_source == 'severus') ? '--exclude-vntr' : '' ].join(' ').trim() } diff --git a/docs/output.md b/docs/output.md index bdf69b86..369f84e8 100644 --- a/docs/output.md +++ b/docs/output.md @@ -469,6 +469,15 @@ steps; tumor-only samples run the combined `savana to` command instead, producin in one step. We strongly recommend matched tumor/normal mode for best performance -- SAVANA's own docs note tumor-only calling is a fallback, best combined with population/panel-of-normals filtering. +SAVANA's PacBio classifier (`--pb`) has no tumour-only mode: it rejects a call only on matched-normal +evidence, so a PacBio tumour-only `sample.classified.somatic.vcf` holds nearly every germline SV in the +genome (~28k on a 30x sample). For these samples the pipeline removes the SVs found in the genome's +Severus germline SV panel of normals (Severus's matching rule) and passes +`sample.pon_filtered.classified.somatic.vcf` to VEP, Padfoot and ReConPlot instead; about 97% of the +calls are removed. What remains is still mostly germline (on a PacBio pair, all 761 survivors had +matched-normal support), so treat PacBio tumour-only SAVANA SVs as a fallback and prefer Severus +or a matched normal. Turn this off with `--savana_pb_pon_filter false`. + ``` ├── savana │ ├── sample.sv_breakpoints.vcf @@ -479,6 +488,8 @@ docs note tumor-only calling is a fallback, best combined with population/panel- │ ├── sample.classified.somatic.vcf │ ├── sample.classified.somatic.bedpe │ ├── sample.classified.germline.vcf +│ ├── sample.pon_filtered.classified.somatic.vcf +│ ├── sample.pon_filtered.pon_filter_summary.tsv │ ├── sample_segmented_absolute_copy_number.tsv │ ├── sample_ranked_solutions.tsv │ ├── sample_fitted_purity_ploidy.tsv @@ -488,23 +499,25 @@ docs note tumor-only calling is a fallback, best combined with population/panel- │ ├── 10kbp_bin_ref_all_sample.bed ``` -| File | Description | -| ----------------------------------------------- | ------------------------------------------------------------------------------ | -| `sample.sv_breakpoints.vcf` | Raw (unclassified) SV breakpoints from `savana run`/`savana to` | -| `sample.sv_breakpoints.bedpe` | Raw SV breakpoints in BEDPE format | -| `sample.sv_breakpoints_read_support.tsv` | Supporting-read evidence per breakpoint | -| `sample.inserted_sequences.fa` | Inserted sequences at breakpoints (insertion SVs) | -| `sample.classified.vcf` | All breakpoints after `savana classify` (somatic + germline) | -| `sample.classified.somatic.vcf` | Classified somatic SV VCF -- fed into VEP for annotation | -| `sample.classified.somatic.bedpe` | Classified somatic SVs in BEDPE format (ONT only: not written with `--pb`) | -| `sample.classified.germline.vcf` | Classified germline SVs | -| `sample_segmented_absolute_copy_number.tsv` | Segmented absolute copy-number calls from `savana cna`/`savana to` | -| `sample_ranked_solutions.tsv` | Candidate purity/ploidy solutions, ranked | -| `sample_fitted_purity_ploidy.tsv` | Selected purity/ploidy fit | -| `sample_raw_read_counts.tsv` | Raw binned read counts used for CN segmentation | -| `sample_read_counts__log2r_segmented.tsv` | Segmented log2 read-count ratios per bin (`` is `mnorm` or `self`) | -| `sample_allele_counts_hetSNPs.bed` | Heterozygous-SNP allele counts (only when SNP/allele-frequency input is given) | -| `10kbp_bin_ref_all_sample.bed` | Genome bins used for read-count binning | +| File | Description | +| ----------------------------------------------- | --------------------------------------------------------------------------------- | +| `sample.sv_breakpoints.vcf` | Raw (unclassified) SV breakpoints from `savana run`/`savana to` | +| `sample.sv_breakpoints.bedpe` | Raw SV breakpoints in BEDPE format | +| `sample.sv_breakpoints_read_support.tsv` | Supporting-read evidence per breakpoint | +| `sample.inserted_sequences.fa` | Inserted sequences at breakpoints (insertion SVs) | +| `sample.classified.vcf` | All breakpoints after `savana classify` (somatic + germline) | +| `sample.classified.somatic.vcf` | Classified somatic SV VCF -- fed into VEP for annotation | +| `sample.classified.somatic.bedpe` | Classified somatic SVs in BEDPE format (ONT only: not written with `--pb`) | +| `sample.classified.germline.vcf` | Classified germline SVs | +| `sample.pon_filtered.classified.somatic.vcf` | PacBio tumour-only: somatic SVs not in the SV panel of normals -- used downstream | +| `sample.pon_filtered.pon_filter_summary.tsv` | PacBio tumour-only: events and records before and after the PoN filter | +| `sample_segmented_absolute_copy_number.tsv` | Segmented absolute copy-number calls from `savana cna`/`savana to` | +| `sample_ranked_solutions.tsv` | Candidate purity/ploidy solutions, ranked | +| `sample_fitted_purity_ploidy.tsv` | Selected purity/ploidy fit | +| `sample_raw_read_counts.tsv` | Raw binned read counts used for CN segmentation | +| `sample_read_counts__log2r_segmented.tsv` | Segmented log2 read-count ratios per bin (`` is `mnorm` or `self`) | +| `sample_allele_counts_hetSNPs.bed` | Heterozygous-SNP allele counts (only when SNP/allele-frequency input is given) | +| `10kbp_bin_ref_all_sample.bed` | Genome bins used for read-count binning | #### `deepvariant` diff --git a/modules/local/reconplot/tests/main.nf.test b/modules/local/reconplot/tests/main.nf.test index 2933ac85..99060d9f 100644 --- a/modules/local/reconplot/tests/main.nf.test +++ b/modules/local/reconplot/tests/main.nf.test @@ -287,4 +287,40 @@ nextflow_process { ) } } + + test("savana single source - real, more SVs than --max-svs") { + + config "./nextflow.maxsvs.config" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'savana', + [ file("${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor_segmented_absolute_copy_number.tsv", checkIfExists: true), + file("${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.vcf", checkIfExists: true), + file("${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor_fitted_purity_ploidy.tsv", checkIfExists: true) ], + 'savana', + [] + ] + input[1] = [ [ id:'reconplot' ], file("${projectDir}/assets/reconplot", checkIfExists: true) ] + input[2] = 'hg38' + """ + } + } + + then { + def reconplot_log = path(process.out.log.get(0).get(1)).text + assertAll( + { assert process.success }, + { assert process.out.per_chromosome.get(0).get(1).collect { f -> file(f).name }.sort() == ['test_chr19.pdf', 'test_chr19.png'] }, + { assert process.out.genome_wide.get(0).get(1).collect { f -> file(f).name }.sort() == ['test_chr19.pdf', 'test_chr19.png'] }, + // both drawing runs (per-chromosome, genome-wide) cap the five junctions at three + { assert reconplot_log.count('WARNING: 5 SV junctions exceed --max-svs 3') == 2 }, + // the harmonised table is written before the cap and keeps every junction + { assert path(process.out.tables.get(0).get(1).find { f -> file(f).name.endsWith('reconplot_sv.tsv') }).readLines().size() == 6 } + ) + } + } } diff --git a/modules/local/reconplot/tests/main.nf.test.snap b/modules/local/reconplot/tests/main.nf.test.snap index 9eff3e6f..0efb957b 100644 --- a/modules/local/reconplot/tests/main.nf.test.snap +++ b/modules/local/reconplot/tests/main.nf.test.snap @@ -134,10 +134,10 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,e6b6a074da8a07cdc818276ac832b53c" ] ], - "timestamp": "2026-10-07T11:27:41.013206989", + "timestamp": "2026-10-29T11:36:47.198264848", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.3" @@ -157,10 +157,10 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,e6b6a074da8a07cdc818276ac832b53c" ] ], - "timestamp": "2026-10-07T11:27:53.7954896", + "timestamp": "2026-10-29T11:36:59.850633017", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.3" @@ -180,10 +180,10 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,e6b6a074da8a07cdc818276ac832b53c" ] ], - "timestamp": "2026-10-07T11:28:15.015932537", + "timestamp": "2026-10-29T11:37:22.752485872", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.3" @@ -203,10 +203,10 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,e6b6a074da8a07cdc818276ac832b53c" ] ], - "timestamp": "2026-10-07T11:28:35.052455573", + "timestamp": "2026-10-29T11:37:43.79629842", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.3" @@ -226,10 +226,10 @@ ] ], [ - "versions.yml:md5,05d093292a6b3bc87f03028726cbab19" + "versions.yml:md5,e6b6a074da8a07cdc818276ac832b53c" ] ], - "timestamp": "2026-10-07T11:28:24.926693238", + "timestamp": "2026-10-29T11:37:32.818963092", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.3" diff --git a/modules/local/reconplot/tests/nextflow.maxsvs.config b/modules/local/reconplot/tests/nextflow.maxsvs.config new file mode 100644 index 00000000..7d7b8f00 --- /dev/null +++ b/modules/local/reconplot/tests/nextflow.maxsvs.config @@ -0,0 +1,7 @@ +// --max-svs below the fixture's five SAVANA junctions: the figures draw the three best supported, the table keeps all five +process { + withName: 'RECONPLOT' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + ext.args = '--max-svs 3' + } +} diff --git a/modules/local/reconplot/tests/nextflow.pipeline.config b/modules/local/reconplot/tests/nextflow.pipeline.config index b01ef582..9c22ed97 100644 --- a/modules/local/reconplot/tests/nextflow.pipeline.config +++ b/modules/local/reconplot/tests/nextflow.pipeline.config @@ -5,6 +5,7 @@ params { reconplot_max_cn = 8 reconplot_format = 'pdf,png' reconplot_min_svlen = 0 + reconplot_max_svs = 1000 reconplot_exclude_vntr = true reconplot_regions = 'chr19:1,000,000-1,300,000' reconplot_genes = 'STK11,NOTAGENE' @@ -19,6 +20,7 @@ process { "--max-cn ${params.reconplot_max_cn}", "--format '${params.reconplot_format}'", params.reconplot_min_svlen ? "--min-svlen ${params.reconplot_min_svlen}" : '', + params.reconplot_max_svs ? "--max-svs ${params.reconplot_max_svs}" : '', (params.reconplot_exclude_vntr && sv_source == 'severus') ? '--exclude-vntr' : '' ].join(' ').trim() } diff --git a/modules/local/savana_pon_filter/environment.yml b/modules/local/savana_pon_filter/environment.yml new file mode 100644 index 00000000..c9b822b2 --- /dev/null +++ b/modules/local/savana_pon_filter/environment.yml @@ -0,0 +1,7 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - "conda-forge::python=3.12" diff --git a/modules/local/savana_pon_filter/main.nf b/modules/local/savana_pon_filter/main.nf new file mode 100644 index 00000000..08f79c40 --- /dev/null +++ b/modules/local/savana_pon_filter/main.nf @@ -0,0 +1,40 @@ +process SAVANA_PON_FILTER { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/python:3.12': + 'biocontainers/python:3.12' }" + + input: + tuple val(meta), path(vcf) + tuple val(meta2), path(pon) + + output: + tuple val(meta), path("${prefix}.classified.somatic.vcf"), emit: vcf + tuple val(meta), path("${prefix}.pon_filter_summary.tsv") , emit: summary + tuple val("${task.process}"), val('python'), eval("python --version | sed 's/Python //'"), topic: versions, emit: versions_python + + when: + task.ext.when == null || task.ext.when + + script: + // The output keeps the .classified.somatic.vcf suffix: the ReConPlot wrapper finds SAVANA's SV file by it + prefix = task.ext.prefix ?: "${meta.id}.pon_filtered" + """ + savana_pon_filter.py \\ + --vcf ${vcf} \\ + --pon ${pon} \\ + --output ${prefix}.classified.somatic.vcf \\ + --summary ${prefix}.pon_filter_summary.tsv \\ + --sample ${meta.id} + """ + + stub: + prefix = task.ext.prefix ?: "${meta.id}.pon_filtered" + """ + touch ${prefix}.classified.somatic.vcf + touch ${prefix}.pon_filter_summary.tsv + """ +} diff --git a/modules/local/savana_pon_filter/meta.yml b/modules/local/savana_pon_filter/meta.yml new file mode 100644 index 00000000..ae133a19 --- /dev/null +++ b/modules/local/savana_pon_filter/meta.yml @@ -0,0 +1,82 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "savana_pon_filter" +description: | + Remove SAVANA SVs that match a germline SV panel of normals, using Severus's PoN matching rule. + SAVANA's PacBio classifier has no tumour-only mode, so without a normal its somatic VCF holds + every germline SV with enough tumour reads. +keywords: + - structural variants + - savana + - panel of normals + - tumour-only +tools: + - "python": + description: "Python, running bin/savana_pon_filter.py from this pipeline." + homepage: "https://www.python.org" + documentation: "https://docs.python.org" + licence: ["PSF-2.0"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing sample information + - vcf: + type: file + description: SAVANA classified somatic SV VCF + pattern: "*.classified.somatic.vcf" + - - meta2: + type: map + description: Groovy Map for the reference + - pon: + type: file + description: | + Severus SV panel of normals: comma-separated chr1,pos1,chr2,pos2,ci1,ci2,svtype,fraction + (pos2 is the length for INS), optionally gzipped and with a header row + pattern: "*.{tsv,tsv.gz}" + +output: + vcf: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}.classified.somatic.vcf": + type: file + description: The input VCF without the events that match the panel of normals + pattern: "*.classified.somatic.vcf" + summary: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}.pon_filter_summary.tsv": + type: file + description: Events and records before and after filtering + pattern: "*.pon_filter_summary.tsv" + versions_python: + - - ${task.process}: + type: string + description: Process name + - python: + type: string + description: Tool name + - version: + type: string + description: Tool version + +topics: + versions: + - - ${task.process}: + type: string + description: Process name + - python: + type: string + description: Tool name + - version: + type: string + description: Tool version + +authors: + - "@ljwharbers" +maintainers: + - "@ljwharbers" diff --git a/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz b/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz new file mode 100644 index 00000000..75831ede Binary files /dev/null and b/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz differ diff --git a/modules/local/savana_pon_filter/tests/main.nf.test b/modules/local/savana_pon_filter/tests/main.nf.test new file mode 100644 index 00000000..5c27ee07 --- /dev/null +++ b/modules/local/savana_pon_filter/tests/main.nf.test @@ -0,0 +1,79 @@ +nextflow_process { + + name "Test Process SAVANA_PON_FILTER" + script "../main.nf" + process "SAVANA_PON_FILTER" + + tag "modules" + tag "modules_local" + tag "savana_pon_filter" + tag "small" + + // Fixture PoN (with the header row the CHM13 PoN has) against the 5-event SAVANA fixture: + // INS at chr19:1100400 len 310 -> matches ID_3 (INS chr19:1100500, 300 bp) + // DEL chr19:1215050-1223040 -> matches ID_1 (chr19:1215001-1223000), both mates dropped + // DEL chr19:1560001-1600000 -> same start as ID_2 but 100 kb shorter: kept + // DEL on chr20 -> other chromosome: no effect + test("removes events in the panel of normals and keeps BND mates together") { + + when { + process { + """ + input[0] = [ + [ id:'test_tumor', platform:'pb' ], + file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.vcf", checkIfExists: true) + ] + input[1] = [ + [:], + file("\${projectDir}/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz", checkIfExists: true) + ] + """ + } + } + + then { + def vcf = path(process.out.vcf[0][1]) + def ids = vcf.readLines().findAll { line -> !line.startsWith('#') }.collect { line -> line.split('\t')[2] } + + assertAll( + { assert process.success }, + { assert file(process.out.vcf[0][1]).name == 'test_tumor.pon_filtered.classified.somatic.vcf' }, + { assert ids == ['ID_4_1', 'ID_4_2', 'ID_5_1', 'ID_2_1', 'ID_2_2', 'ID_5_2'] }, + { assert vcf.readLines().any { line -> line.startsWith('##savana_pon_filter=') } }, + { assert path(process.out.summary[0][1]).readLines()[1] == 'test_tumor\t5\t2\t3\t9\t6' }, + { assert snapshot( + process.out.vcf, + process.out.summary, + process.out.findAll { key, _val -> key.startsWith('versions') } + ).match() } + ) + } + } + + test("stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test_tumor', platform:'pb' ], + file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.vcf", checkIfExists: true) + ] + input[1] = [ + [:], + file("\${projectDir}/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz", checkIfExists: true) + ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/savana_pon_filter/tests/main.nf.test.snap b/modules/local/savana_pon_filter/tests/main.nf.test.snap new file mode 100644 index 00000000..d5be6d95 --- /dev/null +++ b/modules/local/savana_pon_filter/tests/main.nf.test.snap @@ -0,0 +1,99 @@ +{ + "removes events in the panel of normals and keeps BND mates together": { + "content": [ + [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.classified.somatic.vcf:md5,861e715d108dcff56734256d6d0f0be2" + ] + ], + [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.pon_filter_summary.tsv:md5,443d046741ba277c056b58f25ed28d7b" + ] + ], + { + "versions_python": [ + [ + "SAVANA_PON_FILTER", + "python", + "3.12.2" + ] + ] + } + ], + "timestamp": "2026-10-29T11:39:20.243084355", + "meta": { + "nf-test": "0.9.4", + "nextflow": "26.04.3" + } + }, + "stub": { + "content": [ + { + "0": [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.classified.somatic.vcf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.pon_filter_summary.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + "SAVANA_PON_FILTER", + "python", + "3.12.2" + ] + ], + "summary": [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.pon_filter_summary.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "vcf": [ + [ + { + "id": "test_tumor", + "platform": "pb" + }, + "test_tumor.pon_filtered.classified.somatic.vcf:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions_python": [ + [ + "SAVANA_PON_FILTER", + "python", + "3.12.2" + ] + ] + } + ], + "timestamp": "2026-10-29T11:38:05.085736965", + "meta": { + "nf-test": "0.9.4", + "nextflow": "26.04.3" + } + } +} \ No newline at end of file diff --git a/nextflow.config b/nextflow.config index 836884c5..6baf0372 100644 --- a/nextflow.config +++ b/nextflow.config @@ -141,6 +141,7 @@ params { savana_chromosomes = null // e.g. '19' to restrict SAVANA_CNA/SAVANA_TO's allele counting and binning to given chromosome indices savana_allele_min_reads = null // e.g. 2 to lower SAVANA_CNA/SAVANA_TO's het-SNP coverage floor (test profiles only, default is SAVANA's own 10) savana_allele_mapq = null // e.g. 0 to lower SAVANA_CNA/SAVANA_TO's het-SNP mapping-quality floor (test profiles only, default is SAVANA's own 5) + savana_pb_pon_filter = true // drop PacBio tumour-only SAVANA SVs found in the Severus germline SV panel of normals // Padfoot options padfoot_genome = null // hg38 | chm13 | mm10; null = inferred from --genome @@ -152,6 +153,7 @@ params { reconplot_genome = null // hg38 | hg19 | T2T | mm10 | mm39; null = inferred from --genome reconplot_max_cn = 8 reconplot_min_svlen = 0 // drop intra-chromosomal SVs shorter than this (bp); 0 = keep all + reconplot_max_svs = 1000 // draw at most this many SVs per figure, best supported first; 0 = no limit reconplot_exclude_vntr = false // drop Severus SVs flagged inside a VNTR reconplot_regions = null // e.g. "chr8,chr17:30000000-50000000" -> extra focus/ panel reconplot_genes = null // comma-separated HUGO symbols labelled on the focus panel diff --git a/nextflow_schema.json b/nextflow_schema.json index 46c7fe25..b1521850 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -425,6 +425,12 @@ "savana_allele_mapq": { "type": "integer", "description": "Mapping quality threshold for reads counted at a het-SNP site (SAVANA's --allele_mapq, default 5). Intended for test profiles running on downsampled data below that coverage floor." + }, + "savana_pb_pon_filter": { + "type": "boolean", + "default": true, + "description": "Remove SAVANA SVs found in the germline SV panel of normals from PacBio tumour-only samples.", + "help_text": "SAVANA's PacBio classifier (`--pb`) only rejects a call on matched-normal evidence, so without a normal nearly every germline SV is reported as somatic (~28k on a 30x genome). The somatic VCF is then filtered against the genome's Severus panel of normals (`pon_file`) with Severus's matching rule, which removes ~97% of the calls. What remains is still mostly germline (on a PacBio sample with a matched normal, all 761 survivors were germline), so prefer Severus or a matched normal for PacBio tumour-only SVs. The unfiltered VCF is still published. Has no effect on ONT (SAVANA's tumour-only model) or paired samples, or on a genome without a `pon_file`." } } }, @@ -486,6 +492,13 @@ "help_text": "Intra-chromosomal SVs shorter than this many bp are left out of the figures and the harmonised SV table; translocations are always kept. 0 keeps every SV.", "description": "Drop intra-chromosomal SVs shorter than this many bp (translocations are kept). 0 keeps everything." }, + "reconplot_max_svs": { + "type": "integer", + "default": 1000, + "minimum": 0, + "description": "Draw at most this many SVs per ReConPlot figure, the best supported first. 0 draws every SV.", + "help_text": "ReConPlot adds one plot layer per SV, so tens of thousands take hours and the genome-wide panel fails (R node stack overflow) from about 2,000 SVs. Above this limit the figures show the SVs with the most tumour read support (then the longest), and their title says `top N of M SVs`. The harmonised `*.reconplot_sv.tsv` table keeps every SV." + }, "reconplot_exclude_vntr": { "type": "boolean", "help_text": "When true, Severus SVs flagged as lying inside a VNTR are left out of the figures (Severus pairs only).", diff --git a/subworkflows/local/tumor_only/tests/tumoronly_savana.nf.test b/subworkflows/local/tumor_only/tests/tumoronly_savana.nf.test new file mode 100644 index 00000000..45f28193 --- /dev/null +++ b/subworkflows/local/tumor_only/tests/tumoronly_savana.nf.test @@ -0,0 +1,81 @@ +nextflow_workflow { + + name "Test Workflow TUMORONLY_SAVANA" + script "../tumoronly_savana.nf" + workflow "TUMORONLY_SAVANA" + + tag "subworkflows" + tag "subworkflows_local" + tag "tumoronly_savana" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // No test samplesheet has a PacBio tumour-only sample, so the pipeline tests never reach SAVANA_PON_FILTER. + // Under -stub what is tested is the routing: PacBio somatic VCFs go through the PoN filter, ONT ones (SAVANA's + // tumour-only model) do not, and without a PoN every VCF passes unfiltered. + + test("PacBio somatic VCFs are PoN-filtered, ONT ones are not") { + + options "-stub" + + when { + workflow { + """ + // stand-ins: the SAVANA_TO stub does not read its BAM + def bam = file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.vcf", checkIfExists: true) + def bai = file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.bedpe", checkIfExists: true) + input[0] = channel.of( + [ [ id:'sample_pb', paired_data:false, platform:'pb' ], bam, bai, [], [] ], + [ [ id:'sample_ont', paired_data:false, platform:'ont' ], bam, bai, [], [] ] + ) + input[1] = channel.value([ [:], file("\${projectDir}/tests/fixtures/padfoot_reconplot/chr19_1-2Mb.fasta", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/tests/fixtures/padfoot_reconplot/chr19_1-2Mb.fasta.fai", checkIfExists: true) ]) + input[3] = channel.value([ [:], [] ]) + input[4] = channel.value([ [:], [] ]) + input[5] = channel.value([ [:], file("\${projectDir}/modules/local/savana_pon_filter/tests/fixtures/pon.tsv.gz", checkIfExists: true) ]) + """ + } + } + + then { + def names = workflow.out.somatic_vcf.collectEntries { meta, vcf -> [meta.id, file(vcf).name] } + assertAll( + { assert workflow.success }, + { assert names == [ + sample_pb : 'sample_pb.pon_filtered.classified.somatic.vcf', + sample_ont: 'sample_ont.classified.somatic.vcf' + ] } + ) + } + } + + test("without a PoN every somatic VCF passes unfiltered") { + + options "-stub" + + when { + workflow { + """ + // stand-ins: the SAVANA_TO stub does not read its BAM + def bam = file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.vcf", checkIfExists: true) + def bai = file("\${projectDir}/tests/fixtures/padfoot_reconplot/test_tumor.classified.somatic.bedpe", checkIfExists: true) + input[0] = channel.of( + [ [ id:'sample_pb', paired_data:false, platform:'pb' ], bam, bai, [], [] ] + ) + input[1] = channel.value([ [:], file("\${projectDir}/tests/fixtures/padfoot_reconplot/chr19_1-2Mb.fasta", checkIfExists: true) ]) + input[2] = channel.value([ [:], file("\${projectDir}/tests/fixtures/padfoot_reconplot/chr19_1-2Mb.fasta.fai", checkIfExists: true) ]) + input[3] = channel.value([ [:], [] ]) + input[4] = channel.value([ [:], [] ]) + input[5] = channel.value([ [:], [] ]) + """ + } + } + + then { + assertAll( + { assert workflow.success }, + { assert workflow.out.somatic_vcf.collect { meta, vcf -> file(vcf).name } == ['sample_pb.classified.somatic.vcf'] } + ) + } + } +} diff --git a/subworkflows/local/tumor_only/tumoronly_savana.nf b/subworkflows/local/tumor_only/tumoronly_savana.nf index 0267870c..6cf75396 100644 --- a/subworkflows/local/tumor_only/tumoronly_savana.nf +++ b/subworkflows/local/tumor_only/tumoronly_savana.nf @@ -1,5 +1,6 @@ // IMPORT MODULES -include { SAVANA_TO } from '../../../modules/nf-core/savana/to/main' +include { SAVANA_TO } from '../../../modules/nf-core/savana/to/main' +include { SAVANA_PON_FILTER } from '../../../modules/local/savana_pon_filter/main' workflow TUMORONLY_SAVANA { @@ -9,6 +10,7 @@ workflow TUMORONLY_SAVANA { fai // [[:], fai] contigs // [[:], contigs] -- one contig per line, restricts SAVANA to canonical chromosomes g1000_vcf // [[:], '1000g_hg38' | '1000g_t2t'] -- bundled 1000 Genomes population SNP set (genome-specific) + sv_pon // [[:], pon.tsv.gz] or [[:], []] -- Severus germline SV panel of normals; PacBio calls are filtered against it main: // .first() re-establishes a value channel: fasta/fai are singletons reused for every sample, @@ -44,8 +46,32 @@ workflow TUMORONLY_SAVANA { g1000_vcf ) + // + // MODULE: SAVANA_PON_FILTER (label: process_single) + // SAVANA's PacBio classifier only rejects a call on normal-read evidence, so without a normal every germline + // SV with enough tumour reads comes out somatic (~28k on a PacBio genome, against ~2 with its matched normal). + // ONT uses SAVANA's tumour-only model and is left alone. Without a PoN for the genome the VCF passes unfiltered. + // Input: [meta, classified.somatic.vcf] / [[:], pon.tsv.gz] + // Output: .vcf -- [meta, vcf] -- the somatic VCF without events in the panel of normals + // + SAVANA_TO.out.somatic_vcf + .combine(sv_pon) + .branch { meta, vcf, _pon_meta, pon -> + filter: meta.platform == 'pb' && pon + return [meta, vcf] + keep: true + return [meta, vcf] + } + .set { savana_somatic_vcf } + // savana_somatic_vcf.filter / .keep: [meta, vcf] + + SAVANA_PON_FILTER ( + savana_somatic_vcf.filter, + sv_pon + ) + emit: - somatic_vcf = SAVANA_TO.out.somatic_vcf // [meta, vcf] -- classified somatic SV VCF + somatic_vcf = savana_somatic_vcf.keep.mix(SAVANA_PON_FILTER.out.vcf) // [meta, vcf] -- classified somatic SV VCF (PoN-filtered for PacBio) cn_calls = SAVANA_TO.out.cna // [meta, tsv] -- segmented absolute copy number fitted_purity_ploidy = SAVANA_TO.out.fitted_purity_ploidy // [meta, tsv] -- selected purity/ploidy fit (absent when no fit) allele_counts = SAVANA_TO.out.allele_counts // [meta, bed] -- het-SNP allele counts (ReConPlot BAF track) diff --git a/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index e6f56c9e..e9108597 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -1111,13 +1111,17 @@ workflow LRSOMATIC { // Tumor-only has no matched germline control, so allele counting uses the bundled 1000g // population SNP set instead of a (nonexistent) germline VCF -- see TUMORONLY_SAVANA. ch_savana_g1000_vcf = channel.value([[:], params.savana_g1000_vcf ?: []]) + // SAVANA's PacBio classifier has no tumour-only mode, so its PacBio tumour-only calls are filtered against + // the Severus germline SV panel of normals -- see TUMORONLY_SAVANA + ch_savana_sv_pon = channel.value([[:], params.savana_pb_pon_filter && params.pon_file ? file(params.pon_file) : []]) TUMORONLY_SAVANA ( tumoronly_savana_input, ch_fasta, ch_fai, ch_savana_contigs, - ch_savana_g1000_vcf + ch_savana_g1000_vcf, + ch_savana_sv_pon ) PAIRED_SAVANA (