Update Nextclade workflow to handle new clade Ib outbreak sh2023 lineages

nextclade/Snakefile

@@ -106,23 +106,50 @@ rule decompress_metadata:
         """
 
 
+rule filter_known_duplicates:
+    input:
+        known_duplicates="resources/known_duplicates.txt",
+        includes=[
+            "resources/include_accessions.txt",
+            expand(
+                "resources/{build_name}/include_accessions.txt", build_name=BUILD_NAMES
+            ),
+            "resources/all-clades/include_recombinants.txt",
+        ],
+    output:
+        "results/known_duplicates_to_exclude.txt",
+    run:
+        known_duplicates = {
+            line.strip() for line in open(input.known_duplicates) if line.strip()
+        }
+        includes = set()
+        for path in input.includes:
+            for line in open(path):
+                accession = line.split("#", maxsplit=1)[0].strip()
+                if accession:
+                    includes.add(accession)
+        with open(output[0], "w") as file:
+            for accession in sorted(known_duplicates - includes):
+                print(accession, file=file)
+
+
 rule decompress_sequences_and_remove_duplicates:
     "Decompressing sequences and strip version from fasta ID"
     input:
         sequences="data/{file}.fasta.zst",
-        known_duplicates="resources/known_duplicates.txt",
+        known_duplicates="results/known_duplicates_to_exclude.txt",
     output:
         sequences="data/{file}.fasta",
     params:
         min_length=config["min_length"],
     shell:
         """
         zstdcat {input.sequences} \
-        | seqkit replace \
-            -p "^([^.]+).*" -r '$1' \
-        | seqkit grep -w0 -v -f {input.known_duplicates} \
-        | seqkit seq -m={params.min_length} \
-        > {output.sequences}
+            | seqkit replace \
+                -p "^([^.]+).*" -r '$1' \
+            | seqkit grep -w0 -v -f {input.known_duplicates} \
+            | seqkit seq -m={params.min_length} \
+                >{output.sequences}
         """
 
 
@@ -156,18 +183,30 @@ rule premask:
 
 
 rule deduplicate:
-    """
-    Remove identical sequences (even if they have differing Ns)
-    Keep those sequences with fewer Ns
-    Focus for Nextclade is on diversity, not on representativeness
-    """
+    # Remove identical sequences (even if they have differing Ns)
+    # Keep those sequences with fewer Ns
+    # Focus for Nextclade is on diversity, not on representativeness
     input:
         sequences="results/premasked.fasta",
+        includes=[
+            "resources/include_accessions.txt",
+            expand(
+                "resources/{build_name}/include_accessions.txt", build_name=BUILD_NAMES
+            ),
+            "resources/all-clades/include_recombinants.txt",
+        ],
     output:
         "results/duplicates.txt",
+    params:
+        include_options=lambda w, input: " ".join(
+            f"--include {path}" for path in input.includes
+        ),
     shell:
         """
-        python3 scripts/deduplicate.py {input.sequences} {output}
+        python3 scripts/deduplicate.py \
+            {input.sequences} \
+            {output} \
+            {params.include_options}
         """
 
 
@@ -228,7 +267,14 @@ rule reformat_ambiguous:
 
 
         df["date"] = df["date"].apply(replace_ambiguous_date)
-
+        df["outbreakLineage"] = df.apply(
+            lambda row: (
+                f"{row['outbreak']}/{row['lineage']}"
+                if row["outbreak"] and row["lineage"]
+                else ""
+            ),
+            axis=1,
+        )
         df.to_csv(output.metadata, sep="\t", index=False)
 
 
@@ -361,21 +407,21 @@ rule align:
     shell:
         """
         cat {input.root} {input.reference} {input.sequences} \
