diff --git a/bin/plot_wisecondorx_cnv.R b/bin/plot_wisecondorx_cnv.R new file mode 100755 index 0000000..d91b89c --- /dev/null +++ b/bin/plot_wisecondorx_cnv.R @@ -0,0 +1,159 @@ +#!/usr/bin/env Rscript + +suppressPackageStartupMessages({ + library(argparser) + library(readr) + library(dplyr) + library(ggplot2) +}) + +## arguments + +p <- arg_parser("Plot WisecondorX output (bin-level log2 ratio + segments)") +p <- add_argument(p, "--id", help = "Sample ID") +p <- add_argument(p, "--seg_file", help = "Segmented file with calls from WisecondorX", nargs = Inf) +p <- add_argument(p, "--binfile", help = "Bin-level file from WisecondorX") +p <- add_argument(p, "--outdir", help = "Output directory", default = ".") +p <- add_argument(p, "--ratio_limit", help = "Log2ratio limit of additional WisecondorX plots (application-generated plots are unaffected)", default = 1) +argv <- parse_args(p) + +sample_id <- argv$id +bin_path <- argv$binfile +output_dir <- argv$outdir +ratio_limit <- argv$ratio_limit + +# take the right seg_file pattern +seg_candidates <- argv$seg_file +if (length(seg_candidates) > 1) { + seg_candidates <- seg_candidates[!grepl("_gistic\\.seg$", seg_candidates)] +} +seg_path <- seg_candidates[1] + +dir.create(output_dir, showWarnings = FALSE, recursive = TRUE) + + +# Fix column names +read_and_normalize <- function(path) { + df <- read_tsv(path, show_col_types = FALSE) %>% + rename_with(tolower) + + if ("chrom" %in% names(df)) df <- rename(df, chr = chrom) + if ("seg.mean.adj" %in% names(df)) df <- rename(df, ratio = `seg.mean.adj`) + + df +} + +# Remap WisecondorX gain/loss/neut calls to readable labels for the plot + +classify_segments <- function(df) { + call_map <- c(gain = "GAIN", loss = "LOSS", neut = "NEUTRAL") + df %>% mutate(call = recode(tolower(call), !!!call_map, .default = "NEUTRAL")) +} + +# Plot style formatting, genomic coordinates +compute_genomic_coords <- function(bins, segs) { + chr_order <- c(as.character(1:22), "X", "Y") + present_chr <- intersect(chr_order, unique(c(bins$chr, segs$chr))) + + chr_lengths <- bins %>% + filter(chr %in% present_chr) %>% + group_by(chr) %>% + summarise(len = max(end), .groups = "drop") %>% + mutate(chr = factor(chr, levels = present_chr)) %>% + arrange(chr) %>% + mutate(offset = lag(cumsum(len), default = 0)) + + add_offset <- function(df) { + df %>% + filter(chr %in% present_chr) %>% + mutate(chr = factor(chr, levels = present_chr)) %>% + left_join(chr_lengths %>% select(chr, offset), by = "chr") %>% + mutate(start_g = start + offset, end_g = end + offset) + } + + list( + bins = add_offset(bins) %>% mutate(pos_g = (start_g + end_g) / 2), + segs = add_offset(segs), + chr_lengths = chr_lengths %>% mutate(xmax = offset + len) + ) +} + +# Plot wisecondorx plot +plot_wisecondorx_cnv <- function(bins, segs, chr_lengths, sample_id, ratio_limit, output_dir) { + chr_ticks <- chr_lengths %>% mutate(mid = offset + len / 2) + chr_boundaries <- chr_lengths$xmax[-nrow(chr_lengths)] + + color_mapping <- c( + "NEUTRAL" = "#377eb8", + "GAIN" = "#e41a1c", + "LOSS" = "#4daf4a" + ) + + p <- ggplot() + + geom_vline(xintercept = chr_boundaries, color = "grey85", linewidth = 0.4, linetype = "dotted") + + geom_hline(yintercept = 0, color = "grey75", linewidth = 0.4) + + geom_point( + data = bins, + aes(x = pos_g, y = pmin(pmax(ratio, -ratio_limit), ratio_limit)), + color = "grey75", size = 0.25, alpha = 0.35 + ) + + geom_segment( + data = segs, + aes(x = start_g, xend = end_g, y = ratio, yend = ratio, color = call), + linewidth = 1.6, lineend = "round" + ) + + scale_color_manual( + values = color_mapping, + name = "Copy Number Call", + guide = guide_legend(override.aes = list(linewidth = 4)) + ) + + scale_x_continuous(breaks = chr_ticks$mid, labels = chr_ticks$chr, expand = c(0.01, 0.01)) + + scale_y_continuous(limits = c(-ratio_limit, ratio_limit)) + + labs( + x = "Chromosome", y = expression(log[2](ratio)), + title = paste0("Copy Number Profile", if (!is.null(sample_id)) paste0(" - ", sample_id) else "") + ) + + theme_minimal(base_size = 13) + + theme( + panel.grid.minor = element_blank(), + panel.grid.major.x = element_blank(), + panel.grid.major.y = element_line(color = "grey93", linewidth = 0.3), + panel.background = element_rect(fill = "white", color = NA), + plot.background = element_rect(fill = "white", color = NA), + axis.text = element_text(color = "grey30"), + axis.title = element_text(color = "grey20"), + plot.title = element_text(face = "bold"), + plot.subtitle = element_text(color = "grey40", size = 10), + legend.position = "bottom", + legend.title = element_text(face = "bold") + ) + + out_png <- file.path(output_dir, paste0(sample_id, ".copy_number.png")) + out_svg <- file.path(output_dir, paste0(sample_id, ".copy_number.svg")) + ggsave(out_png, plot = p, width = 14, height = 6, dpi = 300) + ggsave(out_svg, plot = p, width = 14, height = 6) + message("Plot (PNG) saved to: ", out_png) + message("Plot (SVG) saved to: ", out_svg) +} + +# ---- main ------------------------------------------------------------- + +message("Processing sample: ", sample_id) + +bins <- read_and_normalize(bin_path) %>% + filter(!is.na(ratio)) +segs <- read_and_normalize(seg_path) %>% + classify_segments() + +message("Data loaded - Bins: ", nrow(bins), " Segments: ", nrow(segs)) + +coords <- compute_genomic_coords(bins, segs) + +plot_wisecondorx_cnv(coords$bins, + coords$segs, + coords$chr_lengths, + sample_id, + ratio_limit, + output_dir) + +message("Plot generation completed!") \ No newline at end of file diff --git a/conf/modules.config b/conf/modules.config index 0245eb6..027fbe2 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -328,11 +328,10 @@ process { withName: PLOT_ICHORCNA { publishDir = [ - path: { "${params.outdir}/ichorcna/corrected_plots" }, + path: { "${params.outdir}/ichorcna/plot_cnv" }, mode: params.publish_dir_mode, saveAs: { filename -> filename.equals('versions.yml') ? null : filename }, ] - ext.when = params.ichorcna_ploidy_aware_plot } withName: FORMAT_ICHORCNA_SEG { @@ -402,6 +401,18 @@ process { ].join(' ') } + withName: PLOT_WISECONDORX_CNV { + publishDir = [ + path: { "${params.outdir}/wisecondorx/plot_cnv" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename }, + ] + ext.args = [ + params.wisecondorx_plot_y_max ? "--ratio_limit ${params.wisecondorx_plot_y_max}" : "" + ].join(' ') + + } + withName: CONVERT_GISTIC_SEG { publishDir = [ path: { "${params.outdir}/wisecondorx/${meta.id}" }, diff --git a/modules/local/plot_wisecondorx_cnv/main.nf b/modules/local/plot_wisecondorx_cnv/main.nf new file mode 100644 index 0000000..d6cb1fd --- /dev/null +++ b/modules/local/plot_wisecondorx_cnv/main.nf @@ -0,0 +1,43 @@ +process PLOT_WISECONDORX_CNV { + tag "Plotting WisecondorX results for $meta.id" + label 'process_low' + container "community.wave.seqera.io/library/procps-ng_r-argparser_r-dplyr_r-ggplot2_pruned:10da72fa04bcba1a" + + input: + tuple val(meta), path(seg_file) + tuple val(meta2), path(bins) + + output: + tuple val(meta), path("*.copy_number.png"), emit: plot_png + tuple val(meta), path("*.copy_number.svg"), emit: plot_svg + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def VERSION = '0.1' + """ + plot_wisecondorx_cnv.R \\ + --id ${prefix} \\ + --seg_file ${seg_file} \\ + --binfile ${bins} \\ + --outdir . \\ + ${args} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + plot_wisecondorx_cnv: $VERSION + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}.copy_number.png + touch ${prefix}.copy_number.svg + touch versions.yml + """ +} diff --git a/nextflow.config b/nextflow.config index d9a9fab..adf5bc4 100644 --- a/nextflow.config +++ b/nextflow.config @@ -97,7 +97,6 @@ params { ichorcna_reptime_wig = null // FIXME: Supply it with the pipeline ichorcna_centromere_file = null - ichorcna_ploidy_aware_plot = false // WisecondorX @@ -107,6 +106,7 @@ params { wisecondorx_ylim = null wisecondorx_zscore = 5 wisecondorx_blacklist = null + wisecondorx_plot_y_max = 1 // ASCAT.sc diff --git a/nextflow_schema.json b/nextflow_schema.json index 559dae3..b76cc33 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -537,10 +537,6 @@ "hidden": true, "description": "Path to a file with centromere data for the specified genome.", "format": "file-path" - }, - "ichorcna_ploidy_aware_plot": { - "type": "boolean", - "description": "Ploidy-aware ichorCNA plots." } } }, @@ -572,6 +568,11 @@ "wisecondorx_ylim": { "type": "string", "description": "Y axis limits for the generated plots" + }, + "wisecondorx_plot_y_max": { + "type": "number", + "default": 1, + "description": "Upper log2ratio limit of additional WisecondorX plots (application-generated plots are unaffected)" } } }, diff --git a/subworkflows/local/liquid_biopsy/main.nf b/subworkflows/local/liquid_biopsy/main.nf index c1753a0..cf5b21e 100644 --- a/subworkflows/local/liquid_biopsy/main.nf +++ b/subworkflows/local/liquid_biopsy/main.nf @@ -6,6 +6,7 @@ include { BAM_CNV_WISECONDORX } from '../../../subworkflows/nf-core/bam include { ICHORCNA } from '../../../subworkflows/local/ichorcna/main' include { ASSEMBLE_WISECONDORX_OUTPUTS } from '../../../modules/local/assemble_wisecondorx_outputs/main' include { CONVERT_WISECONDORX_IMAGES } from '../../../modules/local/convert_wisecondorx_images/main' +include { PLOT_WISECONDORX_CNV } from '../../../modules/local/plot_wisecondorx_cnv/main' workflow LIQUID_BIOPSY { take: @@ -108,6 +109,12 @@ workflow LIQUID_BIOPSY { genome_plot = CONVERT_WISECONDORX_IMAGES.out.genome_plot // For compatibility with workflow output corrected_gistic_file = gistic_file + + PLOT_WISECONDORX_CNV(CONVERT_GISTIC_SEG.out.segfile, + BAM_CNV_WISECONDORX.out.bins_bed) + ch_versions = ch_versions.mix(PLOT_WISECONDORX_CNV.out.versions) + + } else { error("Uknown / unsupported analysis type ${caller}")