Process sample variants in parallel

Also renames the REGION column to SAMPLE, keeping the intended use of REGION in a VCF, and GFF_FEATURE to GB_FEATURE, keeping the intended use of GFF_FEATURE as an output of iVar

Author: ahmig

workflow/rules/vaf.smk

@@ -48,8 +48,6 @@ rule snps_to_ancestor:
                 -m {params.ivar_depth} \
                 -g {input.gff} \
                 -r renamed_reference.fasta
-
-        sed 's/'$ref'/'{wildcards.sample}'/g' {wildcards.sample}.tsv | cat > {output.tsv}
         """
 
 
@@ -69,39 +67,18 @@ rule annotation:
         "../scripts/report/get_annotation.py"
 
 
-rule format_tsv:
-    threads:1
-    shadow: "shallow"
-    input:
-        expand(OUTDIR/"{sample}.tsv", sample = iter_samples())
-    output:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.tsv"
-    log:
-        LOGDIR / "format_tsv" / "log.txt"
-    shell:
-        """
-        path=`echo {input} | awk '{{print $1}}'`
-        grep "^REGION" "$path" > header
-        for tsv in {input}; do
-            tail -n +2 "$tsv"  >> body
-        done
-        cat header body > "{output.tsv}"
-        rm header
-        """
-
-
 rule mask_tsv:
     threads: 1
     conda: "../envs/biopython.yaml"
     params:
          mask_class = ["mask"]
     input: 
-        tsv = OUTDIR/f"{OUTPUT_NAME}.tsv",
+        tsv = OUTDIR/"{sample}.tsv",
         vcf = lambda wildcards: select_problematic_vcf()
     output:
-        masked_tsv = temp(OUTDIR/f"{OUTPUT_NAME}.masked.tsv")
+        masked_tsv = temp(OUTDIR/"{sample}.masked.tsv")
     log:
-        LOGDIR / "mask_tsv" / "log.txt"
+        LOGDIR / "mask_tsv" / "{sample}.log.txt"
     script:
         "../scripts/mask_tsv.py"
 
@@ -114,12 +91,12 @@ rule filter_tsv:
         min_alt_rv = 2,
         min_alt_dp = 2,
     input: 
-        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.tsv",
+        tsv = OUTDIR/"{sample}.masked.tsv",
         annotation = OUTDIR/"annotation.csv"
     output:
-        filtered_tsv = temp(OUTDIR/f"{OUTPUT_NAME}.masked.prefiltered.tsv")
+        filtered_tsv = temp(OUTDIR/"{sample}.masked.prefiltered.tsv")
     log:
-        LOGDIR / "filter_tsv" / "log.txt"
+        LOGDIR / "filter_tsv" / "{sample}.log.txt"
     script:
         "../scripts/filter_tsv.R"
 
@@ -130,11 +107,11 @@ rule tsv_to_vcf:
     params:
         ref_name = config["ALIGNMENT_REFERENCE"],
     input: 
-        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.prefiltered.tsv",
+        tsv = OUTDIR/"{sample}.masked.prefiltered.tsv",
     output:
-        vcf = temp(OUTDIR/f"{OUTPUT_NAME}.vcf")
+        vcf = temp(OUTDIR/"{sample}.vcf")
     log:
-        LOGDIR / "tsv_to_vcf" / "log.txt"
+        LOGDIR / "tsv_to_vcf" / "{sample}.log.txt"
     script:
         "../scripts/tsv_to_vcf.py"
 
@@ -147,11 +124,11 @@ rule variants_effect:
         ref_name = config["ALIGNMENT_REFERENCE"],
         snpeff_data_dir = (BASE_PATH / "config" / "snpeff").resolve()
     input:
-        vcf = OUTDIR/f"{OUTPUT_NAME}.vcf"
+        vcf = OUTDIR/"{sample}.vcf"
     output:
-        ann_vcf = temp(OUTDIR/f"{OUTPUT_NAME}.annotated.vcf")
+        ann_vcf = OUTDIR/"{sample}.annotated.vcf"
     log:
-        LOGDIR / "variants_effect" / "log.txt"
+        LOGDIR / "variants_effect" / "{sample}.log.txt"
     retries: 2
     shell:
         """
@@ -165,7 +142,7 @@ rule variants_effect:
             echo "Local database not found at '{params.snpeff_data_dir}', downloading from repository"
         fi
 
