diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f6e2a92..c51c08f0 100644 --- a/CHANGELOG.md +++ b/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 `/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 `/padfoot//`. 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` diff --git a/CITATIONS.md b/CITATIONS.md index a1403be7..351e0fbb 100644 --- a/CITATIONS.md +++ b/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. diff --git a/README.md b/README.md index d9ffaf9f..4a8c2149 100644 --- a/README.md +++ b/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)) + @@ -108,7 +114,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; 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: @@ -116,6 +122,9 @@ Example output directory structure: ├── Sample 1 │ ├── ascat │ ├── bamfiles +│ ├── padfoot +│ │ ├── savana +│ │ └── severus_wakhan │ ├── qc │ │ ├── tumor │ │ │ ├── cramino_aln @@ -123,6 +132,9 @@ Example output directory structure: │ │ │ ├── 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 diff --git a/assets/reconplot/R/parsers/ascat.R b/assets/reconplot/R/parsers/ascat.R new file mode 100644 index 00000000..08db7a3f --- /dev/null +++ b/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: +## .segments.txt -> CN track (default; ascat.output$segments) +## .segments_raw.txt -> CN track with --ascat-cn-mode raw +## .purityploidy.txt -> purity (aberrant cell fraction) and ploidy +## .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 (".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: Chromosome Position . 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") diff --git a/assets/reconplot/R/parsers/generic.R b/assets/reconplot/R/parsers/generic.R new file mode 100644 index 00000000..8fb334c0 --- /dev/null +++ b/assets/reconplot/R/parsers/generic.R @@ -0,0 +1,108 @@ +## --------------------------------------------------------------------------- +## parsers/generic.R -- fallback parser for any caller not yet given its own file. +## +## Reads an explicit --cn-file and --sv-file and maps common column spellings +## onto the ReConPlot contract. Use it to sanity-check a new caller's output +## before writing a dedicated parser. +## +## CN : chr/chromosome/seqnames, start, end, copyNumber (or nMajor+nMinor, +## total_cn, CN), minorAlleleCopyNumber (or nMinor, minor_cn) +## SV : headerless BEDPE (>=6 cols, orientation in col 9/10 or a "|"-packed +## name column), or a headed table with chr1/pos1/chr2/pos2/strands. +## --------------------------------------------------------------------------- + +.syn <- function(df, target, candidates) { + if (target %in% names(df)) return(target) + hit <- candidates[tolower(candidates) %in% tolower(names(df))] + if (length(hit) == 0) return(NA_character_) + names(df)[match(tolower(hit[1]), tolower(names(df)))] +} + +generic_read_cn <- function(file) { + df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE, + na.strings = c("NA", "nan", "NaN", "")) + names(df) <- sub("^#", "", names(df)) + + cchr <- .syn(df, "chr", c("chromosome", "chrom", "seqnames", "CHR", "Chromosome")) + cstart <- .syn(df, "start", c("startpos", "chromStart", "Start", "begin")) + cend <- .syn(df, "end", c("endpos", "chromEnd", "End", "stop")) + ctot <- .syn(df, "copyNumber", c("total_cn", "totalCN", "cn", "CN", "ntot", "tcn")) + cmin <- .syn(df, "minorAlleleCopyNumber", c("minor_cn", "minorCN", "nMinor", + "nMin", "minor", "mcn")) + cmaj <- .syn(df, "nMajor", c("major_cn", "majorCN", "nMaj", "major")) + + for (nm in c(chr = cchr, start = cstart, end = cend)) { + if (is.na(nm)) stop("Cannot find chr/start/end columns in ", basename(file)) + } + out <- data.frame(chr = df[[cchr]], start = df[[cstart]], end = df[[cend]], + stringsAsFactors = FALSE) + if (!is.na(ctot)) { + out$copyNumber <- df[[ctot]] + } else if (!is.na(cmaj) && !is.na(cmin)) { + out$copyNumber <- as.numeric(df[[cmaj]]) + as.numeric(df[[cmin]]) + } else { + stop("Cannot find a total copy number column in ", basename(file)) + } + out$minorAlleleCopyNumber <- if (!is.na(cmin)) df[[cmin]] else NA_real_ + out +} + +generic_read_sv <- function(file) { + first <- readLines(file, n = 1L, warn = FALSE) + headed <- grepl("chr1|chrom1|CHROM|chr_1", first, ignore.case = TRUE) && + !grepl("^chr[0-9XYM]", first) + + if (headed) { + df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE) + names(df) <- sub("^#", "", names(df)) + c1 <- .syn(df, "chr1", c("chrom1", "chromosome1", "chrA", "CHROM_A")) + p1 <- .syn(df, "pos1", c("start1", "position1", "posA", "START_A")) + c2 <- .syn(df, "chr2", c("chrom2", "chromosome2", "chrB", "CHROM_B")) + p2 <- .syn(df, "pos2", c("start2", "position2", "posB", "START_B")) + st <- .syn(df, "strands", c("strand", "orientation", "bp_notation", "BP_NOTATION")) + if (any(is.na(c(c1, p1, c2, p2)))) { + stop("Cannot find chr1/pos1/chr2/pos2 columns in ", basename(file)) + } + strands <- if (!is.na(st)) df[[st]] else { + s1 <- .syn(df, "strand1", c("str1", "orientation1")) + s2 <- .syn(df, "strand2", c("str2", "orientation2")) + if (any(is.na(c(s1, s2)))) stop("No strand information in ", basename(file)) + paste0(df[[s1]], df[[s2]]) + } + sv <- data.frame(chr1 = df[[c1]], pos1 = df[[p1]], chr2 = df[[c2]], + pos2 = df[[p2]], strands = strands, stringsAsFactors = FALSE) + } else { + df <- data.table::fread(file, sep = "\t", header = FALSE, data.table = FALSE) + if (ncol(df) < 6) stop("BEDPE ", basename(file), " has fewer than 6 columns") + strands <- if (ncol(df) >= 10) { + paste0(df[[9]], df[[10]]) + } else if (ncol(df) >= 7 && any(grepl("\\|", df[[7]]))) { + vapply(strsplit(as.character(df[[7]]), "|", fixed = TRUE), + function(p) p[length(p)], character(1)) + } else { + stop("BEDPE ", basename(file), " carries no orientation (need cols 9-10 or a name field)") + } + sv <- data.frame(chr1 = df[[1]], pos1 = bedpe_point(df[[2]], df[[3]]), + chr2 = df[[4]], pos2 = bedpe_point(df[[5]], df[[6]]), + strands = strands, stringsAsFactors = FALSE) + } + sv$strands <- gsub("[<>]", "", trimws(sv$strands)) + sv$strands[toupper(sv$strands) == "INS"] <- "INS" + sv$strands[toupper(sv$strands) %in% c("SBE", "SBND")] <- "SBE" + sv +} + +parse_generic <- function(args) { + if (is.null(args$cn_file) || is.null(args$sv_file)) { + stop("--source generic requires both --cn-file and --sv-file") + } + log_msg(" CN file: ", basename(args$cn_file)) + log_msg(" SV file: ", basename(args$sv_file)) + list(cn = generic_read_cn(args$cn_file), + sv = generic_read_sv(args$sv_file), + meta = list(sample = args$sample %||% NA_character_), + annotation_fn = NULL) +} + +register_parser("generic", parse_generic, + "Any caller: explicit --cn-file/--sv-file with auto column mapping") diff --git a/assets/reconplot/R/parsers/lrsomatic.R b/assets/reconplot/R/parsers/lrsomatic.R new file mode 100644 index 00000000..8a6f355c --- /dev/null +++ b/assets/reconplot/R/parsers/lrsomatic.R @@ -0,0 +1,77 @@ +## --------------------------------------------------------------------------- +## parsers/lrsomatic.R -- legacy ASCAT+Severus lrsomatic convenience wrapper. +## +## lrsomatic can run different CN/SV callers. Prefer naming those callers +## explicitly with --cn-source/--sv-source. This wrapper remains as the original +## ASCAT+Severus shortcut: +## +## /ascat/ -> parsers/ascat.R (CN) +## /variants/severus/somatic_SVs/ -> parsers/severus.R (SVs) +## +## Equivalent to: +## --cn-source ascat --cn-input /ascat \ +## --sv-source severus --sv-input /variants/severus +## +## For Wakhan CN from the same lrsomatic run, use: +## --cn-source wakhan --sv-source severus --input +## +## --input accepts the per-sample directory (output/P215003155), the pipeline +## outdir (output/, when it holds exactly one sample), or the run directory +## containing output/. +## --------------------------------------------------------------------------- + +#' Resolve --input to an lrsomatic per-sample directory. +lrsomatic_sample_dir <- function(dir, sample = NULL) { + if (is.null(dir)) stop("ASCAT+Severus lrsomatic shortcut needs --input") + is_sample_dir <- function(d) dir.exists(file.path(d, "ascat")) || + dir.exists(file.path(d, "variants", "severus")) + + if (is_sample_dir(dir)) return(dir) + + ## descend through a run directory into output/ + for (d in c(file.path(dir, "output"), dir)) { + if (!dir.exists(d)) next + subs <- list.dirs(d, recursive = FALSE) + subs <- subs[vapply(subs, is_sample_dir, logical(1))] + if (!is.null(sample) && nzchar(sample)) { + pinned <- subs[basename(subs) == sample | startsWith(basename(subs), sample)] + if (length(pinned)) subs <- pinned + } + if (length(subs) == 1) return(subs[1]) + if (length(subs) > 1) { + stop("Several lrsomatic sample directories under ", d, ": ", + paste(basename(subs), collapse = ", "), ". Pass --sample or point ", + "--input at one of them.") + } + } + stop("Could not find an lrsomatic sample directory under ", dir, + " (expected an 'ascat' and/or 'variants/severus' subdirectory)") +} + +parse_lrsomatic <- function(args) { + sample_dir <- lrsomatic_sample_dir(args$input, args$sample) + log_msg(" lrsomatic sample directory: ", sample_dir) + if (is.null(args$sample) || !nzchar(args$sample)) { + args$sample <- basename(sample_dir) + log_msg(" inferred sample: ", args$sample) + } + + args_cn <- args; args_cn$input <- file.path(sample_dir, "ascat") + args_sv <- args; args_sv$input <- file.path(sample_dir, "variants", "severus") + if (!dir.exists(args_cn$input)) stop("No ascat/ directory under ", sample_dir) + if (!dir.exists(args_sv$input)) stop("No variants/severus/ directory under ", sample_dir) + + ## ASCAT names its files ".tumour...", so do not pin on the bare + ## sample id here; let the ascat parser infer its own prefix. + args_cn$sample <- NULL + p_cn <- parse_ascat(args_cn) + p_sv <- parse_severus(args_sv) + + meta <- utils::modifyList(p_sv$meta %||% list(), p_cn$meta %||% list()) + meta$sample <- args$sample + list(cn = p_cn$cn, sv = p_sv$sv, meta = meta, annotation_fn = p_cn$annotation_fn) +} + +register_parser("lrsomatic", parse_lrsomatic, + "Legacy lrsomatic shortcut: ASCAT copy number + Severus somatic SVs", + provides = c("cn", "sv")) diff --git a/assets/reconplot/R/parsers/registry.R b/assets/reconplot/R/parsers/registry.R new file mode 100644 index 00000000..3c24657c --- /dev/null +++ b/assets/reconplot/R/parsers/registry.R @@ -0,0 +1,104 @@ +## --------------------------------------------------------------------------- +## registry.R -- plug-in table mapping a caller name to a parser function. +## +## Adding support for a new SV/CN caller means dropping one file into +## R/parsers/ that ends with a register_parser() call. Nothing else changes. +## +## A parser is a function(args) -> list(...) where `args` is the list of CLI +## options (see run_reconplot.R). What it must return depends on `provides`: +## +## provides = c("cn", "sv") full parser: list(cn = , sv = , ...) +## provides = "cn" CN component: list(cn = , ...) +## provides = "sv" SV component: list(sv = , ...) +## +## Optional extra elements: `meta` (named list; sample/purity/ploidy feed the +## plot title) and `annotation_fn` (function(regions) -> data.frame(chr,pos,y)). +## +## Components exist because a pipeline may call CN and SVs with different tools +## (lrsomatic: ASCAT + Severus). They compose via --cn-source / --sv-source, so +## any CN caller can be paired with any SV caller. +## +## Parsers should return raw-ish tables; validate_cn()/validate_sv() from +## utils_common.R are applied centrally afterwards. +## --------------------------------------------------------------------------- + +.parser_registry <- new.env(parent = emptyenv()) + +register_parser <- function(name, fn, description = "", provides = c("cn", "sv")) { + stopifnot(is.character(name), length(name) == 1, is.function(fn)) + provides <- match.arg(provides, c("cn", "sv"), several.ok = TRUE) + assign(name, list(fn = fn, description = description, provides = provides), + envir = .parser_registry) + invisible(NULL) +} + +get_parser_entry <- function(name) { + if (is.null(name) || !nzchar(name) || + !exists(name, envir = .parser_registry, inherits = FALSE)) { + stop("Unknown source '", name, "'. Available: ", + paste(list_parsers()$name, collapse = ", ")) + } + get(name, envir = .parser_registry, inherits = FALSE) +} + +get_parser <- function(name) get_parser_entry(name)$fn + +parser_provides <- function(name) get_parser_entry(name)$provides + +list_parsers <- function() { + nms <- sort(ls(.parser_registry)) + data.frame( + name = nms, + provides = vapply(nms, function(n) + paste(get(n, envir = .parser_registry)$provides, collapse = "+"), character(1)), + description = vapply(nms, function(n) + get(n, envir = .parser_registry)$description, character(1)), + stringsAsFactors = FALSE, row.names = NULL) +} + +#' Source every parser implementation found in `dir`. +load_parsers <- function(dir) { + files <- setdiff(list.files(dir, pattern = "\\.[Rr]$", full.names = TRUE), + file.path(dir, "registry.R")) + for (f in sort(files)) sys.source(f, envir = globalenv()) + invisible(list_parsers()) +} + +#' Work out which parser supplies CN and which supplies SVs, then run them. +#' +#' Single-source mode calls one full parser. Split mode (--cn-source and/or +#' --sv-source) calls one parser for each half, each with its own --cn-input / +#' --sv-input if given, and merges the results. +run_parsers <- function(args) { + cn_src <- args$cn_source %||% args$source + sv_src <- args$sv_source %||% args$source + split_mode <- !identical(cn_src, sv_src) || + !is.null(args$cn_input) || !is.null(args$sv_input) + + if (!split_mode) { + prov <- parser_provides(cn_src) + if (!all(c("cn", "sv") %in% prov)) { + stop("Source '", cn_src, "' only provides '", paste(prov, collapse = "+"), + "'. Pair it with another via --cn-source/--sv-source, e.g. ", + "--cn-source ascat --sv-source severus") + } + log_msg("source: ", cn_src) + return(get_parser(cn_src)(args)) + } + + log_msg("CN source: ", cn_src, " SV source: ", sv_src) + for (s in c(cn = cn_src, sv = sv_src)) invisible(get_parser_entry(s)) + + args_cn <- args; args_cn$input <- args$cn_input %||% args$input + args_sv <- args; args_sv$input <- args$sv_input %||% args$input + + p_cn <- get_parser(cn_src)(args_cn) + p_sv <- get_parser(sv_src)(args_sv) + if (is.null(p_cn$cn)) stop("Source '", cn_src, "' returned no copy number table") + if (is.null(p_sv$sv)) stop("Source '", sv_src, "' returned no SV table") + + list(cn = p_cn$cn, sv = p_sv$sv, + ## CN metadata wins: purity/ploidy come from the CN fit + meta = utils::modifyList(p_sv$meta %||% list(), p_cn$meta %||% list()), + annotation_fn = p_cn$annotation_fn %||% p_sv$annotation_fn) +} diff --git a/assets/reconplot/R/parsers/savana.R b/assets/reconplot/R/parsers/savana.R new file mode 100644 index 00000000..516c7c40 --- /dev/null +++ b/assets/reconplot/R/parsers/savana.R @@ -0,0 +1,258 @@ +## --------------------------------------------------------------------------- +## parsers/savana.R -- turn a SAVANA output directory into ReConPlot inputs. +## +## Files consumed (all optional except the first two): +## _segmented_absolute_copy_number.tsv -> CN track +## .classified.somatic.bedpe -> SV track (default) +## .classified.somatic.vcf -> SV track (--sv-format vcf) +## _fitted_purity_ploidy.tsv -> purity/ploidy for the title +## _allele_counts_hetSNPs.bed -> optional BAF annotation track +## +## SAVANA writes breakpoint orientation as BP_NOTATION ("+-", "-+", "++", "--", +## ""), which is the same convention ReConPlot uses, so orientations pass +## through unchanged; only ""/"" need renaming. +## --------------------------------------------------------------------------- + +SAVANA_FILE_PATTERNS <- list( + cn = "_segmented_absolute_copy_number\\.tsv$", + bedpe = "\\.classified\\.somatic\\.bedpe$", + vcf = "\\.classified\\.somatic\\.vcf(\\.gz)?$", + purity = "_fitted_purity_ploidy\\.tsv$", + hetsnp = "_allele_counts_hetSNPs\\.bed(\\.gz)?$" +) + +#' Locate a single SAVANA output file inside `dir`, optionally pinned to a sample. +savana_find_file <- function(dir, key, sample = NULL, required = TRUE) { + pat <- SAVANA_FILE_PATTERNS[[key]] + hits <- list.files(dir, pattern = pat, full.names = TRUE) + 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 SAVANA '", key, "' file matching /", pat, "/ in ", dir) + return(NULL) + } + if (length(hits) > 1) { + log_msg(" multiple '", key, "' files found; using ", basename(hits[1]), + " (pass --sample to disambiguate)") + } + hits[1] +} + +#' Infer the sample prefix from the copy number file name. +savana_infer_sample <- function(dir) { + f <- list.files(dir, pattern = SAVANA_FILE_PATTERNS$cn) + if (length(f) == 0) return(NA_character_) + sub("_segmented_absolute_copy_number\\.tsv$", "", f[1]) +} + +## --- copy number ------------------------------------------------------------ + +savana_read_cn <- function(file) { + df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE, + na.strings = c("NA", "nan", "NaN", "")) + names(df)[names(df) == "chromosome"] <- "chr" + need <- c("chr", "start", "end", "copyNumber", "minorAlleleCopyNumber") + missing <- setdiff(need, names(df)) + if (length(missing)) { + stop("SAVANA CN file ", basename(file), " lacks column(s): ", + paste(missing, collapse = ", ")) + } + df[, need, drop = FALSE] +} + +## --- structural variants ---------------------------------------------------- + +#' Normalise a SAVANA BP_NOTATION value to ReConPlot's `strands` vocabulary. +savana_notation_to_strands <- function(x) { + x <- toupper(gsub("[<>]", "", trimws(as.character(x)))) + out <- rep(NA_character_, length(x)) + out[x %in% c("+-", "-+", "++", "--")] <- x[x %in% c("+-", "-+", "++", "--")] + out[x == "INS"] <- "INS" + out[x %in% c("SBE", "SBND", "BND_SINGLE")] <- "SBE" + out +} + +#' Pull the fields SAVANA packs into the BEDPE name column. +#' +#' The name looks like "ID_39590|1608bp|TUMOUR_12|--": an event id, the SV +#' length, the tumour read support and the breakpoint notation. +savana_split_bedpe_name <- function(name) { + parts <- strsplit(as.character(name), "|", fixed = TRUE) + pick <- function(p, rx, default = NA_character_) { + m <- grep(rx, p, value = TRUE) + if (length(m)) m[1] else default + } + data.frame( + sv_id = vapply(parts, function(p) p[1] %||% NA_character_, character(1)), + svlen = vapply(parts, function(p) { + v <- pick(p, "^[0-9]+bp$"); suppressWarnings(as.numeric(sub("bp$", "", v))) + }, numeric(1)), + support = vapply(parts, function(p) { + v <- pick(p, "^TUMOUR_[0-9]+$"); suppressWarnings(as.numeric(sub("^TUMOUR_", "", v))) + }, numeric(1)), + notation = vapply(parts, function(p) p[length(p)] %||% NA_character_, character(1)), + stringsAsFactors = FALSE + ) +} + +#' BEDPE coordinates: SAVANA writes start == end == POS, while a spec-compliant +#' BEDPE writes start == POS-1. Handle both without guessing globally. +bedpe_point <- function(start, end) { + ifelse(end - start == 1, end, start) +} + +savana_read_sv_bedpe <- function(file, min_support = 0) { + df <- data.table::fread(file, sep = "\t", header = FALSE, data.table = FALSE) + if (ncol(df) < 7) stop("SAVANA BEDPE ", basename(file), " has fewer than 7 columns") + names(df)[1:7] <- c("chrom1", "start1", "end1", "chrom2", "start2", "end2", "name") + info <- savana_split_bedpe_name(df$name) + + sv <- data.frame( + sv_id = info$sv_id, + chr1 = df$chrom1, + pos1 = bedpe_point(df$start1, df$end1), + chr2 = df$chrom2, + pos2 = bedpe_point(df$start2, df$end2), + strands = savana_notation_to_strands(info$notation), + svlen = info$svlen, + support = info$support, + stringsAsFactors = FALSE + ) + savana_filter_sv(sv, min_support) +} + +savana_read_sv_vcf <- function(file, min_support = 0, pass_only = TRUE) { + vcf <- vcf_read_records(file, pass_only = pass_only) + if (is.null(vcf) || nrow(vcf) == 0) return(savana_empty_sv()) + + notation <- vcf_info_get(vcf$info, "BP_NOTATION") + svlen <- suppressWarnings(as.numeric(vcf_info_get(vcf$info, "SVLEN"))) + support <- suppressWarnings(as.numeric(vcf_info_get(vcf$info, "TUMOUR_READ_SUPPORT"))) + mate <- vcf_alt_mate(vcf$alt) + + sv <- data.frame( + sv_id = sub("_[12]$", "", vcf$id), + chr1 = vcf$chrom, + pos1 = vcf$pos, + chr2 = ifelse(is.na(mate$chr), vcf$chrom, mate$chr), + pos2 = ifelse(is.na(mate$pos), vcf$pos, mate$pos), + strands = savana_notation_to_strands(notation), + svlen = svlen, + support = support, + stringsAsFactors = FALSE + ) + sv <- dedupe_breakend_pairs(sv) # BND records come in mate pairs + savana_filter_sv(sv, min_support) +} + +savana_empty_sv <- function() { + data.frame(sv_id = character(), chr1 = character(), pos1 = numeric(), + chr2 = character(), pos2 = numeric(), strands = character(), + svlen = numeric(), support = numeric(), stringsAsFactors = FALSE) +} + +savana_filter_sv <- function(sv, min_support = 0) { + n0 <- nrow(sv) + if (min_support > 0) { + keep <- is.na(sv$support) | sv$support >= min_support + if (any(!keep)) { + log_msg(sprintf(" dropped %d SVs with tumour read support < %g", + sum(!keep), min_support)) + } + sv <- sv[keep, , drop = FALSE] + } + log_msg(sprintf(" %d SV junctions read (%d after support filter)", n0, nrow(sv))) + rownames(sv) <- NULL + sv +} + +## --- purity / ploidy -------------------------------------------------------- + +savana_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) + as.list(df[1, , drop = FALSE]) +} + +## --- optional BAF annotation track ----------------------------------------- + +#' Build a `custom_annotation` data frame (chr, pos, y) of het-SNP BAF values. +#' +#' The het-SNP file is large, so it is read once, restricted to the primary +#' contigs and thinned; the per-region subset happens in the returned closure. +savana_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 het-SNP BAF track: ", basename(file)) + d <- data.table::fread(file, sep = "\t", header = FALSE, data.table = FALSE, + select = c(1, 2, 12), + col.names = c("chr", "pos", "y")) + d$chr <- normalize_chrom(d$chr) + d <- d[is_main_chrom(d$chr) & !is.na(d$y), , drop = FALSE] + cache <<- d + } + 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 + } +} + +## --- entry point ------------------------------------------------------------ + +parse_savana <- function(args) { + dir <- args$input + sample <- args$sample + if (is.null(sample) || !nzchar(sample)) { + sample <- if (!is.null(dir)) savana_infer_sample(dir) else NA_character_ + if (!is.na(sample)) log_msg(" inferred sample: ", sample) + } + + cn_file <- args$cn_file %||% savana_find_file(dir, "cn", sample) + log_msg(" CN file: ", basename(cn_file)) + cn <- savana_read_cn(cn_file) + + sv_format <- match.arg(args$sv_format %||% "bedpe", c("bedpe", "vcf")) + sv_file <- args$sv_file %||% savana_find_file(dir, sv_format, sample) + log_msg(" SV file: ", basename(sv_file), " (", sv_format, ")") + sv <- if (sv_format == "bedpe") { + savana_read_sv_bedpe(sv_file, min_support = args$min_support %||% 0) + } else { + savana_read_sv_vcf(sv_file, min_support = args$min_support %||% 0) + } + + meta <- list(sample = sample, cn_file = cn_file, sv_file = sv_file) + pp <- savana_read_purity(args$purity_file %||% + (if (!is.null(dir)) savana_find_file(dir, "purity", sample, required = FALSE))) + if (!is.null(pp)) { + meta$purity <- pp$purity + meta$ploidy <- pp$ploidy + log_msg(sprintf(" purity=%s ploidy=%s", pp$purity, pp$ploidy)) + } + + ann_fn <- NULL + if (isTRUE(args$baf_track)) { + baf_file <- args$annotation_file %||% + (if (!is.null(dir)) savana_find_file(dir, "hetsnp", sample, required = FALSE)) + ann_fn <- savana_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 het-SNP file found; skipping") + } + + list(cn = cn, sv = sv, meta = meta, annotation_fn = ann_fn) +} + +register_parser("savana", parse_savana, + "SAVANA (segmented_absolute_copy_number.tsv + classified.somatic.bedpe/vcf)") diff --git a/assets/reconplot/R/parsers/severus.R b/assets/reconplot/R/parsers/severus.R new file mode 100644 index 00000000..f75a5a3d --- /dev/null +++ b/assets/reconplot/R/parsers/severus.R @@ -0,0 +1,189 @@ +## --------------------------------------------------------------------------- +## parsers/severus.R -- Severus long-read somatic SV calls (SV component). +## +## Files consumed: +## somatic_SVs/severus_somatic.vcf.gz -> SV track (the full callset) +## somatic_SVs/breakpoint_clusters_list.tsv -> optional cluster summary log +## +## Why the VCF and not breakpoint_clusters.tsv: the cluster file lists only the +## junctions Severus assigned to a complex-rearrangement cluster. On sample +## P215003155 that is 213 of the 330 junctions in the VCF (every cluster +## junction is present in the VCF; 117 VCF junctions are in no cluster). It is +## a filter/annotation layer, not the callset. Cluster membership is available +## from the VCF anyway as INFO/CLUSTERID, so --cluster-id filters without ever +## reading the TSV. +## +## Severus INFO/STRANDS uses the same orientation convention as ReConPlot and +## SAVANA (DEL "+-", DUP "-+", INV "++"/"--"), verified against records both +## callers report for this sample, so orientations pass straight through. +## +## Note on defaults: Severus runs at --min-support 3 in lrsomatic and emits +## small VNTR indels alongside rearrangements (99 junctions under 1 kb here, +## versus 19 from SAVANA). Nothing is filtered by default -- use --min-svlen +## and --exclude-vntr to get a rearrangement-only view. +## --------------------------------------------------------------------------- + +SEVERUS_FILE_PATTERNS <- list( + somatic_vcf = "^severus_somatic\\.vcf(\\.gz)?$", + all_vcf = "^severus_all\\.vcf(\\.gz)?$", + clusters = "^breakpoint_clusters\\.tsv$", + cluster_list = "^breakpoint_clusters_list\\.tsv$" +) + +#' Find a Severus file, searching the given directory and the usual subdirs. +#' +#' Accepts the severus/ directory, the somatic_SVs/ directory, or an lrsomatic +#' per-sample directory, so callers do not have to know the layout. +severus_find_file <- function(dir, key, required = TRUE, somatic = TRUE) { + pat <- SEVERUS_FILE_PATTERNS[[key]] + sub <- if (somatic) "somatic_SVs" else "all_SVs" + candidates <- c(dir, + file.path(dir, sub), + file.path(dir, "severus", sub), + file.path(dir, "variants", "severus", sub)) + for (d in candidates) { + if (!dir.exists(d)) next + hits <- list.files(d, pattern = pat, full.names = TRUE) + if (length(hits)) return(hits[1]) + } + if (required) stop("No Severus '", key, "' file matching /", pat, "/ under ", dir) + NULL +} + +#' Log the largest Severus clusters so the user can pick one for --cluster-id. +severus_log_cluster_list <- function(file, top = 5) { + if (is.null(file) || !file.exists(file)) return(invisible(NULL)) + df <- utils::read.table(file, sep = "\t", header = FALSE, skip = 1, + stringsAsFactors = FALSE, quote = "", fill = TRUE) + if (ncol(df) < 5) return(invisible(NULL)) + names(df)[c(1, 2, 5)] <- c("cluster_id", "type", "sv_count") + df <- df[order(-suppressWarnings(as.numeric(df$sv_count))), , drop = FALSE] + n <- min(top, nrow(df)) + log_msg(sprintf(" largest Severus clusters: %s", + paste(sprintf("%s (%s, %s SVs)", df$cluster_id[1:n], + df$type[1:n], df$sv_count[1:n]), collapse = "; "))) + invisible(df) +} + +#' Read severus_somatic.vcf.gz into the ReConPlot SV contract. +severus_read_vcf <- function(file, min_support = 0, min_svlen = 0, + exclude_vntr = FALSE, cluster_id = NULL, + clustered_only = FALSE, pass_only = TRUE) { + vcf <- vcf_read_records(file, pass_only = pass_only) + if (is.null(vcf) || nrow(vcf) == 0) return(severus_empty_sv()) + + svtype <- vcf_info_get(vcf$info, "SVTYPE") + strands <- vcf_info_get(vcf$info, "STRANDS") + end <- suppressWarnings(as.numeric(vcf_info_get(vcf$info, "END"))) + svlen <- suppressWarnings(as.numeric(vcf_info_get(vcf$info, "SVLEN"))) + cluster <- vcf_info_get(vcf$info, "CLUSTERID") + detail <- vcf_info_get(vcf$info, "DETAILED_TYPE") + vntr <- !is.na(vcf_info_get(vcf$info, "INSIDE_VNTR")) + support <- suppressWarnings(as.numeric(vcf_format_get(vcf$format, vcf$sample1, "DV"))) + + ## Mate locus: breakends carry it in the ALT allele, symbolic ALTs + ## (//) carry it as INFO/END on the same chromosome. + mate <- vcf_alt_mate(vcf$alt) + chr2 <- ifelse(!is.na(mate$chr), mate$chr, vcf$chrom) + pos2 <- ifelse(!is.na(mate$pos), mate$pos, + ifelse(!is.na(end), end, vcf$pos)) + + ## STRANDS is absent on insertions; fall back to the SV type. + strands[is.na(strands) & svtype == "INS"] <- "INS" + strands[is.na(strands) & svtype == "DEL"] <- "+-" + strands[is.na(strands) & svtype == "DUP"] <- "-+" + + sv <- data.frame( + sv_id = sub("_[12]$", "", vcf$id), + chr1 = vcf$chrom, pos1 = vcf$pos, + chr2 = chr2, pos2 = pos2, + strands = strands, + svtype = svtype, + detailed_type = detail, + svlen = svlen, + support = support, + cluster_id = cluster, + inside_vntr = vntr, + stringsAsFactors = FALSE) + + n_records <- nrow(sv) + sv <- dedupe_breakend_pairs(sv) + log_msg(sprintf(" %d VCF records -> %d junctions after mate collapsing", + n_records, nrow(sv))) + severus_filter_sv(sv, min_support = min_support, min_svlen = min_svlen, + exclude_vntr = exclude_vntr, cluster_id = cluster_id, + clustered_only = clustered_only) +} + +severus_empty_sv <- function() { + data.frame(sv_id = character(), chr1 = character(), pos1 = numeric(), + chr2 = character(), pos2 = numeric(), strands = character(), + svtype = character(), detailed_type = character(), svlen = numeric(), + support = numeric(), cluster_id = character(), inside_vntr = logical(), + stringsAsFactors = FALSE) +} + +severus_filter_sv <- function(sv, min_support = 0, min_svlen = 0, + exclude_vntr = FALSE, cluster_id = NULL, + clustered_only = FALSE) { + drop <- function(sv, keep, why) { + if (any(!keep)) log_msg(sprintf(" dropped %d junctions: %s", sum(!keep), why)) + sv[keep, , drop = FALSE] + } + if (min_support > 0) { + sv <- drop(sv, is.na(sv$support) | sv$support >= min_support, + sprintf("read support < %g", min_support)) + } + if (min_svlen > 0) { + ## A length only means something within a chromosome; never use it to + ## discard translocations, which carry no SVLEN. + intra <- sv$chr1 == sv$chr2 + span <- ifelse(!is.na(sv$svlen), abs(sv$svlen), abs(sv$pos2 - sv$pos1)) + sv <- drop(sv, !intra | is.na(span) | span >= min_svlen, + sprintf("intra-chromosomal span < %g bp", min_svlen)) + } + if (exclude_vntr) { + sv <- drop(sv, !sv$inside_vntr, "inside a VNTR") + } + if (!is.null(cluster_id) && nzchar(cluster_id)) { + wanted <- trimws(unlist(strsplit(cluster_id, ","))) + sv <- drop(sv, !is.na(sv$cluster_id) & sv$cluster_id %in% wanted, + paste("not in cluster", paste(wanted, collapse = "/"))) + } else if (clustered_only) { + sv <- drop(sv, !is.na(sv$cluster_id), "not assigned to any cluster") + } + rownames(sv) <- NULL + sv +} + +parse_severus <- function(args) { + dir <- args$input + somatic <- !isTRUE(args$severus_all) + sv_file <- args$sv_file %||% + severus_find_file(dir, if (somatic) "somatic_vcf" else "all_vcf", somatic = somatic) + log_msg(" SV file: ", basename(sv_file), + if (somatic) " (somatic)" else " (all SVs, germline included)") + + sv <- severus_read_vcf( + sv_file, + min_support = args$min_support %||% 0, + min_svlen = args$min_svlen %||% 0, + exclude_vntr = isTRUE(args$exclude_vntr), + cluster_id = args$cluster_id, + clustered_only = isTRUE(args$clustered_only)) + + if (nrow(sv) > 0) { + tab <- table(sv$svtype, useNA = "no") + log_msg(" Severus SV types: ", + paste(sprintf("%s=%d", names(tab), as.integer(tab)), collapse = " ")) + } + severus_log_cluster_list(severus_find_file(dir, "cluster_list", required = FALSE, + somatic = somatic)) + + list(sv = sv, meta = list(sample = args$sample, sv_file = sv_file, + sv_caller = "Severus")) +} + +register_parser("severus", parse_severus, + "Severus long-read SVs (somatic_SVs/severus_somatic.vcf.gz)", + provides = "sv") diff --git a/assets/reconplot/R/parsers/wakhan.R b/assets/reconplot/R/parsers/wakhan.R new file mode 100644 index 00000000..a29acccf --- /dev/null +++ b/assets/reconplot/R/parsers/wakhan.R @@ -0,0 +1,163 @@ +## --------------------------------------------------------------------------- +## parsers/wakhan.R -- Wakhan allele-specific copy number (CN component). +## +## Files consumed: +## solutions_ranks.tsv -> best solution +## /bed_output/*_copynumbers_segments_HP_1.bed -> haplotype 1 CN +## /bed_output/*_copynumbers_segments_HP_2.bed -> haplotype 2 CN +## +## The two haplotype BED files cover the full segmentation. ReConPlot wants +## total CN plus minor-allele CN, so join HP1/HP2 by segment coordinates and +## calculate total = HP1 + HP2, minor = min(HP1, HP2). +## --------------------------------------------------------------------------- + +WAKHAN_FILE_PATTERNS <- list( + ranks = "^solutions_ranks\\.tsv$", + hp1 = "_copynumbers_segments_HP_1\\.bed(\\.gz)?$", + hp2 = "_copynumbers_segments_HP_2\\.bed(\\.gz)?$" +) + +wakhan_is_dir <- function(dir) { + if (is.null(dir) || !dir.exists(dir)) return(FALSE) + file.exists(file.path(dir, "solutions_ranks.tsv")) || + dir.exists(file.path(dir, "solution_1", "bed_output")) +} + +#' Resolve --input to a Wakhan output directory. +wakhan_dir <- function(dir, sample = NULL) { + if (is.null(dir)) stop("--cn-source wakhan needs --input or --cn-input") + + candidates <- c(dir, file.path(dir, "wakhan")) + for (d in c(file.path(dir, "output"), dir)) { + if (!dir.exists(d)) next + subs <- list.dirs(d, recursive = FALSE) + if (!is.null(sample) && nzchar(sample)) { + pinned <- subs[basename(subs) == sample | startsWith(basename(subs), sample)] + if (length(pinned)) subs <- pinned + } + candidates <- c(candidates, file.path(subs, "wakhan"), subs) + } + + candidates <- unique(normalizePath(candidates, mustWork = FALSE)) + hits <- candidates[vapply(candidates, wakhan_is_dir, logical(1))] + if (length(hits) == 1) return(hits[1]) + if (length(hits) > 1) { + stop("Several Wakhan output directories found: ", paste(hits, collapse = ", "), + ". Pass --sample or point --cn-input at one of them.") + } + stop("Could not find a Wakhan output directory under ", dir, + " (expected solutions_ranks.tsv and/or solution_1/bed_output)") +} + +wakhan_read_rank <- function(file) { + if (is.null(file) || !file.exists(file)) return(NULL) + df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE, + na.strings = c("NA", "nan", "NaN", "")) + if (nrow(df) == 0) return(NULL) + if ("solution_rank" %in% names(df)) { + df <- df[order(suppressWarnings(as.numeric(df$solution_rank))), , drop = FALSE] + } + df[1, , drop = FALSE] +} + +wakhan_solution_dir <- function(dir, rank = NULL) { + ranks <- wakhan_read_rank(file.path(dir, "solutions_ranks.tsv")) + if (is.null(rank) && !is.null(ranks) && "solution_rank" %in% names(ranks)) { + rank <- suppressWarnings(as.integer(ranks$solution_rank[1])) + } + rank <- rank %||% 1L + sol_candidates <- character() + if (!is.null(ranks) && "repository_name" %in% names(ranks)) { + sol_candidates <- c(sol_candidates, file.path(dir, ranks$repository_name[1])) + } + sol_candidates <- c(sol_candidates, file.path(dir, paste0("solution_", rank))) + subs <- list.dirs(dir, recursive = FALSE, full.names = TRUE) + sol_candidates <- c(sol_candidates, subs[dir.exists(file.path(subs, "bed_output"))]) + sol_candidates <- unique(normalizePath(sol_candidates, mustWork = FALSE)) + sol_dir <- sol_candidates[dir.exists(file.path(sol_candidates, "bed_output"))] + if (length(sol_dir) == 0) { + stop("No Wakhan solution bed_output directory under ", dir, + " (checked ", paste(basename(sol_candidates), collapse = ", "), ")") + } + sol_dir <- sol_dir[1] + list(dir = sol_dir, rank = rank, ranks = ranks) +} + +wakhan_find_file <- function(dir, key, sample = NULL, required = TRUE) { + pat <- WAKHAN_FILE_PATTERNS[[key]] + hits <- list.files(dir, pattern = pat, full.names = TRUE) + 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 Wakhan '", key, "' file matching /", pat, "/ in ", dir) + return(NULL) + } + if (length(hits) > 1) { + log_msg(" multiple Wakhan '", key, "' files; using ", basename(hits[1]), + " (pass --sample to disambiguate)") + } + hits[1] +} + +wakhan_read_hp_bed <- function(file) { + df <- data.table::fread(file, sep = "\t", header = TRUE, data.table = FALSE, + skip = "#chr\tstart\tend", na.strings = c("NA", "nan", "NaN", "")) + names(df) <- sub("^#", "", names(df)) + need <- c("chr", "start", "end", "copynumber_state") + missing <- setdiff(need, names(df)) + if (length(missing)) { + stop("Wakhan BED file ", basename(file), " lacks column(s): ", + paste(missing, collapse = ", ")) + } + data.frame(chr = df$chr, start = df$start, end = df$end, + cn = suppressWarnings(as.numeric(df$copynumber_state)), + stringsAsFactors = FALSE) +} + +wakhan_read_cn <- function(hp1_file, hp2_file) { + hp1 <- wakhan_read_hp_bed(hp1_file) + hp2 <- wakhan_read_hp_bed(hp2_file) + cn <- merge(hp1, hp2, by = c("chr", "start", "end"), suffixes = c("1", "2"), + all = FALSE, sort = FALSE) + if (nrow(cn) == 0) stop("Wakhan HP1/HP2 BED files share no segment coordinates") + if (nrow(cn) < max(nrow(hp1), nrow(hp2))) { + log_msg(sprintf(" kept %d shared Wakhan HP segments from %d/%d rows", + nrow(cn), nrow(hp1), nrow(hp2))) + } + data.frame(chr = cn$chr, start = cn$start, end = cn$end, + copyNumber = cn$cn1 + cn$cn2, + minorAlleleCopyNumber = pmin(cn$cn1, cn$cn2), + stringsAsFactors = FALSE) +} + +parse_wakhan <- function(args) { + dir <- wakhan_dir(args$input, args$sample) + sol <- wakhan_solution_dir(dir) + bed_dir <- file.path(sol$dir, "bed_output") + if (!dir.exists(bed_dir)) stop("No Wakhan bed_output directory under ", sol$dir) + + hp1 <- args$cn_file %||% wakhan_find_file(bed_dir, "hp1", args$sample) + hp2 <- wakhan_find_file(bed_dir, "hp2", args$sample) + log_msg(" Wakhan directory: ", dir) + log_msg(" Wakhan solution: ", basename(sol$dir)) + log_msg(" CN files: ", basename(hp1), " + ", basename(hp2)) + cn <- wakhan_read_cn(hp1, hp2) + + meta <- list(sample = args$sample %||% NA_character_, cn_file = hp1, cn_caller = "Wakhan", + wakhan_solution = basename(sol$dir)) + if (!is.null(sol$ranks)) { + meta$sample <- meta$sample %||% sub("_.*$", "", basename(hp1)) + if ("cell_purity" %in% names(sol$ranks)) meta$purity <- sol$ranks$cell_purity[1] + if ("ploidy" %in% names(sol$ranks)) meta$ploidy <- sol$ranks$ploidy[1] + if ("confidence" %in% names(sol$ranks)) meta$wakhan_confidence <- sol$ranks$confidence[1] + log_msg(sprintf(" purity=%s ploidy=%s", meta$purity %||% NA, meta$ploidy %||% NA)) + } + + list(cn = cn, meta = meta, annotation_fn = NULL) +} + +register_parser("wakhan", parse_wakhan, + "Wakhan allele-specific CN (/bed_output HP BEDs)", + provides = "cn") diff --git a/assets/reconplot/R/plotting.R b/assets/reconplot/R/plotting.R new file mode 100644 index 00000000..88f14be2 --- /dev/null +++ b/assets/reconplot/R/plotting.R @@ -0,0 +1,114 @@ +## --------------------------------------------------------------------------- +## plotting.R -- thin, caller-agnostic wrapper around ReConPlot::ReConPlot(). +## +## Everything above this layer produces the two standard data frames; this file +## only decides sizing, titles, file names and per-panel error handling. +## --------------------------------------------------------------------------- + +suppressPackageStartupMessages({ + library(ggplot2) + library(ReConPlot) +}) + +#' Default canvas size: grows with the number of side-by-side chromosome panels. +auto_plot_size <- function(regions, has_annotation = FALSE) { + n <- nrow(regions) + list(width = min(4 + 2.0 * n, 24), + height = if (has_annotation) 4.6 else 3.4) +} + +build_title <- function(user_title, meta, regions) { + if (!is.null(user_title) && nzchar(user_title)) return(user_title) + bits <- c() + if (!is.null(meta$sample) && !is.na(meta$sample)) bits <- c(bits, meta$sample) + if (!is.null(meta$purity) && !is.null(meta$ploidy)) { + bits <- c(bits, sprintf("purity %.2f / ploidy %.2f", + as.numeric(meta$purity), as.numeric(meta$ploidy))) + } + chrs <- unique(regions$chr) + where <- if (nrow(regions) >= 20 && all(regions$full)) { + "genome-wide" + } else if (length(chrs) > 6) { + paste0(paste(chrs[1:6], collapse = ", "), ", +", length(chrs) - 6, " more") + } else { + paste(chrs, collapse = ", ") + } + bits <- c(bits, where) + paste(bits, collapse = " | ") +} + +#' Render one ReConPlot panel set. +#' +#' @param cn,sv validated data frames +#' @param regions chr_selection data frame +#' @param annotation optional data frame(chr, pos, y) +#' @param extra named list of further arguments forwarded to ReConPlot() +render_reconplot <- function(cn, sv, regions, title = "", genes = NULL, + annotation = NULL, extra = list()) { + ## ReConPlot only reads chr/start/end; drop our bookkeeping columns. + regions <- as.data.frame(regions[, c("chr", "start", "end")], stringsAsFactors = FALSE) + rownames(regions) <- NULL + call_args <- list(sv = sv, cnv = cn, chr_selection = regions, title = title) + if (!is.null(genes) && length(genes)) call_args$genes <- genes + if (!is.null(annotation) && nrow(annotation) > 0) { + call_args$custom_annotation <- annotation + call_args$ann_y_title <- "BAF" + call_args$ann_one_scale <- TRUE + call_args$ann_dot_size <- 0.15 # dense het-SNP track; keep dots small + } + call_args <- utils::modifyList(call_args, extra) + do.call(ReConPlot::ReConPlot, call_args) +} + +save_reconplot <- function(p, outdir, stem, width, height, formats = c("pdf")) { + dir.create(outdir, showWarnings = FALSE, recursive = TRUE) + paths <- character(0) + for (fmt in formats) { + path <- file.path(outdir, paste0(stem, ".", fmt)) + if (fmt == "png") { + ggplot2::ggsave(path, plot = p, width = width, height = height, + units = "in", dpi = 300, limitsize = FALSE) + } else { + ggplot2::ggsave(path, plot = p, width = width, height = height, + units = "in", device = fmt, limitsize = FALSE) + } + paths <- c(paths, path) + } + paths +} + +#' Render and save one plot per element of `region_sets`. +#' +#' @param region_sets list of chr_selection data frames +#' @return character vector of files written +plot_region_sets <- function(region_sets, cn, sv, meta, outdir, prefix, + title = NULL, genes = NULL, annotation_fn = NULL, + formats = c("pdf"), width = NULL, height = NULL, + extra = list()) { + written <- character(0) + for (i in seq_along(region_sets)) { + regions <- region_sets[[i]] + label <- region_label(regions) + stem <- paste(c(prefix, label), collapse = "_") + log_msg(sprintf("plotting %s (%d panel%s)", label, nrow(regions), + if (nrow(regions) == 1) "" else "s")) + + annotation <- if (!is.null(annotation_fn)) annotation_fn(regions) else NULL + size <- auto_plot_size(regions, has_annotation = !is.null(annotation)) + w <- width %||% size$width + h <- height %||% size$height + + p <- tryCatch( + render_reconplot(cn, sv, regions, + title = build_title(title, meta, regions), + genes = genes, annotation = annotation, extra = extra), + error = function(e) { log_msg(" ERROR: ", conditionMessage(e)); NULL }) + if (is.null(p)) next + + out <- tryCatch(save_reconplot(p, outdir, stem, w, h, formats), + error = function(e) { log_msg(" ERROR saving: ", conditionMessage(e)); character(0) }) + for (f in out) log_msg(" wrote ", f) + written <- c(written, out) + } + written +} diff --git a/assets/reconplot/R/utils_common.R b/assets/reconplot/R/utils_common.R new file mode 100644 index 00000000..7cd9cebf --- /dev/null +++ b/assets/reconplot/R/utils_common.R @@ -0,0 +1,267 @@ +## --------------------------------------------------------------------------- +## utils_common.R -- caller-agnostic helpers shared by every ReConPlot parser. +## +## Nothing in this file knows about SAVANA (or any other caller). It defines the +## two data contracts that ReConPlot::ReConPlot() expects, plus the small amount +## of genome bookkeeping needed to turn a region string into a chr_selection +## data frame. +## +## CN contract : chr, start, end, copyNumber, minorAlleleCopyNumber +## SV contract : chr1, pos1, chr2, pos2, strands +## (strands is one of "+-", "-+", "++", "--", "TRA", "INS", "SBE") +## --------------------------------------------------------------------------- + +MAIN_CHROMS <- paste0("chr", c(1:22, "X", "Y")) + +## strand vocabulary understood by ReConPlot's colour mapper. Anything outside +## this set makes the package fail with an uninformative error, so we normalise +## aggressively on the way in. +VALID_STRANDS <- c("+-", "-+", "++", "--", "TRA", "INS", "SBE", + "DEL", "DUP", "h2hINV", "t2tINV") + +CHROM_LENGTHS <- list( + hg38 = c(chr1 = 248956422, chr2 = 242193529, chr3 = 198295559, chr4 = 190214555, + chr5 = 181538259, chr6 = 170805979, chr7 = 159345973, chr8 = 145138636, + chr9 = 138394717, chr10 = 133797422, chr11 = 135086622, chr12 = 133275309, + chr13 = 114364328, chr14 = 107043718, chr15 = 101991189, chr16 = 90338345, + chr17 = 83257441, chr18 = 80373285, chr19 = 58617616, chr20 = 64444167, + chr21 = 46709983, chr22 = 50818468, chrX = 156040895, chrY = 57227415), + hg19 = c(chr1 = 249250621, chr2 = 243199373, chr3 = 198022430, chr4 = 191154276, + chr5 = 180915260, chr6 = 171115067, chr7 = 159138663, chr8 = 146364022, + chr9 = 141213431, chr10 = 135534747, chr11 = 135006516, chr12 = 133851895, + chr13 = 115169878, chr14 = 107349540, chr15 = 102531392, chr16 = 90354753, + chr17 = 81195210, chr18 = 78077248, chr19 = 59128983, chr20 = 63025520, + chr21 = 48129895, chr22 = 51304566, chrX = 155270560, chrY = 59373566) +) + +`%||%` <- function(a, b) if (is.null(a) || length(a) == 0 || all(is.na(a))) b else a + +log_msg <- function(...) { + message(format(Sys.time(), "[%H:%M:%S] "), paste0(..., collapse = "")) +} + +## --- chromosome naming ------------------------------------------------------ + +#' Coerce assorted chromosome spellings to the UCSC "chrN" style ReConPlot wants. +normalize_chrom <- function(x) { + x <- trimws(as.character(x)) + x <- sub("^(chr)?", "chr", x, ignore.case = TRUE) + x <- sub("^chr(chr)+", "chr", x) + x <- sub("^chrMT$", "chrM", x) + x <- sub("^chr23$", "chrX", x) + x <- sub("^chr24$", "chrY", x) + x +} + +is_main_chrom <- function(x) x %in% MAIN_CHROMS + +#' Drop rows touching alt/random/decoy contigs, which ReConPlot rejects outright. +#' +#' @param df data frame +#' @param cols chromosome columns that must all be primary contigs +#' @param what label used in the log line +drop_nonstandard_chroms <- function(df, cols, what = "rows") { + if (nrow(df) == 0) return(df) + keep <- Reduce(`&`, lapply(cols, function(cc) is_main_chrom(df[[cc]]))) + n_drop <- sum(!keep) + if (n_drop > 0) { + bad <- unique(unlist(lapply(cols, function(cc) df[[cc]][!keep]))) + bad <- setdiff(bad, MAIN_CHROMS) + log_msg(sprintf(" dropped %d %s on non-primary contigs (%s%s)", + n_drop, what, paste(utils::head(bad, 5), collapse = ", "), + if (length(bad) > 5) ", ..." else "")) + } + df[keep, , drop = FALSE] +} + +## --- contract validation ---------------------------------------------------- + +validate_cn <- function(cn, drop_na_minor = FALSE) { + req <- c("chr", "start", "end", "copyNumber", "minorAlleleCopyNumber") + missing <- setdiff(req, names(cn)) + if (length(missing)) { + stop("CN table is missing required column(s): ", paste(missing, collapse = ", ")) + } + cn <- as.data.frame(cn, stringsAsFactors = FALSE) + cn$chr <- normalize_chrom(cn$chr) + cn$start <- as.numeric(cn$start) + cn$end <- as.numeric(cn$end) + cn$copyNumber <- suppressWarnings(as.numeric(cn$copyNumber)) + cn$minorAlleleCopyNumber <- suppressWarnings(as.numeric(cn$minorAlleleCopyNumber)) + + cn <- drop_nonstandard_chroms(cn, "chr", "CN segments") + bad <- is.na(cn$start) | is.na(cn$end) | is.na(cn$copyNumber) + if (any(bad)) { + log_msg(sprintf(" dropped %d CN segments with missing coordinates/total CN", sum(bad))) + cn <- cn[!bad, , drop = FALSE] + } + if (drop_na_minor) { + bad <- is.na(cn$minorAlleleCopyNumber) + if (any(bad)) { + log_msg(sprintf(" dropped %d CN segments with missing minor allele CN", sum(bad))) + cn <- cn[!bad, , drop = FALSE] + } + } + ## ReConPlot indexes cnv[, c("chr","start","end","copyNumber","minorAlleleCopyNumber")] + ## positionally in places, so hand it exactly those columns in that order. + cn <- cn[order(cn$chr, cn$start), req, drop = FALSE] + rownames(cn) <- NULL + if (nrow(cn) == 0) stop("No usable copy number segments after filtering.") + cn +} + +validate_sv <- function(sv, interchrom_as_tra = TRUE) { + req <- c("chr1", "pos1", "chr2", "pos2", "strands") + missing <- setdiff(req, names(sv)) + if (length(missing)) { + stop("SV table is missing required column(s): ", paste(missing, collapse = ", ")) + } + sv <- as.data.frame(sv, stringsAsFactors = FALSE) + if (nrow(sv) == 0) return(sv[, union(req, names(sv)), drop = FALSE]) + + sv$chr1 <- normalize_chrom(sv$chr1) + sv$chr2 <- normalize_chrom(sv$chr2) + sv$pos1 <- suppressWarnings(as.integer(round(as.numeric(sv$pos1)))) + sv$pos2 <- suppressWarnings(as.integer(round(as.numeric(sv$pos2)))) + sv$strands <- trimws(as.character(sv$strands)) + + ## single breakends carry no mate; park them on their own locus so the + ## chromosome filter below does not throw them away. + sbe <- sv$strands == "SBE" + if (any(sbe)) { + sv$chr2[sbe] <- sv$chr1[sbe] + sv$pos2[sbe] <- sv$pos1[sbe] + } + + sv <- drop_nonstandard_chroms(sv, c("chr1", "chr2"), "SV junctions") + bad <- is.na(sv$pos1) | is.na(sv$pos2) | is.na(sv$strands) | !(sv$strands %in% VALID_STRANDS) + if (any(bad)) { + log_msg(sprintf(" dropped %d SV junctions with unusable coordinates/orientation", sum(bad))) + sv <- sv[!bad, , drop = FALSE] + } + if (interchrom_as_tra && nrow(sv) > 0) { + tra <- sv$chr1 != sv$chr2 & !(sv$strands %in% c("INS", "SBE")) + if (any(tra)) { + sv$strands[tra] <- "TRA" + log_msg(sprintf(" relabelled %d inter-chromosomal junctions as TRA", sum(tra))) + } + } + rownames(sv) <- NULL + sv +} + +## --- region / chr_selection handling --------------------------------------- + +#' Parse a region specification into a ReConPlot chr_selection data frame. +#' +#' Accepted forms (comma, semicolon or whitespace separated): +#' "all" every primary contig present in the CN table +#' "chr8" whole chromosome +#' "chr8:120000000-130000000" explicit window (commas/underscores allowed) +#' +#' @param spec character scalar, or NULL/"all" +#' @param cn validated CN table, used to bound whole-chromosome requests +#' @param genome one of names(CHROM_LENGTHS); NA lengths fall back to CN extent +parse_regions <- function(spec, cn, genome = "hg38") { + lens <- CHROM_LENGTHS[[genome]] + chrom_end <- function(chr) { + from_ref <- if (!is.null(lens) && chr %in% names(lens)) unname(lens[[chr]]) else NA_real_ + from_cn <- suppressWarnings(max(cn$end[cn$chr == chr], na.rm = TRUE)) + if (!is.finite(from_cn)) from_cn <- NA_real_ + if (is.na(from_ref)) from_cn else from_ref + } + + if (is.null(spec) || length(spec) == 0 || identical(tolower(trimws(spec)), "all")) { + chrs <- MAIN_CHROMS[MAIN_CHROMS %in% unique(cn$chr)] + if (length(chrs) == 0) stop("No primary contigs found in the CN table.") + return(data.frame(chr = chrs, start = 0, + end = vapply(chrs, chrom_end, numeric(1)), + full = TRUE, stringsAsFactors = FALSE)) + } + + tokens <- unlist(strsplit(spec, "[,;[:space:]]+")) + tokens <- tokens[nzchar(tokens)] + out <- lapply(tokens, function(tok) { + parts <- strsplit(tok, ":", fixed = TRUE)[[1]] + chr <- normalize_chrom(parts[1]) + if (length(parts) == 1) { + return(data.frame(chr = chr, start = 0, end = chrom_end(chr), + full = TRUE, stringsAsFactors = FALSE)) + } + rng <- gsub("[,_]", "", parts[2]) + se <- strsplit(rng, "-", fixed = TRUE)[[1]] + if (length(se) != 2) stop("Cannot parse region '", tok, "'. Use chr:start-end.") + data.frame(chr = chr, start = as.numeric(se[1]), end = as.numeric(se[2]), + full = FALSE, stringsAsFactors = FALSE) + }) + do.call(rbind, out) +} + +read_regions_file <- function(path) { + df <- utils::read.table(path, sep = "\t", header = FALSE, stringsAsFactors = FALSE, + comment.char = "#") + if (ncol(df) < 3) stop("Regions BED file needs at least 3 columns: chr, start, end") + data.frame(chr = normalize_chrom(df[[1]]), + start = as.numeric(df[[2]]), + end = as.numeric(df[[3]]), + full = FALSE, stringsAsFactors = FALSE) +} + +#' Keep only regions that ReConPlot can actually draw (primary contig + CN data). +sanitize_regions <- function(regions, cn) { + regions <- as.data.frame(regions, stringsAsFactors = FALSE) + regions$chr <- normalize_chrom(regions$chr) + if (is.null(regions$full)) regions$full <- FALSE + + bad <- !is_main_chrom(regions$chr) + if (any(bad)) { + log_msg(" skipping unsupported contigs in region list: ", + paste(unique(regions$chr[bad]), collapse = ", ")) + regions <- regions[!bad, , drop = FALSE] + } + missing_cn <- !(regions$chr %in% unique(cn$chr)) + if (any(missing_cn)) { + log_msg(" skipping regions without copy number data: ", + paste(unique(regions$chr[missing_cn]), collapse = ", ")) + regions <- regions[!missing_cn, , drop = FALSE] + } + if (nrow(regions) == 0) stop("No plottable regions left after filtering.") + regions$start[is.na(regions$start) | regions$start < 0] <- 0 + regions <- regions[order(match(regions$chr, MAIN_CHROMS), regions$start), , drop = FALSE] + rownames(regions) <- NULL + regions +} + +#' Short, filesystem-safe label for a region set (used in output file names). +region_label <- function(regions) { + ## a full-chromosome sweep gets a name rather than a 24-part file stem + if (nrow(regions) >= 20 && all(regions$full)) return("genome_wide") + parts <- vapply(seq_len(nrow(regions)), function(i) { + r <- regions[i, ] + if (isTRUE(r$full)) r$chr + else sprintf("%s_%.1fMb-%.1fMb", r$chr, r$start / 1e6, r$end / 1e6) + }, character(1)) + if (length(parts) > 6) parts <- c(parts[1:6], sprintf("and%dmore", length(parts) - 6)) + gsub("[^A-Za-z0-9._-]", "_", paste(parts, collapse = "_")) +} + +#' Collapse repeated chromosomes into one spanning window. +#' +#' ReConPlot facets on `factor(chr, levels = unique(chr))`, so a multi-panel +#' selection containing the same chromosome twice fails with +#' "factor level [n] is duplicated". Merging is the useful behaviour: two +#' windows on one chromosome become the interval that spans both. +collapse_duplicate_chroms <- function(regions) { + if (!anyDuplicated(regions$chr)) return(regions) + dup <- unique(regions$chr[duplicated(regions$chr)]) + log_msg(" merging repeated chromosome(s) into a single panel: ", + paste(dup, collapse = ", ")) + parts <- lapply(split(regions, regions$chr), function(g) { + data.frame(chr = g$chr[1], start = min(g$start), end = max(g$end), + full = any(g$full), stringsAsFactors = FALSE) + }) + out <- do.call(rbind, parts) + out <- out[order(match(out$chr, MAIN_CHROMS)), , drop = FALSE] + rownames(out) <- NULL + out +} diff --git a/assets/reconplot/R/utils_vcf.R b/assets/reconplot/R/utils_vcf.R new file mode 100644 index 00000000..a06de0be --- /dev/null +++ b/assets/reconplot/R/utils_vcf.R @@ -0,0 +1,95 @@ +## --------------------------------------------------------------------------- +## utils_vcf.R -- minimal, caller-agnostic VCF reading for SV parsers. +## +## Deliberately not a full VCF library: SV callers only need the fixed columns, +## a few INFO keys and the breakend ALT notation. Shared by savana.R and +## severus.R so both agree on how a BND mate is decoded. +## --------------------------------------------------------------------------- + +#' Read a (optionally bgzipped) VCF into a data frame of the fixed columns. +#' +#' @return data frame with chrom, pos, id, ref, alt, qual, filter, info, +#' format, sample1 (the first sample column, if present) +vcf_read_records <- function(file, pass_only = TRUE) { + con <- if (grepl("\\.gz$", file)) gzfile(file, "rt") else file(file, "rt") + on.exit(close(con), add = TRUE) + lines <- readLines(con, warn = FALSE) + lines <- lines[!startsWith(lines, "##")] + if (length(lines) < 2) return(NULL) + + f <- strsplit(lines[-1], "\t", fixed = TRUE) + ncol_max <- max(lengths(f)) + f <- do.call(rbind, lapply(f, function(x) c(x, rep(NA_character_, ncol_max - length(x))))) + + out <- data.frame( + chrom = f[, 1], pos = suppressWarnings(as.numeric(f[, 2])), id = f[, 3], + ref = f[, 4], alt = f[, 5], qual = f[, 6], filter = f[, 7], info = f[, 8], + format = if (ncol_max >= 9) f[, 9] else NA_character_, + sample1 = if (ncol_max >= 10) f[, 10] else NA_character_, + stringsAsFactors = FALSE) + + if (pass_only) { + keep <- out$filter %in% c("PASS", ".") | is.na(out$filter) + if (any(!keep)) log_msg(sprintf(" dropped %d non-PASS VCF records", sum(!keep))) + out <- out[keep, , drop = FALSE] + } + rownames(out) <- NULL + out +} + +#' Extract one INFO key. Returns NA where the key is absent. +vcf_info_get <- function(info, key) { + rx <- paste0("(^|;)", key, "=([^;]*)") + m <- regexpr(rx, info) + out <- rep(NA_character_, length(info)) + hit <- which(m > 0) + if (length(hit)) { + out[hit] <- sub(paste0("^.*?", key, "="), "", regmatches(info, m), perl = TRUE) + } + out +} + +#' TRUE where a valueless INFO flag (e.g. PRECISE, IMPRECISE) is present. +vcf_info_flag <- function(info, key) { + grepl(paste0("(^|;)", key, "(;|$)"), info) +} + +#' Decode the mate locus from a breakend ALT allele: t[chr:pos[ or ]chr:pos]t +#' +#' @return list(chr, pos); both NA where the ALT is not a breakend +vcf_alt_mate <- function(alt) { + rx <- "[][][^][]+:[0-9]+[][]" + m <- regexpr(rx, alt) + chr <- rep(NA_character_, length(alt)); pos <- rep(NA_real_, length(alt)) + hit <- which(m > 0) + if (length(hit)) { + clean <- gsub("[][]", "", regmatches(alt, m)) + chr[hit] <- sub(":[0-9]+$", "", clean) + pos[hit] <- as.numeric(sub("^.*:", "", clean)) + } + list(chr = chr, pos = pos) +} + +#' Pull one FORMAT subfield (e.g. "DV") out of the sample column. +vcf_format_get <- function(format, sample, key) { + vapply(seq_along(format), function(i) { + if (is.na(format[i]) || is.na(sample[i])) return(NA_character_) + keys <- strsplit(format[i], ":", fixed = TRUE)[[1]] + j <- match(key, keys) + if (is.na(j)) return(NA_character_) + vals <- strsplit(sample[i], ":", fixed = TRUE)[[1]] + if (j > length(vals)) NA_character_ else vals[j] + }, character(1)) +} + +#' Collapse mate-paired breakend records to one row per junction. +#' +#' Keyed on the unordered breakpoint pair so it works whether or not the caller +#' emits MATE_ID. +dedupe_breakend_pairs <- function(sv) { + if (nrow(sv) == 0) return(sv) + a <- paste0(sv$chr1, ":", sv$pos1) + b <- paste0(sv$chr2, ":", sv$pos2) + key <- ifelse(a < b, paste(a, b, sep = "|"), paste(b, a, sep = "|")) + sv[!duplicated(key), , drop = FALSE] +} diff --git a/assets/reconplot/README.md b/assets/reconplot/README.md new file mode 100644 index 00000000..b3ed80a8 --- /dev/null +++ b/assets/reconplot/README.md @@ -0,0 +1,8 @@ +# ReConPlot wrapper + +`run_reconplot.R` and `R/` turn lrsomatic caller output (ASCAT, Wakhan, Severus, SAVANA) into +[ReConPlot](https://github.com/cortes-ciriano-lab/ReConPlot) figures and harmonised CN/SV tables. +The `RECONPLOT` module stages this directory as its wrapper input. + +Vendored from [Tim-Yu/ReConPlot](https://github.com/Tim-Yu/ReConPlot) at the commit recorded in +`VERSION`. To update, copy `run_reconplot.R` and `R/` from that repository and bump `VERSION`. diff --git a/assets/reconplot/VERSION b/assets/reconplot/VERSION new file mode 100644 index 00000000..97614922 --- /dev/null +++ b/assets/reconplot/VERSION @@ -0,0 +1 @@ +2672ccb22fc979909e136f6d2614b5fb763d1de2 diff --git a/assets/reconplot/run_reconplot.R b/assets/reconplot/run_reconplot.R new file mode 100755 index 00000000..28944ccb --- /dev/null +++ b/assets/reconplot/run_reconplot.R @@ -0,0 +1,217 @@ +#!/usr/bin/env Rscript +## --------------------------------------------------------------------------- +## run_reconplot.R -- command line front end for ReConPlot. +## +## Pipeline: --[parser]--> {cn, sv} --[validate]--> ReConPlot +## +## The only caller-specific code lives in R/parsers/. Everything here is +## generic, so supporting a new caller means adding one parser file. +## +## Examples +## ./run_reconplot.R --cn-source savana --sv-source savana \ +## --input /path/to/savana/sample \ +## --outdir /path/to/ReConPlot_output +## +## ./run_reconplot.R --cn-source wakhan --sv-source severus --input DIR --outdir OUT \ +## --regions "chr8,chr17:0-30000000" --layout together --genes MYC,TP53 +## --------------------------------------------------------------------------- + +suppressPackageStartupMessages({ + library(optparse) + library(data.table) +}) + +## --- locate our own installation directory --------------------------------- +script_dir <- function() { + ca <- commandArgs(trailingOnly = FALSE) + f <- sub("^--file=", "", ca[grepl("^--file=", ca)]) + if (length(f)) return(normalizePath(dirname(f[1]))) + normalizePath(".") +} +SCRIPT_DIR <- script_dir() +source(file.path(SCRIPT_DIR, "R", "utils_common.R")) +source(file.path(SCRIPT_DIR, "R", "utils_vcf.R")) +source(file.path(SCRIPT_DIR, "R", "parsers", "registry.R")) +source(file.path(SCRIPT_DIR, "R", "plotting.R")) +load_parsers(file.path(SCRIPT_DIR, "R", "parsers")) + +## --- options ---------------------------------------------------------------- +option_list <- list( + make_option("--source", type = "character", default = "savana", + help = "Single parser for CN+SVs; prefer --cn-source/--sv-source for mixed callers [default %default]"), + make_option("--cn-source", type = "character", default = NULL, + help = "Parser for copy number only (e.g. ascat); overrides --source"), + make_option("--sv-source", type = "character", default = NULL, + help = "Parser for SVs only (e.g. severus); overrides --source"), + make_option("--cn-input", type = "character", default = NULL, + help = "Directory for the CN source [default: --input]"), + make_option("--sv-input", type = "character", default = NULL, + help = "Directory for the SV source [default: --input]"), + make_option("--list-sources", action = "store_true", default = FALSE, + help = "Print the registered parsers and exit"), + make_option("--input", type = "character", default = NULL, + help = "Caller output directory (parser discovers files inside)"), + make_option("--sample", type = "character", default = NULL, + help = "Sample prefix; inferred from file names when omitted"), + make_option("--cn-file", type = "character", default = NULL, + help = "Explicit copy number file (overrides discovery)"), + make_option("--sv-file", type = "character", default = NULL, + help = "Explicit SV file (overrides discovery)"), + make_option("--purity-file", type = "character", default = NULL, + help = "Explicit purity/ploidy file"), + make_option("--sv-format", type = "character", default = "bedpe", + help = "SAVANA SV input: bedpe or vcf [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 = paste("Drop intra-chromosomal SVs shorter than this many bp;", + "translocations are never dropped [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, + help = "Keep only SVs in these Severus cluster IDs (comma-separated)"), + make_option("--clustered-only", action = "store_true", default = FALSE, + help = "Keep only SVs assigned to some Severus cluster"), + make_option("--severus-all", action = "store_true", default = FALSE, + help = "Use severus_all.vcf.gz (germline included) instead of the somatic set"), + make_option("--ascat-cn-mode", type = "character", default = "segments", + help = paste("ASCAT CN table: segments (fitted integer calls),", + "raw (pre-fit segmentation), raw-fractional", + "(pre-fit, unrounded nAraw/nBraw) [default %default]")), + + make_option("--outdir", type = "character", default = NULL, + help = "Output directory (required)"), + make_option("--prefix", type = "character", default = NULL, + help = "Output file name prefix [default: sample name]"), + + make_option("--regions", type = "character", default = "all", + help = "'all', or e.g. 'chr8,chr17:0-30000000' [default %default]"), + make_option("--regions-file", type = "character", default = NULL, + help = "BED file of regions (overrides --regions)"), + make_option("--layout", type = "character", default = "separate", + help = "separate | together | both [default %default]"), + + make_option("--genes", type = "character", default = NULL, + help = "Comma-separated HUGO gene symbols to annotate"), + make_option("--title", type = "character", default = NULL, + help = "Plot title [default: sample, purity/ploidy, region]"), + make_option("--genome", type = "character", default = "hg38", + help = "Genome build: hg38, hg19, T2T, mm10, mm39 [default %default]"), + make_option("--max-cn", type = "double", default = 8, + help = "Copy number axis ceiling [default %default]"), + make_option("--format", type = "character", default = "pdf,png", + help = "Comma-separated output formats [default %default]"), + make_option("--width", type = "double", default = NULL, help = "Figure width (in)"), + make_option("--height", type = "double", default = NULL, help = "Figure height (in)"), + + make_option("--baf-track", action = "store_true", default = FALSE, + help = "Add a het-SNP BAF annotation panel (SAVANA only)"), + make_option("--annotation-file", type = "character", default = NULL, + help = "Explicit annotation source (het-SNP BED for --baf-track)"), + make_option("--baf-max-points", type = "integer", default = 5000, + help = "Het-SNPs drawn per panel after thinning [default %default]"), + make_option("--extra", type = "character", default = NULL, + help = paste("Escape hatch for any other ReConPlot() argument, e.g.", + "'size_text=6,curvature_intrachr_SVs=-0.2'")), + make_option("--drop-na-minor", action = "store_true", default = FALSE, + help = "Drop CN segments lacking a minor allele CN"), + make_option("--keep-interchrom-strands", action = "store_true", default = FALSE, + help = "Colour inter-chromosomal SVs by orientation instead of as TRA"), + make_option("--write-tables", action = "store_true", default = FALSE, + help = "Also write the harmonised CN/SV tables as TSV"), + make_option("--seed", type = "integer", default = 1, + help = "RNG seed (BAF thinning) [default %default]") +) + +opt <- parse_args(OptionParser( + usage = "%prog --source SOURCE --input DIR --outdir DIR [options]", + option_list = option_list)) + +## optparse turns --cn-file into opt$`cn-file`; give everything snake_case names +## so parsers can use args$cn_file. +args <- opt +names(args) <- gsub("-", "_", names(args)) + +if (isTRUE(args$list_sources)) { + print(list_parsers(), right = FALSE) + quit(status = 0) +} +if (is.null(args$outdir)) stop("--outdir is required") +if (is.null(args$input) && is.null(args$cn_input) && is.null(args$cn_file)) { + stop("Provide --input (a caller output directory) or explicit --cn-file/--sv-file") +} +set.seed(args$seed) + +## --- 1. parse --------------------------------------------------------------- +parsed <- run_parsers(args) +meta <- parsed$meta %||% list() + +## --- 2. harmonise / validate ----------------------------------------------- +log_msg("validating inputs") +cn <- validate_cn(parsed$cn, drop_na_minor = isTRUE(args$drop_na_minor)) +sv <- validate_sv(parsed$sv, interchrom_as_tra = !isTRUE(args$keep_interchrom_strands)) +log_msg(sprintf(" %d CN segments, %d SV junctions retained", nrow(cn), nrow(sv))) +if (nrow(sv) > 0) { + tab <- table(sv$strands) + log_msg(" SV types: ", paste(sprintf("%s=%d", names(tab), as.integer(tab)), collapse = " ")) +} + +dir.create(args$outdir, showWarnings = FALSE, recursive = TRUE) +prefix <- args$prefix %||% meta$sample %||% args$source +if (isTRUE(args$write_tables)) { + data.table::fwrite(cn, file.path(args$outdir, paste0(prefix, ".reconplot_cn.tsv")), sep = "\t") + data.table::fwrite(sv, file.path(args$outdir, paste0(prefix, ".reconplot_sv.tsv")), sep = "\t") + log_msg(" wrote harmonised tables to ", args$outdir) +} + +## --- 3. regions ------------------------------------------------------------- +regions <- if (!is.null(args$regions_file)) { + read_regions_file(args$regions_file) +} else { + parse_regions(args$regions, cn, genome = args$genome) +} +regions <- sanitize_regions(regions, cn) +log_msg(sprintf("%d region(s) to plot", nrow(regions))) + +combined <- collapse_duplicate_chroms(regions) +region_sets <- switch( + args$layout, + separate = split(regions, seq_len(nrow(regions))), + together = list(combined), + both = c(list(combined), split(regions, seq_len(nrow(regions)))), + stop("--layout must be one of: separate, together, both") +) +region_sets <- lapply(region_sets, function(x) { rownames(x) <- NULL; x }) + +## --- 4. plot ---------------------------------------------------------------- +genes <- if (!is.null(args$genes)) trimws(unlist(strsplit(args$genes, ","))) else NULL +formats <- trimws(unlist(strsplit(args$format, ","))) + +## --extra lets any remaining ReConPlot() argument through without a dedicated flag +extra <- list(max.cn = args$max_cn, genome_version = args$genome) +if (!is.null(args$extra)) { + for (kv in trimws(unlist(strsplit(args$extra, ",")))) { + if (!nzchar(kv)) next + parts <- strsplit(kv, "=", fixed = TRUE)[[1]] + if (length(parts) != 2) stop("--extra entries must look like key=value: ", kv) + val <- suppressWarnings(as.numeric(parts[2])) + if (is.na(val)) { + val <- switch(parts[2], "TRUE" = TRUE, "FALSE" = TRUE, parts[2]) + if (parts[2] == "FALSE") val <- FALSE + } + extra[[trimws(parts[1])]] <- val + } + log_msg(" extra ReConPlot args: ", paste(names(extra), unlist(extra), sep = "=", collapse = " ")) +} + +written <- plot_region_sets( + region_sets, cn = cn, sv = sv, meta = meta, + outdir = args$outdir, prefix = prefix, + title = args$title, genes = genes, + annotation_fn = parsed$annotation_fn, + formats = formats, width = args$width, height = args$height, + extra = extra +) + +log_msg(sprintf("done: %d file(s) written to %s", length(written), args$outdir)) +if (length(written) == 0) quit(status = 1) diff --git a/conf/modules.config b/conf/modules.config index 8b20e4fd..0acbfc14 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -727,6 +727,75 @@ process { ] } + // + // Annotation of SVs / CNAs + // + + withName : '.*:PADFOOT_WGET' { + ext.suffix = 'tar.gz' + publishDir = [ enabled: false ] + } + + withName : '.*:PADFOOT_UNTAR' { + publishDir = [ enabled: false ] + } + + withName : '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { + ext.args = { + [ + params.padfoot_run_repeatmasker ? '' : '--skip_RepeatMasker', + genome == 'mm10' ? '--specie mouse' : '' // RepeatMasker -species; Padfoot defaults to human + ].join(' ').trim() + } + ext.prefix = { sv_caller == cna_caller ? "${sv_caller}" : "${sv_caller}_${cna_caller}" } + publishDir = [ + path: { "${params.outdir}/${meta.id}/padfoot" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + // + // Visualisation of SVs / CNAs + // + + withName : '.*:RECONPLOT_PKG_WGET' { + ext.suffix = 'tar.gz' + publishDir = [ enabled: false ] + } + + withName : '.*:RECONPLOT_PKG_UNTAR' { + publishDir = [ enabled: false ] + } + + withName : '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + // shared filters / rendering + ext.args = { + [ + "--max-cn ${params.reconplot_max_cn}", + "--format ${params.reconplot_format}", + params.reconplot_min_svlen ? "--min-svlen ${params.reconplot_min_svlen}" : '', + (params.reconplot_exclude_vntr && sv_source == 'severus') ? '--exclude-vntr' : '' + ].join(' ').trim() + } + // genome-wide strip + ext.args2 = '--width 34 --height 4 --extra size_text=4,size_chr_labels=5,scale_ticks=50000000' + // optional focus panel; empty => not produced + ext.args3 = { + params.reconplot_regions ? [ + "--regions '${params.reconplot_regions}'", + params.reconplot_genes ? "--genes ${params.reconplot_genes}" : '', + (params.reconplot_baf_track && cn_source in ['ascat', 'savana']) ? '--baf-track' : '' + ].join(' ').trim() : '' + } + publishDir = [ + path: { "${params.outdir}/${meta.id}/reconplot" }, + 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 1829a02e..1baaf6dd 100644 --- a/conf/test.config +++ b/conf/test.config @@ -70,6 +70,9 @@ params { skip_wakhan = true skip_ascat = true skip_modkit = true + // Padfoot/ReConPlot need whole-genome caller output; covered by the module- and subworkflow-level nf-tests + skip_padfoot = true + skip_reconplot = 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 d0bf023c..96617892 100644 --- a/conf/test_full.config +++ b/conf/test_full.config @@ -36,6 +36,9 @@ params { skip_wakhan = true skip_ascat = true skip_vep_plugins = true + // Padfoot/ReConPlot need whole-genome caller output; covered by the module- and subworkflow-level nf-tests + skip_padfoot = true + skip_reconplot = true // Both callers, union combine mode germline_var_keep = 'clair, deepvariant' diff --git a/containers/padfoot/Dockerfile b/containers/padfoot/Dockerfile new file mode 100644 index 00000000..1169baa2 --- /dev/null +++ b/containers/padfoot/Dockerfile @@ -0,0 +1,30 @@ +FROM mambaorg/micromamba:2.3.2 + +COPY --chown=$MAMBA_USER modules/local/padfoot/environment.yml /tmp/environment.yml + +ARG DFAM_URL_BASE=https://www.dfam.org/releases/current/families/FamDB +ARG DFAM_ROOT_SHA256=4cdb3d571126de3061d61c9e79d0ffcf0b1d9dc1d3c555bf9bcb8db695efa91a +ARG DFAM_CONSENSUS_SHA256=ce50a75fe90a5481864cd87792a80564667a4e7f373322730a42a7b737b8bbdd + +RUN micromamba create --yes --name padfoot --file /tmp/environment.yml \ + && micromamba clean --all --yes \ + && export PATH=/opt/conda/envs/padfoot/bin:${PATH} \ + && mkdir -p /home/mambauser/dfam /tmp/repeatmasker-smoke \ + && wget --quiet --tries=3 --output-document=/tmp/dfam40.0.h5.gz "${DFAM_URL_BASE}/dfam40.0.h5.gz" \ + && wget --quiet --tries=3 --output-document=/tmp/dfam40.curated.consensus.0.h5.gz "${DFAM_URL_BASE}/dfam40.curated.consensus.0.h5.gz" \ + && gunzip /tmp/dfam40.0.h5.gz /tmp/dfam40.curated.consensus.0.h5.gz \ + && echo "${DFAM_ROOT_SHA256} /tmp/dfam40.0.h5" | sha256sum --check \ + && echo "${DFAM_CONSENSUS_SHA256} /tmp/dfam40.curated.consensus.0.h5" | sha256sum --check \ + && mv /tmp/dfam40.0.h5 /tmp/dfam40.curated.consensus.0.h5 /home/mambauser/dfam/ \ + && FAMDB_CONFIG="$(find /opt/conda/envs/padfoot/share -path '*/famdb-*/famdb.conf' -print -quit)" \ + && test -n "${FAMDB_CONFIG}" \ + && printf '[famdb]\nFAMDB_DATA_DIR = /home/mambauser/dfam\n' > "${FAMDB_CONFIG}" \ + && export PATH=/opt/conda/envs/padfoot/bin:${PATH} FAMDB_DATA_DIR=/home/mambauser/dfam \ + && famdb.py info \ + && printf '>smoke\nACGTACGTACGTACGT\n' > /tmp/repeatmasker-smoke/input.fa \ + && RepeatMasker -species human -dir /tmp/repeatmasker-smoke /tmp/repeatmasker-smoke/input.fa \ + && test -s /tmp/repeatmasker-smoke/input.fa.out \ + && rm -rf /tmp/dfam40.* /tmp/repeatmasker-smoke + +ENV PATH=/opt/conda/envs/padfoot/bin:${PATH} +ENV FAMDB_DATA_DIR=/home/mambauser/dfam diff --git a/containers/padfoot/README.md b/containers/padfoot/README.md new file mode 100644 index 00000000..3cb9aa6f --- /dev/null +++ b/containers/padfoot/README.md @@ -0,0 +1,21 @@ +# Padfoot RepeatMasker image + +This is the image the `PADFOOT` module runs in (pinned by digest in `modules/local/padfoot/main.nf`). It embeds the Dfam 4.0 root and curated-consensus FamDB partitions, configures FamDB at build time, and validates human RepeatMasker during the build. + +Build and publish from the pipeline root: + +```bash +export PADFOOT_IMAGE=/padfoot-repeatmasker:4.2.4-dfam4 +docker build -f containers/padfoot/Dockerfile -t "$PADFOOT_IMAGE" . +docker push "$PADFOOT_IMAGE" +``` + +Resolve the pushed digest and pin it in the `container` directive of `modules/local/padfoot/main.nf`, or override per site: + +```groovy +process { withName: '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { container = '/padfoot-repeatmasker@sha256:' } } +``` + +The Dockerfile verifies the decompressed Dfam partition checksums. A changed Dfam `current` release therefore fails the build rather than silently changing the annotation database. + +The module currently pins `ghcr.io/tim-yu/padfoot-repeatmasker@sha256:f98b0d352ec47cd9fa015f321959dcba7f6c6cd4da05beaea8d811b97dff70f5`; RepeatMasker runs by default (`--padfoot_run_repeatmasker`). diff --git a/containers/reconplot/Dockerfile b/containers/reconplot/Dockerfile new file mode 100644 index 00000000..9c9286c6 --- /dev/null +++ b/containers/reconplot/Dockerfile @@ -0,0 +1,19 @@ +FROM mambaorg/micromamba:2.3.2 + +# R runtime for the ReConPlot wrapper, with the ReConPlot package (not on conda) pre-installed +# from a pinned upstream commit. The wrapper itself is staged by the pipeline (params.reconplot_url). +COPY --chown=$MAMBA_USER modules/local/reconplot/environment.yml /tmp/environment.yml + +ARG RECONPLOT_COMMIT=88bb76011703d7908aaa2e5bdecfe18ba5737db8 +ADD --chown=$MAMBA_USER https://github.com/cortes-ciriano-lab/ReConPlot/archive/${RECONPLOT_COMMIT}.tar.gz /tmp/reconplot.tar.gz + +RUN micromamba create --yes --name reconplot --file /tmp/environment.yml \ + && micromamba clean --all --yes \ + && export PATH=/opt/conda/envs/reconplot/bin:${PATH} \ + && mkdir -p /tmp/reconplot-src \ + && tar -xzf /tmp/reconplot.tar.gz -C /tmp/reconplot-src --strip-components 1 \ + && R CMD INSTALL --no-docs --no-html /tmp/reconplot-src \ + && Rscript -e 'suppressPackageStartupMessages(library(ReConPlot)); cat("ReConPlot", as.character(packageVersion("ReConPlot")), "\n")' \ + && rm -rf /tmp/reconplot.tar.gz /tmp/reconplot-src + +ENV PATH=/opt/conda/envs/reconplot/bin:${PATH} diff --git a/containers/reconplot/README.md b/containers/reconplot/README.md new file mode 100644 index 00000000..eaec0a6c --- /dev/null +++ b/containers/reconplot/README.md @@ -0,0 +1,21 @@ +# ReConPlot image + +R runtime for the [ReConPlot wrapper](https://github.com/Tim-Yu/ReConPlot) used by the +`RECONPLOT` module, with the upstream [ReConPlot](https://github.com/cortes-ciriano-lab/ReConPlot) +R package (not distributed on conda) installed from a pinned commit. The wrapper scripts live in +`assets/reconplot/` and are staged by the pipeline. + +Build and publish from the pipeline root: + +```bash +export RECONPLOT_IMAGE=/reconplot:0.2-r4.4 +docker build -f containers/reconplot/Dockerfile -t "$RECONPLOT_IMAGE" . +docker push "$RECONPLOT_IMAGE" +``` + +Pin the pushed digest in the `container` directive of `modules/local/reconplot/main.nf` (or override per site via +`process { withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = ... } }`). The module currently pins +`ghcr.io/tim-yu/reconplot@sha256:1145fc5aebe0227bec371f4c59b08b9a09871498e403c01b83f83973149ae9e7`. + +Under `-profile conda` the module builds `modules/local/reconplot/environment.yml` and installs +ReConPlot at run time from the source tree staged via `--reconplot_pkg_url` / `--reconplot_pkg_dir`. diff --git a/docs/output.md b/docs/output.md index 495ab48c..47ef5ae8 100644 --- a/docs/output.md +++ b/docs/output.md @@ -19,6 +19,9 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ ├── methylation │ │ └── tumor │ │ └── modkit_pileup +│ ├── padfoot +│ │ ├── severus_wakhan +│ │ └── savana │ ├── qc │ │ ├── tumor │ │ │ ├── cramino_aln @@ -29,6 +32,9 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ │ ├── nanoplot_ubam_rep1 │ │ │ └── samtools │ │ └── whatshap_stats +│ ├── reconplot +│ │ ├── severus_wakhan +│ │ └── savana │ ├── signatures │ │ ├── assignment │ │ └── matrices @@ -37,6 +43,7 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ ├── deepsomatic │ │ ├── deepvariant │ │ ├── phased +│ │ ├── savana │ │ └── severus │ ├── vep │ │ ├── somatic @@ -56,6 +63,9 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ │ └── modkit_pileup │ │ └── normal │ │ └── modkit_pileup +│ ├── padfoot +│ │ ├── severus_wakhan +│ │ └── savana │ ├── qc │ │ ├── tumor │ │ │ ├── cramino_aln @@ -74,6 +84,10 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ │ ├── nanoplot_ubam_rep1 │ │ │ └── samtools │ │ └── whatshap_stats +│ ├── reconplot +│ │ ├── severus_ascat +│ │ ├── severus_wakhan +│ │ └── savana │ ├── signatures │ │ ├── assignment │ │ └── matrices @@ -83,6 +97,7 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ ├── deepsomatic │ │ ├── deepvariant │ │ ├── phased +│ │ ├── savana │ │ └── severus │ ├── vep │ │ ├── germline @@ -94,6 +109,8 @@ 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` and a matched normal (ASCAT is not run for tumour-only samples), 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`
@@ -594,6 +611,66 @@ Mutational signature analysis of the PASS SNVs and indels in the phased somatic
+### `padfoot` + +
+Output files + +``` +├── padfoot +│ ├── severus_wakhan +│ │ ├── annotated_svs.tsv +│ │ ├── by_gene.tsv +│ │ └── padfoot.log +│ └── savana +│ ├── annotated_svs.tsv +│ ├── by_gene.tsv +│ └── padfoot.log +``` + +| File | Description | +| ------------------- | ---------------------------------------------------------------------------------------------------------------------------------------- | +| `annotated_svs.tsv` | One row per somatic SV: breakpoints, support/VAF, overlapping genes and exons per breakend, repeat annotation, microhomology, VNTR, type | +| `by_gene.tsv` | One row per gene: SV and copy-number impact per haplotype | +| `padfoot.log` | Padfoot log | + +`severus_wakhan/` combines Severus somatic SVs with the top-ranked Wakhan copy-number solution; `savana/` combines SAVANA classified somatic SVs with SAVANA absolute copy number (only present when SAVANA CNA was produced). + +
+ +### `reconplot` + +
+Output files + +``` +├── reconplot +│ ├── severus_ascat +│ │ ├── per_chromosome/sample_chr{1..22,X,Y}.{pdf,png} +│ │ ├── genome_wide/sample_genome_wide.{pdf,png} +│ │ ├── focus/sample_.{pdf,png} +│ │ ├── sample.reconplot_cn.tsv +│ │ ├── sample.reconplot_sv.tsv +│ │ └── reconplot.log +│ ├── severus_wakhan +│ │ └── (same layout) +│ └── savana +│ └── (same layout) +``` + +| File | Description | +| ------------------------- | ----------------------------------------------------------------------------------------------------- | +| `per_chromosome/*` | One ReConPlot figure per chromosome: copy number (total + minor allele) with SV arcs coloured by type | +| `genome_wide/*` | All chromosomes side by side in one strip | +| `focus/*` | Multi-panel figure for `--reconplot_regions`, with gene labels / BAF track if requested (optional) | +| `sample.reconplot_cn.tsv` | Harmonised CN table (`chr,start,end,copyNumber,minorAlleleCopyNumber`) as passed to ReConPlot | +| `sample.reconplot_sv.tsv` | Harmonised SV table (`chr1,pos1,chr2,pos2,strands`) as passed to ReConPlot | +| `reconplot.log` | Wrapper log (parser choices, purity/ploidy read, filters applied) | + +`severus_ascat/` and `severus_wakhan/` pair Severus somatic SVs with ASCAT or the top-ranked Wakhan copy-number solution; `savana/` uses SAVANA's own SVs and absolute copy number. A pair is only produced when both callers ran for the sample. + +
+ ### `wakhan`
diff --git a/docs/usage.md b/docs/usage.md index c9fe4d26..d6697d5e 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -169,6 +169,8 @@ If the loci cannot belong to the reference, ClairS-TO disables Verdict with a wa | `--skip_mosdepth` | A boolean to skip `mosdepth`. Default = `false` | | `--skip_ascat` | A boolean to skip `ascat`. ClairS-TO's Verdict germline tagging then falls back to Verdict's own purity and copy number estimate, which is still up to 0.14 from ASCAT's on the samples it was measured on — see [Verdict tags](output.md#clairs-to). Default = `false` | | `--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_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` | @@ -378,6 +380,64 @@ Running with neither, and without `--skip_signatures`, stops the pipeline at sta Both tools run from `ghcr.io/ljwharbers/sigprofiler`, which adds CHM13 support not yet in a SigProfiler release: SigProfilerMatrixGenerator from the branch behind [SigProfilerSuite/SigProfilerMatrixGenerator#250](https://github.com/SigProfilerSuite/SigProfilerMatrixGenerator/pull/250) and SigProfilerAssignment from [ljwharbers/SigProfilerAssignment](https://github.com/ljwharbers/SigProfilerAssignment/tree/chm13-t2t-support), with COSMIC SBS/DBS signatures renormalised to CHM13. ID83 signatures always use the GRCh37 set. Conda is not supported for this step. +#### Padfoot Options + +[Padfoot](https://github.com/KolmogorovLab/Padfoot) annotates somatic SVs and CNAs with gene/exon overlap, repeat context and complex-SV grouping. It is run for every available SV/CNA caller pair, for paired and tumour-only samples alike: + +- `padfoot/severus_wakhan/` -- Severus somatic SVs + the top-ranked (`solution_1`) Wakhan integer copy-number VCF (requires Wakhan not skipped) +- `padfoot/savana/` -- SAVANA classified somatic SVs + SAVANA segmented absolute copy number (requires SAVANA CNA, i.e. an SNP source: the phased germline VCF for paired samples, or the bundled 1000G panel for tumour-only samples on GRCh38/CHM13). Samples without SAVANA CNA are silently skipped. + +Padfoot is not distributed on bioconda. The pipeline downloads the source tree from `--padfoot_url` (GitHub archive, once per run; the default is a pinned commit of the [Tim-Yu/Padfoot](https://github.com/Tim-Yu/Padfoot) fork, which adds SAVANA input support not yet in `KolmogorovLab/Padfoot`) and runs it inside a container / conda environment that provides its dependencies (python, pysam, pandas, biopython, samtools, minimap2, bedtools). On systems without internet access on compute nodes, clone Padfoot once and pass the checkout with `--padfoot_dir`. + +RepeatMasker (used only to classify the sequence of novel insertions) runs by default. Under Docker/Singularity/Apptainer the module image (`ghcr.io/tim-yu/padfoot-repeatmasker`, recipe in `containers/padfoot/`) ships the full Dfam 4.0 database. Under `-profile conda` it uses the small curated Dfam subset bundled with bioconda RepeatMasker (sufficient for common human repeats such as Alu/L1/SVA). Use `--padfoot_run_repeatmasker false` to disable it; all other Padfoot annotations are unaffected. + +Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. For other genomes (e.g. CHM13) provide `--padfoot_gff` and `--padfoot_rm`, otherwise Padfoot is skipped with a warning. + +| Parameter | Description | +| ---------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `--padfoot_url` | URL of a Padfoot source tarball (GitHub archive). Default = pinned commit of `Tim-Yu/Padfoot` with Savana support | +| `--padfoot_dir` | Local Padfoot checkout (directory with `padfoot.py` and `beds/`); overrides `--padfoot_url`. Default = `null` | +| `--padfoot_genome` | Padfoot genome preset (`hg38`, `chm13`, `mm10`). Default = `null` (inferred from `--genome`) | +| `--padfoot_gff` | Custom GFF3 gene annotation. Default = `null` (bundled) | +| `--padfoot_rm` | Custom RepeatMasker annotation. Default = `null` (bundled) | +| `--padfoot_run_repeatmasker` | Run RepeatMasker on inserted sequences (repeat class of novel insertions). Containers ship the full Dfam 4.0 database; conda uses the curated Dfam subset bundled with bioconda RepeatMasker. Default = `true` | + +#### ReConPlot Options + +[ReConPlot](https://github.com/cortes-ciriano-lab/ReConPlot) rearrangement + copy-number figures are generated through the wrapper shipped in `assets/reconplot/` (vendored from [Tim-Yu/ReConPlot](https://github.com/Tim-Yu/ReConPlot)) for every CN/SV caller pair available for a sample, into `reconplot//`: + +- `severus_ascat/` -- ASCAT allele-specific CN + Severus somatic SVs +- `severus_wakhan/` -- Wakhan top-ranked solution CN + Severus somatic SVs +- `savana/` -- SAVANA absolute CN + SAVANA classified somatic SVs + +Each pair produces `per_chromosome/` (one figure per chromosome), `genome_wide/` (all chromosomes in one strip), the harmonised CN/SV tables, and, when `--reconplot_regions` is set, a `focus/` multi-panel figure with optional gene labels and BAF track. Both the wrapper and the ReConPlot R package (neither on conda) are staged as source from GitHub (or local checkouts for offline systems); the default container ships the package pre-installed, while `-profile conda` installs it at run time. + +| Parameter | Description | +| -------------------------- | ------------------------------------------------------------------------------------------------------------- | +| `--reconplot_pkg_url` | ReConPlot R package source tarball. Default = pinned commit of `cortes-ciriano-lab/ReConPlot` | +| `--reconplot_pkg_dir` | Local ReConPlot package checkout; overrides `--reconplot_pkg_url`. Default = `null` | +| `--reconplot_genome` | ReConPlot genome preset (`hg38`, `hg19`, `T2T`, `mm10`, `mm39`). Default = `null` (inferred from `--genome`) | +| `--reconplot_max_cn` | Copy-number axis ceiling. Default = `8` | +| `--reconplot_min_svlen` | Drop intra-chromosomal SVs shorter than this (bp); translocations kept. Default = `0` | +| `--reconplot_exclude_vntr` | Drop Severus SVs flagged inside a VNTR. Default = `false` | +| `--reconplot_regions` | Regions for an extra `focus/` panel, e.g. `"chr8,chr17:30000000-50000000"`. Default = `null` (no focus panel) | +| `--reconplot_genes` | Comma-separated HUGO symbols labelled on the focus panel. Default = `null` | +| `--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` | + +##### Offline and air-gapped systems + +Both images are pinned by digest directly in the module `container` directives, 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: + +```groovy +process { + withName: '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/padfoot-repeatmasker.sif' } + withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/reconplot.sif' } +} +``` + +The Padfoot source tree and the ReConPlot R package are downloaded from GitHub at run time; use `--padfoot_dir` and `--reconplot_pkg_dir` to point at local checkouts instead. + #### Variant Filtering and Combining Options These options control how variants from multiple callers are filtered and merged. @@ -597,8 +657,7 @@ Two of these predictors get there anyway, because they score _proteins_ rather t AlphaMissense covers the same class of variant and is CC BY 4.0. - **SpliceAI** — not currently wired up on either assembly. -> [!IMPORTANT] -> **REVEL is enabled by default and is free for non-commercial use only**; CADD and EVE, if you +> [!IMPORTANT] > **REVEL is enabled by default and is free for non-commercial use only**; CADD and EVE, if you > enable them, are the same. The pipeline cannot accept those terms on your behalf: if your work is > commercial, pass `--skip_vep_plugins`, or set only the resources you are licensed for. diff --git a/modules/local/padfoot/environment.yml b/modules/local/padfoot/environment.yml new file mode 100644 index 00000000..1fcba5d0 --- /dev/null +++ b/modules/local/padfoot/environment.yml @@ -0,0 +1,18 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +# Runtime dependencies for Padfoot (https://github.com/KolmogorovLab/Padfoot); Padfoot itself is not on +# bioconda and is staged as source (see params.padfoot_url / params.padfoot_dir). +channels: + - conda-forge + - bioconda +dependencies: + - "conda-forge::python=3.12" + - "conda-forge::numpy=2.2.6" + - "conda-forge::pandas=2.3.1" + - "bioconda::pysam=0.23.3" + - "conda-forge::biopython=1.85" + - "bioconda::samtools=1.22.1" + - "bioconda::minimap2=2.30" + - "bioconda::bedtools=2.31.1" + - "bioconda::repeatmasker=4.2.4" + - "conda-forge::procps-ng" # `ps` required by Nextflow task metrics diff --git a/modules/local/padfoot/main.nf b/modules/local/padfoot/main.nf new file mode 100644 index 00000000..ea5d1216 --- /dev/null +++ b/modules/local/padfoot/main.nf @@ -0,0 +1,71 @@ +process PADFOOT { + tag "${meta.id}:${sv_caller}+${cna_caller}" + label 'process_medium' + + conda "${moduleDir}/environment.yml" + // Built from containers/padfoot/Dockerfile: Padfoot deps + RepeatMasker 4.2.4 + Dfam 4.0 (Padfoot itself is staged as source). + // Override per site with `process { withName: '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { container = ... } }`. + container "ghcr.io/tim-yu/padfoot-repeatmasker@sha256:f98b0d352ec47cd9fa015f321959dcba7f6c6cd4da05beaea8d811b97dff70f5" + + input: + tuple val(meta), path(sv_vcf), val(sv_caller), path(cna_file), val(cna_caller) + tuple val(meta2), path(fasta) + tuple val(meta3), path(fai) + tuple val(meta4), path(padfoot_src) // Padfoot source tree (contains padfoot.py + beds/) + tuple val(meta5), val(genome), path(gff), path(rm) // gff/rm may be [] -> Padfoot bundled annotations for `genome` + + output: + tuple val(meta), path("${prefix}/annotated_svs.tsv"), emit: annotated_svs + tuple val(meta), path("${prefix}/by_gene.tsv") , emit: by_gene + tuple val(meta), path("${prefix}/padfoot.log") , emit: log + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: "${sv_caller}_${cna_caller}" + def gff_arg = gff ? "--gff ${gff}" : '' + def rm_arg = rm ? "--rm ${rm}" : '' + + """ + python3 ${padfoot_src}/padfoot.py \\ + --sv-vcf ${sv_vcf} \\ + --sv-caller ${sv_caller} \\ + --cna-file ${cna_file} \\ + --cna-caller ${cna_caller} \\ + --ref ${fasta} \\ + --genome ${genome} \\ + ${gff_arg} \\ + ${rm_arg} \\ + --threads ${task.cpus} \\ + --out-dir ${prefix} \\ + ${args} + + rm -rf ${prefix}/temp + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + padfoot: \$(python3 ${padfoot_src}/padfoot.py --version 2>&1 | tail -1) + minimap2: \$(minimap2 --version 2>&1) + samtools: \$(samtools --version | head -1 | sed 's/samtools //') + repeatmasker: \$(command -v RepeatMasker >/dev/null && RepeatMasker -v 2>&1 | sed -n 's/^RepeatMasker version //p' || echo 'not available') + END_VERSIONS + """ + + stub: + prefix = task.ext.prefix ?: "${sv_caller}_${cna_caller}" + """ + mkdir -p ${prefix} + touch ${prefix}/annotated_svs.tsv ${prefix}/by_gene.tsv ${prefix}/padfoot.log + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + padfoot: stub + minimap2: stub + samtools: stub + repeatmasker: stub + END_VERSIONS + """ +} diff --git a/modules/local/padfoot/meta.yml b/modules/local/padfoot/meta.yml new file mode 100644 index 00000000..b183e2f6 --- /dev/null +++ b/modules/local/padfoot/meta.yml @@ -0,0 +1,107 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "padfoot" +description: Functional annotation of somatic structural variants and copy number alterations (gene/exon overlap, repeats, microhomology, complex SV grouping) from Severus/Wakhan or SAVANA output. RepeatMasker insertion annotation requires the dedicated Dfam image. +keywords: + - structural-variation + - copy-number + - annotation + - somatic + - long-read +tools: + - "padfoot": + description: "SV and CNA annotation tool for Severus, Wakhan and Savana output" + homepage: "https://github.com/KolmogorovLab/Padfoot" + documentation: "https://github.com/KolmogorovLab/Padfoot" + tool_dev_url: "https://github.com/Tim-Yu/Padfoot" + licence: ["BSD-3-Clause"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - sv_vcf: + type: file + description: Somatic SV VCF (Severus `severus_somatic.vcf(.gz)` or SAVANA `*.classified.somatic.vcf(.gz)`) + pattern: "*.vcf{,.gz}" + - sv_caller: + type: string + description: SV input format, `severus` or `savana` + - cna_file: + type: file + description: CNA file (Wakhan `*_wakhan_cna_integers.vcf` or SAVANA `*_segmented_absolute_copy_number.tsv`) + pattern: "*.{vcf,tsv}" + - cna_caller: + type: string + description: CNA input format, `wakhan` or `savana` + - - meta2: + type: map + description: Groovy Map for the reference + - fasta: + type: file + description: Reference genome FASTA + pattern: "*.{fa,fasta}" + - - meta3: + type: map + description: Groovy Map for the reference index + - fai: + type: file + description: Reference genome FASTA index (must sit next to the FASTA) + pattern: "*.fai" + - - meta4: + type: map + description: Groovy Map for the Padfoot source + - padfoot_src: + type: directory + description: Padfoot source tree containing `padfoot.py`, `padfoot/` and `beds/` (Padfoot is not on bioconda) + - - meta5: + type: map + description: Groovy Map for annotations + - genome: + type: string + description: Padfoot genome preset (`hg38`, `chm13`, `mm10`); bundled annotations exist for hg38 and mm10 only + - gff: + type: file + description: Optional custom GFF3 gene annotation (overrides bundled) + pattern: "*.gff{,3}{,.gz}" + - rm: + type: file + description: Optional custom RepeatMasker annotation (overrides bundled) + pattern: "*.{out,bed}{,.gz}" + +output: + - annotated_svs: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/annotated_svs.tsv": + type: file + description: Per-SV annotation table (genes, exons, repeats, microhomology, VNTR, complex SV cluster) + pattern: "annotated_svs.tsv" + - by_gene: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/by_gene.tsv": + type: file + description: Per-gene table of SV and copy-number impact + pattern: "by_gene.tsv" + - log: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/padfoot.log": + type: file + description: Padfoot log + pattern: "padfoot.log" + - versions: + - versions.yml: + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@Tim-Yu" +maintainers: + - "@Tim-Yu" diff --git a/modules/local/padfoot/tests/main.nf.test b/modules/local/padfoot/tests/main.nf.test new file mode 100644 index 00000000..5ebf5396 --- /dev/null +++ b/modules/local/padfoot/tests/main.nf.test @@ -0,0 +1,74 @@ +nextflow_process { + + name "Test Process PADFOOT" + script "../main.nf" + process "PADFOOT" + config "./nextflow.config" + + tag "modules" + tag "modules_local" + tag "padfoot" + + // Real Severus/Wakhan/SAVANA outputs are not available as public test data; only the I/O contract is tested here. + // An end-to-end run on real data is exercised through the pipeline (see docs/usage.md, Padfoot options). + + test("severus + wakhan - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test', paired_data:'test' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/vcf/test.rnaseq.vcf', checkIfExists: true), + 'severus', + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gvcf/test.genome.vcf', checkIfExists: true), + 'wakhan' + ] + input[1] = [ [ id:'ref' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'ref' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'padfoot' ], file("${moduleDir}", checkIfExists: true) ] + input[4] = [ [:], 'hg38', [], [] ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("savana + savana - tumour only - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test', paired_data:null ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/vcf/test.rnaseq.vcf', checkIfExists: true), + 'savana', + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.bed', checkIfExists: true), + 'savana' + ] + input[1] = [ [ id:'ref' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'ref' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'padfoot' ], file("${moduleDir}", checkIfExists: true) ] + input[4] = [ [:], 'hg38', [], [] ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/padfoot/tests/main.nf.test.snap b/modules/local/padfoot/tests/main.nf.test.snap new file mode 100644 index 00000000..56884a0e --- /dev/null +++ b/modules/local/padfoot/tests/main.nf.test.snap @@ -0,0 +1,144 @@ +{ + "savana + savana - tumour only - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test", + "paired_data": null + }, + "annotated_svs.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "test", + "paired_data": null + }, + "by_gene.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + { + "id": "test", + "paired_data": null + }, + "padfoot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + "versions.yml:md5,a39dd9bafbb792d343dd3d5f7b148606" + ], + "annotated_svs": [ + [ + { + "id": "test", + "paired_data": null + }, + "annotated_svs.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "by_gene": [ + [ + { + "id": "test", + "paired_data": null + }, + "by_gene.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "log": [ + [ + { + "id": "test", + "paired_data": null + }, + "padfoot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions": [ + "versions.yml:md5,a39dd9bafbb792d343dd3d5f7b148606" + ] + } + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-03T11:06:52.264251811" + }, + "severus + wakhan - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "annotated_svs.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "by_gene.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "padfoot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + "versions.yml:md5,a39dd9bafbb792d343dd3d5f7b148606" + ], + "annotated_svs": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "annotated_svs.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "by_gene": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "by_gene.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "log": [ + [ + { + "id": "test", + "paired_data": "test" + }, + "padfoot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions": [ + "versions.yml:md5,a39dd9bafbb792d343dd3d5f7b148606" + ] + } + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-03T11:06:46.571491429" + } +} \ No newline at end of file diff --git a/modules/local/padfoot/tests/nextflow.config b/modules/local/padfoot/tests/nextflow.config new file mode 100644 index 00000000..42fcd305 --- /dev/null +++ b/modules/local/padfoot/tests/nextflow.config @@ -0,0 +1,6 @@ +process { + withName: 'PADFOOT' { + ext.args = '--skip_RepeatMasker' + ext.prefix = { sv_caller == cna_caller ? "${sv_caller}" : "${sv_caller}_${cna_caller}" } + } +} diff --git a/modules/local/reconplot/environment.yml b/modules/local/reconplot/environment.yml new file mode 100644 index 00000000..4b9011f1 --- /dev/null +++ b/modules/local/reconplot/environment.yml @@ -0,0 +1,18 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +# R runtime for the ReConPlot wrapper (https://github.com/Tim-Yu/ReConPlot). The ReConPlot package itself +# (https://github.com/cortes-ciriano-lab/ReConPlot) is not on conda; the module installs it from the staged +# source tree when it is not already present in the environment. +channels: + - conda-forge + - bioconda +dependencies: + - "conda-forge::r-base=4.4.3" + - "conda-forge::r-ggplot2=4.0.3" + - "conda-forge::r-data.table=1.18.4" + - "conda-forge::r-optparse=1.8.2" + - "conda-forge::r-cowplot=1.2.0" + - "conda-forge::r-dplyr=1.2.1" + - "conda-forge::r-scales=1.4.0" + - "conda-forge::r-magrittr=2.0.5" + - "conda-forge::procps-ng" # `ps` required by Nextflow task metrics diff --git a/modules/local/reconplot/main.nf b/modules/local/reconplot/main.nf new file mode 100644 index 00000000..467c54a6 --- /dev/null +++ b/modules/local/reconplot/main.nf @@ -0,0 +1,95 @@ +process RECONPLOT { + tag "${meta.id}:${cn_source}+${sv_source}" + label 'process_low' + + conda "${moduleDir}/environment.yml" + // Built from containers/reconplot/Dockerfile: R deps + ReConPlot package (not on conda; the wrapper is staged as source). + // Override per site with `process { withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = ... } }`. + container "ghcr.io/tim-yu/reconplot@sha256:1145fc5aebe0227bec371f4c59b08b9a09871498e403c01b83f83973149ae9e7" + + input: + // cn_files / sv_files are the caller output files the wrapper's parsers discover by name. + // When cn_source == sv_source (e.g. 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(reconplot_src) // ReConPlot wrapper (contains run_reconplot.R + R/) + tuple val(meta3), path(reconplot_pkg) // ReConPlot R package source; installed only if the env lacks it + val(genome) // hg38 | hg19 | T2T | mm10 | mm39 + + output: + tuple val(meta), path("${prefix}/per_chromosome/*.{pdf,png}"), emit: per_chromosome + tuple val(meta), path("${prefix}/genome_wide/*.{pdf,png}") , emit: genome_wide + tuple val(meta), path("${prefix}/focus/*.{pdf,png}") , emit: focus , optional: true + tuple val(meta), path("${prefix}/*.reconplot_{cn,sv}.tsv") , emit: tables + tuple val(meta), path("${prefix}/reconplot.log") , emit: log + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' // shared filters (e.g. --min-svlen, --max-cn) + def args2 = task.ext.args2 ?: '' // genome-wide strip extras + def args3 = task.ext.args3 ?: '' // focus panel extras (--regions/--genes/--baf-track); focus skipped if empty + prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}") + def sample = meta.id + def source_args = cn_source == sv_source + ? "--source ${cn_source} --input cn_input" + : "--cn-source ${cn_source} --cn-input cn_input --sv-source ${sv_source} --sv-input sv_input" + // Wakhan's parser expects /solutions_ranks.tsv + //bed_output/*.bed + def layout_cmd = cn_source == 'wakhan' + ? "mkdir -p cn_input/solution_1/bed_output && mv cn_input/*.bed cn_input/solution_1/bed_output/" + : "" + def focus_cmd = args3 + ? """ + Rscript ${reconplot_src}/run_reconplot.R ${source_args} --sample ${sample} --prefix ${sample} \\ + --genome ${genome} --outdir ${prefix}/focus --layout together ${args} ${args3} 2>&1 | tee -a ${prefix}/reconplot.log + """ + : "" + + """ + # ReConPlot is pre-installed in the container; conda envs get it from the staged source tree + if ! Rscript -e 'suppressMessages(library(ReConPlot))' 2>/dev/null; then + mkdir -p rlib + R CMD INSTALL --no-docs --no-html -l rlib ${reconplot_pkg} > rlib_install.log 2>&1 + export R_LIBS=\$PWD/rlib\${R_LIBS:+:\$R_LIBS} + fi + + ${layout_cmd} + mkdir -p ${prefix} + + Rscript ${reconplot_src}/run_reconplot.R ${source_args} --sample ${sample} --prefix ${sample} \\ + --genome ${genome} --outdir ${prefix}/per_chromosome --regions all --layout separate \\ + --write-tables ${args} 2>&1 | tee ${prefix}/reconplot.log + mv ${prefix}/per_chromosome/*.reconplot_{cn,sv}.tsv ${prefix}/ + + Rscript ${reconplot_src}/run_reconplot.R ${source_args} --sample ${sample} --prefix ${sample} \\ + --genome ${genome} --outdir ${prefix}/genome_wide --regions all --layout together \\ + ${args} ${args2} 2>&1 | tee -a ${prefix}/reconplot.log + ${focus_cmd} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + reconplot: \$(Rscript -e 'cat(as.character(packageVersion("ReConPlot")))' 2>/dev/null) + reconplot_wrapper: \$(cat ${reconplot_src}/VERSION 2>/dev/null || echo unknown) + r-base: \$(Rscript -e 'cat(R.version\$major, R.version\$minor, sep=".")' 2>/dev/null) + ggplot2: \$(Rscript -e 'cat(as.character(packageVersion("ggplot2")))' 2>/dev/null) + END_VERSIONS + """ + + stub: + prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}") + """ + mkdir -p ${prefix}/per_chromosome ${prefix}/genome_wide + touch ${prefix}/per_chromosome/${meta.id}_chr1.pdf ${prefix}/per_chromosome/${meta.id}_chr1.png + touch ${prefix}/genome_wide/${meta.id}_genome_wide.pdf ${prefix}/genome_wide/${meta.id}_genome_wide.png + touch ${prefix}/${meta.id}.reconplot_cn.tsv ${prefix}/${meta.id}.reconplot_sv.tsv ${prefix}/reconplot.log + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + reconplot: stub + reconplot_wrapper: stub + r-base: stub + ggplot2: stub + END_VERSIONS + """ +} diff --git a/modules/local/reconplot/meta.yml b/modules/local/reconplot/meta.yml new file mode 100644 index 00000000..7c080ed2 --- /dev/null +++ b/modules/local/reconplot/meta.yml @@ -0,0 +1,103 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "reconplot" +description: Publication-style rearrangement + copy-number figures with ReConPlot, via the Tim-Yu/ReConPlot wrapper that reshapes ASCAT/Wakhan/Severus/SAVANA output into ReConPlot's input tables. +keywords: + - structural-variation + - copy-number + - visualisation + - somatic + - long-read +tools: + - "reconplot": + description: "ReConPlot: R package for visualising complex genome rearrangements together with copy number" + homepage: "https://github.com/cortes-ciriano-lab/ReConPlot" + documentation: "https://github.com/Tim-Yu/ReConPlot" + tool_dev_url: "https://github.com/cortes-ciriano-lab/ReConPlot" + doi: "10.1093/bioinformatics/btad719" + licence: ["MIT"] + identifier: "" + +input: + - - meta: + type: map + description: Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - cn_source: + type: string + description: CN parser name understood by the wrapper (`ascat`, `wakhan`, `savana`) + - cn_files: + type: file + description: Caller output files the CN parser discovers by name (staged into `cn_input/`). For `wakhan`, the two `*_copynumbers_segments_HP_[12].bed` plus `solutions_ranks.tsv`. When `cn_source == sv_source`, also holds the SV files. + - sv_source: + type: string + description: SV parser name understood by the wrapper (`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 + - reconplot_src: + type: directory + description: ReConPlot wrapper directory (`assets/reconplot`, contains `run_reconplot.R`, `R/` and `VERSION`) + - - meta3: + type: map + description: Groovy Map for the R package source + - reconplot_pkg: + type: directory + description: cortes-ciriano-lab/ReConPlot R package source; installed at run time only if the environment lacks it (containers ship it pre-installed) + - - genome: + type: string + description: ReConPlot genome preset (`hg38`, `hg19`, `T2T`, `mm10`, `mm39`) + +output: + - per_chromosome: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/per_chromosome/*.{pdf,png}": + type: file + description: One figure per chromosome + pattern: "*.{pdf,png}" + - genome_wide: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/genome_wide/*.{pdf,png}": + type: file + description: Genome-wide strip with all chromosomes side by side + pattern: "*.{pdf,png}" + - focus: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/focus/*.{pdf,png}": + type: file + description: Optional multi-panel figure for `--reconplot_regions` (with gene labels / BAF track) + pattern: "*.{pdf,png}" + - tables: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/*.reconplot_{cn,sv}.tsv": + type: file + description: Harmonised CN and SV tables as passed to ReConPlot + pattern: "*.reconplot_{cn,sv}.tsv" + - log: + - meta: + type: map + description: Groovy Map containing sample information + - "${prefix}/reconplot.log": + type: file + description: Wrapper log + pattern: "reconplot.log" + - versions: + - versions.yml: + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@Tim-Yu" +maintainers: + - "@Tim-Yu" diff --git a/modules/local/reconplot/tests/main.nf.test b/modules/local/reconplot/tests/main.nf.test new file mode 100644 index 00000000..19e3e4c6 --- /dev/null +++ b/modules/local/reconplot/tests/main.nf.test @@ -0,0 +1,72 @@ +nextflow_process { + + name "Test Process RECONPLOT" + script "../main.nf" + process "RECONPLOT" + config "./nextflow.config" + + tag "modules" + tag "modules_local" + tag "reconplot" + + // Real ASCAT/Wakhan/Severus/SAVANA outputs are not available as public test data; only the I/O contract is + // tested here. The three combinations were run end to end on real data through the pipeline. + + test("ascat + severus - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + 'ascat', + [ file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.bed', checkIfExists: true) ], + 'severus', + [ file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gvcf/test.genome.vcf', checkIfExists: true) ] + ] + input[1] = [ [ id:'reconplot' ], file("${projectDir}/modules/local/reconplot", checkIfExists: true) ] + input[2] = [ [ id:'reconplot_pkg' ], file("${projectDir}/modules/local/padfoot", 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(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.bed', checkIfExists: true) ], + 'savana', + [] + ] + input[1] = [ [ id:'reconplot' ], file("${projectDir}/modules/local/reconplot", checkIfExists: true) ] + input[2] = [ [ id:'reconplot_pkg' ], file("${projectDir}/modules/local/padfoot", checkIfExists: true) ] + input[3] = 'hg38' + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/local/reconplot/tests/main.nf.test.snap b/modules/local/reconplot/tests/main.nf.test.snap new file mode 100644 index 00000000..b49c6a5d --- /dev/null +++ b/modules/local/reconplot/tests/main.nf.test.snap @@ -0,0 +1,212 @@ +{ + "ascat + severus - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + [ + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_chr1.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "1": [ + [ + { + "id": "test" + }, + [ + "test_genome_wide.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_genome_wide.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "2": [ + + ], + "3": [ + [ + { + "id": "test" + }, + [ + "test.reconplot_cn.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.reconplot_sv.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "4": [ + [ + { + "id": "test" + }, + "reconplot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "5": [ + "versions.yml:md5,e5066ab640c75ba5b7dffcc4a0c48dd3" + ], + "focus": [ + + ], + "genome_wide": [ + [ + { + "id": "test" + }, + [ + "test_genome_wide.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_genome_wide.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "log": [ + [ + { + "id": "test" + }, + "reconplot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "per_chromosome": [ + [ + { + "id": "test" + }, + [ + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_chr1.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "tables": [ + [ + { + "id": "test" + }, + [ + "test.reconplot_cn.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.reconplot_sv.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "versions": [ + "versions.yml:md5,e5066ab640c75ba5b7dffcc4a0c48dd3" + ] + } + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-03T15:19:50.861183455" + }, + "savana single source - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + [ + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_chr1.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "1": [ + [ + { + "id": "test" + }, + [ + "test_genome_wide.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_genome_wide.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "2": [ + + ], + "3": [ + [ + { + "id": "test" + }, + [ + "test.reconplot_cn.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.reconplot_sv.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "4": [ + [ + { + "id": "test" + }, + "reconplot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "5": [ + "versions.yml:md5,e5066ab640c75ba5b7dffcc4a0c48dd3" + ], + "focus": [ + + ], + "genome_wide": [ + [ + { + "id": "test" + }, + [ + "test_genome_wide.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_genome_wide.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "log": [ + [ + { + "id": "test" + }, + "reconplot.log:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "per_chromosome": [ + [ + { + "id": "test" + }, + [ + "test_chr1.pdf:md5,d41d8cd98f00b204e9800998ecf8427e", + "test_chr1.png:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "tables": [ + [ + { + "id": "test" + }, + [ + "test.reconplot_cn.tsv:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.reconplot_sv.tsv:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ] + ], + "versions": [ + "versions.yml:md5,e5066ab640c75ba5b7dffcc4a0c48dd3" + ] + } + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.10.4" + }, + "timestamp": "2026-09-03T15:19:56.660161308" + } +} \ No newline at end of file diff --git a/modules/local/reconplot/tests/nextflow.config b/modules/local/reconplot/tests/nextflow.config new file mode 100644 index 00000000..355e1561 --- /dev/null +++ b/modules/local/reconplot/tests/nextflow.config @@ -0,0 +1,5 @@ +process { + withName: 'RECONPLOT' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } + } +} diff --git a/nextflow.config b/nextflow.config index 696d5161..a507df8c 100644 --- a/nextflow.config +++ b/nextflow.config @@ -83,6 +83,8 @@ params { skip_bamstats = false skip_ascat = false skip_savana = false + skip_padfoot = false + skip_reconplot = false skip_wakhan = false skip_fiber = false skip_normalfiber = false @@ -113,6 +115,26 @@ params { 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) + // Padfoot options + padfoot_url = 'https://github.com/Tim-Yu/Padfoot/archive/bf942055d4a0c710aed1178869998350823b15e0.tar.gz' + padfoot_dir = null // local Padfoot checkout (offline); overrides padfoot_url + padfoot_genome = null // hg38 | chm13 | mm10; null = inferred from --genome + padfoot_gff = null // custom GFF (required for chm13) + padfoot_rm = null // custom RepeatMasker BED (required for chm13) + padfoot_run_repeatmasker = true // containers ship full Dfam 4.0; conda uses RepeatMasker's bundled curated Dfam subset + + // ReConPlot options + reconplot_pkg_url = 'https://github.com/cortes-ciriano-lab/ReConPlot/archive/88bb76011703d7908aaa2e5bdecfe18ba5737db8.tar.gz' // R package source + reconplot_pkg_dir = null // local ReConPlot package checkout; overrides reconplot_pkg_url + 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_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 + reconplot_baf_track = false // het-SNP BAF track on the focus panel (ASCAT/SAVANA) + reconplot_format = 'pdf,png' + // ASCAT options ascat_ploidy = null ascat_min_base_qual = 20 diff --git a/nextflow_schema.json b/nextflow_schema.json index 713e2baa..9902bd4d 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -376,6 +376,103 @@ } } }, + "padfoot_options": { + "title": "Padfoot options", + "type": "object", + "description": "Options for Padfoot SV/CNA functional annotation", + "default": "", + "properties": { + "padfoot_url": { + "type": "string", + "default": "https://github.com/Tim-Yu/Padfoot/archive/bf942055d4a0c710aed1178869998350823b15e0.tar.gz", + "description": "URL of a Padfoot source tarball (GitHub archive). Padfoot is not distributed on bioconda; the tarball is downloaded once per run." + }, + "padfoot_dir": { + "type": "string", + "format": "directory-path", + "exists": true, + "description": "Local Padfoot checkout (directory containing padfoot.py and beds/). Overrides --padfoot_url; use for offline systems." + }, + "padfoot_genome": { + "type": "string", + "enum": ["hg38", "chm13", "mm10"], + "description": "Padfoot genome preset (`hg38`, `chm13`, `mm10`). Default: inferred from `--genome`.", + "help_text": "GRCh38 -> hg38, CHM13 -> chm13. Padfoot bundles gene and repeat annotations for hg38 and mm10 only; chm13 and other presets need `--padfoot_gff` and `--padfoot_rm`, otherwise Padfoot is skipped with a warning." + }, + "padfoot_gff": { + "type": "string", + "format": "file-path", + "exists": true, + "description": "Custom GFF3 gene annotation for Padfoot" + }, + "padfoot_rm": { + "type": "string", + "format": "file-path", + "exists": true, + "description": "Custom RepeatMasker annotation (.fa.out / BED) for Padfoot" + }, + "padfoot_run_repeatmasker": { + "type": "boolean", + "default": true, + "description": "Run RepeatMasker on inserted sequences. Containers ship the full Dfam 4.0 database; conda uses the curated Dfam subset bundled with bioconda RepeatMasker." + } + } + }, + "reconplot_options": { + "title": "ReConPlot options", + "type": "object", + "description": "Options for ReConPlot rearrangement/copy-number figures", + "default": "", + "properties": { + "reconplot_pkg_url": { + "type": "string", + "default": "https://github.com/cortes-ciriano-lab/ReConPlot/archive/88bb76011703d7908aaa2e5bdecfe18ba5737db8.tar.gz", + "description": "URL of the ReConPlot R package source (installed at run time when the environment lacks it, i.e. under conda)." + }, + "reconplot_pkg_dir": { + "type": "string", + "format": "directory-path", + "exists": true, + "description": "Local checkout of the ReConPlot R package. Overrides --reconplot_pkg_url." + }, + "reconplot_genome": { + "type": "string", + "enum": ["hg38", "hg19", "T2T", "mm10", "mm39"], + "description": "ReConPlot genome build for gene/chromosome annotations. Default: inferred from --genome (GRCh38 -> hg38, CHM13 -> T2T)." + }, + "reconplot_max_cn": { + "type": "integer", + "default": 8, + "description": "Copy number axis ceiling" + }, + "reconplot_min_svlen": { + "type": "integer", + "default": 0, + "description": "Drop intra-chromosomal SVs shorter than this many bp (translocations are kept). 0 keeps everything." + }, + "reconplot_exclude_vntr": { + "type": "boolean", + "description": "Drop Severus SVs flagged as inside a VNTR" + }, + "reconplot_regions": { + "type": "string", + "description": "Comma-separated regions (e.g. 'chr8,chr17:30000000-50000000') for an additional multi-panel focus figure" + }, + "reconplot_genes": { + "type": "string", + "description": "Comma-separated HUGO gene symbols to label on the focus figure" + }, + "reconplot_baf_track": { + "type": "boolean", + "description": "Add a het-SNP BAF track to the focus figure (ASCAT and SAVANA copy number only)" + }, + "reconplot_format": { + "type": "string", + "default": "pdf,png", + "description": "Comma-separated output formats: pdf, png" + } + } + }, "ascat_parameters": { "title": "ASCAT parameters", "type": "object", @@ -547,6 +644,14 @@ "type": "boolean", "description": "Skip SAVANA (SV + copy-number calling)" }, + "skip_padfoot": { + "type": "boolean", + "description": "Skips Padfoot SV/CNA annotation" + }, + "skip_reconplot": { + "type": "boolean", + "description": "Skips ReConPlot SV/CNA figures" + }, "skip_m6a": { "type": "boolean", "description": "Skip m6a calling by Fibertools" @@ -771,6 +876,12 @@ { "$ref": "#/$defs/savana_options" }, + { + "$ref": "#/$defs/padfoot_options" + }, + { + "$ref": "#/$defs/reconplot_options" + }, { "$ref": "#/$defs/ascat_parameters" }, diff --git a/subworkflows/local/padfoot_annotation.nf b/subworkflows/local/padfoot_annotation.nf new file mode 100644 index 00000000..de46a99d --- /dev/null +++ b/subworkflows/local/padfoot_annotation.nf @@ -0,0 +1,75 @@ +// IMPORT MODULES +include { PADFOOT as PADFOOT_SEVERUS_WAKHAN } from '../../modules/local/padfoot/main' +include { PADFOOT as PADFOOT_SAVANA } from '../../modules/local/padfoot/main' +include { WGET as PADFOOT_WGET } from '../../modules/nf-core/wget/main' +include { UNTAR as PADFOOT_UNTAR } from '../../modules/nf-core/untar/main' + +// +// Padfoot annotation of somatic SVs + CNAs, once per caller pair that produced output for a sample: +// Severus SVs + the top-ranked Wakhan integer-CN VCF, and SAVANA SVs + SAVANA absolute CN. +// Pass channel.empty() for a caller that did not run. Padfoot is not on bioconda: its source tree +// comes from params.padfoot_url (GitHub archive) or a local checkout in params.padfoot_dir. +// +workflow PADFOOT_ANNOTATION { + + take: + severus_vcf // [meta, severus_somatic.vcf.gz] + wakhan_vcf_files // [meta, [wakhan_cna_*.vcf, ...]] -- every fitted solution + savana_vcf // [meta, classified.somatic.vcf] + savana_cna // [meta, segmented_absolute_copy_number.tsv] + fasta // [[:], fasta] + fai // [[:], fai] + annot // [[:], padfoot_genome, gff | [], rm | []] + + main: + ch_versions = channel.empty() + + if (params.padfoot_dir) { + padfoot_src = channel.value([[id: 'padfoot'], file(params.padfoot_dir, type: 'dir', checkIfExists: true)]) + } + else { + PADFOOT_WGET( channel.value([[id: 'padfoot'], params.padfoot_url]) ) + PADFOOT_UNTAR( PADFOOT_WGET.out.outfile ) + padfoot_src = PADFOOT_UNTAR.out.untar + ch_versions = ch_versions.mix(PADFOOT_WGET.out.versions) + } + // padfoot_src: [meta, dir] -- padfoot.py + beds/ + + // + // MODULE: PADFOOT_SEVERUS_WAKHAN (label: process_medium) + // Input: [meta, severus_somatic.vcf.gz, 'severus', wakhan_cna_integers.vcf, 'wakhan'] + // + // Wakhan writes every fitted solution; solution_1/ holds the top-ranked one + wakhan_vcf_files + .map { meta, vcfs -> + def integers = [vcfs].flatten().findAll { vcf -> vcf.name.endsWith('_wakhan_cna_integers.vcf') } + return [meta, integers.find { vcf -> vcf.toString().contains('/solution_1/') } ?: integers[0]] + } + .filter { _meta, vcf -> vcf != null } + .set { wakhan_best_cna } + // wakhan_best_cna: [meta, wakhan_cna_integers.vcf] + + severus_vcf + .join(wakhan_best_cna) + .map { meta, sv, cna -> [meta, sv, 'severus', cna, 'wakhan'] } + .set { severus_wakhan_input } + + PADFOOT_SEVERUS_WAKHAN( severus_wakhan_input, fasta, fai, padfoot_src, annot ) + ch_versions = ch_versions.mix(PADFOOT_SEVERUS_WAKHAN.out.versions) + + // + // MODULE: PADFOOT_SAVANA (label: process_medium) + // Input: [meta, classified.somatic.vcf, 'savana', segmented_absolute_copy_number.tsv, 'savana'] + // SAVANA CN is only present when a fit was found, so the join drops unfitted samples + // + savana_vcf + .join(savana_cna) + .map { meta, sv, cna -> [meta, sv, 'savana', cna, 'savana'] } + .set { savana_input } + + PADFOOT_SAVANA( savana_input, fasta, fai, padfoot_src, annot ) + ch_versions = ch_versions.mix(PADFOOT_SAVANA.out.versions) + + emit: + versions = ch_versions // [versions.yml] +} diff --git a/subworkflows/local/paired/paired_savana.nf b/subworkflows/local/paired/paired_savana.nf index e719be03..b7295e13 100644 --- a/subworkflows/local/paired/paired_savana.nf +++ b/subworkflows/local/paired/paired_savana.nf @@ -76,6 +76,9 @@ workflow PAIRED_SAVANA { ) emit: - somatic_vcf = SAVANA_CLASSIFY.out.somatic_vcf // [meta, vcf] -- classified somatic SV VCF - cn_calls = SAVANA_CNA.out.cna // [meta, tsv] -- segmented absolute copy number + somatic_vcf = SAVANA_CLASSIFY.out.somatic_vcf // [meta, vcf] -- classified somatic SV VCF + somatic_bedpe = SAVANA_CLASSIFY.out.somatic_bedpe // [meta, bedpe] -- classified somatic SVs in BEDPE (ReConPlot) + cn_calls = SAVANA_CNA.out.cna // [meta, tsv] -- segmented absolute copy number + fitted_purity_ploidy = SAVANA_CNA.out.fitted_purity_ploidy // [meta, tsv] -- selected purity/ploidy fit (absent when no fit) + allele_counts = SAVANA_CNA.out.allele_counts // [meta, bed] -- het-SNP allele counts (ReConPlot BAF track) } diff --git a/subworkflows/local/reconplot_figures.nf b/subworkflows/local/reconplot_figures.nf new file mode 100644 index 00000000..0801db6e --- /dev/null +++ b/subworkflows/local/reconplot_figures.nf @@ -0,0 +1,103 @@ +// IMPORT MODULES +include { RECONPLOT as RECONPLOT_SEVERUS_ASCAT } from '../../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_SEVERUS_WAKHAN } from '../../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_SAVANA } from '../../modules/local/reconplot/main' +include { WGET as RECONPLOT_PKG_WGET } from '../../modules/nf-core/wget/main' +include { UNTAR as RECONPLOT_PKG_UNTAR } from '../../modules/nf-core/untar/main' + +// +// ReConPlot rearrangement + copy-number figures 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 wrapper (assets/reconplot) is shipped with the pipeline; the ReConPlot R +// package is staged as source from params.reconplot_pkg_url or a local checkout in params.reconplot_pkg_dir. +// +workflow RECONPLOT_FIGURES { + + take: + severus_vcf // [meta, severus_somatic.vcf.gz] + ascat_segments // [meta, segments.txt] + ascat_purityploidy // [meta, purityploidy.txt] + ascat_bafs // [meta, [*BAF.txt]] + wakhan_bed_files // [meta, [bed_output/*.bed, ...]] -- every fitted solution + wakhan_solutions_ranks // [meta, solutions_ranks.tsv] + savana_cna // [meta, segmented_absolute_copy_number.tsv] + savana_bedpe // [meta, classified.somatic.bedpe] + savana_purity_ploidy // [meta, fitted_purity_ploidy.tsv] + savana_allele_counts // [meta, allele_counts_hetSNPs.bed] -- absent without an SNP source + genome // ReConPlot genome preset: hg38, hg19, T2T, mm10 or mm39 + + main: + ch_versions = channel.empty() + + reconplot_src = channel.value([[id: 'reconplot'], file("${projectDir}/assets/reconplot", type: 'dir', checkIfExists: true)]) + if (params.reconplot_pkg_dir) { + reconplot_pkg = channel.value([[id: 'reconplot_pkg'], file(params.reconplot_pkg_dir, type: 'dir', checkIfExists: true)]) + } + else { + RECONPLOT_PKG_WGET( channel.value([[id: 'reconplot_pkg'], params.reconplot_pkg_url]) ) + RECONPLOT_PKG_UNTAR( RECONPLOT_PKG_WGET.out.outfile ) + reconplot_pkg = RECONPLOT_PKG_UNTAR.out.untar + ch_versions = ch_versions.mix(RECONPLOT_PKG_WGET.out.versions) + } + // reconplot_src: [meta, dir] -- wrapper; reconplot_pkg: [meta, dir] -- R package source (conda installs it at run time) + + severus_sv_files = severus_vcf.map { meta, vcf -> [meta, [vcf]] } + // severus_sv_files: [meta, [severus_somatic.vcf.gz]] + + // + // MODULE: RECONPLOT_SEVERUS_ASCAT (label: process_low) + // Input: [meta, 'ascat', [segments.txt, purityploidy.txt, *BAF.txt], 'severus', [vcf]] + // the wrapper picks .tumour_tumourBAF.txt from the BAF tables by name + // + ascat_segments + .join(ascat_purityploidy) + .join(ascat_bafs) + .map { meta, seg, pp, bafs -> [meta, [seg, pp, bafs].flatten()] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'ascat', cn, 'severus', sv] } + .set { ascat_input } + + RECONPLOT_SEVERUS_ASCAT( ascat_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_SEVERUS_ASCAT.out.versions) + + // + // MODULE: RECONPLOT_SEVERUS_WAKHAN (label: process_low) + // Input: [meta, 'wakhan', [HP_1.bed, HP_2.bed, solutions_ranks.tsv], 'severus', [vcf]] + // the two allele-specific segment BEDs of the top-ranked solution (solution_1/) + // + wakhan_bed_files + .map { meta, beds -> + def hp = [beds].flatten().findAll { bed -> bed.name ==~ /.*_copynumbers_segments_HP_[12]\.bed/ } + def best = hp.findAll { bed -> bed.toString().contains('/solution_1/') } ?: hp + return [meta, best.unique { bed -> bed.name }] + } + .filter { _meta, beds -> beds.size() == 2 } + .join(wakhan_solutions_ranks) + .map { meta, beds, ranks -> [meta, beds + [ranks]] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'wakhan', cn, 'severus', sv] } + .set { wakhan_input } + + RECONPLOT_SEVERUS_WAKHAN( wakhan_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_SEVERUS_WAKHAN.out.versions) + + // + // MODULE: RECONPLOT_SAVANA (label: process_low) + // Input: [meta, 'savana', [cna.tsv, somatic.bedpe, fitted_purity_ploidy.tsv, allele_counts.bed], 'savana', []] + // single-source mode; allele_counts is optional, and a sample with allele counts but no CN fit + // only exists on the right of the remainder join ([meta, null, bed]) and is dropped + // + savana_cna + .join(savana_bedpe) + .join(savana_purity_ploidy) + .join(savana_allele_counts, remainder: true) + .filter { row -> row[1] != null } + .map { meta, cna, bedpe, pp, hetsnp -> [meta, 'savana', [cna, bedpe, pp, hetsnp].findAll { f -> f != null }, 'savana', []] } + .set { savana_input } + + RECONPLOT_SAVANA( savana_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_SAVANA.out.versions) + + emit: + versions = ch_versions // [versions.yml] +} diff --git a/subworkflows/local/tests/padfoot_annotation.nf.test b/subworkflows/local/tests/padfoot_annotation.nf.test new file mode 100644 index 00000000..bede62c7 --- /dev/null +++ b/subworkflows/local/tests/padfoot_annotation.nf.test @@ -0,0 +1,73 @@ +nextflow_workflow { + + name "Test Workflow PADFOOT_ANNOTATION" + script "../padfoot_annotation.nf" + workflow "PADFOOT_ANNOTATION" + + tag "subworkflows" + tag "subworkflows_local" + tag "padfoot_annotation" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // Padfoot needs whole-genome caller output, so the pipeline tests skip it. Under -stub what is + // tested is the pairing: Severus SVs + the top-ranked Wakhan solution, SAVANA SVs + SAVANA copy + // number, and that a sample missing one side of a pair is dropped rather than failing the run. + + test("Severus + Wakhan and SAVANA pairs form per sample") { + + options "-stub" + + when { + params { + outdir = "$outputDir" + } + workflow { + """ + // Empty stand-ins carrying the real callers' file names; Padfoot's stub does not read them + def mk = { String rel -> def f = file("\${launchDir}/fixtures/\${rel}"); f.parent.mkdirs(); f.text = ''; return f } + 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' ] + + // Severus somatic SVs for both paired samples + input[0] = channel.of( + [ s1, mk('severus/sample1/severus_somatic.vcf.gz') ], + [ s2, mk('severus/sample2/severus_somatic.vcf.gz') ] + ) + // Wakhan: sample1 has two solutions (solution_1 must be picked); sample2 has no integer-CN VCF + input[1] = channel.of( + [ s1, [ mk('wakhan/sample1/solution_1/sample1_2.0_1.0_wakhan_cna_integers.vcf'), + mk('wakhan/sample1/solution_1/sample1_2.0_1.0_wakhan_cna_subclonals.vcf'), + mk('wakhan/sample1/solution_2/sample1_3.0_0.5_wakhan_cna_integers.vcf') ] ], + [ s2, [ mk('wakhan/sample2/solution_1/sample2_2.0_1.0_wakhan_cna_subclonals.vcf') ] ] + ) + // SAVANA: SVs for sample1 and sample3, a copy-number fit for sample1 only + input[2] = channel.of( + [ s1, mk('savana/sample1/sample1.classified.somatic.vcf') ], + [ s3, mk('savana/sample3/sample3.classified.somatic.vcf') ] + ) + input[3] = channel.of([ s1, mk('savana/sample1/sample1_segmented_absolute_copy_number.tsv') ]) + input[4] = channel.value([ [:], mk('ref/genome.fasta') ]) + input[5] = channel.value([ [:], mk('ref/genome.fasta.fai') ]) + input[6] = channel.value([ [:], 'hg38', [], [] ]) + """ + } + } + + then { + def tasks = workflow.trace.tasks().collect { task -> task.name } + + assertAll( + { assert workflow.success }, + // sample1 only: sample2 has no integer-CN VCF to pair with + { assert tasks.count { name -> name.contains('PADFOOT_SEVERUS_WAKHAN') } == 1 }, + { assert tasks.any { name -> name.contains('PADFOOT_SEVERUS_WAKHAN (sample1:') } }, + // sample1 only: sample3 has SAVANA SVs but no copy-number fit + { assert tasks.count { name -> name.contains('PADFOOT_SAVANA') } == 1 }, + { assert tasks.any { name -> name.contains('PADFOOT_SAVANA (sample1:') } }, + { assert workflow.out.versions.size() >= 1 } + ) + } + } +} diff --git a/subworkflows/local/tests/reconplot_figures.nf.test b/subworkflows/local/tests/reconplot_figures.nf.test new file mode 100644 index 00000000..49745db5 --- /dev/null +++ b/subworkflows/local/tests/reconplot_figures.nf.test @@ -0,0 +1,98 @@ +nextflow_workflow { + + name "Test Workflow RECONPLOT_FIGURES" + script "../reconplot_figures.nf" + workflow "RECONPLOT_FIGURES" + + tag "subworkflows" + tag "subworkflows_local" + tag "reconplot_figures" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // ReConPlot 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 + // (with or without allele counts), and that a SAVANA sample with allele counts but no fit is dropped. + + test("each CN/SV caller pair forms per sample, no-fit SAVANA samples are dropped") { + + options "-stub" + + when { + params { + outdir = "$outputDir" + } + workflow { + """ + // Empty stand-ins carrying the real callers' file names; ReConPlot's stub does not read them + def mk = { String rel -> def f = file("\${launchDir}/fixtures/\${rel}"); f.parent.mkdirs(); f.text = ''; return f } + 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' ] + def s4 = [ id:'sample4', paired_data:false, platform:'ont' ] + + // Severus somatic SVs for both paired samples + input[0] = channel.of( + [ s1, mk('severus/sample1/severus_somatic.vcf.gz') ], + [ s2, mk('severus/sample2/severus_somatic.vcf.gz') ] + ) + // ASCAT for sample1 only + input[1] = channel.of([ s1, mk('ascat/sample1.segments.txt') ]) + input[2] = channel.of([ s1, mk('ascat/sample1.purityploidy.txt') ]) + input[3] = channel.of([ s1, [ mk('ascat/sample1.tumour_tumourBAF.txt'), mk('ascat/sample1.tumour_normalBAF.txt') ] ]) + // Wakhan: sample1 has two solutions (solution_1's two HP beds must be picked); sample2 has one HP bed only + input[4] = channel.of( + [ s1, [ mk('wakhan/sample1/solution_1/sample1_2.0_1.0_copynumbers_segments_HP_1.bed'), + mk('wakhan/sample1/solution_1/sample1_2.0_1.0_copynumbers_segments_HP_2.bed'), + mk('wakhan/sample1/solution_1/sample1_2.0_1.0_copynumbers_subclonal_segments_HP_1.bed'), + mk('wakhan/sample1/solution_2/sample1_3.0_0.5_copynumbers_segments_HP_1.bed'), + mk('wakhan/sample1/solution_2/sample1_3.0_0.5_copynumbers_segments_HP_2.bed') ] ], + [ s2, [ mk('wakhan/sample2/solution_1/sample2_2.0_1.0_copynumbers_segments_HP_1.bed') ] ] + ) + input[5] = channel.of( + [ s1, mk('wakhan/sample1/solutions_ranks.tsv') ], + [ s2, mk('wakhan/sample2/solutions_ranks.tsv') ] + ) + // SAVANA: fits for sample1 and sample3; sample4 has allele counts but no fit + input[6] = channel.of( + [ s1, mk('savana/sample1/sample1_segmented_absolute_copy_number.tsv') ], + [ s3, mk('savana/sample3/sample3_segmented_absolute_copy_number.tsv') ] + ) + input[7] = channel.of( + [ s1, mk('savana/sample1/sample1.classified.somatic.bedpe') ], + [ s3, mk('savana/sample3/sample3.classified.somatic.bedpe') ] + ) + input[8] = channel.of( + [ s1, mk('savana/sample1/sample1_fitted_purity_ploidy.tsv') ], + [ s3, mk('savana/sample3/sample3_fitted_purity_ploidy.tsv') ] + ) + input[9] = channel.of( + [ s1, mk('savana/sample1/sample1_allele_counts_hetSNPs.bed') ], + [ s4, mk('savana/sample4/sample4_allele_counts_hetSNPs.bed') ] + ) + input[10] = 'hg38' + """ + } + } + + then { + def tasks = workflow.trace.tasks().collect { task -> task.name } + + assertAll( + { assert workflow.success }, + // ASCAT + Severus: sample1 only + { assert tasks.count { name -> name.contains('RECONPLOT_SEVERUS_ASCAT') } == 1 }, + { assert tasks.any { name -> name.contains('RECONPLOT_SEVERUS_ASCAT (sample1:') } }, + // Wakhan + Severus: sample1 only, sample2 lacks the second HP bed + { assert tasks.count { name -> name.contains('RECONPLOT_SEVERUS_WAKHAN') } == 1 }, + { assert tasks.any { name -> name.contains('RECONPLOT_SEVERUS_WAKHAN (sample1:') } }, + // SAVANA: sample1 (with allele counts) and sample3 (without); sample4 has no fit and is dropped + { assert tasks.count { name -> name.contains('RECONPLOT_SAVANA') } == 2 }, + { assert tasks.any { name -> name.contains('RECONPLOT_SAVANA (sample1:') } }, + { assert tasks.any { name -> name.contains('RECONPLOT_SAVANA (sample3:') } }, + { assert !tasks.any { name -> name.contains('(sample4:') } }, + { assert workflow.out.versions.size() >= 1 } + ) + } + } +} diff --git a/subworkflows/local/tumor_only/tumoronly_savana.nf b/subworkflows/local/tumor_only/tumoronly_savana.nf index ba4b32b7..6a0ba350 100644 --- a/subworkflows/local/tumor_only/tumoronly_savana.nf +++ b/subworkflows/local/tumor_only/tumoronly_savana.nf @@ -45,6 +45,9 @@ workflow TUMORONLY_SAVANA { ) emit: - somatic_vcf = SAVANA_TO.out.somatic_vcf // [meta, vcf] -- classified somatic SV VCF - cn_calls = SAVANA_TO.out.cna // [meta, tsv] -- segmented absolute copy number + somatic_vcf = SAVANA_TO.out.somatic_vcf // [meta, vcf] -- classified somatic SV VCF + somatic_bedpe = SAVANA_TO.out.somatic_bedpe // [meta, bedpe] -- classified somatic SVs in BEDPE (ReConPlot) + 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/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf b/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf index a64c405c..a2396b5a 100644 --- a/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf +++ b/subworkflows/local/utils_nfcore_lrsomatic_pipeline/main.nf @@ -239,6 +239,7 @@ workflow PIPELINE_COMPLETION { def validateInputParameters() { genomeExistsError() validateReportGenePanels() + validateSvAnnotationParams() } // @@ -316,6 +317,28 @@ def validateReportGenePanels() { } } +// +// Warn on SV/CNA annotation and plotting parameter combinations that cannot produce output +// +def validateSvAnnotationParams() { + if (!params.skip_padfoot) { + def padfoot_genome = params.padfoot_genome ?: + (params.genome == 'GRCh38' ? 'hg38' : params.genome == 'CHM13' ? 'chm13' : null) + // Must mirror `padfoot_annot_ok` in workflows/lrsomatic.nf: a null genome disables + // Padfoot even when --padfoot_gff/--padfoot_rm are supplied. + def padfoot_annot_ok = padfoot_genome && ((padfoot_genome in ['hg38', 'mm10']) || (params.padfoot_gff && params.padfoot_rm)) + if (!padfoot_annot_ok) { + log.warn "Padfoot will be skipped: no annotations for genome '${params.genome}' (padfoot_genome=${padfoot_genome}). " + + "Set --padfoot_genome hg38|mm10, or set --padfoot_genome together with --padfoot_gff and --padfoot_rm." + } + } + + 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." + } +} + // // Validate channels from input samplesheet // diff --git a/tests/.nftignore b/tests/.nftignore index dfc63dec..5808bc53 100644 --- a/tests/.nftignore +++ b/tests/.nftignore @@ -34,3 +34,6 @@ pipeline_info/*.{html,json,txt,yml} # samtools merge gives sample4's colliding @PG IDs a random suffix, so this BAM's md5 differs every run (reads asserted in tests/clair_only.nf.test) sample4/bamfiles/sample4_tumor.bam sample4/bamfiles/sample4_tumor.bam.bai +*/padfoot/*/padfoot.log +*/reconplot/*/reconplot.log +*/reconplot/*/{per_chromosome,genome_wide,focus}/*.{pdf,png} diff --git a/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index 37bd1c4c..04c3e558 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -60,6 +60,8 @@ include { PAIRED_SMALLVAR_GERMLINE } from '../subworkflows/local/paired/p include { PHASING_HAPLOTYPING } from '../subworkflows/local/phasing_haplotyping' include { TUMORONLY_SAVANA } from '../subworkflows/local/tumor_only/tumoronly_savana' include { PAIRED_SAVANA } from '../subworkflows/local/paired/paired_savana' +include { PADFOOT_ANNOTATION } from '../subworkflows/local/padfoot_annotation' +include { RECONPLOT_FIGURES } from '../subworkflows/local/reconplot_figures' @@ -1112,7 +1114,11 @@ workflow LRSOMATIC { // SV_VEP below, alongside Severus's SVs. // - savana_somatic_vcf = channel.empty() + savana_somatic_vcf = channel.empty() + savana_cna = channel.empty() + savana_somatic_bedpe = channel.empty() + savana_fitted_purity_ploidy = channel.empty() + savana_allele_counts = channel.empty() if (!params.skip_savana) { // SAVANA reads the HP (haplotype) tag per read and its README recommends phased BAMs, @@ -1197,6 +1203,15 @@ workflow LRSOMATIC { .set { savana_somatic_vcf } // savana_somatic_vcf: [meta, vcf] + // Copy-number products consumed by Padfoot / ReConPlot below. All optional: absent without + // an SNP source, and cna/fitted_purity_ploidy absent when SAVANA finds no acceptable fit. + TUMORONLY_SAVANA.out.cn_calls.mix(PAIRED_SAVANA.out.cn_calls).set { savana_cna } + TUMORONLY_SAVANA.out.somatic_bedpe.mix(PAIRED_SAVANA.out.somatic_bedpe).set { savana_somatic_bedpe } + TUMORONLY_SAVANA.out.fitted_purity_ploidy.mix(PAIRED_SAVANA.out.fitted_purity_ploidy).set { savana_fitted_purity_ploidy } + TUMORONLY_SAVANA.out.allele_counts.mix(PAIRED_SAVANA.out.allele_counts).set { savana_allele_counts } + // savana_cna: [meta, segmented_absolute_copy_number.tsv] savana_somatic_bedpe: [meta, classified.somatic.bedpe] + // savana_fitted_purity_ploidy: [meta, tsv] savana_allele_counts: [meta, allele_counts_hetSNPs.bed] + if (!params.skip_vep) { // // MODULE: VEP_SAVANA (ENSEMBLVEP_VEP alias; label: process_medium) @@ -1262,6 +1277,50 @@ workflow LRSOMATIC { // ch_wakhan_files: [meta, [file_or_dir, ...]] } + // + // SUBWORKFLOW: PADFOOT_ANNOTATION -- Padfoot SV/CNA annotation per caller pair (Severus + Wakhan, SAVANA) + // Padfoot bundles annotations for hg38 and mm10 only; other genomes need --padfoot_gff and --padfoot_rm + // (validateInputParameters() warns when this gate is not met). + // + def padfoot_genome = params.padfoot_genome ?: + (params.genome == 'GRCh38' ? 'hg38' : params.genome == 'CHM13' ? 'chm13' : null) + def padfoot_annot_ok = padfoot_genome && ((padfoot_genome in ['hg38', 'mm10']) || (params.padfoot_gff && params.padfoot_rm)) + + if (!params.skip_padfoot && padfoot_annot_ok) { + PADFOOT_ANNOTATION ( + SEVERUS.out.somatic_vcf, + params.skip_wakhan ? channel.empty() : WAKHAN.out.vcf_files, + savana_somatic_vcf, + savana_cna, + ch_fasta, + ch_fai, + [[:], padfoot_genome, + params.padfoot_gff ? file(params.padfoot_gff, checkIfExists: true) : [], + params.padfoot_rm ? file(params.padfoot_rm, checkIfExists: true) : []] + ) + ch_versions = ch_versions.mix(PADFOOT_ANNOTATION.out.versions) + } + + // + // SUBWORKFLOW: RECONPLOT_FIGURES -- ReConPlot figures per CN/SV caller pair (ASCAT + Severus, Wakhan + Severus, SAVANA) + // + if (!params.skip_reconplot) { + RECONPLOT_FIGURES ( + SEVERUS.out.somatic_vcf, + params.skip_ascat ? channel.empty() : ASCAT.out.segments, + params.skip_ascat ? channel.empty() : ASCAT.out.purityploidy, + params.skip_ascat ? channel.empty() : ASCAT.out.bafs, + params.skip_wakhan ? channel.empty() : WAKHAN.out.bed_files, + params.skip_wakhan ? channel.empty() : WAKHAN.out.solutions_ranks, + savana_cna, + savana_somatic_bedpe, + savana_fitted_purity_ploidy, + savana_allele_counts, + params.reconplot_genome ?: (params.genome == 'CHM13' ? 'T2T' : 'hg38') + ) + ch_versions = ch_versions.mix(RECONPLOT_FIGURES.out.versions) + } + // // MODULE: LRSOMATICREPORT -- per-sample HTML report; all inputs optional, so joins use remainder: true on the tumor id //