Add E gene trees

phylogenetic/README.md

@@ -62,7 +62,7 @@ example sequences into the `data/` directory like so:
 
 With access to AWS, this can be more quickly run as:
 
-    nextstrain build --aws-batch --aws-batch-cpus 4 --aws-batch-memory 7200 . --jobs 4
+    nextstrain build --aws-batch --aws-batch-cpus 10 --aws-batch-memory 7200 . --jobs 10
 
 [Nextstrain]: https://nextstrain.org
 [augur]: https://docs.nextstrain.org/projects/augur/en/stable/

phylogenetic/Snakefile

@@ -1,12 +1,18 @@
 configfile: "config/config_dengue.yaml"
 
 serotypes = ['all', 'denv1', 'denv2', 'denv3', 'denv4']
+genes = ['E', 'genome']
+
+wildcard_constraints:
+    serotype = "|".join(serotypes),
+    gene = "|".join(genes)
 
 rule all:
     input:
-        auspice_json = expand("auspice/dengue_{serotype}_genome.json", serotype=serotypes)
+        auspice_json = expand("auspice/dengue_{serotype}_{gene}.json", serotype=serotypes, gene=genes),
 
 include: "rules/prepare_sequences.smk"
+include: "rules/prepare_sequences_E.smk"
 include: "rules/construct_phylogeny.smk"
 include: "rules/annotate_phylogeny.smk"
 include: "rules/export.smk"

phylogenetic/bin/newreference.py

@@ -0,0 +1,46 @@
+from Bio import SeqIO
+from Bio.SeqRecord import SeqRecord
+from Bio.SeqFeature import SeqFeature, FeatureLocation, Seq
+import shutil
+import argparse
+
+def new_reference(referencefile, outgenbank, outfasta, gene):
+    ref = SeqIO.read(referencefile, "genbank")
+    startofgene = None
+    endofgene = None
+    for feature in ref.features:
+        if feature.type == 'source':
+            ref_source_feature = feature
+        if feature.type =='gene' or feature.type == 'CDS':
+            a = list(feature.qualifiers.items())[0][-1][0]
+            if a == gene:
+                startofgene = int(list(feature.location)[0])
+                endofgene =  int(list(feature.location)[-1])+1
+
+    record = ref[startofgene:endofgene]
+    source_feature =  SeqFeature(FeatureLocation(start=0, end=len(record)), type='source',
+                                 qualifiers=ref_source_feature.qualifiers)
+    record.features.append(source_feature)
+
+    SeqIO.write(record, outgenbank, 'genbank')
+    SeqIO.write(record, outfasta, "fasta")
+
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        description="make new reference depending on whether the entire genome or only part is to be used for the tree",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+    parser.add_argument("--reference", required=True, help="GenBank file with reference sequences")
+    parser.add_argument("--output-fasta", required=True, help="GenBank new reference file")
+    parser.add_argument("--output-genbank", required=True, help="GenBank new reference file")
+    parser.add_argument("--gene", help="gene name or genome for entire genome")
+    args = parser.parse_args()
+
+    if args.gene=='genome':
+        shutil.copy(args.reference, args.output_genbank)
+        SeqIO.write(SeqIO.read(args.reference, 'genbank'), args.output_fasta, 'fasta')
+    else:
+        new_reference(args.reference, args.output_genbank, args.output_fasta, args.gene)
+

phylogenetic/config/config_dengue.yaml

@@ -5,21 +5,34 @@ filter:
   exclude: "config/exclude.txt"
   group_by: "year region"
   min_length: 5000
+  # Length of the E gene is approximately 1400
+  E_min_length: 1000
   sequences_per_group:
     all: '10'
     denv1: '36'
     denv2: '36'
     denv3: '36'
     denv4: '36'
+  E_root_sequence:
+    all: 'NC_002640'
+    denv1: 'NC_001477'
+    denv2: 'NC_001474'
+    denv3: 'NC_001475'
+    denv4: 'NC_002640'
 
 traits:
   sampling_bias_correction: '3'
   traits_columns:
-    all: 'region nextclade_subtype'
-    denv1: 'country region nextclade_subtype'
-    denv2: 'country region nextclade_subtype'
-    denv3: 'country region nextclade_subtype'
-    denv4: 'country region nextclade_subtype'
+    all_genome: 'region nextclade_subtype'
+    denv1_genome: 'country region nextclade_subtype'
+    denv2_genome: 'country region nextclade_subtype'
+    denv3_genome: 'country region nextclade_subtype'
+    denv4_genome: 'country region nextclade_subtype'
+    all_E: 'region nextclade_type'
+    denv1_E: 'country region nextclade_subtype'
+    denv2_E: 'country region nextclade_subtype'
+    denv3_E: 'country region nextclade_subtype'
+    denv4_E: 'country region nextclade_subtype'
 
 clades:
   clade_definitions:

