Filter PacBio tumour-only SAVANA SVs against the SV PoN; cap ReConPlot SVs (#225)

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

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

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

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.

assets/reconplot/VERSION

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

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)

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') == '<INS>' or info.get('SVTYPE') == 'INS' or alt == '<INS>':
+        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=<Description="Events matching the germline SV panel of normals '
+                          f'{args.pon.split("/")[-1]} removed (Severus PoN matching rule)">\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()

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

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_<mode>_log2r_segmented.tsv` | Segmented log2 read-count ratios per bin (`<mode>` 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_<mode>_log2r_segmented.tsv` | Segmented log2 read-count ratios per bin (`<mode>` 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`