Feature snpeff

workflow/envs/snpeff.yaml

@@ -0,0 +1,4 @@
+channels:
+  - bioconda
+dependencies:
+  - snpeff==5.1d

workflow/rules/vaf.smk

@@ -110,8 +110,55 @@ rule filter_tsv:
         tsv = OUTDIR/f"{OUTPUT_NAME}.masked.tsv",
         annotation = OUTDIR/"annotation.csv"
     output:
-        filtered_tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv"
+        filtered_tsv = OUTDIR/f"{OUTPUT_NAME}.masked.prefiltered.tsv"
     log:
         LOGDIR / "filter_tsv" / "log.txt"
     script:
         "../scripts/filter_tsv.R"
+
+rule tsv_to_vcf:
+    threads: 1
+    conda: "../envs/biopython.yaml"
+    input: 
+        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.prefiltered.tsv",
+    output:
+        vcf = OUTDIR/f"{OUTPUT_NAME}.vcf"
+    log:
+        LOGDIR / "tsv_to_vcf" / "log.txt"
+    script:
+        "../scripts/tsv_to_vcf.py"
+
+rule variants_effect:
+    threads: 1
+    conda: "../envs/snpeff.yaml"
+    params:
+        ref_name = config["ALIGNMENT_REFERENCE"]
+    input:
+        vcf = OUTDIR/f"{OUTPUT_NAME}.vcf"
+    output:
+        ann_vcf = OUTDIR/f"{OUTPUT_NAME}.annotated.vcf"
+    log:
+        LOGDIR / "variants_effect" / "log.txt"
+    shell:
+        """
+        exec >{log}                                                                    
+        exec 2>&1
+
+        snpEff eff {params.ref_name} {input.vcf} > {output.ann_vcf} || true
+        rm snpEff_genes.txt snpEff_summary.html
+        """
+
+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"
+    output:
+        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv"
+    log:
+        LOGDIR / "vcf_to_tsv" / "log.txt"
+    script:
+        "../scripts/vcf_to_tsv.R"
+
+

workflow/scripts/report/NV_description.R

@@ -57,21 +57,17 @@ date_order <- read_csv(snakemake@params[["metadata"]]) %>%
   unique()
 
 # Create SNP variable and select useful variables
+
 vcf <- vcf %>%