phylogenetic/config/exclude.txt

@@ -44,7 +44,6 @@ ON046272 # 01488/09 # too divergent
 ON046271 # 01542/04 # too divergent
 MZ284953 # dev1 # too divergent
 MZ215848 # DKE_121 # too divergent
-MW946564 # SENDAK_HD_10674 # sylvatic
 OK605757 # DENV2_1_DAK_HD_76395 # sylvatic
 MW945427 # DENV3/PUERTORICO/1963/PRS_228762_AC27 # too divergent
 OM258630 # PR_6 # too divergent
@@ -891,3 +890,35 @@ GU131859 # duplicate of GU131847
 GU131867 # duplicate of GU131862
 HQ332171 # duplicate of HQ332170
 HQ332178 # duplicate of HQ332177
+OP909734 # DENV2-16681-rcCap-NS1 # too divergent E gene
+KR919820 # Brun2014 # highly divergent denv1
+MZ285058 # denv1 # This record was removed from NCBI
+OP741244 # D20 # Too divergent in Denv1 tree
+OP741243 # 05 # Too divergent in Denv1 tree
+ON127553 # BCSIR-IBS-D19 # Too divergent in Denv1 tree
+ON127552 # NC_001475.2 # Too divergent in Denv1 tree
+ON127551 # BCSIR-IBS-D18 # Too divergent in Denv1 tree
+ON127550 # NC_001475.2 # Too divergent in Denv1 tree
+ON127538 # BCSIR-IBS-D17 # Too divergent in Denv1 tree
+ON127537 # NC_001475.2 # Too divergent in Denv1 tree
+ON127534 # strain_1 # Too divergent in Denv1 tree
+ON127533 # NC_001475.2 # Too divergent in Denv1 tree
+ON127419 # BCSIR-IBS-D14 # Too divergent in Denv1 tree
+ON116130 # BCSIR-IBS-D13 # Too divergent in Denv1 tree
+ON116129 # NC_001475.2 # Too divergent in Denv1 tree
+ON115818 # NC_001475.2 # Too divergent in Denv1 tree
+ON115817 # BCSIR-IBS-D11 # Too divergent in Denv1 tree
+ON115814 # NC_001475.2 # Too divergent in Denv1 tree
+ON115812 # BCSIR-IBS-D9 # Too divergent in Denv1 tree
+ON115811 # NC_001475.2 # Too divergent in Denv1 tree
+ON115183 # BCSIR-IBS-D8 # Too divergent in Denv1 tree
+ON115030 # BCSIR-IBS-D7 # Too divergent in Denv1 tree
+ON115026 # BCSIR-IBS-D6 # Too divergent in Denv1 tree
+ON115025 # NC_001475.2 # Too divergent in Denv1 tree
+ON115024 # BCSIR-IBS-D5 # Too divergent in Denv1 tree
+ON115023 # NC_001475.2 # Too divergent in Denv1 tree
+ON115022 # BCSIR-IBS-D4 # Too divergent in Denv1 tree
+ON115021 # NC_001475.2 # Too divergent in Denv1 tree
+ON111445 # NC_001475.2 # Too divergent in Denv1 tree
+ON103303 # BCSIR-IBS-D1 # Too divergent in Denv1 tree
+OK079104 # BCSIR-05 # Too divergent in Denv1 tree

phylogenetic/config/reference_dengue_denv2.gb → phylogenetic/config/reference_denv2_genome.gb

@@ -1,3 +1,4 @@
+
 LOCUS       DENV2/THAILAND/REFERENCE/1964                  10723 bp    DNA     VRL 15-SEP-2015
 DEFINITION  Dengue virus 2, complete genome.
 ACCESSION   NC_001474
@@ -49,13 +50,13 @@ FEATURES             Location/Qualifiers
                      /db_xref="GeneID:1494449"
                      /gene="flavivirus polyprotein gene"
      CDS             97..438
-                     /db_xref="VBRC:35917"
                      /gene="C"
+                     /db_xref="VBRC:35917"
                      /product="anchored capsid protein C"
                      /protein_id="NP_739581.2"
      CDS             439..936
-                     /db_xref="VBRC:35919"
                      /gene="M"
+                     /db_xref="VBRC:35919"
                      /product="membrane glycoprotein precursor M"
                      /protein_id="NP_739582.2"
      CDS             439..711
@@ -64,50 +65,50 @@ FEATURES             Location/Qualifiers
                      /product="protein pr"
                      /protein_id="YP_009164954.1"
      CDS             937..2421
