temp: consider filtered sites before imputing zero alt frequency

Author: ahmig

workflow/rules/report.smk

@@ -347,6 +347,7 @@ rule af_time_correlation_data:
     input:
         variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"],
+        sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
     output:
         fmt_variants = temp(REPORT_DIR_TABLES/"variants.filled.dated.tsv"),
         correlations = report(REPORT_DIR_TABLES/"af_time_correlation.csv"),

workflow/rules/vaf.smk

@@ -252,23 +252,31 @@ rule filter_mpileup_all_sites:
     input:
         OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.query.tsv",
     output:
-        OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.filtered.tsv",
+        temp(OUTDIR / "bcftools_mpileup_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)
 
 
+use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
+    input:
+        expand(OUTDIR / "bcftools_mpileup_all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()),
+    output:
+        OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
+
+
 rule pairwise_trajectory_correlation:
     conda: "../envs/renv.yaml"
     input:

workflow/scripts/report/af_time_correlation_data.R

@@ -13,6 +13,13 @@ library(logger)
 log_threshold(INFO)
 
 # Read data
+log_info("Reading filtered sites")
+sites <- read_tsv(
+  snakemake@input[["sites"]],
+  col_select = c("SAMPLE", "POS")
+) %>%
+  mutate(SELECTED_SITE = TRUE)
+
 log_info("Reading variants table")
 variants <- read_delim(
   snakemake@input[["variants"]],
@@ -24,11 +31,20 @@ variants <- read_delim(
     "POS"
   )
 ) %>%
+  # Annotate sites selected after filter
+  left_join(sites, by = c("SAMPLE", "POS")) %>%
+  replace_na(list(SELECTED_SITE = FALSE)) %>%
   # Fill positions without alt frequency with 0
+  # TODO: add a flag to choose behavior
   complete(
     nesting(REGION, VARIANT_NAME, POS),
     SAMPLE,
-    fill = list(ALT_FREQ = 0)
+    fill = list(ALT_FREQ = NA)
+  ) %>%
+  mutate(
+    ALT_FREQ = case_when(
+      is.na(ALT_FREQ) & SELECTED_POS ~ 0
+    )
   )
 
 log_info("Reading metadata")