-  mutate(
-    SNP = paste(REF, POS, ALT, sep = "-")
-    ) %>%
   dplyr::select(
-    SNP,
+    variant,
     REGION,
     ALT_FREQ,
     GFF_FEATURE,
-    synonimous
-    ) %>%
-  rowwise() %>%
-  mutate(POS = strsplit(SNP, "-")[[1]][2]) %>%
-  ungroup()
-
+    synonimous,
+    POS,
+    ALT
+      )
 
 # Df with gene length for scheme
 notation_empty <- data.frame(
@@ -98,8 +94,8 @@ log_info("Classifying variants")
 vcf <- vcf %>%
   mutate(
     NV_class = case_when(
-      str_detect(SNP, fixed("--")) |
-      str_detect(SNP, fixed("+")) ~ "INDEL",
+      str_detect(ALT, fixed("-")) |
+      str_detect(ALT, fixed("+")) ~ "INDEL",
       TRUE ~ "SNP"
       ),
     Class = case_when(
@@ -111,12 +107,8 @@ vcf <- vcf %>%
   rowwise() %>%
   mutate(
     indel_len = case_when(
-      NV_class == "INDEL" &
-        str_detect(SNP, fixed("--")) ~
-        str_length(strsplit(SNP, "--")[[1]][2]),
-      NV_class == "INDEL" &
-        str_detect(SNP, fixed("-+")) ~
-        str_length(strsplit(SNP, "-+")[[1]][2])
+      NV_class == "INDEL"  ~
+        str_length(ALT) - 1
       ),
     indel_class = case_when(
       GFF_FEATURE == "Intergenic" ~ "Intergenic",
@@ -362,8 +354,6 @@ ggsave(
 log_info("Plotting nº of heterozygus sites for each sample")
 figur_SNP_time <- vcf_snp %>%
   filter(ALT_FREQ <= 0.95) %>%
-  select(!GFF_FEATURE) %>%
-  unique() %>%
   left_join(
     read_csv(snakemake@params[["metadata"]]),
     by = c("REGION" = "ID")
@@ -409,14 +399,14 @@ vcf %>%
   select(
     REGION,
     POS,
-    SNP,
+    variant,
     ALT_FREQ,
     NV_class,
     group
   ) %>%
   rename(
     sample = REGION,
-    NV = SNP,
+    Variant = variant,
     Class = group
   ) %>%
   filter(ALT_FREQ > 0) %>%
@@ -462,11 +452,11 @@ vcf_snp %>%
 
 n_indels <- vcf %>%
   filter(NV_class == "INDEL") %>%
-  pull(SNP) %>%
+  pull(variant) %>%
   unique() %>%
   length()
 
-n_snv <- length(unique(vcf$SNP)) - n_indels
+n_snv <- length(unique(vcf$variant)) - n_indels
 
 list(
   "INDELS" = n_indels,

workflow/scripts/report/evo_plots.R

@@ -22,18 +22,14 @@ N_S_position <- read_delim(snakemake@input[["N_S"]])
 
 # Create SNP variable and select useful variables
 vcf <- vcf %>%
-  mutate(
-    SNP = paste(REF, POS, ALT, sep = "-")) %>%
   dplyr::select(
-    SNP,
+    variant,
     REGION,
     ALT_FREQ,
     GFF_FEATURE,
-    synonimous
-  ) %>%
-  rowwise() %>%
-  mutate(POS = strsplit(SNP, "-")[[1]][2]) %>%
-  ungroup()
+    synonimous,
+    POS
+  ) 
 
 # Create variable for days sins first sample in metadata
 metadata <- metadata %>%

workflow/scripts/report/heatmap.R

@@ -16,27 +16,14 @@ date_order <- read_csv(snakemake@params[["metadata"]]) %>%
 
 # Create SNP variable and select useful variables from vcf
 vcf <- vcf %>%
-  mutate(
-    SNP = case_when(
-      !is.na(REF_AA) ~ paste(
-        GFF_FEATURE,
-        ":",
-        REF_AA,
-        POS_AA,
-        ALT_AA,
-        sep = ""
-      ),
-      TRUE ~ paste(REF, POS, ALT, sep = "")
-    )
-  ) %>%
-  unique() %>%
-  dplyr::select(SNP, REGION, ALT_FREQ)
+  dplyr::select(variant, REGION, ALT_FREQ)
 
 vcf <- vcf %>%
   pivot_wider(
-    names_from = SNP,
+    names_from = variant,
     values_from = ALT_FREQ,
-    values_fill = 0
+    values_fill = 0,
+    values_fn = sum
   ) %>%
   arrange(factor(REGION, levels = date_order)) %>%
   column_to_rownames(var = "REGION")

workflow/scripts/report/snp_plots.R

@@ -36,34 +36,8 @@ date_order <- read_csv(snakemake@params[["metadata"]]) %>%
 
 # Simplify features names and create SNP variable
 vcf <- vcf %>%
-  mutate(
-    SNP = case_when(
-      !is.na(REF_AA) ~ paste(
-        GFF_FEATURE,
-        ":",
-        REF_AA,
-        POS_AA,
-        ALT_AA,
-        sep = ""
-        ),
-      GFF_FEATURE != "Intergenic" ~ paste(
-        GFF_FEATURE,
-        ":",
-        POS - 1,
-        "-",
-        ALT
-      ),
-      TRUE ~ paste(
-        REF,
-        POS,
-        ALT,
-        sep = ""
-        )
-      )
-    ) %>%
-  unique() %>%
   dplyr::select(
-    SNP,
+    variant,
     REGION,
     ALT_FREQ,
     POS
@@ -107,7 +81,7 @@ vcf <- arrange(
 # Get list with all different polymorphisms
 SNPs <- pull(
     vcf,
-    SNP
+    variant
     ) %>%
   unique()
 
@@ -124,7 +98,7 @@ cor.df.fill <- lapply(
     function(snp) {
       df <- filter(
         vcf,
-        SNP == snp
+        variant == snp
       )
 
       test <- cor.test(
@@ -196,14 +170,14 @@ sign <- filter(
 # SNPs which are in positions with more than one alternative allele
 dup <- vcf %>%
   select(
-    SNP,
+    variant,
     POS
   ) %>%
   unique() %>%
   group_by(POS) %>%
   filter(n() > 1) %>%
   ungroup() %>%
-  pull(SNP) %>%
+  pull(variant) %>%
   unique()
 
 subset <- c(sign, dup) %>%
@@ -214,12 +188,12 @@ plot.height <- ceiling(length(subset) / 4) * 42
 
 log_info("PLotting SNPs trends in time")
 panel <- vcf %>%
-  filter(SNP %in% subset) %>%
+  filter(variant %in% subset) %>%
   ggplot() +
   aes(
     x = interval,
     y = ALT_FREQ,
-    color = SNP
+    color = variant
   ) +
   scale_color_viridis_d() +
   geom_point() +
@@ -269,13 +243,14 @@ cor.df.fill %>%
 
 log_info("Saving SNPs trends table")
 vcf %>%
-  filter(SNP %in% subset) %>%
+  filter(variant %in% subset) %>%
   transmute(
     sample = REGION,
     POS = POS,
-    NV = SNP,
+    NV = variant,
     ALT_FREQ = ALT_FREQ,
     DaysSinceFirst = interval
     ) %>%
     write.csv(snakemake@output[["table_2"]], row.names = FALSE)
 
+

workflow/scripts/tsv_to_vcf.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python3
+
+import pandas as pd
+import sys
+
+def tsv_to_vcf(tsv_file, vcf_file):
+    # Read the TSV file
+    tsv_df = pd.read_csv(tsv_file, sep='\t')
+
+    # Open a new VCF file for writing
+    with open(vcf_file, 'w') as vcf:
+        # Write the VCF header
+        vcf.write('##fileformat=VCFv4.2\n')
+        vcf.write('#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n')
+
+        # Process each row in the TSV file
+        for index, row in tsv_df.iterrows():
+            # Extract fields from the TSV row
+            chrom = 'NC_045512.2'
+            pos = row['POS']
+            ref = row['REF']
+            alt = row['ALT']
+            # Add other necessary fields for the VCF file
+
+            # Handle INDELs
+            if alt[0] == "+":
+                alt = ref + alt[1:]
+            elif alt[0] == "-":
+                alt2 = ref
+                ref += alt[1:]
+                alt = alt2  # Adjust this line to use the correct value for ALT after a deletion
+
+            # Write the VCF row
+            vcf.write(f'{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\t.\n')
+
+def main():
+    input_tsv_file = snakemake.input.tsv
+    output_vcf_file = snakemake.output.vcf
+    tsv_to_vcf(input_tsv_file, output_vcf_file)
+
+
+if __name__ == '__main__':
+    main()

workflow/scripts/vcf_to_tsv.R

@@ -0,0 +1,32 @@
+#!/usr/bin/env Rscript
+
+library(tidyverse)
+library(logger)
+log_threshold(INFO)
+
+# Write stdout and stderr to log file
+log <- file(snakemake@log[[1]], open = "wt")
+sink(log, type = "message")
+sink(log, type = "output")
+
+# read data
+log_info("Reading data")
+vcf <- read_tsv(snakemake@input[["ann_vcf"]], comment = "##")
+tsv <- read_tsv(snakemake@input[["pre_tsv"]])
+
+tsv["variant"] <- str_extract(vcf$INFO, "p\\.([^|]*)", group = 1)
+tsv["nuc_variant"] <- str_extract(vcf$INFO, "c\\.([^|]*)", group = 1)
+
+tsv <- tsv %>%
+    mutate(
+        variant = case_when(
+            is.na(variant) ~ paste(POS, REF, ">", ALT, sep = ""),
+            TRUE ~ paste(GFF_FEATURE, ":", variant, sep = "")
+        )
+    )
+
+log_info("Saving results")
+write_tsv(
+    tsv,
+    snakemake@output[["tsv"]]
+)