-                     /db_xref="VBRC:35921"
                      /gene="E"
+                     /db_xref="VBRC:35921"
                      /product="envelope protein E"
                      /protein_id="NP_739583.2"
      CDS             2422..3477
-                     /db_xref="VBRC:35922"
                      /gene="NS1"
+                     /db_xref="VBRC:35922"
                      /product="nonstructural protein NS1"
                      /protein_id="NP_739584.2"
      CDS             3478..4131
-                     /db_xref="VBRC:35923"
                      /gene="NS2A"
+                     /db_xref="VBRC:35923"
                      /product="nonstructural protein NS2A"
                      /protein_id="NP_739585.2"
      CDS             4132..4521
-                     /db_xref="VBRC:35924"
                      /gene="NS2B"
+                     /db_xref="VBRC:35924"
                      /product="nonstructural protein NS2B"
                      /protein_id="NP_739586.2"
      CDS             4522..6375
-                     /db_xref="VBRC:35925"
                      /gene="NS3"
+                     /db_xref="VBRC:35925"
                      /note="ATPase; component of capping enzyme (RNA
                      thriphosphatase); protease; RNA-helicase"
                      /product="nonstructural protein NS3"
                      /protein_id="NP_739587.2"
      CDS             6376..6756
-                     /db_xref="VBRC:35926"
                      /gene="NS4A"
+                     /db_xref="VBRC:35926"
                      /product="nonstructural protein NS4A"
                      /protein_id="NP_739588.2"
      CDS             6757..6825
-                     /db_xref="VBRC:35927"
                      /gene="2K"
+                     /db_xref="VBRC:35927"
                      /product="protein 2K"
                      /protein_id="NP_739593.2"
      CDS             6826..7569
-                     /db_xref="VBRC:35928"
                      /gene="NS4B"
+                     /db_xref="VBRC:35928"
                      /product="nonstructural protein NS4B"
                      /protein_id="NP_739589.2"
      CDS             7570..10269
-                     /db_xref="VBRC:35929"
                      /gene="NS5"
+                     /db_xref="VBRC:35929"
                      /note="methyltransferase component of capping enzyme;
                      nonstructural protein NS5"
                      /product="RNA-dependent RNA polymerase NS5"

phylogenetic/rules/annotate_phylogeny.smk

@@ -24,10 +24,10 @@ to the ones produced by Augur commands.
 rule ancestral:
     """Reconstructing ancestral sequences and mutations"""
     input:
-        tree = "results/tree_{serotype}.nwk",
-        alignment = "results/aligned_{serotype}.fasta"
+        tree = "results/tree_{serotype}_{gene}.nwk",
+        alignment = "results/aligned_{serotype}_{gene}.fasta"
     output:
-        node_data = "results/nt-muts_{serotype}.json"
+        node_data = "results/nt-muts_{serotype}_{gene}.json"
     params:
         inference = "joint"
     shell:
@@ -42,11 +42,12 @@ rule ancestral:
 rule translate:
     """Translating amino acid sequences"""
     input:
-        tree = "results/tree_{serotype}.nwk",
-        node_data = "results/nt-muts_{serotype}.json",
-        reference = "config/reference_dengue_{serotype}.gb"
+        tree = "results/tree_{serotype}_{gene}.nwk",
+        node_data = "results/nt-muts_{serotype}_{gene}.json",
+        # The genbank references for the E gene are dynamically generated files located within the results folder.
+        reference = lambda wildcard: "config/reference_{serotype}_{gene}.gb" if wildcard.gene in ["genome"] else "results/config/reference_{serotype}_{gene}.gb"
     output:
-        node_data = "results/aa-muts_{serotype}.json"
+        node_data = "results/aa-muts_{serotype}_{gene}.json"
     shell:
         """
         augur translate \
@@ -62,12 +63,12 @@ rule traits:
       - increase uncertainty of reconstruction by {params.sampling_bias_correction} to partially account for sampling bias
     """
     input:
-        tree = "results/tree_{serotype}.nwk",
+        tree = "results/tree_{serotype}_{gene}.nwk",
         metadata = "data/metadata_{serotype}.tsv"
     output:
-        node_data = "results/traits_{serotype}.json",
+        node_data = "results/traits_{serotype}_{gene}.json",
     params:
-        columns = lambda wildcards: config['traits']['traits_columns'][wildcards.serotype],
+        columns = lambda wildcards: config['traits']['traits_columns'][wildcards.serotype + "_" + wildcards.gene],
         sampling_bias_correction = config['traits']['sampling_bias_correction'],
         strain_id = config.get("strain_id_field", "strain"),
     shell:
