From 1f2ae56332d1362e1ed19f4232317513ed2b5aad Mon Sep 17 00:00:00 2001 From: m-jahn Date: Mon, 5 Oct 2026 14:10:41 +0200 Subject: [PATCH 1/5] fix: #48 --- workflow/notebooks/report.Rmd | 75 +++++++++++++++++++++++++---------- 1 file changed, 55 insertions(+), 20 deletions(-) diff --git a/workflow/notebooks/report.Rmd b/workflow/notebooks/report.Rmd index 783feab..d169086 100644 --- a/workflow/notebooks/report.Rmd +++ b/workflow/notebooks/report.Rmd @@ -36,6 +36,9 @@ Loaded required R packages: - `GenomicFeatures` - `Biostrings` +```{r, echo = FALSE} +knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE) +``` ```{r} suppressPackageStartupMessages({ @@ -49,10 +52,6 @@ suppressPackageStartupMessages({ }) ``` -```{r} -knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE) -``` - ```{r} # messages messages <- c() @@ -81,8 +80,9 @@ custom_range <- function(n = 5, quantitative = FALSE) { } custom_theme <- function( - base_size = 12, base_line_size = 1.0, - base_rect_size = 1.0, ...) { + base_size = 12, base_line_size = 1.0, + base_rect_size = 1.0, ... +) { theme_light( base_size = base_size, base_line_size = base_line_size, @@ -129,8 +129,9 @@ figheight <- 4.5 # function to export an image as svg and png save_plot <- function( - pl, path = figure_dir, - width = figwidth, height = figheight, index = NA) { + pl, path = figure_dir, + width = figwidth, height = figheight, index = NA +) { pl_name <- deparse(substitute(pl)) if (!is.na(index)) pl_name <- paste0(pl_name, "_", index) svg( @@ -163,10 +164,6 @@ list_variants <- unlist(snakemake@input$variants) print(list_variants) ``` -The imported variant data has the following layout (first 10 rows shown): - -- *NOTE: tables are only rendered in HTML output but not PDF.* - ```{r} df <- lapply(list_variants, function(file) { properties <- strsplit(file, "/")[[1]] @@ -197,11 +194,25 @@ if ("CSQ" %in% colnames(df)) { dplyr::select(-matches("CSQ_")) %>% mutate(VARTYPE = gsub("coding_sequence_variant&", "", VARTYPE)) } - -head(df, n = 10) ``` +**Consensus variants** + +- All variants which are detected in a number of samples exceeding the predefined threshold are listed in the following consensus table. +- The table can also be found in `report/_report.tsv`. +- *NOTE: Tables are only rendered in HTML output but not PDF.* + + ```{r} +# import sample sheet +df_samples <- read_tsv(snakemake@config$samplesheet, show_col_types = FALSE) + +# annotation columns to select +cols_select <- c( + "CHROM", "POS", "REF", "ALT", "VARTYPE", "IMPACT", + "GENE", "AA_CHANGE", "CODON_CHANGE", "INDEL" +) + # import consensus file list_consensus <- snakemake@input$consensus df_consensus <- read_tsv( @@ -209,7 +220,32 @@ df_consensus <- read_tsv( show_col_types = FALSE, col_names = c("CHROM", "POS", "REF", "ALT", "APPEARANCE") ) %>% - mutate(POS = as.numeric(POS)) + mutate( + POS = as.numeric(POS), + APPEARANCE = as.character(APPEARANCE) + ) %>% + separate_wider_position(APPEARANCE, width = setNames(rep_along(df_samples$sample, 1), df_samples$sample)) %>% + rowwise() %>% + mutate( + across(df_samples$sample, as.numeric), + N_SAMPLES_ABS = sum(c_across(df_samples$sample)), + N_SAMPLES_REL = N_SAMPLES_ABS / length(df_samples$sample) + ) + +# add variant information and export annotated table +df_consensus <- df_consensus %>% + left_join( + df %>% dplyr::select(any_of(cols_select)) %>% + distinct(), + by = join_by(CHROM, POS, REF, ALT) + ) %>% + dplyr::select(any_of(c(cols_select, "N_SAMPLES_ABS", "N_SAMPLES_REL")), everything()) %>% + arrange(CHROM, POS) + +output_table <- str_replace(snakemake@output$html, "\\.html$", ".tsv") +write_tsv(df_consensus, file = output_table) + +df_consensus ``` ```{r} @@ -332,7 +368,10 @@ for (chrom in chromosomes) { ```{r} df_filtered <- df %>% - inner_join(df_consensus, by = join_by(CHROM, POS, REF, ALT)) + inner_join( + dplyr::select(df_consensus, CHROM, POS, REF, ALT), + by = join_by(CHROM, POS, REF, ALT) + ) no_var_message <- paste0( "No variants detected that appear more than ", vcf_count, " times" @@ -558,10 +597,6 @@ if (nrow(df_filtered)) { ```{r} filter_variants <- function(df, vartype, topn) { - cols_select <- c( - "CHROM", "POS", "REF", "ALT", "VARTYPE", "IMPACT", - "GENE", "AA_CHANGE", "CODON_CHANGE", "INDEL" - ) if ("VARTYPE" %in% colnames(df)) { return( df %>% From ba84424ffe8d3064a7992406f9a5904dd73c0979 Mon Sep 17 00:00:00 2001 From: jahn Date: Tue, 6 Oct 2026 07:58:01 +0200 Subject: [PATCH 2/5] fix: #51 --- workflow/notebooks/report.Rmd | 21 ++++++++++++++++----- workflow/rules/report.smk | 11 +++++++++-- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/workflow/notebooks/report.Rmd b/workflow/notebooks/report.Rmd index d169086..e57d10c 100644 --- a/workflow/notebooks/report.Rmd +++ b/workflow/notebooks/report.Rmd @@ -213,8 +213,19 @@ cols_select <- c( "GENE", "AA_CHANGE", "CODON_CHANGE", "INDEL" ) -# import consensus file +# import consensus files list_consensus <- snakemake@input$consensus +index_consensus <- snakemake@input$consensus_index + +names_consensus <- read_tsv( + index_consensus, + show_col_types = FALSE, + skip = 4, + col_names = c("col1", "col2", "col3") +) %>% + mutate(name = str_split_i(col3, "\\/", i = 4) %>% str_remove("\\_(vep|snpeff)\\.vcf\\.gz$")) %>% + pull(name) + df_consensus <- read_tsv( list_consensus, show_col_types = FALSE, @@ -224,12 +235,12 @@ df_consensus <- read_tsv( POS = as.numeric(POS), APPEARANCE = as.character(APPEARANCE) ) %>% - separate_wider_position(APPEARANCE, width = setNames(rep_along(df_samples$sample, 1), df_samples$sample)) %>% + separate_wider_position(APPEARANCE, width = setNames(rep_along(names_consensus, 1), names_consensus)) %>% rowwise() %>% mutate( - across(df_samples$sample, as.numeric), - N_SAMPLES_ABS = sum(c_across(df_samples$sample)), - N_SAMPLES_REL = N_SAMPLES_ABS / length(df_samples$sample) + across(names_consensus, as.numeric), + N_SAMPLES_ABS = sum(c_across(names_consensus)), + N_SAMPLES_REL = N_SAMPLES_ABS / length(names_consensus) ) # add variant information and export annotated table diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index e154a1e..293ecde 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -18,7 +18,8 @@ rule bcftools_intersection: input: get_variants, output: - "results/{caller}/consensus/variants.vcf", + consensus="results/{caller}/consensus/variants.vcf", + consensus_index="results/{caller}/consensus/README.txt", log: "results/{caller}/consensus/variants.log", conda: @@ -29,7 +30,12 @@ rule bcftools_intersection: message: "intersect variant calls from different callers" shell: - "bcftools isec -n+{params.vcf_count} {params.extra} -o {output} {input}" + """ + output_dir="$(dirname {output.consensus})" + bcftools isec -p ${{output_dir}} -n+{params.vcf_count} {params.extra} {input} >{log} 2>&1 + mv ${{output_dir}}/sites.txt {output.consensus} + rm -rf ${{output_dir}}/0*.vcf + """ rule report_html: @@ -38,6 +44,7 @@ rule report_html: gff=rules.get_genome.output.gff, variants=get_variants, consensus="results/{caller}/consensus/variants.vcf", + consensus_index="results/{caller}/consensus/README.txt", output: html="results/report/{caller}_report.html", log: From da30a256b18fe8a6e3cf37e5a21c420e1db4f585 Mon Sep 17 00:00:00 2001 From: jahn Date: Tue, 6 Oct 2026 11:35:12 +0200 Subject: [PATCH 3/5] fix: review comments --- workflow/notebooks/report.Rmd | 8 +++----- workflow/rules/report.smk | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/workflow/notebooks/report.Rmd b/workflow/notebooks/report.Rmd index e57d10c..1e0f5d1 100644 --- a/workflow/notebooks/report.Rmd +++ b/workflow/notebooks/report.Rmd @@ -199,14 +199,11 @@ if ("CSQ" %in% colnames(df)) { **Consensus variants** - All variants which are detected in a number of samples exceeding the predefined threshold are listed in the following consensus table. -- The table can also be found in `report/_report.tsv`. +- The table can also be found in `results/report/_report.tsv`. - *NOTE: Tables are only rendered in HTML output but not PDF.* ```{r} -# import sample sheet -df_samples <- read_tsv(snakemake@config$samplesheet, show_col_types = FALSE) - # annotation columns to select cols_select <- c( "CHROM", "POS", "REF", "ALT", "VARTYPE", "IMPACT", @@ -241,7 +238,8 @@ df_consensus <- read_tsv( across(names_consensus, as.numeric), N_SAMPLES_ABS = sum(c_across(names_consensus)), N_SAMPLES_REL = N_SAMPLES_ABS / length(names_consensus) - ) + ) %>% + ungroup() # add variant information and export annotated table df_consensus <- df_consensus %>% diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index 293ecde..0e3bca1 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -18,7 +18,7 @@ rule bcftools_intersection: input: get_variants, output: - consensus="results/{caller}/consensus/variants.vcf", + consensus="results/{caller}/consensus/variants.tsv", consensus_index="results/{caller}/consensus/README.txt", log: "results/{caller}/consensus/variants.log", @@ -43,7 +43,7 @@ rule report_html: fasta=rules.get_genome.output.fasta, gff=rules.get_genome.output.gff, variants=get_variants, - consensus="results/{caller}/consensus/variants.vcf", + consensus="results/{caller}/consensus/variants.tsv", consensus_index="results/{caller}/consensus/README.txt", output: html="results/report/{caller}_report.html", From 0d770490312b1a45ee89e0c77ce07a8f022c343a Mon Sep 17 00:00:00 2001 From: jahn Date: Tue, 6 Oct 2026 13:13:38 +0200 Subject: [PATCH 4/5] fix: allow to switch off variant calling --- config/schemas/config.schema.yml | 5 ++-- workflow/Snakefile | 6 +---- workflow/rules/common.smk | 45 ++++++++++++++++++++------------ 3 files changed, 33 insertions(+), 23 deletions(-) diff --git a/config/schemas/config.schema.yml b/config/schemas/config.schema.yml index aa96695..1ec047e 100644 --- a/config/schemas/config.schema.yml +++ b/config/schemas/config.schema.yml @@ -180,8 +180,9 @@ properties: type: object properties: tool: - type: [string, array] - description: the variant caller to use, one of 'freebayes', 'bcftools', or both + type: [string, array, "null"] + description: the variant caller to use, one of 'freebayes', 'bcftools', both, or none (empty list, 'null') + default: freebayes bcftools_pileup: type: object properties: diff --git a/workflow/Snakefile b/workflow/Snakefile index c70666e..8ab3932 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -43,8 +43,4 @@ onerror: rule all: default_target: True input: - "results/multiqc/multiqc_report.html", - expand( - "results/report/{caller}_report.pdf", - caller=config["variant_calling"]["tool"], - ), + get_targets(), diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index e80dc8d..0c513a7 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -105,22 +105,35 @@ def get_multiqc_input(wildcards): sample=samples.index, ext=["png", "raw", "metrics"], ) - result += expand( - "results/{caller}/call/{sample}{ext}", - sample=samples.index, - caller=config["variant_calling"]["tool"], - ext=["_stats.txt", "_all.vcf", "_variants.vcf"], - ) - result += expand( - "results/{caller}/effect/{sample}{ext}", - sample=samples.index, - caller=config["variant_calling"]["tool"], - ext=( - ["_vep.vcf", "_vep.html"] - if config["variant_annotation"]["tool"] == "vep" - else ["_snpeff.vcf", "_snpeff.csv"] - ), - ) + if config["variant_calling"]["tool"] not in [None, [], "null"]: + result += expand( + "results/{caller}/call/{sample}{ext}", + sample=samples.index, + caller=config["variant_calling"]["tool"], + ext=["_stats.txt", "_all.vcf", "_variants.vcf"], + ) + result += expand( + "results/{caller}/effect/{sample}{ext}", + sample=samples.index, + caller=config["variant_calling"]["tool"], + ext=( + ["_vep.vcf", "_vep.html"] + if config["variant_annotation"]["tool"] == "vep" + else ["_snpeff.vcf", "_snpeff.csv"] + ), + ) + return result + + +# get targets for entire workflow +def get_targets(): + result = [] + result += ["results/multiqc/multiqc_report.html"] + if config["variant_calling"]["tool"] not in [None, [], "null"]: + result += expand( + "results/report/{caller}_report.pdf", + caller=config["variant_calling"]["tool"], + ) return result From 34034eb51118ca04476c7ee3f30b4a35f5cfcdc1 Mon Sep 17 00:00:00 2001 From: jahn Date: Tue, 6 Oct 2026 14:25:55 +0200 Subject: [PATCH 5/5] fix: report snpeff variant effects --- workflow/notebooks/report.Rmd | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/workflow/notebooks/report.Rmd b/workflow/notebooks/report.Rmd index 1e0f5d1..a77b2ee 100644 --- a/workflow/notebooks/report.Rmd +++ b/workflow/notebooks/report.Rmd @@ -184,15 +184,29 @@ df <- lapply(list_variants, function(file) { # parse the variant effect column if ("CSQ" %in% colnames(df)) { + # FROM VEP df <- df %>% mutate(CSQ = str_split(CSQ, "\\|")) %>% unnest_wider(CSQ, names_sep = "_") %>% dplyr::rename( VARTYPE = `CSQ_2`, IMPACT = `CSQ_3`, GENE = `CSQ_4`, - AA_CHANGE = `CSQ_11`, CODON_CHANGE = `CSQ_16` + DNA_CHANGE = `CSQ_11`, CODON_CHANGE = `CSQ_16` ) %>% dplyr::select(-matches("CSQ_")) %>% mutate(VARTYPE = gsub("coding_sequence_variant&", "", VARTYPE)) +} else if ("ANN" %in% colnames(df)) { + # FROM SNPEFF + df <- df %>% + mutate(ANN = str_split(ANN, "\\|")) %>% + unnest_wider(ANN, names_sep = "_") %>% + dplyr::rename( + VARTYPE = `ANN_2`, IMPACT = `ANN_3`, GENE = `ANN_4`, + CODON_CHANGE = `ANN_11` + ) %>% + unite("DNA_CHANGE", matches("ANN_(7|8|9|10)"), sep = "|", remove = TRUE) %>% + dplyr::select(-matches("ANN_")) +} else { + stop("No recognized variant effect columns found ('CSQ' from VEP or 'ANN' from SNPEFF).") } ``` @@ -207,7 +221,7 @@ if ("CSQ" %in% colnames(df)) { # annotation columns to select cols_select <- c( "CHROM", "POS", "REF", "ALT", "VARTYPE", "IMPACT", - "GENE", "AA_CHANGE", "CODON_CHANGE", "INDEL" + "GENE", "DNA_CHANGE", "CODON_CHANGE", "INDEL" ) # import consensus files @@ -235,7 +249,7 @@ df_consensus <- read_tsv( separate_wider_position(APPEARANCE, width = setNames(rep_along(names_consensus, 1), names_consensus)) %>% rowwise() %>% mutate( - across(names_consensus, as.numeric), + across(all_of(names_consensus), as.numeric), N_SAMPLES_ABS = sum(c_across(names_consensus)), N_SAMPLES_REL = N_SAMPLES_ABS / length(names_consensus) ) %>%