Extract SARS-CoV-2 feature coordinates to JSON file

Author: ahmig

.test/targets.yaml

@@ -42,5 +42,7 @@ REPORT_QMD:
   "template.qmd"
 FEATURES_JSON:
   "config/sarscov2_features.json"
+COORDINATES_JSON:
+  "config/sarscov2_coordinates.json"
 GENETIC_CODE_JSON:
   "config/standard_genetic_code.json"

config/README.md

@@ -131,6 +131,7 @@ All of the following variables are pre-defined in [config.yaml](/config/config.y
 - `ALIGNMENT_REFERENCE`: NCBI accession number of the reference record for sequence alignment.
 - `PROBLEMATIC_VCF`: URL or path of a VCF file containing problematic genome positions for masking.
 - `FEATURES_JSON`: path of a JSON file containing name equivalences of genome features for data visualization.
+- `COORDINATES_JSON`: path of a JSON file containing genome feature coordinates for data visualization.
 - `GENETIC_CODE_JSON`: path of a JSON file containing a genetic code for gene translation.
 - `TREE_MODEL`: substitution model used by IQTREE (see [docs](http://www.iqtree.org/doc/Substitution-Models)).
 - `UFBOOT_REPS`: ultrafast bootstrap replicates for IQTREE (see [UFBoot](https://doi.org/10.1093/molbev/msx281)).

config/config.yaml

@@ -4,6 +4,8 @@ PROBLEMATIC_VCF:
   "https://raw.githubusercontent.com/W-L/ProblematicSites_SARS-CoV2/da322c32004f7b16bdaa6a8ee7fd24d56e79d9dc/problematic_sites_sarsCov2.vcf"
 FEATURES_JSON:
   "config/sarscov2_features.json"
+COORDINATES_JSON:
+  "config/sarscov2_coordinates.json"
 GENETIC_CODE_JSON:
   "config/standard_genetic_code.json"
 TREE_MODEL:

config/sarscov2_coordinates.json

@@ -0,0 +1,23 @@
+{
+    "five_prime_UTR": [1, 265],
+    "orf1ab": [266, 21555],
+    "Intergenic_1": [21556, 21562],
+    "S": [21563, 25384],
+    "Intergenic_2": [25385, 25392],
+    "ORF3a": [25393, 26220],
+    "Intergenic_3": [26221, 26244],
+    "E": [26245, 26472],
+    "Intergenic_4": [26473, 26522],
+    "M": [26523, 27191],
+    "Intergenic_5": [27192, 27201],
+    "ORF6": [27202, 27387],
+    "Intergenic_6": [27388, 27393],
+    "ORF7": [27394, 27759],
+    "Intergenic_7": [27760, 27893],
+    "ORF8": [27894, 28259],
+    "Intergenic_8": [28260, 28273],
+    "N": [28274, 29533],
+    "Intergenic_9": [29534, 29557],
+    "ORF10": [29558, 29674],
+    "three_prime_UTR": [29675, 29903]
+}

workflow/Snakefile

@@ -6,7 +6,7 @@ import subprocess
 
 min_version("7.19")
 
-__version__ = "1.2.2"
+__version__ = "1.3.0"
 
 containerized: "docker://ahmig/vipera:v" + __version__
 

workflow/rules/report.smk

@@ -99,6 +99,7 @@ rule general_NV_description:
         samples = expand("{sample}", sample = iter_samples()),
         design = config["PLOTS"],
         regions = config["PLOT_GENOME_REGIONS"],
+        coordinates = config["COORDINATES_JSON"],
         window = config["WINDOW"]["WIDTH"],
         step = config["WINDOW"]["STEP"],
         max_alt_freq = 1.0 - config["VC"]["IVAR_FREQ"]

workflow/scripts/report/NV_description.R

@@ -16,30 +16,10 @@ log_threshold(INFO)
 # Import file with plots style
 source(snakemake@params[["design"]])
 
-# SARS-CoV-2 anotation for genome scheme
-SCov2_annotation <- list(
-  "five_prime_UTR" = c(1:265),
-  "orf1ab" = c(266:21555),
-  "Intergenic_1" = c(21556:21562),
-  "S" = c(21563:25384),
-  "Intergenic_2" = c(25385:25392),
-  "ORF3a" = c(25393:26220),
-  "Intergenic_3" = c(26221:26244),
-  "E" = c(26245:26472),
-  "Intergenic_4" = c(26473:26522),
-  "M" = c(26523:27191),
-  "Intergenic_5" = c(27192:27201),
-  "ORF6" = c(27202:27387),
-  "Intergenic_6" = c(27388:27393),
-  "ORF7" = c(27394:27759),
-  "Intergenic_7" = c(27760:27893),
-  "ORF8" = c(27894:28259),
-  "Intergenic_8" = c(28260:28273),
-  "N" = c(28274:29533),
-  "Intergenic_9" = c(29534:29557),
-  "ORF10" = c(29558:29674),
-  "three_prime_UTR" = c(29675:29903)
-)
+log_info("Reading inputs")
+
+# Anotation for genome scheme
+coordinates <- read_json(snakemake@input$coordinates)
 
 vcf <- read_delim(snakemake@input[["vcf"]])
 vcf_snp <- vcf
@@ -78,20 +58,20 @@ vcf <- vcf %>%
     ALT
   )
 
-# Df with gene length for scheme
+# Create dataframe with gene length for scheme
 notation_empty <- data.frame(
   gene = "",
   len = 0
 ) %>%
   filter(len != 0)
 
 notation <- lapply(
-  names(SCov2_annotation),
+  names(coordinates),
   function(x) {
     add_row(
       notation_empty,
       gene = x,
-      len = length(SCov2_annotation[[x]])
+      len = length(coordinates[[x]])
     )
   }
 ) %>%
@@ -184,7 +164,7 @@ variants <- vcf %>%
     aes(
       x = len,
       y = 0.3,
-      fill = factor(gene, rev(names(SCov2_annotation)))
+      fill = factor(gene, rev(names(coordinates)))
     ),
     inherit.aes = FALSE,
     width = 0.3