-        | seqkit seq -w0 -i \
-        | seqkit rmdup -n \
-        | nextclade3 run \
-            --retry-reverse-complement \
-            --input-ref {input.reference} \
-            --excess-bandwidth 100 \
-            --terminal-bandwidth 300 \
-            --allowed-mismatches 8 \
-            --window-size 40 \
-            --min-seed-cover 0.1 \
-            --input-annotation {input.annotation} \
-            --gap-alignment-side left \
-            --output-fasta /dev/stdout \
-            --output-translations {params.translation_template} \
-        | seqkit seq -w0 -i > {output.alignment}
+            | seqkit seq -w0 -i \
+            | seqkit rmdup -n \
+            | nextclade3 run \
+                --retry-reverse-complement \
+                --input-ref {input.reference} \
+                --excess-bandwidth 100 \
+                --terminal-bandwidth 300 \
+                --allowed-mismatches 8 \
+                --window-size 40 \
+                --min-seed-cover 0.1 \
+                --input-annotation {input.annotation} \
+                --gap-alignment-side left \
+                --output-fasta /dev/stdout \
+                --output-translations {params.translation_template} \
+            | seqkit seq -w0 -i >{output.alignment}
         """
 
 
@@ -395,18 +441,25 @@ rule mask:
 
 
 rule deduplicate_2:
-    """
-    Remove identical sequences (even if they have differing Ns)
-    Keep those sequences with fewer Ns
-    Focus for Nextclade is on diversity, not on representativeness
-    """
+    # Remove identical sequences (even if they have differing Ns)
+    # Keep those sequences with fewer Ns
+    # Focus for Nextclade is on diversity, not on representativeness
     input:
         sequences="results/{build_name}/masked_with_dups.fasta",
+        include=lambda w: config[w.build_name]["general_include"],
+        specific_include="resources/{build_name}/include_accessions.txt",
     output:
         duplicates="results/{build_name}/duplicates.txt",
+    params:
+        include_options=lambda w, input: " ".join(
+            f"--include {path}" for path in [input.include, input.specific_include]
+        ),
     shell:
         """
-        python3 scripts/deduplicate.py {input.sequences} {output.duplicates}
+        python3 scripts/deduplicate.py \
+            {input.sequences} \
+            {output.duplicates} \
+            {params.include_options}
         """
 
 
@@ -422,7 +475,7 @@ rule pre_tree_exclude:
         seqkit grep -w0 -v -f \
             <(cat {input.duplicates} {input.recombinants}) \
             {input.sequences} \
-        > {output.sequences}
+            >{output.sequences}
         """
 
 
@@ -451,23 +504,22 @@ else:
     rule tree:
         input:
             alignment="results/{build_name}/masked.fasta",
+            constraint_tree=lambda w: config[w.build_name].get("constraint_tree", ""),
         output:
             tree="results/{build_name}/masked.fasta.treefile",
+        params:
+            constraint_option=lambda w: (
+                f"-g {config[w.build_name]['constraint_tree']}"
+                if "constraint_tree" in config[w.build_name]
+                else ""
+            ),
         shell:
-            # """
-            # augur tree \
-            #     --alignment {input.alignment} \
-            #     --tree-builder-args "--polytomy --ninit 2 -n 2 --epsilon 0.05 -T 4 --redo" \
-            #     --nthreads 4 \
-            #     --output {output.tree}
-            # """
             """
-            ~/code/pree/rust/target/release/rust_parsimony build \
-                {input.alignment} \
-                --output-file {output.tree} \
-                --pruning-threshold 30 \
-                --num-trees 20 \
-                --spr-trees 10
+            augur tree \
+                --alignment {input.alignment} \
+                --tree-builder-args "--polytomy --ninit 2 -n 2 --epsilon 0.05 -T 4 --redo {params.constraint_option}" \
+                --nthreads 4 \
+                --output {output.tree}
             """
 
 
@@ -565,7 +617,7 @@ rule ancestral:
             --genes {params.genes} \
             --output-node-data {output.node_data} \
             --inference joint \
-            2>&1 > {output.log}
+            2>&1 >{output.log}
         """
 
 
@@ -590,14 +642,21 @@ rule clades:
 rule rename_clades:
     input:
         clades="results/{build_name}/clades_raw.json",
+        tree="results/{build_name}/tree.nwk",
         script="scripts/clades_renaming.py",
     output:
         node_data="results/{build_name}/clades.json",
+    params:
+        outgroup_clade_name=lambda w: config[w.build_name].get(
+            "outgroup_clade_name", "unassigned"
+        ),
     shell:
         """
         python {input.script} \
-        --input-node-data {input.clades} \
-        --output-node-data {output.node_data}
+            --input-node-data {input.clades} \
+            --output-node-data {output.node_data} \
+            --outgroup-clade-name {params.outgroup_clade_name} \
+            --tree-file {input.tree}
         """
 
 
@@ -642,7 +701,7 @@ rule patch_auspice_config:
         auspice_config="results/{build_name}/auspice_config.json",
     shell:
         """
