From 5e0f888984287dcbcb4bd622d6b0309a97b99f17 Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Thu, 10 Sep 2026 09:48:26 +0100 Subject: [PATCH 01/10] Expose SAVANA BEDPE, purity/ploidy and allele counts for downstream annotation PAIRED_SAVANA and TUMORONLY_SAVANA only emitted the classified somatic VCF and the segmented copy number. Padfoot and ReConPlot also need the somatic BEDPE, the fitted purity/ploidy and the het-SNP allele counts, all of which the nf-core savana/classify, savana/cna and savana/to modules already produce. Add those three emits to both subworkflows and mix the paired and tumour-only streams into per-product channels in the main workflow. All are optional: absent without an SNP source, and cna/fitted_purity_ploidy absent when SAVANA finds no acceptable fit (No_fit_found_PARAMS.tsv). Co-Authored-By: Claude Fable 5.1 --- subworkflows/local/paired/paired_savana.nf | 7 +++++-- subworkflows/local/tumor_only/tumoronly_savana.nf | 7 +++++-- workflows/lrsomatic.nf | 15 ++++++++++++++- 3 files changed, 24 insertions(+), 5 deletions(-) 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/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/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index 37bd1c4c..e78b3f00 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -1112,7 +1112,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 +1201,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) From 95b1e1700a34fc0327dcf7584a9a0a56c99b918e Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Thu, 10 Sep 2026 09:48:26 +0100 Subject: [PATCH 02/10] Add Padfoot annotation and ReConPlot figures on top of the SAVANA subworkflows Padfoot (functional annotation of somatic SVs + CNAs) runs once per available caller pair: Severus + top-ranked Wakhan integer CN VCF, and SAVANA classified somatic SVs + SAVANA segmented absolute CN. ReConPlot (rearrangement + copy-number figures) runs per CN/SV pair: ASCAT + Severus, Wakhan + Severus, and SAVANA alone, producing per-chromosome and genome-wide figures, an optional region focus panel and the harmonised CN/SV tables. SAVANA inputs come from PAIRED_SAVANA / TUMORONLY_SAVANA (plain VCF from the nf-core modules; Padfoot reads it via pysam, the ReConPlot wrapper uses the BEDPE). Samples without a copy-number fit are dropped before the map that builds the ReConPlot input, since the remainder join yields a shorter tuple for them. Neither tool is on Bioconda: the pinned source trees are fetched with WGET+UNTAR (or taken from --padfoot_dir / --reconplot_dir / --reconplot_pkg_dir) and run in digest-pinned public images (ghcr.io/tim-yu/padfoot-repeatmasker with RepeatMasker 4.2.4 + Dfam 4.0, ghcr.io/tim-yu/reconplot); recipes under containers/. Both are skipped in the test profiles (whole-genome inputs, multi-GB images) and covered by module-level nf-tests instead. Parameter checks for unsupported genome presets live in validateSvAnnotationParams(). Co-Authored-By: Claude Fable 5.1 --- conf/modules.config | 69 ++++++ conf/test.config | 4 + conf/test_full.config | 4 + containers/padfoot/Dockerfile | 30 +++ containers/padfoot/README.md | 21 ++ containers/reconplot/Dockerfile | 19 ++ containers/reconplot/README.md | 21 ++ modules/local/padfoot/environment.yml | 18 ++ modules/local/padfoot/main.nf | 71 ++++++ modules/local/padfoot/meta.yml | 107 +++++++++ modules/local/padfoot/tests/main.nf.test | 74 ++++++ modules/local/padfoot/tests/main.nf.test.snap | 144 ++++++++++++ modules/local/padfoot/tests/nextflow.config | 6 + modules/local/reconplot/environment.yml | 18 ++ modules/local/reconplot/main.nf | 95 ++++++++ modules/local/reconplot/meta.yml | 103 +++++++++ modules/local/reconplot/tests/main.nf.test | 72 ++++++ .../local/reconplot/tests/main.nf.test.snap | 212 ++++++++++++++++++ modules/local/reconplot/tests/nextflow.config | 5 + nextflow.config | 24 ++ nextflow_schema.json | 121 ++++++++++ .../utils_nfcore_lrsomatic_pipeline/main.nf | 23 ++ tests/.nftignore | 3 + workflows/lrsomatic.nf | 194 ++++++++++++++++ 24 files changed, 1458 insertions(+) create mode 100644 containers/padfoot/Dockerfile create mode 100644 containers/padfoot/README.md create mode 100644 containers/reconplot/Dockerfile create mode 100644 containers/reconplot/README.md create mode 100644 modules/local/padfoot/environment.yml create mode 100644 modules/local/padfoot/main.nf create mode 100644 modules/local/padfoot/meta.yml create mode 100644 modules/local/padfoot/tests/main.nf.test create mode 100644 modules/local/padfoot/tests/main.nf.test.snap create mode 100644 modules/local/padfoot/tests/nextflow.config create mode 100644 modules/local/reconplot/environment.yml create mode 100644 modules/local/reconplot/main.nf create mode 100644 modules/local/reconplot/meta.yml create mode 100644 modules/local/reconplot/tests/main.nf.test create mode 100644 modules/local/reconplot/tests/main.nf.test.snap create mode 100644 modules/local/reconplot/tests/nextflow.config diff --git a/conf/modules.config b/conf/modules.config index 8b20e4fd..7043d271 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_(WGET|PKG_WGET)' { + ext.suffix = 'tar.gz' + publishDir = [ enabled: false ] + } + + withName : '.*:RECONPLOT_(UNTAR|PKG_UNTAR)' { + publishDir = [ enabled: false ] + } + + withName : '.*:RECONPLOT_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${cn_source}_${sv_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..adb63b3e 100644 --- a/conf/test.config +++ b/conf/test.config @@ -70,6 +70,10 @@ params { skip_wakhan = true skip_ascat = true skip_modkit = true + // Padfoot/ReConPlot need whole-genome caller output and fetch multi-GB images/annotations; + // covered by the module-level nf-tests instead. + 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..dab4ef7c 100644 --- a/conf/test_full.config +++ b/conf/test_full.config @@ -36,6 +36,10 @@ params { skip_wakhan = true skip_ascat = true skip_vep_plugins = true + // Padfoot/ReConPlot need whole-genome caller output and fetch multi-GB images/annotations; + // covered by the module-level nf-tests instead. + 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..776543e8 --- /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 are not +baked in; the pipeline stages them from `--reconplot_url` / `--reconplot_dir`. + +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_(ASCAT_SEVERUS|WAKHAN_SEVERUS|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/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..20624cb4 --- /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_(ASCAT_SEVERUS|WAKHAN_SEVERUS|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}" : "${cn_source}_${sv_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: ${params.reconplot_dir ? 'local checkout' : params.reconplot_url} + 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}" : "${cn_source}_${sv_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..3a5c3e4d --- /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: Tim-Yu/ReConPlot wrapper checkout (contains `run_reconplot.R` and `R/`) + - - 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..25fc5a08 --- /dev/null +++ b/modules/local/reconplot/tests/nextflow.config @@ -0,0 +1,5 @@ +process { + withName: 'RECONPLOT' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${cn_source}_${sv_source}" } + } +} diff --git a/nextflow.config b/nextflow.config index 696d5161..b19a4b8a 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,28 @@ 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_url = 'https://github.com/Tim-Yu/ReConPlot/archive/2672ccb22fc979909e136f6d2614b5fb763d1de2.tar.gz' // wrapper scripts + reconplot_dir = null // local wrapper checkout (offline); overrides reconplot_url + 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..872cd03f 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -376,6 +376,113 @@ } } }, + "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. Default: inferred from --genome (GRCh38 -> hg38, CHM13 -> chm13). Padfoot only bundles hg38/mm10 annotations; chm13 requires --padfoot_gff and --padfoot_rm." + }, + "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_url": { + "type": "string", + "default": "https://github.com/Tim-Yu/ReConPlot/archive/2672ccb22fc979909e136f6d2614b5fb763d1de2.tar.gz", + "description": "URL of the ReConPlot wrapper source tarball (GitHub archive of Tim-Yu/ReConPlot)." + }, + "reconplot_dir": { + "type": "string", + "format": "directory-path", + "exists": true, + "description": "Local checkout of the ReConPlot wrapper (contains run_reconplot.R). Overrides --reconplot_url." + }, + "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 +654,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 +886,12 @@ { "$ref": "#/$defs/savana_options" }, + { + "$ref": "#/$defs/padfoot_options" + }, + { + "$ref": "#/$defs/reconplot_options" + }, { "$ref": "#/$defs/ascat_parameters" }, 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 e78b3f00..eebbba24 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -46,6 +46,17 @@ include { MODKIT_PILEUP } from '../modules/nf-core/modkit/pi include { BCFTOOLS_VIEW as SIGNATURES_BCFTOOLS_VIEW } from '../modules/nf-core/bcftools/view/main' include { SIGPROFILER_MATRIXGENERATOR } from '../modules/local/sigprofiler/matrixgenerator/main' include { SIGPROFILER_ASSIGNMENT } from '../modules/local/sigprofiler/assignment/main' +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' +include { RECONPLOT as RECONPLOT_ASCAT_SEVERUS } from '../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_WAKHAN_SEVERUS } from '../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_SAVANA } from '../modules/local/reconplot/main' +include { WGET as RECONPLOT_WGET } from '../modules/nf-core/wget/main' +include { UNTAR as RECONPLOT_UNTAR } from '../modules/nf-core/untar/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' // // IMPORT SUBWORKFLOWS @@ -1275,6 +1286,189 @@ workflow LRSOMATIC { // ch_wakhan_files: [meta, [file_or_dir, ...]] } + // + // MODULE: PADFOOT (label: process_medium) + // Functional annotation of somatic SVs + CNAs. Run once per available SV/CNA caller pair: + // PADFOOT_SEVERUS_WAKHAN -- Severus somatic SVs + Wakhan best-solution integer CNA VCF + // PADFOOT_SAVANA -- SAVANA classified somatic SVs + SAVANA segmented absolute CN TSV + // Both paired and tumour-only samples are annotated. SAVANA CNA only exists when an SNP source was + // available (phased germline VCF for paired; 1000G panel for tumour-only on GRCh38/CHM13), so the + // join silently drops samples without CNA. + // SAVANA products come from PAIRED_SAVANA/TUMORONLY_SAVANA (nf-core savana modules, plain VCF). + // Padfoot is not on bioconda: the source tree is downloaded (params.padfoot_url) or taken from a + // local checkout (params.padfoot_dir); the container/conda env provide only its dependencies. + // Unsupported genomes (no bundled or user-supplied annotations) skip Padfoot with a warning. + // + + def padfoot_genome = params.padfoot_genome ?: + (params.genome == 'GRCh38' ? 'hg38' : params.genome == 'CHM13' ? 'chm13' : null) + // Padfoot only bundles hg38/mm10 annotations; unsupported genomes are reported by validateInputParameters() + 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_annot = [ + [:], + padfoot_genome, + params.padfoot_gff ? file(params.padfoot_gff, checkIfExists: true) : [], + params.padfoot_rm ? file(params.padfoot_rm, checkIfExists: true) : [] + ] + + if (params.padfoot_dir) { + padfoot_src = channel.value([[id: 'padfoot'], file(params.padfoot_dir, type: 'dir', checkIfExists: true)]) + } + else { + // MODULE: PADFOOT_WGET + PADFOOT_UNTAR -- fetch Padfoot source tarball (GitHub archive) + 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] -- directory containing padfoot.py and beds/ + + if (!params.skip_wakhan) { + // Wakhan writes every fitted solution; solution_1/ holds the top-ranked one + WAKHAN.out.vcf_files + .map { meta, vcfs -> + def files = [vcfs].flatten() + def integers = files.findAll { vcf -> vcf.name.endsWith('_wakhan_cna_integers.vcf') } + def best = integers.find { vcf -> vcf.toString().contains('/solution_1/') } ?: integers[0] + return [meta, best] + } + .filter { _meta, vcf -> vcf != null } + .set { wakhan_best_cna } + // wakhan_best_cna: [meta, wakhan_cna_integers.vcf] + + SEVERUS.out.somatic_vcf + .join(wakhan_best_cna) + .map { meta, sv, cna -> [meta, sv, 'severus', cna, 'wakhan'] } + .set { padfoot_severus_wakhan_input } + // padfoot_severus_wakhan_input: [meta, severus_somatic.vcf.gz, 'severus', wakhan_cna_integers.vcf, 'wakhan'] + + PADFOOT_SEVERUS_WAKHAN ( + padfoot_severus_wakhan_input, + ch_fasta, + ch_fai, + padfoot_src, + padfoot_annot + ) + ch_versions = ch_versions.mix(PADFOOT_SEVERUS_WAKHAN.out.versions) + } + + if (!params.skip_savana) { + savana_somatic_vcf + .join(savana_cna) + .map { meta, sv, cna -> [meta, sv, 'savana', cna, 'savana'] } + .set { padfoot_savana_input } + // padfoot_savana_input: [meta, classified.somatic.vcf.gz, 'savana', segmented_absolute_copy_number.tsv, 'savana'] + + PADFOOT_SAVANA ( + padfoot_savana_input, + ch_fasta, + ch_fai, + padfoot_src, + padfoot_annot + ) + ch_versions = ch_versions.mix(PADFOOT_SAVANA.out.versions) + } + } + + // + // MODULE: RECONPLOT (label: process_low) + // Rearrangement + copy-number figures (ReConPlot) for each available CN/SV caller pair, written to + // {outdir}/{sample}/reconplot/{ascat_severus,wakhan_severus,savana}/ with per-chromosome figures, a + // genome-wide strip and the harmonised CN/SV tables. Paired and tumour-only samples alike; a pair is + // only produced when both callers emitted output for the sample (join semantics). + // The wrapper (Tim-Yu/ReConPlot) and the ReConPlot R package are staged as source (WGET+UNTAR or + // local dirs); the container ships the package pre-installed, conda installs it at run time. + // + + if (!params.skip_reconplot) { + + def reconplot_genome = params.reconplot_genome ?: + (params.genome == 'CHM13' ? 'T2T' : 'hg38') + + if (params.reconplot_dir) { + reconplot_src = channel.value([[id: 'reconplot'], file(params.reconplot_dir, type: 'dir', checkIfExists: true)]) + } + else { + RECONPLOT_WGET( channel.value([[id: 'reconplot'], params.reconplot_url]) ) + RECONPLOT_UNTAR( RECONPLOT_WGET.out.outfile ) + reconplot_src = RECONPLOT_UNTAR.out.untar + ch_versions = ch_versions.mix(RECONPLOT_WGET.out.versions) + } + 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 (run_reconplot.R + R/) + // reconplot_pkg: [meta, dir] -- ReConPlot R package source + + // Severus somatic SVs are the SV component for the two lrsomatic CN callers + severus_sv_files = SEVERUS.out.somatic_vcf.map { meta, vcf -> [meta, [vcf]] } + // severus_sv_files: [meta, [severus_somatic.vcf.gz]] + + if (!params.skip_ascat) { + ASCAT.out.segments + .join(ASCAT.out.purityploidy) + .join(ASCAT.out.bafs) + // segments = fitted segments.txt (segments_raw.txt is a separate emit); bafs = every *BAF.txt, + // the wrapper picks .tumour_tumourBAF.txt by name + .map { meta, seg, pp, bafs -> [meta, [seg, pp, bafs].flatten()] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'ascat', cn, 'severus', sv] } + .set { reconplot_ascat_input } + // reconplot_ascat_input: [meta, 'ascat', [segments.txt, purityploidy.txt, *BAF.txt], 'severus', [vcf]] + + RECONPLOT_ASCAT_SEVERUS( reconplot_ascat_input, reconplot_src, reconplot_pkg, reconplot_genome ) + ch_versions = ch_versions.mix(RECONPLOT_ASCAT_SEVERUS.out.versions) + } + + if (!params.skip_wakhan) { + // Top-ranked Wakhan solution: allele-specific segment BEDs (HP1 + HP2) plus the ranking table + WAKHAN.out.bed_files + .map { meta, beds -> + def files = [beds].flatten() + def hp = files.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.out.solutions_ranks) + .map { meta, beds, ranks -> [meta, beds + [ranks]] } + .join(severus_sv_files) + .map { meta, cn, sv -> [meta, 'wakhan', cn, 'severus', sv] } + .set { reconplot_wakhan_input } + // reconplot_wakhan_input: [meta, 'wakhan', [HP_1.bed, HP_2.bed, solutions_ranks.tsv], 'severus', [vcf]] + + RECONPLOT_WAKHAN_SEVERUS( reconplot_wakhan_input, reconplot_src, reconplot_pkg, reconplot_genome ) + ch_versions = ch_versions.mix(RECONPLOT_WAKHAN_SEVERUS.out.versions) + } + + if (!params.skip_savana) { + // Single-source mode: all SAVANA files in cn_files, sv_files empty. + // allele_counts is optional (absent without an SNP source), so join with remainder and drop nulls. + // A sample with allele counts but no CN fit (SAVANA "No_fit_found") only exists on the right-hand + // side and surfaces as [meta, null, bed]; there is nothing to plot for it, so drop it before the map. + savana_cna + .join(savana_somatic_bedpe) + .join(savana_fitted_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 { reconplot_savana_input } + // reconplot_savana_input: [meta, 'savana', [segmented_absolute_copy_number.tsv, classified.somatic.bedpe, fitted_purity_ploidy.tsv, allele_counts_hetSNPs.bed], 'savana', []] + + RECONPLOT_SAVANA( reconplot_savana_input, reconplot_src, reconplot_pkg, reconplot_genome ) + ch_versions = ch_versions.mix(RECONPLOT_SAVANA.out.versions) + } + } + // // MODULE: LRSOMATICREPORT -- per-sample HTML report; all inputs optional, so joins use remainder: true on the tumor id // From 5bc9dc5dcf6e7dd133b43327d4243dd12bfaf64e Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Thu, 10 Sep 2026 09:48:26 +0100 Subject: [PATCH 03/10] Document Padfoot and ReConPlot across README, CITATIONS, CHANGELOG and docs Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 2 ++ CITATIONS.md | 24 ++++++++++++++++ README.md | 6 ++++ docs/output.md | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++ docs/usage.md | 68 +++++++++++++++++++++++++++++++++++++++++--- 5 files changed, 173 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f6e2a92..af8a32db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,8 @@ 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). +- Added Padfoot annotation of somatic SVs + CNAs for both Severus/Wakhan and SAVANA (#189) outputs (paired and tumour-only). Padfoot source is fetched from GitHub (or `--padfoot_dir`) and run in a public image bundling its dependencies with RepeatMasker 4.2.4 + Dfam 4.0 (`ghcr.io/tim-yu/padfoot-repeatmasker`, recipe in `containers/padfoot/`) or in a conda env; RepeatMasker annotation of inserted sequences runs by default. New params `skip_padfoot`, `padfoot_url`, `padfoot_dir`, `padfoot_genome`, `padfoot_gff`, `padfoot_rm`, `padfoot_run_repeatmasker` (@Tim-Yu). +- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{ascat_severus,wakhan_severus,savana}/`): per-chromosome, genome-wide and optional region-focus panels. Wrapper (`Tim-Yu/ReConPlot`) and ReConPlot R package are staged from GitHub or local checkouts; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_url`, `reconplot_dir`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@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..a071d6bb 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)) + diff --git a/docs/output.md b/docs/output.md index 495ab48c..077b6965 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 +│ │ ├── wakhan_severus +│ │ └── 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 +│ │ ├── ascat_severus +│ │ ├── wakhan_severus +│ │ └── 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`/`wakhan_severus` need `--skip_wakhan false`, `ascat_severus` 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 +│ ├── ascat_severus +│ │ ├── 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 +│ ├── wakhan_severus +│ │ └── (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) | + +`ascat_severus/` and `wakhan_severus/` 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..a914ae3a 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,66 @@ 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) 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 [Tim-Yu/ReConPlot](https://github.com/Tim-Yu/ReConPlot) wrapper for every CN/SV caller pair available for a sample, into `reconplot//`: + +- `ascat_severus/` -- ASCAT allele-specific CN + Severus somatic SVs +- `wakhan_severus/` -- 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_url` | Wrapper source tarball (GitHub archive). Default = pinned commit of `Tim-Yu/ReConPlot` | +| `--reconplot_dir` | Local wrapper checkout (contains `run_reconplot.R`); overrides `--reconplot_url`. Default = `null` | +| `--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_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { container = '/path/to/reconplot.sif' } +} +``` + +The Padfoot and ReConPlot **source trees** are downloaded from GitHub at run time; use `--padfoot_dir`, `--reconplot_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. @@ -493,8 +555,7 @@ Prepared REVEL and EVE files are published to `/vep_plugins/`, alongside a later run can point `--vep_revel` / `--vep_revel_tbi` (or the `--vep_eve` pair) at them and skip both the download and the reshaping. The two AlphaMissense tables are fetched already indexed. -> [!IMPORTANT] -> `/vep_plugins/` holds data that is free for non-commercial use only. Exclude it when you +> [!IMPORTANT] > `/vep_plugins/` holds data that is free for non-commercial use only. Exclude it when you > share or archive a results directory — passing it on is redistribution, which those licences do > not grant you. @@ -597,8 +658,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. From e667aae9072dc8b0ea5a3527f356c37e49653a9b Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Tue, 22 Sep 2026 19:34:35 +0100 Subject: [PATCH 04/10] Move the Padfoot and ReConPlot wiring into local subworkflows Nextflow stores each workflow body as a single string constant, which the JVM caps at 65,535 characters. After merging dev (report, signatures and Verdict blocks) the inline Padfoot and ReConPlot wiring pushed LRSOMATIC's body to 66,936 characters and the pipeline failed to compile ("String too long"). PADFOOT_ANNOTATION and RECONPLOT_FIGURES now hold the source staging, the caller-pair joins and the module calls, mirroring PAIRED_SAVANA; the main workflow passes channel.empty() for skipped callers and is back to 59,214 characters. No change to process inputs, published paths or the conf/modules.config selectors. Co-Authored-By: Claude Fable 5.1 --- subworkflows/local/padfoot_annotation.nf | 75 ++++++++ subworkflows/local/reconplot_figures.nf | 113 ++++++++++++ workflows/lrsomatic.nf | 212 ++++------------------- 3 files changed, 220 insertions(+), 180 deletions(-) create mode 100644 subworkflows/local/padfoot_annotation.nf create mode 100644 subworkflows/local/reconplot_figures.nf 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/reconplot_figures.nf b/subworkflows/local/reconplot_figures.nf new file mode 100644 index 00000000..ed2e003d --- /dev/null +++ b/subworkflows/local/reconplot_figures.nf @@ -0,0 +1,113 @@ +// IMPORT MODULES +include { RECONPLOT as RECONPLOT_ASCAT_SEVERUS } from '../../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_WAKHAN_SEVERUS } from '../../modules/local/reconplot/main' +include { RECONPLOT as RECONPLOT_SAVANA } from '../../modules/local/reconplot/main' +include { WGET as RECONPLOT_WGET } from '../../modules/nf-core/wget/main' +include { UNTAR as RECONPLOT_UNTAR } from '../../modules/nf-core/untar/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 (run_reconplot.R + R/) and the ReConPlot R package are staged +// as source from params.reconplot_url / params.reconplot_pkg_url, or local checkouts via the *_dir params. +// +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() + + if (params.reconplot_dir) { + reconplot_src = channel.value([[id: 'reconplot'], file(params.reconplot_dir, type: 'dir', checkIfExists: true)]) + } + else { + RECONPLOT_WGET( channel.value([[id: 'reconplot'], params.reconplot_url]) ) + RECONPLOT_UNTAR( RECONPLOT_WGET.out.outfile ) + reconplot_src = RECONPLOT_UNTAR.out.untar + ch_versions = ch_versions.mix(RECONPLOT_WGET.out.versions) + } + 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_ASCAT_SEVERUS (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_ASCAT_SEVERUS( ascat_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_ASCAT_SEVERUS.out.versions) + + // + // MODULE: RECONPLOT_WAKHAN_SEVERUS (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_WAKHAN_SEVERUS( wakhan_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_WAKHAN_SEVERUS.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/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index eebbba24..04c3e558 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -46,17 +46,6 @@ include { MODKIT_PILEUP } from '../modules/nf-core/modkit/pi include { BCFTOOLS_VIEW as SIGNATURES_BCFTOOLS_VIEW } from '../modules/nf-core/bcftools/view/main' include { SIGPROFILER_MATRIXGENERATOR } from '../modules/local/sigprofiler/matrixgenerator/main' include { SIGPROFILER_ASSIGNMENT } from '../modules/local/sigprofiler/assignment/main' -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' -include { RECONPLOT as RECONPLOT_ASCAT_SEVERUS } from '../modules/local/reconplot/main' -include { RECONPLOT as RECONPLOT_WAKHAN_SEVERUS } from '../modules/local/reconplot/main' -include { RECONPLOT as RECONPLOT_SAVANA } from '../modules/local/reconplot/main' -include { WGET as RECONPLOT_WGET } from '../modules/nf-core/wget/main' -include { UNTAR as RECONPLOT_UNTAR } from '../modules/nf-core/untar/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' // // IMPORT SUBWORKFLOWS @@ -71,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' @@ -1287,186 +1278,47 @@ workflow LRSOMATIC { } // - // MODULE: PADFOOT (label: process_medium) - // Functional annotation of somatic SVs + CNAs. Run once per available SV/CNA caller pair: - // PADFOOT_SEVERUS_WAKHAN -- Severus somatic SVs + Wakhan best-solution integer CNA VCF - // PADFOOT_SAVANA -- SAVANA classified somatic SVs + SAVANA segmented absolute CN TSV - // Both paired and tumour-only samples are annotated. SAVANA CNA only exists when an SNP source was - // available (phased germline VCF for paired; 1000G panel for tumour-only on GRCh38/CHM13), so the - // join silently drops samples without CNA. - // SAVANA products come from PAIRED_SAVANA/TUMORONLY_SAVANA (nf-core savana modules, plain VCF). - // Padfoot is not on bioconda: the source tree is downloaded (params.padfoot_url) or taken from a - // local checkout (params.padfoot_dir); the container/conda env provide only its dependencies. - // Unsupported genomes (no bundled or user-supplied annotations) skip Padfoot with a warning. + // 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) - // Padfoot only bundles hg38/mm10 annotations; unsupported genomes are reported by validateInputParameters() 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_annot = [ - [:], - padfoot_genome, - params.padfoot_gff ? file(params.padfoot_gff, checkIfExists: true) : [], - params.padfoot_rm ? file(params.padfoot_rm, checkIfExists: true) : [] - ] - - if (params.padfoot_dir) { - padfoot_src = channel.value([[id: 'padfoot'], file(params.padfoot_dir, type: 'dir', checkIfExists: true)]) - } - else { - // MODULE: PADFOOT_WGET + PADFOOT_UNTAR -- fetch Padfoot source tarball (GitHub archive) - 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] -- directory containing padfoot.py and beds/ - - if (!params.skip_wakhan) { - // Wakhan writes every fitted solution; solution_1/ holds the top-ranked one - WAKHAN.out.vcf_files - .map { meta, vcfs -> - def files = [vcfs].flatten() - def integers = files.findAll { vcf -> vcf.name.endsWith('_wakhan_cna_integers.vcf') } - def best = integers.find { vcf -> vcf.toString().contains('/solution_1/') } ?: integers[0] - return [meta, best] - } - .filter { _meta, vcf -> vcf != null } - .set { wakhan_best_cna } - // wakhan_best_cna: [meta, wakhan_cna_integers.vcf] - - SEVERUS.out.somatic_vcf - .join(wakhan_best_cna) - .map { meta, sv, cna -> [meta, sv, 'severus', cna, 'wakhan'] } - .set { padfoot_severus_wakhan_input } - // padfoot_severus_wakhan_input: [meta, severus_somatic.vcf.gz, 'severus', wakhan_cna_integers.vcf, 'wakhan'] - - PADFOOT_SEVERUS_WAKHAN ( - padfoot_severus_wakhan_input, - ch_fasta, - ch_fai, - padfoot_src, - padfoot_annot - ) - ch_versions = ch_versions.mix(PADFOOT_SEVERUS_WAKHAN.out.versions) - } - - if (!params.skip_savana) { - savana_somatic_vcf - .join(savana_cna) - .map { meta, sv, cna -> [meta, sv, 'savana', cna, 'savana'] } - .set { padfoot_savana_input } - // padfoot_savana_input: [meta, classified.somatic.vcf.gz, 'savana', segmented_absolute_copy_number.tsv, 'savana'] - - PADFOOT_SAVANA ( - padfoot_savana_input, - ch_fasta, - ch_fai, - padfoot_src, - padfoot_annot - ) - ch_versions = ch_versions.mix(PADFOOT_SAVANA.out.versions) - } + 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) } // - // MODULE: RECONPLOT (label: process_low) - // Rearrangement + copy-number figures (ReConPlot) for each available CN/SV caller pair, written to - // {outdir}/{sample}/reconplot/{ascat_severus,wakhan_severus,savana}/ with per-chromosome figures, a - // genome-wide strip and the harmonised CN/SV tables. Paired and tumour-only samples alike; a pair is - // only produced when both callers emitted output for the sample (join semantics). - // The wrapper (Tim-Yu/ReConPlot) and the ReConPlot R package are staged as source (WGET+UNTAR or - // local dirs); the container ships the package pre-installed, conda installs it at run time. + // SUBWORKFLOW: RECONPLOT_FIGURES -- ReConPlot figures per CN/SV caller pair (ASCAT + Severus, Wakhan + Severus, SAVANA) // - if (!params.skip_reconplot) { - - def reconplot_genome = params.reconplot_genome ?: - (params.genome == 'CHM13' ? 'T2T' : 'hg38') - - if (params.reconplot_dir) { - reconplot_src = channel.value([[id: 'reconplot'], file(params.reconplot_dir, type: 'dir', checkIfExists: true)]) - } - else { - RECONPLOT_WGET( channel.value([[id: 'reconplot'], params.reconplot_url]) ) - RECONPLOT_UNTAR( RECONPLOT_WGET.out.outfile ) - reconplot_src = RECONPLOT_UNTAR.out.untar - ch_versions = ch_versions.mix(RECONPLOT_WGET.out.versions) - } - 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 (run_reconplot.R + R/) - // reconplot_pkg: [meta, dir] -- ReConPlot R package source - - // Severus somatic SVs are the SV component for the two lrsomatic CN callers - severus_sv_files = SEVERUS.out.somatic_vcf.map { meta, vcf -> [meta, [vcf]] } - // severus_sv_files: [meta, [severus_somatic.vcf.gz]] - - if (!params.skip_ascat) { - ASCAT.out.segments - .join(ASCAT.out.purityploidy) - .join(ASCAT.out.bafs) - // segments = fitted segments.txt (segments_raw.txt is a separate emit); bafs = every *BAF.txt, - // the wrapper picks .tumour_tumourBAF.txt by name - .map { meta, seg, pp, bafs -> [meta, [seg, pp, bafs].flatten()] } - .join(severus_sv_files) - .map { meta, cn, sv -> [meta, 'ascat', cn, 'severus', sv] } - .set { reconplot_ascat_input } - // reconplot_ascat_input: [meta, 'ascat', [segments.txt, purityploidy.txt, *BAF.txt], 'severus', [vcf]] - - RECONPLOT_ASCAT_SEVERUS( reconplot_ascat_input, reconplot_src, reconplot_pkg, reconplot_genome ) - ch_versions = ch_versions.mix(RECONPLOT_ASCAT_SEVERUS.out.versions) - } - - if (!params.skip_wakhan) { - // Top-ranked Wakhan solution: allele-specific segment BEDs (HP1 + HP2) plus the ranking table - WAKHAN.out.bed_files - .map { meta, beds -> - def files = [beds].flatten() - def hp = files.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.out.solutions_ranks) - .map { meta, beds, ranks -> [meta, beds + [ranks]] } - .join(severus_sv_files) - .map { meta, cn, sv -> [meta, 'wakhan', cn, 'severus', sv] } - .set { reconplot_wakhan_input } - // reconplot_wakhan_input: [meta, 'wakhan', [HP_1.bed, HP_2.bed, solutions_ranks.tsv], 'severus', [vcf]] - - RECONPLOT_WAKHAN_SEVERUS( reconplot_wakhan_input, reconplot_src, reconplot_pkg, reconplot_genome ) - ch_versions = ch_versions.mix(RECONPLOT_WAKHAN_SEVERUS.out.versions) - } - - if (!params.skip_savana) { - // Single-source mode: all SAVANA files in cn_files, sv_files empty. - // allele_counts is optional (absent without an SNP source), so join with remainder and drop nulls. - // A sample with allele counts but no CN fit (SAVANA "No_fit_found") only exists on the right-hand - // side and surfaces as [meta, null, bed]; there is nothing to plot for it, so drop it before the map. - savana_cna - .join(savana_somatic_bedpe) - .join(savana_fitted_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 { reconplot_savana_input } - // reconplot_savana_input: [meta, 'savana', [segmented_absolute_copy_number.tsv, classified.somatic.bedpe, fitted_purity_ploidy.tsv, allele_counts_hetSNPs.bed], 'savana', []] - - RECONPLOT_SAVANA( reconplot_savana_input, reconplot_src, reconplot_pkg, reconplot_genome ) - ch_versions = ch_versions.mix(RECONPLOT_SAVANA.out.versions) - } + 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) } // From 5bbcc20d79874ff6e92c4ad83d5a9c8590a00889 Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Tue, 22 Sep 2026 19:38:51 +0100 Subject: [PATCH 05/10] Ship the ReConPlot wrapper under assets/ instead of downloading it at run time run_reconplot.R and R/ are pipeline glue (parsers for ASCAT, Wakhan, Severus and SAVANA output), so they now live in assets/reconplot, vendored from Tim-Yu/ReConPlot at the commit recorded in assets/reconplot/VERSION. RECONPLOT_FIGURES stages that directory directly; RECONPLOT_WGET and RECONPLOT_UNTAR and the reconplot_url / reconplot_dir parameters are removed. The ReConPlot R package itself is still fetched from its pinned upstream commit (reconplot_pkg_url / reconplot_pkg_dir) because -profile conda installs it at run time. versions.yml records the wrapper commit from VERSION. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 2 +- assets/reconplot/R/parsers/ascat.R | 192 +++++++++++++++++ assets/reconplot/R/parsers/generic.R | 108 ++++++++++ assets/reconplot/R/parsers/lrsomatic.R | 77 +++++++ assets/reconplot/R/parsers/registry.R | 104 +++++++++ assets/reconplot/R/parsers/savana.R | 258 +++++++++++++++++++++++ assets/reconplot/R/parsers/severus.R | 189 +++++++++++++++++ assets/reconplot/R/parsers/wakhan.R | 163 +++++++++++++++ assets/reconplot/R/plotting.R | 114 ++++++++++ assets/reconplot/R/utils_common.R | 267 ++++++++++++++++++++++++ assets/reconplot/R/utils_vcf.R | 95 +++++++++ assets/reconplot/README.md | 8 + assets/reconplot/VERSION | 1 + assets/reconplot/run_reconplot.R | 217 +++++++++++++++++++ conf/modules.config | 4 +- containers/reconplot/README.md | 4 +- docs/usage.md | 6 +- modules/local/reconplot/main.nf | 2 +- modules/local/reconplot/meta.yml | 2 +- nextflow.config | 2 - nextflow_schema.json | 11 - subworkflows/local/reconplot_figures.nf | 16 +- 22 files changed, 1805 insertions(+), 37 deletions(-) create mode 100644 assets/reconplot/R/parsers/ascat.R create mode 100644 assets/reconplot/R/parsers/generic.R create mode 100644 assets/reconplot/R/parsers/lrsomatic.R create mode 100644 assets/reconplot/R/parsers/registry.R create mode 100644 assets/reconplot/R/parsers/savana.R create mode 100644 assets/reconplot/R/parsers/severus.R create mode 100644 assets/reconplot/R/parsers/wakhan.R create mode 100644 assets/reconplot/R/plotting.R create mode 100644 assets/reconplot/R/utils_common.R create mode 100644 assets/reconplot/R/utils_vcf.R create mode 100644 assets/reconplot/README.md create mode 100644 assets/reconplot/VERSION create mode 100755 assets/reconplot/run_reconplot.R diff --git a/CHANGELOG.md b/CHANGELOG.md index af8a32db..c1fc7215 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,7 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - [#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). - Added Padfoot annotation of somatic SVs + CNAs for both Severus/Wakhan and SAVANA (#189) outputs (paired and tumour-only). Padfoot source is fetched from GitHub (or `--padfoot_dir`) and run in a public image bundling its dependencies with RepeatMasker 4.2.4 + Dfam 4.0 (`ghcr.io/tim-yu/padfoot-repeatmasker`, recipe in `containers/padfoot/`) or in a conda env; RepeatMasker annotation of inserted sequences runs by default. New params `skip_padfoot`, `padfoot_url`, `padfoot_dir`, `padfoot_genome`, `padfoot_gff`, `padfoot_rm`, `padfoot_run_repeatmasker` (@Tim-Yu). -- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{ascat_severus,wakhan_severus,savana}/`): per-chromosome, genome-wide and optional region-focus panels. Wrapper (`Tim-Yu/ReConPlot`) and ReConPlot R package are staged from GitHub or local checkouts; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_url`, `reconplot_dir`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@Tim-Yu). +- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{ascat_severus,wakhan_severus,savana}/`): per-chromosome, genome-wide and optional region-focus panels. The wrapper ships in `assets/reconplot/`; the ReConPlot R package is staged from GitHub or a local checkout; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@Tim-Yu). ### `Changed` 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 7043d271..2b46a27b 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -759,12 +759,12 @@ process { // Visualisation of SVs / CNAs // - withName : '.*:RECONPLOT_(WGET|PKG_WGET)' { + withName : '.*:RECONPLOT_PKG_WGET' { ext.suffix = 'tar.gz' publishDir = [ enabled: false ] } - withName : '.*:RECONPLOT_(UNTAR|PKG_UNTAR)' { + withName : '.*:RECONPLOT_PKG_UNTAR' { publishDir = [ enabled: false ] } diff --git a/containers/reconplot/README.md b/containers/reconplot/README.md index 776543e8..febc7cef 100644 --- a/containers/reconplot/README.md +++ b/containers/reconplot/README.md @@ -2,8 +2,8 @@ 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 are not -baked in; the pipeline stages them from `--reconplot_url` / `--reconplot_dir`. +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: diff --git a/docs/usage.md b/docs/usage.md index a914ae3a..f24d997c 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -404,7 +404,7 @@ Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. For othe #### ReConPlot Options -[ReConPlot](https://github.com/cortes-ciriano-lab/ReConPlot) rearrangement + copy-number figures are generated through the [Tim-Yu/ReConPlot](https://github.com/Tim-Yu/ReConPlot) wrapper for every CN/SV caller pair available for a sample, into `reconplot//`: +[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//`: - `ascat_severus/` -- ASCAT allele-specific CN + Severus somatic SVs - `wakhan_severus/` -- Wakhan top-ranked solution CN + Severus somatic SVs @@ -414,8 +414,6 @@ Each pair produces `per_chromosome/` (one figure per chromosome), `genome_wide/` | Parameter | Description | | -------------------------- | ------------------------------------------------------------------------------------------------------------- | -| `--reconplot_url` | Wrapper source tarball (GitHub archive). Default = pinned commit of `Tim-Yu/ReConPlot` | -| `--reconplot_dir` | Local wrapper checkout (contains `run_reconplot.R`); overrides `--reconplot_url`. Default = `null` | | `--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`) | @@ -438,7 +436,7 @@ process { } ``` -The Padfoot and ReConPlot **source trees** are downloaded from GitHub at run time; use `--padfoot_dir`, `--reconplot_dir` and `--reconplot_pkg_dir` to point at local checkouts instead. +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 diff --git a/modules/local/reconplot/main.nf b/modules/local/reconplot/main.nf index 20624cb4..7abf0c51 100644 --- a/modules/local/reconplot/main.nf +++ b/modules/local/reconplot/main.nf @@ -70,7 +70,7 @@ process RECONPLOT { cat <<-END_VERSIONS > versions.yml "${task.process}": reconplot: \$(Rscript -e 'cat(as.character(packageVersion("ReConPlot")))' 2>/dev/null) - reconplot_wrapper: ${params.reconplot_dir ? 'local checkout' : params.reconplot_url} + 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 diff --git a/modules/local/reconplot/meta.yml b/modules/local/reconplot/meta.yml index 3a5c3e4d..7c080ed2 100644 --- a/modules/local/reconplot/meta.yml +++ b/modules/local/reconplot/meta.yml @@ -39,7 +39,7 @@ input: description: Groovy Map for the wrapper source - reconplot_src: type: directory - description: Tim-Yu/ReConPlot wrapper checkout (contains `run_reconplot.R` and `R/`) + description: ReConPlot wrapper directory (`assets/reconplot`, contains `run_reconplot.R`, `R/` and `VERSION`) - - meta3: type: map description: Groovy Map for the R package source diff --git a/nextflow.config b/nextflow.config index b19a4b8a..a507df8c 100644 --- a/nextflow.config +++ b/nextflow.config @@ -124,8 +124,6 @@ params { padfoot_run_repeatmasker = true // containers ship full Dfam 4.0; conda uses RepeatMasker's bundled curated Dfam subset // ReConPlot options - reconplot_url = 'https://github.com/Tim-Yu/ReConPlot/archive/2672ccb22fc979909e136f6d2614b5fb763d1de2.tar.gz' // wrapper scripts - reconplot_dir = null // local wrapper checkout (offline); overrides reconplot_url 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 diff --git a/nextflow_schema.json b/nextflow_schema.json index 872cd03f..b6f29b90 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -423,17 +423,6 @@ "description": "Options for ReConPlot rearrangement/copy-number figures", "default": "", "properties": { - "reconplot_url": { - "type": "string", - "default": "https://github.com/Tim-Yu/ReConPlot/archive/2672ccb22fc979909e136f6d2614b5fb763d1de2.tar.gz", - "description": "URL of the ReConPlot wrapper source tarball (GitHub archive of Tim-Yu/ReConPlot)." - }, - "reconplot_dir": { - "type": "string", - "format": "directory-path", - "exists": true, - "description": "Local checkout of the ReConPlot wrapper (contains run_reconplot.R). Overrides --reconplot_url." - }, "reconplot_pkg_url": { "type": "string", "default": "https://github.com/cortes-ciriano-lab/ReConPlot/archive/88bb76011703d7908aaa2e5bdecfe18ba5737db8.tar.gz", diff --git a/subworkflows/local/reconplot_figures.nf b/subworkflows/local/reconplot_figures.nf index ed2e003d..cb72c6bc 100644 --- a/subworkflows/local/reconplot_figures.nf +++ b/subworkflows/local/reconplot_figures.nf @@ -2,16 +2,14 @@ include { RECONPLOT as RECONPLOT_ASCAT_SEVERUS } from '../../modules/local/reconplot/main' include { RECONPLOT as RECONPLOT_WAKHAN_SEVERUS } from '../../modules/local/reconplot/main' include { RECONPLOT as RECONPLOT_SAVANA } from '../../modules/local/reconplot/main' -include { WGET as RECONPLOT_WGET } from '../../modules/nf-core/wget/main' -include { UNTAR as RECONPLOT_UNTAR } from '../../modules/nf-core/untar/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 (run_reconplot.R + R/) and the ReConPlot R package are staged -// as source from params.reconplot_url / params.reconplot_pkg_url, or local checkouts via the *_dir params. +// 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 { @@ -31,15 +29,7 @@ workflow RECONPLOT_FIGURES { main: ch_versions = channel.empty() - if (params.reconplot_dir) { - reconplot_src = channel.value([[id: 'reconplot'], file(params.reconplot_dir, type: 'dir', checkIfExists: true)]) - } - else { - RECONPLOT_WGET( channel.value([[id: 'reconplot'], params.reconplot_url]) ) - RECONPLOT_UNTAR( RECONPLOT_WGET.out.outfile ) - reconplot_src = RECONPLOT_UNTAR.out.untar - ch_versions = ch_versions.mix(RECONPLOT_WGET.out.versions) - } + 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)]) } From e19467d4ca1ef2a9c50dabd9e0e193c281a70086 Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Tue, 22 Sep 2026 19:41:10 +0100 Subject: [PATCH 06/10] Name ReConPlot caller-pair folders SV caller first, matching Padfoot Padfoot published severus_wakhan/ while ReConPlot published wakhan_severus/ for the same pair. Both now use _: ReConPlot writes reconplot/{severus_ascat,severus_wakhan,savana}/ and the process aliases are RECONPLOT_SEVERUS_ASCAT and RECONPLOT_SEVERUS_WAKHAN. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 2 +- conf/modules.config | 4 ++-- containers/reconplot/README.md | 2 +- docs/output.md | 14 +++++++------- docs/usage.md | 6 +++--- modules/local/reconplot/main.nf | 6 +++--- modules/local/reconplot/tests/nextflow.config | 2 +- subworkflows/local/reconplot_figures.nf | 16 ++++++++-------- 8 files changed, 26 insertions(+), 26 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c1fc7215..ff387c33 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,7 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - [#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). - Added Padfoot annotation of somatic SVs + CNAs for both Severus/Wakhan and SAVANA (#189) outputs (paired and tumour-only). Padfoot source is fetched from GitHub (or `--padfoot_dir`) and run in a public image bundling its dependencies with RepeatMasker 4.2.4 + Dfam 4.0 (`ghcr.io/tim-yu/padfoot-repeatmasker`, recipe in `containers/padfoot/`) or in a conda env; RepeatMasker annotation of inserted sequences runs by default. New params `skip_padfoot`, `padfoot_url`, `padfoot_dir`, `padfoot_genome`, `padfoot_gff`, `padfoot_rm`, `padfoot_run_repeatmasker` (@Tim-Yu). -- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{ascat_severus,wakhan_severus,savana}/`): per-chromosome, genome-wide and optional region-focus panels. The wrapper ships in `assets/reconplot/`; the ReConPlot R package is staged from GitHub or a local checkout; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@Tim-Yu). +- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{severus_ascat,severus_wakhan,savana}/`): per-chromosome, genome-wide and optional region-focus panels. The wrapper ships in `assets/reconplot/`; the ReConPlot R package is staged from GitHub or a local checkout; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@Tim-Yu). ### `Changed` diff --git a/conf/modules.config b/conf/modules.config index 2b46a27b..0acbfc14 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -768,8 +768,8 @@ process { publishDir = [ enabled: false ] } - withName : '.*:RECONPLOT_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { - ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${cn_source}_${sv_source}" } + withName : '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } // shared filters / rendering ext.args = { [ diff --git a/containers/reconplot/README.md b/containers/reconplot/README.md index febc7cef..eaec0a6c 100644 --- a/containers/reconplot/README.md +++ b/containers/reconplot/README.md @@ -14,7 +14,7 @@ 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_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { container = ... } }`). The module currently pins +`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 diff --git a/docs/output.md b/docs/output.md index 077b6965..0dc949d4 100644 --- a/docs/output.md +++ b/docs/output.md @@ -33,7 +33,7 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ │ └── samtools │ │ └── whatshap_stats │ ├── reconplot -│ │ ├── wakhan_severus +│ │ ├── severus_wakhan │ │ └── savana │ ├── signatures │ │ ├── assignment @@ -85,8 +85,8 @@ The pipeline produces per-sample output directories. Two modes exist depending o │ │ │ └── samtools │ │ └── whatshap_stats │ ├── reconplot -│ │ ├── ascat_severus -│ │ ├── wakhan_severus +│ │ ├── severus_ascat +│ │ ├── severus_wakhan │ │ └── savana │ ├── signatures │ │ ├── assignment @@ -109,7 +109,7 @@ 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`/`wakhan_severus` need `--skip_wakhan false`, `ascat_severus` 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. +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`/`severus_wakhan` need `--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` @@ -645,14 +645,14 @@ Mutational signature analysis of the PASS SNVs and indels in the phased somatic ``` ├── reconplot -│ ├── ascat_severus +│ ├── 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 -│ ├── wakhan_severus +│ ├── severus_wakhan │ │ └── (same layout) │ └── savana │ └── (same layout) @@ -667,7 +667,7 @@ Mutational signature analysis of the PASS SNVs and indels in the phased somatic | `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) | -`ascat_severus/` and `wakhan_severus/` 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. +`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.
diff --git a/docs/usage.md b/docs/usage.md index f24d997c..8b6f1b7c 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -406,8 +406,8 @@ Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. For othe [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//`: -- `ascat_severus/` -- ASCAT allele-specific CN + Severus somatic SVs -- `wakhan_severus/` -- Wakhan top-ranked solution CN + Severus somatic SVs +- `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. @@ -432,7 +432,7 @@ Both images are pinned by digest directly in the module `container` directives, ```groovy process { withName: '.*:PADFOOT_(SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/padfoot-repeatmasker.sif' } - withName: '.*:RECONPLOT_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { container = '/path/to/reconplot.sif' } + withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = '/path/to/reconplot.sif' } } ``` diff --git a/modules/local/reconplot/main.nf b/modules/local/reconplot/main.nf index 7abf0c51..467c54a6 100644 --- a/modules/local/reconplot/main.nf +++ b/modules/local/reconplot/main.nf @@ -4,7 +4,7 @@ process RECONPLOT { 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_(ASCAT_SEVERUS|WAKHAN_SEVERUS|SAVANA)' { container = ... } }`. + // Override per site with `process { withName: '.*:RECONPLOT_(SEVERUS_ASCAT|SEVERUS_WAKHAN|SAVANA)' { container = ... } }`. container "ghcr.io/tim-yu/reconplot@sha256:1145fc5aebe0227bec371f4c59b08b9a09871498e403c01b83f83973149ae9e7" input: @@ -30,7 +30,7 @@ process RECONPLOT { 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}" : "${cn_source}_${sv_source}") + 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" @@ -77,7 +77,7 @@ process RECONPLOT { """ stub: - prefix = task.ext.prefix ?: (cn_source == sv_source ? "${cn_source}" : "${cn_source}_${sv_source}") + 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 diff --git a/modules/local/reconplot/tests/nextflow.config b/modules/local/reconplot/tests/nextflow.config index 25fc5a08..355e1561 100644 --- a/modules/local/reconplot/tests/nextflow.config +++ b/modules/local/reconplot/tests/nextflow.config @@ -1,5 +1,5 @@ process { withName: 'RECONPLOT' { - ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${cn_source}_${sv_source}" } + ext.prefix = { cn_source == sv_source ? "${cn_source}" : "${sv_source}_${cn_source}" } } } diff --git a/subworkflows/local/reconplot_figures.nf b/subworkflows/local/reconplot_figures.nf index cb72c6bc..0801db6e 100644 --- a/subworkflows/local/reconplot_figures.nf +++ b/subworkflows/local/reconplot_figures.nf @@ -1,6 +1,6 @@ // IMPORT MODULES -include { RECONPLOT as RECONPLOT_ASCAT_SEVERUS } from '../../modules/local/reconplot/main' -include { RECONPLOT as RECONPLOT_WAKHAN_SEVERUS } from '../../modules/local/reconplot/main' +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' @@ -45,7 +45,7 @@ workflow RECONPLOT_FIGURES { // severus_sv_files: [meta, [severus_somatic.vcf.gz]] // - // MODULE: RECONPLOT_ASCAT_SEVERUS (label: process_low) + // 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 // @@ -57,11 +57,11 @@ workflow RECONPLOT_FIGURES { .map { meta, cn, sv -> [meta, 'ascat', cn, 'severus', sv] } .set { ascat_input } - RECONPLOT_ASCAT_SEVERUS( ascat_input, reconplot_src, reconplot_pkg, genome ) - ch_versions = ch_versions.mix(RECONPLOT_ASCAT_SEVERUS.out.versions) + RECONPLOT_SEVERUS_ASCAT( ascat_input, reconplot_src, reconplot_pkg, genome ) + ch_versions = ch_versions.mix(RECONPLOT_SEVERUS_ASCAT.out.versions) // - // MODULE: RECONPLOT_WAKHAN_SEVERUS (label: process_low) + // 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/) // @@ -78,8 +78,8 @@ workflow RECONPLOT_FIGURES { .map { meta, cn, sv -> [meta, 'wakhan', cn, 'severus', sv] } .set { wakhan_input } - RECONPLOT_WAKHAN_SEVERUS( wakhan_input, reconplot_src, reconplot_pkg, genome ) - ch_versions = ch_versions.mix(RECONPLOT_WAKHAN_SEVERUS.out.versions) + 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) From 931d2d119c62f5ed04f69c7899af238152322d78 Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Tue, 22 Sep 2026 19:51:34 +0100 Subject: [PATCH 07/10] Add stub nf-tests for the Padfoot and ReConPlot subworkflows Padfoot and ReConPlot need whole-genome caller output, so the pipeline tests skip them. These -stub tests feed PADFOOT_ANNOTATION and RECONPLOT_FIGURES stand-in files named like the real caller output and check the pairing: Severus + the top-ranked Wakhan solution, ASCAT + Severus, SAVANA with and without allele counts, and that a sample missing one side of a pair (no integer-CN VCF, one HP bed, no SAVANA fit) is dropped rather than failing the run. Tagged small so they run on every pull request. Co-Authored-By: Claude Fable 5.1 --- .../local/tests/padfoot_annotation.nf.test | 73 ++++++++++++++ .../local/tests/reconplot_figures.nf.test | 98 +++++++++++++++++++ 2 files changed, 171 insertions(+) create mode 100644 subworkflows/local/tests/padfoot_annotation.nf.test create mode 100644 subworkflows/local/tests/reconplot_figures.nf.test 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 } + ) + } + } +} From edab7a2b927d81c4c7789793c6ae86ba209ab7f1 Mon Sep 17 00:00:00 2001 From: Tim-Yu Date: Tue, 22 Sep 2026 19:51:34 +0100 Subject: [PATCH 08/10] Bring the Padfoot and ReConPlot docs in line with dev's conventions CHANGELOG entries take dev's linked one-to-three-sentence form (PR number to fill in once the pull request exists). README lists padfoot/ and reconplot/ in the output tree and the output paragraph. docs/usage.md says why Padfoot is pinned to the Tim-Yu fork (SAVANA input support not yet upstream). The padfoot_genome schema description moves its detail into help_text, the duplicated caller-pair name in docs/output.md is fixed, and the test-profile skip comments shrink to one line. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 6 ++++-- README.md | 15 ++++++++++++++- conf/test.config | 3 +-- conf/test_full.config | 3 +-- docs/output.md | 2 +- docs/usage.md | 2 +- nextflow_schema.json | 3 ++- 7 files changed, 24 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ff387c33..72049021 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,8 +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). -- Added Padfoot annotation of somatic SVs + CNAs for both Severus/Wakhan and SAVANA (#189) outputs (paired and tumour-only). Padfoot source is fetched from GitHub (or `--padfoot_dir`) and run in a public image bundling its dependencies with RepeatMasker 4.2.4 + Dfam 4.0 (`ghcr.io/tim-yu/padfoot-repeatmasker`, recipe in `containers/padfoot/`) or in a conda env; RepeatMasker annotation of inserted sequences runs by default. New params `skip_padfoot`, `padfoot_url`, `padfoot_dir`, `padfoot_genome`, `padfoot_gff`, `padfoot_rm`, `padfoot_run_repeatmasker` (@Tim-Yu). -- Added ReConPlot rearrangement + copy-number figures for each available CN/SV caller pair (`reconplot/{severus_ascat,severus_wakhan,savana}/`): per-chromosome, genome-wide and optional region-focus panels. The wrapper ships in `assets/reconplot/`; the ReConPlot R package is staged from GitHub or a local checkout; runs in a public image (`ghcr.io/tim-yu/reconplot`, recipe in `containers/reconplot/`) or conda. New params `skip_reconplot`, `reconplot_pkg_url`, `reconplot_pkg_dir`, `reconplot_genome`, `reconplot_max_cn`, `reconplot_min_svlen`, `reconplot_exclude_vntr`, `reconplot_regions`, `reconplot_genes`, `reconplot_baf_track`, `reconplot_format` (@Tim-Yu). +- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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). +- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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). +- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - `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). +- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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/README.md b/README.md index a071d6bb..4a8c2149 100644 --- a/README.md +++ b/README.md @@ -114,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: @@ -122,6 +122,9 @@ Example output directory structure: ├── Sample 1 │ ├── ascat │ ├── bamfiles +│ ├── padfoot +│ │ ├── savana +│ │ └── severus_wakhan │ ├── qc │ │ ├── tumor │ │ │ ├── cramino_aln @@ -129,6 +132,9 @@ Example output directory structure: │ │ │ ├── fibertoolsrs │ │ │ ├── mosdepth │ │ │ ├── samtools +│ ├── reconplot +│ │ ├── savana +│ │ └── severus_wakhan │ ├── variants │ │ ├──clairS-TO │ │ ├──severus @@ -142,6 +148,9 @@ Example output directory structure: ├── Sample 2 │ ├── ascat │ ├── bamfiles +│ ├── padfoot +│ │ ├── savana +│ │ └── severus_wakhan │ ├── qc │ │ ├── tumor │ │ │ ├── cramino_aln @@ -155,6 +164,10 @@ Example output directory structure: │ │ │ ├── fibertoolsrs │ │ │ ├── mosdepth │ │ │ ├── samtools +│ ├── reconplot +│ │ ├── savana +│ │ ├── severus_ascat +│ │ └── severus_wakhan │ ├── variants │ │ ├── clair3 │ │ ├── clairS diff --git a/conf/test.config b/conf/test.config index adb63b3e..1baaf6dd 100644 --- a/conf/test.config +++ b/conf/test.config @@ -70,8 +70,7 @@ params { skip_wakhan = true skip_ascat = true skip_modkit = true - // Padfoot/ReConPlot need whole-genome caller output and fetch multi-GB images/annotations; - // covered by the module-level nf-tests instead. + // 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" diff --git a/conf/test_full.config b/conf/test_full.config index dab4ef7c..96617892 100644 --- a/conf/test_full.config +++ b/conf/test_full.config @@ -36,8 +36,7 @@ params { skip_wakhan = true skip_ascat = true skip_vep_plugins = true - // Padfoot/ReConPlot need whole-genome caller output and fetch multi-GB images/annotations; - // covered by the module-level nf-tests instead. + // Padfoot/ReConPlot need whole-genome caller output; covered by the module- and subworkflow-level nf-tests skip_padfoot = true skip_reconplot = true diff --git a/docs/output.md b/docs/output.md index 0dc949d4..47ef5ae8 100644 --- a/docs/output.md +++ b/docs/output.md @@ -109,7 +109,7 @@ 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`/`severus_wakhan` need `--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. +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` diff --git a/docs/usage.md b/docs/usage.md index 8b6f1b7c..240b3baf 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -387,7 +387,7 @@ Both tools run from `ghcr.io/ljwharbers/sigprofiler`, which adds CHM13 support n - `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) 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`. +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. diff --git a/nextflow_schema.json b/nextflow_schema.json index b6f29b90..9902bd4d 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -396,7 +396,8 @@ "padfoot_genome": { "type": "string", "enum": ["hg38", "chm13", "mm10"], - "description": "Padfoot genome preset. Default: inferred from --genome (GRCh38 -> hg38, CHM13 -> chm13). Padfoot only bundles hg38/mm10 annotations; chm13 requires --padfoot_gff and --padfoot_rm." + "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", From e793b614217d8680e555a1d29fe7c1fde3efe51c Mon Sep 17 00:00:00 2001 From: Tim Yu Date: Wed, 23 Sep 2026 22:24:11 +0100 Subject: [PATCH 09/10] Update usage instructions for VEP plugins Clarify usage of `--vep_revel` and `--vep_eve` options and licensing information for `vep_plugins`. Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- docs/usage.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/usage.md b/docs/usage.md index 240b3baf..d6697d5e 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -553,7 +553,8 @@ Prepared REVEL and EVE files are published to `/vep_plugins/`, alongside a later run can point `--vep_revel` / `--vep_revel_tbi` (or the `--vep_eve` pair) at them and skip both the download and the reshaping. The two AlphaMissense tables are fetched already indexed. -> [!IMPORTANT] > `/vep_plugins/` holds data that is free for non-commercial use only. Exclude it when you +> [!IMPORTANT] +> `/vep_plugins/` holds data that is free for non-commercial use only. Exclude it when you > share or archive a results directory — passing it on is redistribution, which those licences do > not grant you. From 8049f36cbc6c774bbca93155cb49b9008ad1949f Mon Sep 17 00:00:00 2001 From: Boyu Yu Date: Wed, 23 Sep 2026 15:00:18 +0100 Subject: [PATCH 10/10] Point the CHANGELOG entries at PR #202 --- CHANGELOG.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 72049021..c51c08f0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,10 +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). -- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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). -- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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). -- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - `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). -- [#XXX](https://github.com/IntGenomicsLab/lrsomatic/pull/XXX) - 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). +- [#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`