Extract NV panel regions from GenBank file

Removes the need of a JSON file with coordinates, but is flexible enough to allow its usage

Author: ahmig

.test/targets.yaml

@@ -40,7 +40,5 @@ PLOT_GENOME_REGIONS:
   "config/nsp_annotation.csv"
 REPORT_QMD:
   "template.qmd"
-COORDINATES_JSON:
-  "config/sarscov2_coordinates.json"
 GENETIC_CODE_JSON:
   "config/standard_genetic_code.json"

config/README.md

@@ -130,7 +130,6 @@ 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.
-- `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` & `SHALRT`: settings for ultrafast bootstrap (see [UFBoot](https://doi.org/10.1093/molbev/msx281))

config/config.yaml

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

config/sarscov2_coordinates.json

@@ -95,6 +95,20 @@ rule diversity_plot:
         "../scripts/report/diversity_plot.R"
 
 
+rule extract_genbank_regions:
+    conda: "../envs/biopython.yaml"
+    params:
+        gb_qualifier = "gene",
+    input:
+        gb = OUTDIR/"reference.cds.gb",
+    output:
+        regions = temp(REPORT_DIR_TABLES/"genbank_regions.json"),
+    log:
+        LOGDIR / "extract_genbank_regions" / "log.txt"
+    script:
+        "../scripts/report/extract_genbank_regions.py"
+
+
 rule general_NV_description:
     conda: "../envs/renv.yaml"
     params:
@@ -104,15 +118,15 @@ rule general_NV_description:
         step = config["WINDOW"]["STEP"],
         max_alt_freq = 1.0 - config["VC"]["IVAR_FREQ"]
     input:
-        coordinates = config["COORDINATES_JSON"],
+        coordinates = REPORT_DIR_TABLES/"genbank_regions.json",
         regions = config["PLOT_GENOME_REGIONS"],
         window = OUTDIR/f"{OUTPUT_NAME}.window.csv",
         vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"]
     output:
-        fig = report(REPORT_DIR_PLOTS/"figure_5a.png"),
-        fig_s = report(REPORT_DIR_PLOTS/"figure_5b.png"),
-        fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
+        fig = report(REPORT_DIR_PLOTS/"figure_5a.png"),      # panel
+        fig_s = report(REPORT_DIR_PLOTS/"figure_5b.png"),    # panel spike
+        fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),   # SNPs vs time
         json = temp(OUTDIR/"summary_nv.json"),
         table_1 = report(REPORT_DIR_TABLES/"figure_5a.csv"),
         table_2 = report(REPORT_DIR_TABLES/"figure_5b.csv"),

workflow/rules/report.smk

@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+
+import logging
+import json
+
+from Bio import SeqIO
+from Bio.SeqFeature import ExactPosition
+
+
+if __name__ == "__main__":
+    logging.basicConfig(
+        filename=snakemake.log[0], format=snakemake.config["LOG_PY_FMT"],
+        level=logging.INFO
+    )
+
+    logging.info("Reading GenBank file")
+    gb = SeqIO.read(snakemake.input.gb, format="gb")
+
+    logging.info("Calculating CDS regions")
+    regions = {}
+    for feature in gb.features:
+        identifier = "|".join(feature.qualifiers.get(snakemake.params.gb_qualifier, []))
+        if identifier == "":
+            logging.error(f"Feature at {feature.location} has no qualifier '{snakemake.params.gb_qualifier_display}'")
+        elif identifier in regions:
+            logging.warning(f"Identifier '{identifier}' is already among coding records and will not be replaced by feature at {feature.location}")
+        else:
+            logging.debug(f"Adding feature")
+            if type(feature.location.start) is not ExactPosition or type(feature.location.start) is not ExactPosition:
+                logging.warning(f"Feature '{identifier}' location is not exact but will be treated as such: {feature.location}")
+            regions[identifier] = (
+                int(feature.location.start.real + 1),
+                int(feature.location.end.real)
+            )
+    
+    logging.info("Calculating intergenic regions")
+    cds_names = tuple(regions.keys())
+    intergenic_count = 0
+    for i in range(1, len(cds_names)):
+        start1, end1 = regions[cds_names[i-1]]
+        start2, end2 = regions[cds_names[i]]
+        if (start2 - end1) > 1:
+            intergenic_count += 1
+            regions[f"Intergenic_{intergenic_count}"] = (end1 + 1, start2 - 1)
+    
+    logging.info("Writing regions to JSON file")
+    with open(snakemake.output.regions, "w") as fw:
+        json.dump(regions, fw, indent=2)