-        jsonpatch {input.auspice_config} {input.patch_json} > {output.auspice_config}
+        jsonpatch {input.auspice_config} {input.patch_json} >{output.auspice_config}
         """
 
 

nextclade/config/config.yaml

@@ -4,18 +4,25 @@ clade-iib:
   exclude_where:
     b1: "--exclude-where clade!=IIb"
     nonb1: >-
-      --exclude-where clade!=IIb lineage=B.1 lineage=B.1.1 lineage=B.1.2 lineage=B.1.3 lineage=B.1.4 lineage=B.1.5 lineage=B.1.6
-      lineage=B.1.7 lineage=B.1.8 lineage=B.1.9 lineage=B.1.10 lineage=B.1.11 lineage=B.1.12 lineage=B.1.13 lineage=B.1.14
-      lineage=B.1.15
-      lineage=B.1.16 lineage=B.1.17 lineage=B.1.18 lineage=B.1.19 lineage=B.1.20 lineage=B.1.21 lineage=B.1.22 lineage=B.1.23
-      lineage=C.1 lineage=C.2 lineage=C.3 lineage=C.4 lineage=C.5 lineage=C.6 lineage=D.1 lineage=D.2 lineage=D.3 lineage=D.4
-      lineage=D.5 lineage=D.6 lineage=E.1 lineage=E.2 lineage=E.3 lineage=E.4 lineage=E.5 lineage=E.6 lineage=F.1 lineage=F.2
-      lineage=F.3 lineage=F.4 lineage=F.5 lineage=F.6 lineage=C.1.2 lineage=C.1.3 lineage=E.4 lineage=E.1.1 lineage=E.2.1
-      lineage=E.3.1 lineage=F.2.1 lineage=F.4.1 lineage=J.1
+      --exclude-where clade!=IIb outbreakLineage=sh2017/B.1 outbreakLineage=sh2017/B.1.1 outbreakLineage=sh2017/B.1.2 outbreakLineage=sh2017/B.1.3
+      outbreakLineage=sh2017/B.1.4 outbreakLineage=sh2017/B.1.5 outbreakLineage=sh2017/B.1.6
+      outbreakLineage=sh2017/B.1.7 outbreakLineage=sh2017/B.1.8 outbreakLineage=sh2017/B.1.9 outbreakLineage=sh2017/B.1.10
+      outbreakLineage=sh2017/B.1.11 outbreakLineage=sh2017/B.1.12 outbreakLineage=sh2017/B.1.13 outbreakLineage=sh2017/B.1.14
+      outbreakLineage=sh2017/B.1.15
+      outbreakLineage=sh2017/B.1.16 outbreakLineage=sh2017/B.1.17 outbreakLineage=sh2017/B.1.18 outbreakLineage=sh2017/B.1.19
+      outbreakLineage=sh2017/B.1.20 outbreakLineage=sh2017/B.1.21 outbreakLineage=sh2017/B.1.22 outbreakLineage=sh2017/B.1.23
+      outbreakLineage=sh2017/C.1 outbreakLineage=sh2017/C.2 outbreakLineage=sh2017/C.3 outbreakLineage=sh2017/C.4 outbreakLineage=sh2017/C.5
+      outbreakLineage=sh2017/C.6 outbreakLineage=sh2017/D.1 outbreakLineage=sh2017/D.2 outbreakLineage=sh2017/D.3 outbreakLineage=sh2017/D.4
+      outbreakLineage=sh2017/D.5 outbreakLineage=sh2017/D.6 outbreakLineage=sh2017/E.1 outbreakLineage=sh2017/E.2 outbreakLineage=sh2017/E.3
+      outbreakLineage=sh2017/E.4 outbreakLineage=sh2017/E.5 outbreakLineage=sh2017/E.6 outbreakLineage=sh2017/F.1 outbreakLineage=sh2017/F.2
+      outbreakLineage=sh2017/F.3 outbreakLineage=sh2017/F.4 outbreakLineage=sh2017/F.5 outbreakLineage=sh2017/F.6 outbreakLineage=sh2017/C.1.2
+      outbreakLineage=sh2017/C.1.3 outbreakLineage=sh2017/E.4 outbreakLineage=sh2017/E.1.1 outbreakLineage=sh2017/E.2.1
+      outbreakLineage=sh2017/E.3.1 outbreakLineage=sh2017/F.2.1 outbreakLineage=sh2017/F.4.1 outbreakLineage=sh2017/J.1
   subsampling:
-    b1: "--subsample-max-sequences 1000 --group-by clade lineage year month country --probabilistic-sampling"
-    nonb1: "--sequences-per-group 10 --group-by clade lineage year month country --probabilistic-sampling"
+    b1: "--subsample-max-sequences 1000 --group-by clade outbreakLineage year month country --probabilistic-sampling"
+    nonb1: "--sequences-per-group 10 --group-by clade outbreakLineage year month country --probabilistic-sampling"
   general_include: "resources/include_accessions.txt"
+  outgroup_clade_name: "outgroup"
 
 clade-i:
   root: "reconstructed_ancestral_mpox"
@@ -24,8 +31,8 @@ clade-i:
     b1: "--exclude-where clade=II clade=IIa clade=IIb clade=outgroup clade=Ia"
     nonb1: "--exclude-where clade=II clade=IIa clade=IIb clade=outgroup clade=Ib"
   subsampling:
-    b1: "--subsample-max-sequences 1000 --group-by clade lineage year month country --probabilistic-sampling"
-    nonb1: "--subsample-max-sequences 1000 --group-by clade lineage year month country --probabilistic-sampling"
+    b1: "--subsample-max-sequences 1000 --group-by clade outbreakLineage year month country --probabilistic-sampling"
+    nonb1: "--subsample-max-sequences 1000 --group-by clade outbreakLineage year month country --probabilistic-sampling"
   annotation: "resources/clade-i/genome_annotation.gff3"
   genes: "resources/clade-i/genes.txt"
   general_include: "resources/clade-i/include_accessions.txt"
@@ -36,37 +43,45 @@ clade-i:
     - OPG005
     - OPG015_dup
     - OPG016
+  outgroup_clade_name: "outgroup"
+  constraint_tree: "resources/clade-i/constraint_tree.nwk"
 
 all-clades:
   root: "reconstructed_ancestral_mpox"
   min_date: ""
   exclude_where:
-    b1: "--exclude-where outbreak=hMPXV-1"
-    nonb1: "--exclude-where outbreak!=hMPXV-1"
+    b1: "--exclude-where outbreak=sh2017"
+    nonb1: "--exclude-where outbreak!=sh2017"
   subsampling:
-    b1: "--subsample-max-sequences 500 --group-by clade lineage year country --probabilistic-sampling"
-    nonb1: "--subsample-max-sequences 800 --group-by clade lineage year country --probabilistic-sampling"
+    b1: "--subsample-max-sequences 600 --group-by clade outbreakLineage year country --probabilistic-sampling"
+    nonb1: "--subsample-max-sequences 1000 --group-by clade outbreakLineage year country --probabilistic-sampling"
   general_include: "resources/include_accessions.txt"
   genes_to_exclude:
     - OPG001_dup
     - OPG002_dup
     - OPG003_dup
     - OPG015_dup
+  outgroup_clade_name: "unassigned"
 
 lineage-b.1:
   root: "ref_in_coord"
   min_date: ""
   exclude_where:
     b1: >-
-      --min-date 2022 --exclude-where outbreak!=hMPXV-1 lineage=A.2 lineage=A.2.1 lineage=A.2.2 lineage=A.2.3 lineage=A.2.4
-      lineage=A.2.5 lineage=G.1 lineage=H.1 lineage=H.2 lineage=A lineage=A.1
+      --min-date 2022 --exclude-where outbreak!=sh2017 outbreakLineage=sh2017/A.2 outbreakLineage=sh2017/A.2.1 outbreakLineage=sh2017/A.2.2
+      outbreakLineage=sh2017/A.2.3 outbreakLineage=sh2017/A.2.4
+      outbreakLineage=sh2017/A.2.5 outbreakLineage=sh2017/G.1 outbreakLineage=sh2017/H.1 outbreakLineage=sh2017/H.2 outbreakLineage=sh2017/A
+      outbreakLineage=sh2017/A.1
     nonb1: >-
-      --min-date 2024 --exclude-where outbreak!=hMPXV-1 lineage=A.2 lineage=A.2.1 lineage=A.2.2 lineage=A.2.3 lineage=A.2.4
-      lineage=A.2.5 lineage=G.1 lineage=H.1 lineage=H.2 lineage=A lineage=A.1
+      --min-date 2024 --exclude-where outbreak!=sh2017 outbreakLineage=sh2017/A.2 outbreakLineage=sh2017/A.2.1 outbreakLineage=sh2017/A.2.2
+      outbreakLineage=sh2017/A.2.3 outbreakLineage=sh2017/A.2.4
+      outbreakLineage=sh2017/A.2.5 outbreakLineage=sh2017/G.1 outbreakLineage=sh2017/H.1 outbreakLineage=sh2017/H.2 outbreakLineage=sh2017/A
+      outbreakLineage=sh2017/A.1
   subsampling:
-    b1: "--subsample-max-sequences 2000 --group-by clade lineage year month country --probabilistic-sampling"
-    nonb1: "--subsample-max-sequences 2000 --group-by clade lineage year month country --probabilistic-sampling"
+    b1: "--subsample-max-sequences 2000 --group-by clade outbreakLineage year month country --probabilistic-sampling"
+    nonb1: "--subsample-max-sequences 2000 --group-by clade outbreakLineage year month country --probabilistic-sampling"
   general_include: "resources/include_accessions.txt"
+  outgroup_clade_name: "outgroup"
 
 ## filter
 min_length: 160000