-        snpEff eff -dataDir {params.snpeff_data_dir} -noStats {params.ref_name} {input.vcf} > {output.ann_vcf}
+        snpEff eff -dataDir {params.snpeff_data_dir} -noStats {params.ref_name} {input.vcf} >{output.ann_vcf}
         """
 
 
@@ -180,11 +157,11 @@ rule extract_vcf_fields:
         ],
         sep = ","
     input:
-        vcf = OUTDIR/f"{OUTPUT_NAME}.annotated.vcf"
+        vcf = OUTDIR/"{sample}.annotated.vcf"
     output:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.vcf_fields.tsv"
+        tsv = OUTDIR/"{sample}.vcf_fields.tsv"
     log:
-        LOGDIR / "tsv_to_vcf" / "log.txt"
+        LOGDIR / "tsv_to_vcf" / "{sample}.log.txt"
     shell:
         'SnpSift extractFields -s {params.sep:q} {input.vcf:q} {params.extract_columns} >{output.tsv:q} 2>{log:q}'
 
@@ -194,11 +171,11 @@ rule format_vcf_fields_longer:
     params:
         sep = ","
     input:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.vcf_fields.tsv"
+        tsv = OUTDIR/"{sample}.vcf_fields.tsv"
     output:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.vcf_fields.longer.tsv"
+        tsv = OUTDIR/"{sample}.vcf_fields.longer.tsv"
     log:
-        LOGDIR / "format_vcf_fields_longer" / "log.txt"
+        LOGDIR / "format_vcf_fields_longer" / "{sample}.log.txt"
     script:
         "../scripts/format_vcf_fields_longer.R"
 
@@ -207,11 +184,22 @@ rule vcf_to_tsv:
     threads: 1
     conda: "../envs/renv.yaml"
     input:
-        ann_vcf = OUTDIR/f"{OUTPUT_NAME}.annotated.vcf",
-        pre_tsv = OUTDIR/f"{OUTPUT_NAME}.masked.prefiltered.tsv"
+        ann_vcf = OUTDIR/"{sample}.annotated.vcf",
+        pre_tsv = OUTDIR/"{sample}.masked.prefiltered.tsv"
     output:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv"
+        tsv = OUTDIR/"{sample}.masked.filtered.tsv"
     log:
-        LOGDIR / "vcf_to_tsv" / "log.txt"
+        LOGDIR / "vcf_to_tsv" / "{sample}.log.txt"
     script:
         "../scripts/vcf_to_tsv.R"
+
+
+rule compile_variants:
+    threads: 1
+    input: expand(OUTDIR/"{sample}.masked.filtered.tsv", sample=iter_samples())
+    output:
+        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv"
+    log:
+        LOGDIR / "compile_variants" / "log.txt"
+    script:
+        "../scripts/compile_variants.R"

workflow/scripts/compile_variants.R

@@ -0,0 +1,23 @@
+#!/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(tidyverse)
+library(logger)
+log_threshold(INFO)
+
+log_info("Reading input partial variant tables and adding sample column")
+tables <- lapply(
+    snakemake@input,
+    function(path) {
+        read_tsv(path) %>%
+            mutate(SAMPLE = sub('\\.tsv$', '', path))
+    }
+)
+
+log_info("Binding and writing table")
+bind_rows(tables) %>%
+    write_tsv(snakemake@output$tsv)

workflow/scripts/filter_tsv.R

@@ -47,9 +47,8 @@ data <- distinct(data, pick(!GFF_FEATURE), .keep_all = TRUE)
 features <- read_csv(snakemake@input[["annotation"]])
 
 data <- data %>%
-  select(!GFF_FEATURE) %>%
   left_join(features) %>%
-  rename(GFF_FEATURE = GEN)
+  rename(GB_FEATURE = GEN)
 
 log_info("Saving results")
 write_tsv(

workflow/scripts/report/NV_description.R

@@ -38,10 +38,11 @@ date_order <- metadata %>%
   unique()
 
 empty_vcf <- tibble(
-  REGION = date_order,
+  SAMPLE = date_order,
+  REGION = as.character(NA),
   variant = as.character(NA),
   ALT_FREQ = as.numeric(NA),
-  GFF_FEATURE = as.character(NA),
+  GB_FEATURE = as.character(NA),
   synonimous = as.character(NA),
   POS = as.numeric(NA),
   ALT = as.character(NA),
@@ -52,10 +53,11 @@ empty_vcf <- tibble(
 # Create SNP variable and select useful variables
 vcf <- vcf %>%
   dplyr::select(
-    variant,
+    SAMPLE,
     REGION,
+    variant,
     ALT_FREQ,
-    GFF_FEATURE,
+    GB_FEATURE,
     synonimous,
     POS,
     ALT
@@ -91,7 +93,7 @@ vcf <- vcf %>%
       TRUE ~ "SNP"
     ),
     Class = case_when(
-      GFF_FEATURE == "Intergenic" ~ "Intergenic",
+      GB_FEATURE == "Intergenic" ~ "Intergenic",
       TRUE ~ synonimous
     ),
     POS = as.numeric(POS)
@@ -102,7 +104,7 @@ vcf <- vcf %>%
       NV_class == "INDEL" ~ str_length(ALT) - 1
     ),
     indel_class = case_when(
-      GFF_FEATURE == "Intergenic" ~ "Intergenic",
+      GB_FEATURE == "Intergenic" ~ "Intergenic",
       NV_class == "INDEL" &
         indel_len %% 3 == 0 ~
         "In frame",
@@ -114,7 +116,7 @@ vcf <- vcf %>%
   ungroup() %>%
   mutate(
     group = case_when(
-      GFF_FEATURE == "Intergenic" ~ "Intergenic",
+      GB_FEATURE == "Intergenic" ~ "Intergenic",
       NV_class == "SNP" ~ Class,
       NV_class == "INDEL" ~ indel_class
     )
@@ -156,7 +158,7 @@ variants <- vcf %>%
   ggplot() +
   aes(
     x = POS,
-    y = factor(REGION, date_order),
+    y = factor(SAMPLE, date_order),
     shape = factor(NV_class, c("SNP", "INDEL")),
     color = group,
     alpha = ALT_FREQ
@@ -312,7 +314,7 @@ variants_spike <- vcf_spike %>%
   ggplot() +
   aes(
     x = POS,
-    y = factor(REGION, date_order),
+    y = factor(SAMPLE, date_order),
     shape = factor(NV_class, c("SNP", "INDEL")),
     color = group,
     alpha = ALT_FREQ
@@ -363,9 +365,9 @@ figur_SNP_table <- vcf_snp %>%
   filter(ALT_FREQ <= snakemake@params$max_alt_freq) %>%
   left_join(
     metadata,
-    by = c("REGION" = "ID")
+    by = c("SAMPLE" = "ID")
   ) %>%
-  group_by(REGION) %>%
+  group_by(SAMPLE) %>%
   summarise(
     CollectionDate = min(as.Date(CollectionDate)),
     n = n_distinct(POS)
@@ -410,15 +412,15 @@ log_info("Saving plot tables")
 # Variants summary table
 vcf %>%
   select(
-    REGION,
+    SAMPLE,
     POS,
     variant,
     ALT_FREQ,
     NV_class,
     group
   ) %>%
   rename(
-    sample = REGION,
+    sample = SAMPLE,
     Variant = variant,
     Class = group
   ) %>%
@@ -444,18 +446,18 @@ window %>%
 # Heterozygous sites per sample table
 vcf_snp %>%
   filter(ALT_FREQ <= snakemake@params$max_alt_freq) %>%
-  select(!GFF_FEATURE) %>%
+  select(!GB_FEATURE) %>%
   left_join(
     metadata,
-    by = c("REGION" = "ID")
+    by = c("SAMPLE" = "ID")
   ) %>%
-  group_by(REGION) %>%
+  group_by(SAMPLE) %>%
   summarise(
     CollectionDate = min(as.Date(CollectionDate)),
     n_PolymorphicSites = n_distinct(POS)
   ) %>%
   ungroup() %>%
-  rename(sample = REGION) %>%
+  rename(sample = SAMPLE) %>%
   unique() %>%
   write.csv(snakemake@output[["table_3"]], row.names = FALSE)
 

workflow/scripts/report/evo_plots.R

@@ -21,10 +21,11 @@ N_S_position <- read_delim(snakemake@input[["N_S"]])
 # Create SNP variable and select useful variables
 vcf <- vcf %>%
   dplyr::select(
+    SAMPLE,
     variant,
     REGION,
     ALT_FREQ,
-    GFF_FEATURE,
+    GB_FEATURE,
     synonimous,
     POS
   )
@@ -37,14 +38,14 @@ metadata <- metadata %>%
     )
   ) %>%
   select(ID, interval) %>%
-  rename(REGION = ID)
+  rename(SAMPLE = ID)
 
 vcf <- left_join(vcf, metadata)
 
 # PLOT
 log_info("Ploting dN and dS over time")
 plot_df <- vcf %>%
-  group_by(REGION, synonimous) %>%
+  group_by(SAMPLE, synonimous) %>%
   summarise(
     Freq = sum(ALT_FREQ, na.rm = TRUE)
   ) %>%
@@ -67,7 +68,7 @@ plot_df <- vcf %>%
     values_to = "value",
     names_to = "d"
   ) %>%
-  left_join(unique(select(vcf, REGION, interval)))
+  left_join(unique(select(vcf, SAMPLE, interval)))
 
 plot <- plot_df %>%
   filter(d != "w") %>%
@@ -134,7 +135,7 @@ ggsave(
 # PLOT TABLES
 log_info("Saving plot table")
 vcf %>%
-  group_by(REGION, synonimous) %>%
+  group_by(SAMPLE, synonimous) %>%
   summarise(
     Freq = sum(ALT_FREQ, na.rm = TRUE)
   ) %>%
@@ -148,9 +149,9 @@ vcf %>%
     ds = Yes / sum(N_S_position$S)
   ) %>%
   ungroup() %>%
-  left_join(unique(select(vcf, REGION, interval))) %>%
+  left_join(unique(select(vcf, SAMPLE, interval))) %>%
   transmute(
-    sample = REGION,
+    sample = SAMPLE,
     DaysSinceFirst = interval,
     dN = dn,
     dS = ds,