@@ -85,12 +86,12 @@ rule traits:
 rule clades:
     """Annotating serotypes / genotypes"""
     input:
-        tree = "results/tree_{serotype}.nwk",
-        nt_muts = "results/nt-muts_{serotype}.json",
-        aa_muts = "results/aa-muts_{serotype}.json",
+        tree = "results/tree_{serotype}_genome.nwk",
+        nt_muts = "results/nt-muts_{serotype}_genome.json",
+        aa_muts = "results/aa-muts_{serotype}_genome.json",
         clade_defs = lambda wildcards: config['clades']['clade_definitions'][wildcards.serotype],
     output:
-        clades = "results/clades_{serotype}.json"
+        clades = "results/clades_{serotype}_genome.json"
     shell:
         """
         augur clades \

phylogenetic/rules/construct_phylogeny.smk

@@ -15,9 +15,9 @@ See Augur's usage docs for these commands for more details.
 rule tree:
     """Building tree"""
     input:
-        alignment = "results/aligned_{serotype}.fasta"
+        alignment = "results/aligned_{serotype}_{gene}.fasta"
     output:
-        tree = "results/tree-raw_{serotype}.nwk"
+        tree = "results/tree-raw_{serotype}_{gene}.nwk"
     shell:
         """
         augur tree \
@@ -35,12 +35,12 @@ rule refine:
       - filter tips more than {params.clock_filter_iqd} IQDs from clock expectation
     """
     input:
-        tree = "results/tree-raw_{serotype}.nwk",
-        alignment = "results/aligned_{serotype}.fasta",
+        tree = "results/tree-raw_{serotype}_{gene}.nwk",
+        alignment = "results/aligned_{serotype}_{gene}.fasta",
         metadata = "data/metadata_{serotype}.tsv"
     output:
-        tree = "results/tree_{serotype}.nwk",
-        node_data = "results/branch-lengths_{serotype}.json",
+        tree = "results/tree_{serotype}_{gene}.nwk",
+        node_data = "results/branch-lengths_{serotype}_{gene}.json",
     params:
         coalescent = "const",
         date_inference = "marginal",

phylogenetic/rules/export.smk

@@ -40,7 +40,7 @@ rule colors:
 rule prepare_auspice_config:
     """Prepare the auspice config file for each serotypes"""
     output:
-        auspice_config="results/config/auspice_config_{serotype}.json",
+        auspice_config="results/config/auspice_config_{serotype}_{gene}.json",
     params:
         replace_clade_key="clade_membership",
         replace_clade_title=lambda wildcard: r"Serotype" if wildcard.serotype in ['all'] else r"DENV genotype",
@@ -109,17 +109,17 @@ rule prepare_auspice_config:
 rule export:
     """Exporting data files for auspice"""
     input:
-        tree = "results/tree_{serotype}.nwk",
+        tree = "results/tree_{serotype}_{gene}.nwk",
         metadata = "data/metadata_{serotype}.tsv",
-        branch_lengths = "results/branch-lengths_{serotype}.json",
-        traits = "results/traits_{serotype}.json",
-        clades = "results/clades_{serotype}.json",
-        nt_muts = "results/nt-muts_{serotype}.json",
-        aa_muts = "results/aa-muts_{serotype}.json",
-        auspice_config = "results/config/auspice_config_{serotype}.json",
+        branch_lengths = "results/branch-lengths_{serotype}_{gene}.json",
+        traits = "results/traits_{serotype}_{gene}.json",
+        clades = lambda wildcard: "results/clades_{serotype}_{gene}.json" if wildcard.gene in ['genome'] else [],
+        nt_muts = "results/nt-muts_{serotype}_{gene}.json",
+        aa_muts = "results/aa-muts_{serotype}_{gene}.json",
+        auspice_config = "results/config/auspice_config_{serotype}_{gene}.json",
         colors = "results/colors_{serotype}.tsv",
     output:
-        auspice_json = "results/raw_dengue_{serotype}.json"
+        auspice_json = "results/raw_dengue_{serotype}_{gene}.json",
     params:
         strain_id = config.get("strain_id_field", "strain"),
     shell:
@@ -137,10 +137,10 @@ rule export:
 
 rule final_strain_name:
     input:
-        auspice_json="results/raw_dengue_{serotype}.json",
+        auspice_json="results/raw_dengue_{serotype}_{gene}.json",
         metadata="data/metadata_{serotype}.tsv",
     output:
-        auspice_json="auspice/dengue_{serotype}_genome.json"
+        auspice_json="auspice/dengue_{serotype}_{gene}.json"
     params:
         strain_id=config.get("strain_id_field", "strain"),
         display_strain_field=config.get("display_strain_field", "strain"),