WIP: SH tree

Author: j23414

phylogenetic/Snakefile

@@ -8,7 +8,7 @@ workdir: workflow.current_basedir
 # Use default configuration values. Override with Snakemake's --configfile/--config options.
 configfile: "defaults/config.yaml"
 
-builds = ["north-america","global"]
+builds = ["north-america","global","sh"]
 
 wildcard_constraints:
     build = "|".join(builds)

phylogenetic/defaults/color_orderings.tsv

@@ -298,11 +298,17 @@ MuV_genotype	A
 MuV_genotype	B
 MuV_genotype	C
 MuV_genotype	D
+MuV_genotype	D1
 MuV_genotype	F
 MuV_genotype	G
+MuV_genotype	G1
+MuV_genotype	G2
 MuV_genotype	H
+MuV_genotype	H1
+MuV_genotype	H2
 MuV_genotype	I
 MuV_genotype	J
 MuV_genotype	K
+MuV_genotype	K/M
 MuV_genotype	L
 MuV_genotype	N

phylogenetic/defaults/config.yaml

@@ -4,25 +4,30 @@ sequences_url: "https://data.nextstrain.org/files/workflows/mumps/sequences.fast
 metadata_url: "https://data.nextstrain.org/files/workflows/mumps/metadata.tsv.zst"
 
 strain_id_field: "accession"
-reference: "defaults/reference.gb"
+reference:
+  north-america: "defaults/reference.gb"
+  global: "defaults/reference.gb"
+  sh: "defaults/sh/reference.gb"
 
 filter:
-  min_length: 8000
   group_by: country year month MuV_genotype division
   specific:
-    north-america: --subsample-max-sequences 4000 --min-date 2006 --query "region=='North America' & (MuV_genotype=='G')"
-    global: --subsample-max-sequences 4000 --min-date 1950
+    north-america: --min-length 8000 --subsample-max-sequences 4000 --min-date 2006 --query "region=='North America' & (MuV_genotype=='G')"
+    global: --min-length 8000 --subsample-max-sequences 4000 --min-date 1950
+    sh: --exclude-all
 
 refine:
-  north-america: "--clock-filter-iqd 4"
-  global: ""
+  north-america: "--timetree --clock-filter-iqd 4 --coalescent opt --date-confidence --date-inference marginal"
+  global: "--timetree --coalescent opt --date-confidence --date-inference marginal"
+  sh: "--timetree --clock-rate 0.0001"
 
 ancestral:
   inference: "joint"
 
 traits:
   north-america: country division MuV_genotype
   global: region MuV_genotype
+  sh: region
   sampling_bias_correction: 3
 
 colors:

phylogenetic/rules/annotate_phylogeny.smk

@@ -59,7 +59,7 @@ rule translate:
     input:
         tree = "results/{build}/tree.nwk",
         node_data = "results/{build}/nt_muts.json",
-        reference = config['reference'],
+        reference = lambda wildcard: config['reference'][wildcard.build],
     output:
         node_data = "results/{build}/aa_muts.json",
     log:

phylogenetic/rules/construct_phylogeny.smk

@@ -56,8 +56,6 @@ rule refine:
     benchmark:
         "benchmarks/{build}/refine.txt",
     params:
-        coalescent = "opt",
-        date_inference = "marginal",
         clock_filter_iqd = lambda wildcard: config['refine'][wildcard.build],
         strain_id = config.get("strain_id_field", "strain"),
     shell:
@@ -69,9 +67,5 @@ rule refine:
             --metadata-id-columns {params.strain_id:q} \
             --output-tree {output.tree:q} \
             --output-node-data {output.node_data:q} \
