Add Padfoot annotation and ReConPlot figures for the SV/CNA callers

CHANGELOG.md

@@ -17,6 +17,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - [#176](https://github.com/IntGenomicsLab/lrsomatic/pull/176) - Added a `solution_dirs` output to the WAKHAN module so its per-solution copy-number plots can be staged downstream (@ljwharbers).
 - [#193](https://github.com/IntGenomicsLab/lrsomatic/pull/193) - Added VEP plugins: AlphaMissense, SIFT/PolyPhen, ClinVar and REVEL on GRCh38, and AlphaMissense plus SIFT/PolyPhen on CHM13 via protein-space lookup. Enabled by default with `--genome GRCh38` or `--genome CHM13` (a first GRCh38 run downloads around 1.4 GB); any `--vep_*` path overrides its default and `--skip_vep_plugins` turns the set off. CADD and EVE are opt-in behind `--vep_cadd_snv`/`--vep_cadd_indel` and `--vep_eve` because of their size (81 GB and 9.6 GB); prepared REVEL and EVE files are published to `<outdir>/vep_plugins/` for reuse. Two lab-hosted AlphaMissense files are CC BY 4.0 with attribution in `CITATIONS.md` (@AmberVerhasselt).
 - [#189](https://github.com/IntGenomicsLab/lrsomatic/pull/189) - Added SAVANA structural variant and copy-number calling, running alongside Severus/ASCAT (@yannvrb).
+- [#202](https://github.com/IntGenomicsLab/lrsomatic/pull/202) - Added Padfoot annotation of somatic SVs and CNAs for every caller pair with output for a sample, Severus + Wakhan and SAVANA, paired and tumour-only, published to `<sample>/padfoot/<pair>/`. Padfoot is fetched from a pinned commit of [Tim-Yu/Padfoot](https://github.com/Tim-Yu/Padfoot) (adds SAVANA input) and runs in `ghcr.io/tim-yu/padfoot-repeatmasker` with RepeatMasker on by default; `--skip_padfoot` and the `--padfoot_*` parameters control it (@Tim-Yu).
+- [#202](https://github.com/IntGenomicsLab/lrsomatic/pull/202) - Added ReConPlot rearrangement and copy-number figures per CN/SV caller pair, `reconplot/{severus_ascat,severus_wakhan,savana}/`, with per-chromosome, genome-wide and optional focus panels plus the harmonised CN/SV tables. The wrapper ships in `assets/reconplot/`, the ReConPlot R package is fetched from a pinned commit, and the module runs in `ghcr.io/tim-yu/reconplot`; `--skip_reconplot` and the `--reconplot_*` parameters control it (@Tim-Yu).
+- [#202](https://github.com/IntGenomicsLab/lrsomatic/pull/202) - `PAIRED_SAVANA` and `TUMORONLY_SAVANA` also emit the somatic BEDPE, the fitted purity/ploidy and the het-SNP allele counts, so Padfoot and ReConPlot can consume SAVANA output (@Tim-Yu).
+- [#202](https://github.com/IntGenomicsLab/lrsomatic/pull/202) - Added stub nf-tests for `PADFOOT_ANNOTATION` and `RECONPLOT_FIGURES` that check every caller pair is formed and unfitted samples are dropped (tag `small`) (@Tim-Yu).
 
 ### `Changed`
 

CITATIONS.md

@@ -50,6 +50,10 @@
 
   > De Coster W. cramino: A fast and simple tool for quality control of long read sequencing data [Software]. GitHub. https://github.com/wdecoster/cramino
 
+- [Dfam](https://pubmed.ncbi.nlm.nih.gov/33436076/)
+
+  > Storer J, Hubley R, Rosen J, Wheeler TJ, Smit AF. The Dfam community resource of transposable element families, sequence models, and genome annotations. Mob DNA. 2021 Jan 11;12(1):2. doi: 10.1186/s13100-020-00230-y. PubMed PMID: 33436076; PubMed Central PMCID: PMC7805199.
+
 - [Ensembl VEP](https://pubmed.ncbi.nlm.nih.gov/27268795/)
 
   > McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GR, Thormann A, Flicek P, Cunningham F. The Ensembl Variant Effect Predictor. Genome Biol. 2016 Jun 6;17(1):122. doi: 10.1186/s13059-016-0974-4. PubMed PMID: 27268795; PubMed Central PMCID: PMC4893825.
@@ -62,6 +66,10 @@
 
   > Jha, A., Bohaczuk, S. C., Mao, Y., Ranchalis, J., Mallory, B. J., Min, A. T., Hamm, M. O., Swanson, E., Dubocanin, D., Finkbeiner, C., Li, T., Whittington, D., Noble, W. S., Stergachis, A. B., & Vollger, M. R. (2024). DNA-m6A calling and integrated long-read epigenetic and genetic analysis with fibertools. Genome Research.
 
+- [ggplot2](https://ggplot2.tidyverse.org)
+
+  > Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016. ISBN 978-3-319-24277-4. https://ggplot2.tidyverse.org
+
 - [LongPhase](https://pubmed.ncbi.nlm.nih.gov/35104333/)
 
   > Lin JH, Chen LC, Yu SC, Huang YT. LongPhase: an ultra-fast chromosome-scale phasing algorithm for small and large variants. Bioinformatics. 2022 Apr 28;38(9):2452-2455. doi: 10.1093/bioinformatics/btac126. PubMed PMID: 35253834; PubMed Central PMCID: PMC9048675.
@@ -90,10 +98,26 @@
 
   > De Coster W, Rademakers R. NanoPack2: population-scale evaluation of long-read sequencing data. Bioinformatics. 2023 May 4;39(5):btad311. doi: 10.1093/bioinformatics/btad311. PubMed PMID: 37171891; PubMed Central PMCID: PMC10196664.
 
+- [Padfoot](https://github.com/KolmogorovLab/Padfoot)
+
+  > Functional annotation of somatic structural variants and copy number alterations from long-read cancer genomes [Software]. GitHub. https://github.com/KolmogorovLab/Padfoot
+
 - [PolyPhen-2](https://doi.org/10.1038/nmeth0410-248)
 
   > Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, Kondrashov AS, Sunyaev SR. A method and server for predicting damaging missense mutations. Nat Methods. 2010 Apr;7(4):248-9. doi: 10.1038/nmeth0410-248.
 
+- [R](https://www.R-project.org/)
+
+  > R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
+
+- [ReConPlot](https://pubmed.ncbi.nlm.nih.gov/38058190/)
+
+  > 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.
+
+- [RepeatMasker](https://www.repeatmasker.org/)
+
+  > Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0 [Software]. https://www.repeatmasker.org/
+
 - [REVEL](https://doi.org/10.1016/j.ajhg.2016.08.016)
 
   > Ioannidis NM, Rothstein JH, Pejaver V, Middha S, McDonnell SK, Baheti S, Musolf A, Li Q, Holzinger E, Karyadi D, Cannon-Albright LA, Teerlink CC, Stanford JL, Isaacs WB, Xu J, Cooney KA, Lange EM, Schleutker J, Carpten JD, Powell IJ, Cussenot O, Cancel-Tassin G, Giles GG, MacInnis RJ, Maier C, Hsieh CL, Wiklund F, Catalona WJ, Foulkes WD, Mandal D, Eeles RA, Kote-Jarai Z, Bustamante CD, Schaid DJ, Hastie T, Ostrander EA, Bailey-Wilson JE, Radivojac P, Thibodeau SN, Whittemore AS, Sieh W. REVEL: An Ensemble Method for Predicting the Pathogenicity of Rare Missense Variants. Am J Hum Genet. 2016 Oct 6;99(4):877-885. doi: 10.1016/j.ajhg.2016.08.016.

README.md

@@ -66,6 +66,12 @@ a. Small variant annotation ([`VEP`](https://github.com/Ensembl/ensembl-vep)), w
 
 b. Structural variant annotation ([`VEP`](https://github.com/Ensembl/ensembl-vep))
 
+c. Somatic SV and CNA functional annotation ([`Padfoot`](https://github.com/KolmogorovLab/Padfoot))
+
+**5) Visualisation:**
+
+a. Rearrangement and copy-number figures per CN/SV caller pair ([`ReConPlot`](https://github.com/cortes-ciriano-lab/ReConPlot))
+
 <!-- TODO nf-core: Include a figure that guides the user through the major workflow steps. Many nf-core
      workflows use the "tube map" design for that. See https://nf-co.re/docs/guidelines/graphic_design/workflow_diagrams#examples for examples.   -->
 
@@ -108,21 +114,27 @@ IntGenomicsLab/lr_somatic was originally written by Luuk Harbers, Robert Forsyth
 
 The main output is an aligned and phased tumour BAM, per-sample QC from `cramino`, `mosdepth`, `samtools` and optionally `fibertools`, a MultiQC report, and a self-contained per-sample HTML report (`<sample>/report/<sample>_report.html`; disable it with `--skip_report`).
 
-Variant and copy number callers (`clairS`, `clairS-TO`, `severus`, `savana`, `ascat`) write to their own folders; 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 and ReConPlot figures go to `padfoot/` and `reconplot/`, one subfolder per caller pair. See the [output documentation](/docs/output.md).
 
 Example output directory structure:
 
 ```
 ├── Sample 1
 │    ├── ascat
 │    ├── bamfiles
+│    ├── padfoot
+│    │   ├── savana
+│    │   └── severus_wakhan
 │    ├── qc
 │    │    ├── tumor
 │    │    │   ├── cramino_aln
 │    │    │   ├── cramino_ubam
 │    │    │   ├── fibertoolsrs
 │    │    │   ├── mosdepth
 │    │    │   ├── samtools
+│    ├── reconplot
+│    │   ├── savana
+│    │   └── severus_wakhan
 │    ├── variants
 │    │   ├──clairS-TO
 │    │   ├──severus
@@ -136,6 +148,9 @@ Example output directory structure:
 ├── Sample 2
 │    ├── ascat
 │    ├── bamfiles
+│    ├── padfoot
+│    │   ├── savana
+│    │   └── severus_wakhan
 │    ├── qc
 │    │    ├── tumor
 │    │    │   ├── cramino_aln
@@ -149,6 +164,10 @@ Example output directory structure:
 │    │    │   ├── fibertoolsrs
 │    │    │   ├── mosdepth
 │    │    │   ├── samtools
+│    ├── reconplot
+│    │   ├── savana
+│    │   ├── severus_ascat
+│    │   └── severus_wakhan
 │    ├── variants
 │    │   ├── clair3
 │    │   ├── clairS

assets/reconplot/R/parsers/ascat.R

@@ -0,0 +1,192 @@
+## ---------------------------------------------------------------------------
+## parsers/ascat.R -- ASCAT allele-specific copy number (CN component).
+##
+## Used by lrsomatic's long-read ASCAT step, but works for any ASCAT run
+## (nf-core/sarek included) since the file names come from ASCAT itself.
+##
+## Files consumed:
+##   <sample>.segments.txt      -> CN track (default; ascat.output$segments)
+##   <sample>.segments_raw.txt  -> CN track with --ascat-cn-mode raw
+##   <sample>.purityploidy.txt  -> purity (aberrant cell fraction) and ploidy
+##   <sample>.tumour_tumourBAF.txt -> optional per-SNP BAF annotation track
+##
+## segments.txt holds the final fitted integer allele-specific calls
+## (nMajor/nMinor), which is what ReConPlot wants: copyNumber = nMajor + nMinor
+## and minorAlleleCopyNumber = nMinor. cnvs.txt is the same table minus the
+## sample column, so it is only used as a fallback. segments_raw.txt is the
+## pre-fit ASPCF segmentation and additionally carries the unrounded nAraw /
+## nBraw, exposed via --ascat-cn-mode raw / raw-fractional.
+##
+## ASCAT writes chromosomes without a "chr" prefix; normalize_chrom() fixes that.
+## ---------------------------------------------------------------------------
+
+ASCAT_FILE_PATTERNS <- list(
+  segments     = "\\.segments\\.txt$",
+  segments_raw = "\\.segments_raw\\.txt$",
+  cnvs         = "\\.cnvs\\.txt$",
+  purityploidy = "\\.purityploidy\\.txt$",
+  metrics      = "\\.metrics\\.txt$",
+  baf          = "_tumourBAF\\.txt$"      # not _tumourBAF_rawBAF.txt
+)
+
+ascat_find_file <- function(dir, key, sample = NULL, required = TRUE) {
+  pat <- ASCAT_FILE_PATTERNS[[key]]
+  hits <- list.files(dir, pattern = pat, full.names = TRUE)
+  if (key == "segments") hits <- hits[!grepl("segments_raw", hits)]
+  if (!is.null(sample) && nzchar(sample)) {
+    pinned <- hits[startsWith(basename(hits), sample)]
+    if (length(pinned)) hits <- pinned
+  }
+  if (length(hits) == 0) {
+    if (required) stop("No ASCAT '", key, "' file matching /", pat, "/ in ", dir)
+    return(NULL)
+  }
+  if (length(hits) > 1) {
+    log_msg("  multiple ASCAT '", key, "' files; using ", basename(hits[1]),
+            " (pass --sample to disambiguate)")
+  }
+  hits[1]
+}
+
+ascat_infer_sample <- function(dir) {
+  f <- list.files(dir, pattern = ASCAT_FILE_PATTERNS$segments)
+  f <- f[!grepl("segments_raw", f)]
+  if (length(f) == 0) return(NA_character_)
+  sub("\\.segments\\.txt$", "", f[1])
+}
+
+#' Read ASCAT segments into the ReConPlot CN contract.
+#'
+#' @param mode "segments" (fitted integer calls), "raw" (pre-fit segmentation,
+#'   still using the rounded nMajor/nMinor) or "raw-fractional" (pre-fit
+#'   segmentation using the unrounded nAraw/nBraw)
+ascat_read_cn <- function(file, mode = "segments", sample = NULL) {
+  df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE,
+                          na.strings = c("NA", "nan", "NaN", ""))
+  need <- c("chr", "startpos", "endpos", "nMajor", "nMinor")
+  missing <- setdiff(need, names(df))
+  if (length(missing)) {
+    stop("ASCAT file ", basename(file), " lacks column(s): ", paste(missing, collapse = ", "))
+  }
+  ## segments.txt keeps a sample column ("<prefix>.tumour"); cnvs.txt does not.
+  if ("sample" %in% names(df) && !is.null(sample) && nzchar(sample)) {
+    hit <- grepl(sample, df$sample, fixed = TRUE)
+    if (any(hit)) {
+      if (any(!hit)) log_msg(sprintf("  kept %d/%d segments for sample '%s'",
+                                     sum(hit), nrow(df), sample))
+      df <- df[hit, , drop = FALSE]
+    }
+  }
+
+  if (mode == "raw-fractional") {
+    if (!all(c("nAraw", "nBraw") %in% names(df))) {
+      stop("--ascat-cn-mode raw-fractional needs nAraw/nBraw (segments_raw.txt)")
+    }
+    total <- as.numeric(df$nAraw) + as.numeric(df$nBraw)
+    minor <- as.numeric(df$nBraw)
+  } else {
+    total <- as.numeric(df$nMajor) + as.numeric(df$nMinor)
+    minor <- as.numeric(df$nMinor)
+  }
+
+  data.frame(chr = df$chr, start = df$startpos, end = df$endpos,
+             copyNumber = total, minorAlleleCopyNumber = minor,
+             stringsAsFactors = FALSE)
+}
+
+ascat_read_purity <- function(file) {
+  if (is.null(file) || !file.exists(file)) return(NULL)
+  df <- utils::read.table(file, sep = "\t", header = TRUE, stringsAsFactors = FALSE)
+  if (nrow(df) == 0) return(NULL)
+  purity <- if ("AberrantCellFraction" %in% names(df)) df$AberrantCellFraction[1] else NA_real_
+  ploidy <- if ("Ploidy" %in% names(df)) df$Ploidy[1] else NA_real_
+  if (is.na(purity) && is.na(ploidy)) return(NULL)
+  list(purity = purity, ploidy = ploidy)
+}
+
+#' BAF annotation track from ASCAT's per-SNP tumour BAF table.
+#'
+#' Layout is: <rowname> Chromosome Position <sample>. Values are already
+#' restricted to germline het SNPs, so no extra filtering is needed.
+ascat_baf_annotation_fn <- function(file, max_points_per_region = 5000) {
+  if (is.null(file) || !file.exists(file)) return(NULL)
+  cache <- NULL
+  function(regions) {
+    if (is.null(cache)) {
+      log_msg("  reading ASCAT BAF track: ", basename(file))
+      d <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE)
+      cols <- names(d)
+      chr_col <- cols[tolower(cols) == "chromosome"][1]
+      pos_col <- cols[tolower(cols) == "position"][1]
+      val_col <- setdiff(cols, c(chr_col, pos_col, cols[1]))
+      if (is.na(chr_col) || is.na(pos_col) || length(val_col) == 0) {
+        log_msg("  BAF file has an unexpected layout; skipping annotation track")
+        return(NULL)
+      }
+      d <- data.frame(chr = normalize_chrom(d[[chr_col]]),
+                      pos = as.numeric(d[[pos_col]]),
+                      y = suppressWarnings(as.numeric(d[[val_col[1]]])),
+                      stringsAsFactors = FALSE)
+      cache <<- d[is_main_chrom(d$chr) & !is.na(d$y), , drop = FALSE]
+    }
+    out <- do.call(rbind, lapply(seq_len(nrow(regions)), function(i) {
+      r <- regions[i, ]
+      sub <- cache[cache$chr == r$chr & cache$pos >= r$start & cache$pos <= r$end, , drop = FALSE]
+      if (nrow(sub) > max_points_per_region) {
+        sub <- sub[sort(sample.int(nrow(sub), max_points_per_region)), , drop = FALSE]
+      }
+      sub
+    }))
+    if (is.null(out) || nrow(out) == 0) return(NULL)
+    rownames(out) <- NULL
+    out
+  }
+}
+
+parse_ascat <- function(args) {
+  dir <- args$input
+  sample <- args$sample
+  if ((is.null(sample) || !nzchar(sample)) && !is.null(dir)) {
+    sample <- ascat_infer_sample(dir)
+    if (!is.na(sample)) log_msg("  inferred ASCAT sample: ", sample)
+  }
+
+  mode <- match.arg(args$ascat_cn_mode %||% "segments",
+                    c("segments", "raw", "raw-fractional"))
+  cn_file <- args$cn_file
+  if (is.null(cn_file)) {
+    key <- if (mode == "segments") "segments" else "segments_raw"
+    cn_file <- ascat_find_file(dir, key, sample, required = FALSE)
+    if (is.null(cn_file) && mode == "segments") {
+      cn_file <- ascat_find_file(dir, "cnvs", sample, required = FALSE)
+      if (!is.null(cn_file)) log_msg("  segments.txt absent; falling back to cnvs.txt")
+    }
+    if (is.null(cn_file)) stop("No ASCAT copy number file found in ", dir)
+  }
+  log_msg("  CN file: ", basename(cn_file), " (mode: ", mode, ")")
+  cn <- ascat_read_cn(cn_file, mode = mode, sample = sample)
+
+  meta <- list(sample = sample, cn_file = cn_file, cn_caller = "ASCAT")
+  pp <- ascat_read_purity(args$purity_file %||%
+                            (if (!is.null(dir)) ascat_find_file(dir, "purityploidy", sample, FALSE)))
+  if (!is.null(pp)) {
+    meta$purity <- pp$purity
+    meta$ploidy <- pp$ploidy
+    log_msg(sprintf("  purity=%s ploidy=%s", pp$purity, round(as.numeric(pp$ploidy), 3)))
+  }
+
+  ann_fn <- NULL
+  if (isTRUE(args$baf_track)) {
+    baf_file <- args$annotation_file %||%
+      (if (!is.null(dir)) ascat_find_file(dir, "baf", sample, required = FALSE))
+    ann_fn <- ascat_baf_annotation_fn(baf_file,
+                                      max_points_per_region = args$baf_max_points %||% 5000)
+    if (is.null(ann_fn)) log_msg("  --baf-track requested but no ASCAT BAF file found; skipping")
+  }
+
+  list(cn = cn, meta = meta, annotation_fn = ann_fn)
+}
+
+register_parser("ascat", parse_ascat,
+                "ASCAT allele-specific CN (segments.txt / segments_raw.txt)",
+                provides = "cn")