Create E gene phylogenetic trees

phylogenetic/Snakefile

@@ -1,6 +1,6 @@
 configfile: "defaults/config_zika.yaml"
 
-genes = ['genome']
+genes = ['genome', 'ENV']
 
 wildcard_constraints:
     gene = "|".join(genes)
@@ -10,6 +10,7 @@ rule all:
         auspice_json = expand("auspice/zika_{gene}.json", gene=genes)
 
 include: "rules/prepare_sequences.smk"
+include: "rules/prepare_gene_sequences.smk"
 include: "rules/merge_sequences_usvi.smk"
 include: "rules/construct_phylogeny.smk"
 include: "rules/annotate_phylogeny.smk"

phylogenetic/defaults/config_zika.yaml

@@ -10,7 +10,9 @@ filter:
   group_by: "country year month"
   sequences_per_group: 40
   min_date: 2012
-  min_length: 5385
+  min_length:
+    genome: 5385
+    ENV: 1000
 
 refine:
   coalescent: "opt"

phylogenetic/rules/prepare_gene_sequences.smk

@@ -0,0 +1,57 @@
+"""
+This part of the workflow prepares reference files and sequences for constructing the gene phylogenetic trees.
+REQUIRED INPUTS:
+    reference   = path to reference sequence or genbank
+    sequences   = path to all sequences from which gene sequences can be extracted
+
+OUTPUTS:
+    gene_fasta = reference fasta for the gene (e.g. E gene)
+    gene_genbank = reference genbank for the gene (e.g. E gene)
+    sequences = sequences with gene sequences extracted and aligned to the reference gene sequence
+This part of the workflow usually includes the following steps:
+    - newreference.py: Creates new gene genbank and gene reference FASTA from the whole genome reference genbank
+    - nextclade: Aligns sequences to the reference gene sequence and extracts the gene sequences to ensure the reference files are valid
+See Nextclade or script usage docs for these commands for more details.
+"""
+
+rule generate_gene_reference_files:
+    """
+    Generating reference files for gene builds
+    """
+    input:
+        reference = "defaults/zika_reference.gb",
+    output:
+        fasta = "results/config/reference_{gene}.fasta",
+        genbank = "results/config/reference_{gene}.gb",
+    shell:
+        """
+        python3 scripts/newreference.py \
+            --reference {input.reference} \
+            --output-fasta {output.fasta} \
+            --output-genbank {output.genbank} \
+            --gene {wildcards.gene}
+        """
+
+
+rule align_and_extract_gene:
+    """
+    Cutting sequences to the length of the gene reference sequence
+    """
+    input:
+        sequences = "data/sequences_all.fasta",
+        reference = "results/config/reference_{gene}.fasta"
+    output:
+        sequences = "results/{gene}/sequences.fasta"
+    params:
+        min_length = lambda wildcard: config["filter"]["min_length"][wildcard.gene],
+    shell:
+        """
+        nextclade run \
+           -j 1 \
+           --input-ref {input.reference} \
+           --output-fasta {output.sequences} \
+           --min-seed-cover 0.01 \
+           --min-length {params.min_length} \
+           --silent \
+           {input.sequences}
+        """

phylogenetic/rules/prepare_sequences.smk

@@ -67,7 +67,7 @@ rule filter:
         group_by = config["filter"]["group_by"],
         sequences_per_group = config["filter"]["sequences_per_group"],
         min_date = config["filter"]["min_date"],
-        min_length = config["filter"]["min_length"],
+        min_length = lambda wildcard: config["filter"]["min_length"][wildcard.gene],
         strain_id = config.get("strain_id_field", "strain"),
     shell:
         """