-            --timetree \
-            --coalescent {params.coalescent:q} \
-            --date-confidence \
-            --date-inference {params.date_inference:q} \
             {params.clock_filter_iqd} 2>&1 | tee {log:q}
         """

phylogenetic/rules/export.smk

@@ -51,7 +51,7 @@ rule export:
     """Exporting data files for for auspice"""
     input:
         tree = "results/{build}/tree.nwk",
-        metadata = "data/metadata.tsv",
+        metadata = "results/sh/metadata_merged.tsv",
         branch_lengths = "results/{build}/branch_lengths.json",
         traits = "results/{build}/traits.json",
         nt_muts = "results/{build}/nt_muts.json",

phylogenetic/rules/prepare_sequences.smk

@@ -50,14 +50,31 @@ rule decompress:
         zstd -d -c {input.metadata} > {output.metadata}
         """
 
+rule merge_annotations:
+    """Merge identical sequence annotations"""
+    input:
+        metadata = "data/metadata.tsv",
+        identical = "defaults/sh/metadata_identical.tsv",
+    output:
+        merged_metadata = "data/metadata_merged.tsv",
+    params:
+        id_column = "accession",
+    shell:
+        r"""
+        augur merge --metadata \
+          a={input.metadata:q} \
+          b={input.identical:q} \
+          --metadata-id-columns {params.id_column} \
+          --output-metadata {output.merged_metadata}
+        """
+
 rule filter:
     """
     Filtering to
       - various criteria based on the auspice JSON target
       - from {params.min_date} onwards
       - excluding strains in {input.exclude}
       - including strains in {input.include}
-      - minimum genome length of {params.min_length} (50% of Zika virus genome)
     """
     input:
         sequences = "data/sequences.fasta",
@@ -72,7 +89,6 @@ rule filter:
     benchmark:
         "benchmarks/{build}/filtered.txt",
     params:
-        min_length = config['filter']['min_length'],
         group_by = config['filter']['group_by'],
         filter_params = lambda wildcard: config['filter']['specific'][wildcard.build],
         strain_id = config.get("strain_id_field", "strain"),
@@ -86,19 +102,63 @@ rule filter:
             --include {input.include:q} \
             --output {output.sequences:q} \
             --output-metadata {output.metadata:q} \
-            --min-length {params.min_length:q} \
             --group-by {params.group_by} \
             {params.filter_params} 2>&1 | tee {log:q}
         """
 
+ruleorder: filter_sh > filter
+
+rule filter_sh:
+    """
+    Filtering to
+      - various criteria based on the auspice JSON target
+      - from {params.min_date} onwards
+      - excluding strains in {input.exclude}
+      - including strains in {input.include}
+    """
+    input:
+        sequences = "data/sequences.fasta",
+        metadata = "data/metadata.tsv",
+        clade_membership = "defaults/sh/metadata_duplicate.txt",
+        exclude = "defaults/sh/exclude.txt",
+        include = "defaults/sh/include.txt"
+    output:
+        merged_metadata = "results/sh/metadata_merged.tsv",
+        sequences = "results/sh/filtered.fasta",
+        metadata = "results/sh/metadata.tsv",
+    log:
+        "logs/sh/filtered.txt",
+    benchmark:
+        "benchmarks/sh/filtered.txt",
+    params:
+        group_by = config['filter']['group_by'],
+        filter_params = config['filter']['specific']['sh'],
+        strain_id = config.get("strain_id_field", "strain"),
+    shell:
+        r"""
+        augur merge \
+        --metadata a={input.metadata:q} b={input.clade_membership:q} \
+        --metadata-id-columns a={params.strain_id:q} b={params.strain_id:q} \
+        --output-metadata {output.merged_metadata:q}
+
+        augur filter \
+            --sequences {input.sequences:q} \
+            --metadata {output.merged_metadata:q} \
+            --metadata-id-columns {params.strain_id:q} \
+            --include {input.include:q} \
+            --output-sequences {output.sequences:q} \
+            --output-metadata {output.metadata:q} \
+            {params.filter_params} 2>&1 | tee {log:q}
+        """
+
 rule align:
     """
     Aligning sequences to {input.reference}
       - filling gaps with N
     """
     input:
         sequences = "results/{build}/filtered.fasta",
-        reference = config['reference'],
+        reference = lambda wildcard: config['reference'][wildcard.build],
     output:
         alignment = "results/{build}/aligned.fasta",
     log: