Make tree for 450bp of the N gene ("N450")

phylogenetic/Snakefile

@@ -5,6 +5,7 @@ rule all:
         auspice_json = "auspice/measles.json",
 
 include: "rules/prepare_sequences.smk"
+include: "rules/prepare_sequences_N450.smk"
 include: "rules/construct_phylogeny.smk"
 include: "rules/annotate_phylogeny.smk"
 include: "rules/export.smk"

phylogenetic/defaults/config.yaml

@@ -2,13 +2,20 @@ strain_id_field: "accession"
 files:
     exclude: "defaults/dropped_strains.txt"
     reference: "defaults/measles_reference.gb"
+    reference_N450: "defaults/measles_reference_N450.gb"
+    reference_N450_fasta: "defaults/measles_reference_N450.fasta"
     colors: "defaults/colors.tsv"
     auspice_config: "defaults/auspice_config.json"
 filter: 
     group_by: "country year month"
     sequences_per_group: 20
     min_date: 1950
     min_length: 5000
+filter_N450:
+    group_by: "country year"
+    subsample_max_sequences: 3000
+    min_date: 1950
+    min_length: 400
 refine:
     coalescent: "opt"
     date_inference: "marginal"

phylogenetic/rules/prepare_sequences.smk

@@ -74,7 +74,7 @@ rule align:
         sequences = "results/filtered.fasta",
         reference = config["files"]["reference"]
     output:
-        alignment = "results/aligned.fasta"
+        alignment = "results/aligned_genome.fasta"
     shell:
         """
         augur align \

phylogenetic/rules/prepare_sequences_N450.smk

@@ -0,0 +1,58 @@
+"""
+This part of the workflow prepares sequences for constructing the phylogenetic tree for 450bp of the N gene.
+
+See Augur's usage docs for these commands for more details.
+"""
+
+rule align_and_extract_N450:
+    input:
+        sequences = "data/sequences.fasta",
+        reference = config["files"]["reference_N450_fasta"]
+    output:
+        sequences = "results/sequences_N450.fasta"
+    params:
+        min_length = config['filter_N450']['min_length']
+    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}
+        """
+rule filter_N450:
+    """
+    Filtering to
+      - {params.sequences_per_group} sequence(s) per {params.group_by!s}
+      - excluding strains in {input.exclude}
+      - minimum genome length of {params.min_length}
+      - excluding strains with missing region, country or date metadata
+    """
+    input:
+        sequences = "results/sequences_N450.fasta",
+        metadata = "data/metadata.tsv",
+        exclude = config["files"]["exclude"]
+    output:
+        sequences = "results/aligned_N450.fasta"
+    params:
+        group_by = config['filter_N450']['group_by'],
+        subsample_max_sequences = config["filter_N450"]["subsample_max_sequences"],
+        min_date = config["filter_N450"]["min_date"],
+        min_length = config['filter_N450']['min_length'],
+        strain_id = config["strain_id_field"]
+    shell:
+        """
+        augur filter \
+            --sequences {input.sequences} \
+            --metadata {input.metadata} \
+            --metadata-id-columns {params.strain_id} \
+            --exclude {input.exclude} \
+            --output {output.sequences} \
+            --group-by {params.group_by} \
+            --subsample-max-sequences {params.subsample_max_sequences} \
+            --min-date {params.min_date} \
+            --min-length {params.min_length}
+        """