feat: add low-confidence site filtering

Ensures missing variants are distinguished from zero frequency values by applying depth and quality filters. Prevents low-confidence sites from biasing allele frequency analyses.

Author: ahmig

template.qmd

@@ -92,8 +92,7 @@ n.samples <- table %>% pull(Sample) %>% unique() %>% length()
 heatmap.df <- read.csv(params$heat_tab)
 row.names(heatmap.df) <- heatmap.df$`...1`
 heatmap.df[1] <- NULL
-cor.mat <- cor(heatmap.df, method = params$cor_method)
-cor.mat[is.na(cor.mat)] <- 0
+cor.mat <- cor(heatmap.df, method = params$cor_method, use = "pairwise.complete.obs")
 if (params$cor_method == "pearson") {
   cor.method.name <- "Pearson's"
 } else if (params$cor_method == "spearman") {
@@ -247,10 +246,11 @@ Labeled dots represent nucleotide variants whose allele frequency is correlated
 
 Significantly correlated nucleotide variants are described in more detail in @fig-panel.
 
-![Time series of relative allele frequencies. The shown positions include
-nucleotide variants with a significant correlation with time and sites with more
-than two possible states. Each subplot depicts the progression of the allele
-frequencies in time for a given genome position.](`r params$panel`){#fig-panel}
+![Time series of relative allele frequencies. Nucleotide variants with a significant 
+correlation with time and sites with more than two possible states are included. Each 
+subplot depicts the progression of the allele frequencies in time for a given genome 
+position; lines connecting points indicate data from contiguous time points, while 
+unconnected points denote intervals with missing data.](`r params$panel`){#fig-panel}
 
 A `r dist.tree.algo` tree has been constructed using pairwise distances
 between target samples (@fig-tree), based on the allele frequencies measured from
@@ -277,8 +277,18 @@ The heatmap is interactive and allows zooming in on specific regions.
 ```{r fig-heatmap, echo = F, message = F, warning = F, fig.align = 'center'}
 #| fig-cap: "Interactive heatmap with hierarchical clustering of the pairwise correlation coefficients between the time series of allele frequencies in the case study."
 
+# Calculate dendrograms using data with no NAs
+cor.mat.clustering <- cor.mat
+cor.mat.clustering[is.na(cor.mat.clustering)] <- 0 
+dist.rows <- dist(cor.mat.clustering)
+dend.rows <- hclust(dist.rows)
+dend.cols <- hclust(dist(t(cor.mat.clustering)))
+
+# Display custom dendrograms over real data
 heatmaply_cor(
   cor.mat,
+  Rowv = as.dendrogram(dend.rows),
+  Colv = as.dendrogram(dend.cols),
   grid_gap = 0.25,
   fontsize_row = 7, 
   fontsize_col = 7,

workflow/rules/report.smk

@@ -345,7 +345,7 @@ rule af_time_correlation_data:
         cor_method = config["COR"]["METHOD"],
         cor_exact = config["COR"]["EXACT"],
     input:
-        variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
+        variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
         metadata = config["METADATA"],
     output:
         fmt_variants = temp(REPORT_DIR_TABLES/"variants.filled.dated.tsv"),
@@ -392,6 +392,19 @@ rule af_trajectory_panel_plot:
         "../scripts/report/af_trajectory_panel_plot.R"
 
 
+rule pairwise_trajectory_correlation_data:
+    conda: "../envs/renv.yaml"
+    input:
+        variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
+        metadata = config["METADATA"],
+    output:
+        table = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation.csv"),
+    log:
+        LOGDIR / "pairwise_trajectory_correlation_data" / "log.txt"
+    script:
+        "../scripts/pairwise_trajectory_correlation_data.R"
+
+
 rule summary_table:
     conda: "../envs/renv.yaml"
     input:
@@ -423,7 +436,7 @@ rule report:
         temest     = report(REPORT_DIR_PLOTS/"time_signal.png"),
         evo        = report(REPORT_DIR_PLOTS/"dn_and_ds.png"),
         omega_plot = report(REPORT_DIR_PLOTS/"dnds.png"),
-        heat_table = report(OUTDIR/"vaf"/"pairwise_trajectory_correlation.csv"),
+        heat_table = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation.csv"),
         freyja_ts  = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
         value      = REPORT_DIR_TABLES/"diversity.json",
         stats_lm   = REPORT_DIR_TABLES/"time_signal.json",

workflow/rules/vaf.smk

@@ -233,8 +233,8 @@ rule bcftools_mpileup_all_sites:
         bam = get_input_bam,
         reference = OUTDIR/"vaf"/"{sample}.reference.fasta",
     output:
-        mpileup = OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.vcf",
-        query = OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.query.tsv",
+        mpileup = temp(OUTDIR / "all_sites" / "{sample}.mpileup.vcf"),
+        query = temp(OUTDIR / "all_sites" / "{sample}.query.tsv"),
     log:
         mpileup = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.txt",
         query = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.query.txt",
@@ -251,33 +251,41 @@ rule filter_mpileup_all_sites:
         min_total_ADF = 0,
         min_total_ADR = 0,
     input:
-        OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.query.tsv",
+        OUTDIR / "all_sites" / "{sample}.query.tsv",
     output:
-        OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.filtered.tsv",
+        temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
     log:
         LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt"
     run:
         import pandas as pd
         df = pd.read_csv(input[0], sep="\t")
-        df["ref_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0]))
-        df["total_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
-        df["total_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
-        df["total_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
+        df["SAMPLE"] = wildcards.sample
+        df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0]))
+        df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
+        df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
+        df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
         df[
-            (df.total_AD >= params.min_total_AD) &
-            (df.total_ADF >= params.min_total_ADF) &
-            (df.total_ADR >= params.min_total_ADR)
+            (df.TOTAL_AD >= params.min_total_AD) &
+            (df.TOTAL_ADF >= params.min_total_ADF) &
+            (df.TOTAL_ADR >= params.min_total_ADR)
         ].to_csv(output[0], sep="\t", index=False)
 
 
-rule pairwise_trajectory_correlation:
+use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
+    input:
+        expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()),
+    output:
+        OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
+
+
+rule fill_all_sites:
     conda: "../envs/renv.yaml"
     input:
-        variants =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
-        metadata = config["METADATA"],
+        variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
+        sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
     output:
-        table = report(OUTDIR/"vaf"/"pairwise_trajectory_correlation.csv"),
+        variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
     log:
-        LOGDIR / "pairwise_trajectory_correlation" / "log.txt"
+        LOGDIR / "fill_all_sites" / "log.txt"
     script:
-        "../scripts/pairwise_trajectory_correlation.R"
+        "../scripts/fill_all_sites.R"

workflow/scripts/fill_all_sites.R

@@ -0,0 +1,54 @@
+#!/usr/bin/env Rscript
+
+# Write stdout and stderr to log file
+log <- file(snakemake@log[[1]], open = "wt")
+sink(log, type = "message")
+sink(log, type = "output")
+
+library(dplyr)
+library(readr)
+library(tidyr)
+library(logger)
+
+log_threshold(INFO)
+
+log_info("Reading variants")
+variants <- read_tsv(snakemake@input[["variants"]])
+
+# Create a mapping of variant names to their genomic position
+variant_coords <- variants %>%
+  select(VARIANT_NAME, REGION, POS) %>%
+  distinct()
+
+log_info("Reading filtered sites")
+sites <- read_tsv(snakemake@input[["sites"]]) %>%
+  select(SAMPLE, POS) %>%  # TODO: consider region/chrom
+  distinct() %>%
+  mutate(FILTER_PASS = TRUE)
+
+log_info("Processing variants")
+all_variants <- variants %>%
+  # Select minimal columns
+  select(VARIANT_NAME, REGION, SAMPLE, ALT_FREQ) %>%
+  # Handle duplicates
+  group_by(SAMPLE, VARIANT_NAME, REGION) %>%
+  summarise(ALT_FREQ = sum(ALT_FREQ, na.rm = TRUE), .groups = "drop") %>%
+  # Complete with NA
+  complete(SAMPLE, VARIANT_NAME, REGION) %>%
+  # Assign genomic positions for all combinations
+  left_join(variant_coords, by = c("REGION", "VARIANT_NAME")) %>%
+  # Merge filtered sites
+  # TODO: consider region/chrom
+  left_join(sites, by = c("SAMPLE", "POS")) %>%
+  replace_na(list(FILTER_PASS = FALSE)) %>%
+  # Fill missing frequencies conditionally
+  mutate(
+    ALT_FREQ = case_when(
+      !is.na(ALT_FREQ) ~ ALT_FREQ,  # variant was called, keep the frequency
+      FILTER_PASS ~ 0,              # missing but site is covered: 0 (reference)
+      TRUE ~ NA                     # missing and not covered: NA (unknown)
+    )
+  )
+
+log_info("Saving table")
+write_tsv(all_variants, snakemake@output[["variants"]])

…cripts/pairwise_trajectory_correlation.R …s/pairwise_trajectory_correlation_data.Rworkflow/scripts/pairwise_trajectory_correlation.R renamed to workflow/scripts/pairwise_trajectory_correlation_data.R workflow/scripts/pairwise_trajectory_correlation.R renamed to workflow/scripts/pairwise_trajectory_correlation_data.R

@@ -14,6 +14,7 @@ library(logger)
 
 log_threshold(INFO)
 
+log_info("Reading variants")
 variants <- read_tsv(snakemake@input[["variants"]])
 
 # Obtain sample names ordered by CollectionDate
@@ -22,20 +23,17 @@ date_order <- read_csv(snakemake@input[["metadata"]]) %>%
   pull(ID) %>%
   unique()
 
-# Create SNP variable and select useful variables
-variants <- variants %>% select(VARIANT_NAME, SAMPLE, ALT_FREQ)
-
-variants <- variants %>%
+all_variants_wider <- variants %>%
+  select(SAMPLE, VARIANT_NAME, ALT_FREQ) %>%
   pivot_wider(
     names_from = VARIANT_NAME,
-    values_from = ALT_FREQ,
-    values_fill = 0,
-    values_fn = sum
+    values_from = ALT_FREQ
   ) %>%
+  # Apply chronological ordering
   arrange(factor(SAMPLE, levels = date_order)) %>%
-  # Removes "|"-separated annotations, keeping the first one + ellipsis (clarifies heatmap)
+  # Removes "|"-separated annotations, keeping the first one + ellipsis
   rename_with(~ str_replace(., "^([^|]+)\\|.*$", "\\1(...)"), -SAMPLE) %>%
   column_to_rownames(var = "SAMPLE")
 
-log_info("Saving table to create heatmap")
-write.csv(variants, snakemake@output[["table"]])
+log_info("Saving table")
+write.csv(all_variants_wider, snakemake@output[["table"]])

workflow/scripts/report/af_time_correlation_data.R

@@ -24,11 +24,11 @@ variants <- read_delim(
     "POS"
   )
 ) %>%
-  # Fill positions without alt frequency with 0
+  # Fill positions without alt frequency with NA
   complete(
     nesting(REGION, VARIANT_NAME, POS),
     SAMPLE,
-    fill = list(ALT_FREQ = 0)
+    fill = list(ALT_FREQ = NA)
   )
 
 log_info("Reading metadata")