From 3304ffef9689eb09c8450cf61eedfcabdaf02028 Mon Sep 17 00:00:00 2001 From: Trevor Bedford Date: Sun, 2 Aug 2026 11:12:45 +0200 Subject: [PATCH 1/2] Add cumulative root-to-node S1 mutation count S1_mutations uses augur distance --compare-to root, which counts S1 sites (spike aa 14-685) that currently differ from the root. As S1 sites now mutate recurrently, a second hit at an already-changed site adds nothing and a reversion subtracts, so the metric saturates and looks like it is slowing down. Add scripts/root_to_node_distance.py, which sums the per-branch S1 distance along the path from the root to each node (reusing augur's read_distance_map, get_distance_between_nodes and load_alignments), so every substitution event accrues -- a molecular-clock-like tally rather than a net difference. Gaps and ambiguities are excluded via a new ignored_characters key in S1.json so artifactual gaps in submitted sequences cannot inflate divergence. Wire this in as a new rule cumulative_distances producing the attribute S1_cumulative_mutations, alongside the existing vs-root S1_mutations, with a new "S1 mutations (cumulative)" coloring, so the two can be compared head-to-head before S1_mutations itself is switched over. Co-Authored-By: Claude Opus 4.8 --- defaults/auspice_config.json | 5 ++ defaults/distance_maps/S1.json | 1 + scripts/root_to_node_distance.py | 63 ++++++++++++++++++++++ workflow/snakemake_rules/main_workflow.smk | 26 +++++++++ 4 files changed, 95 insertions(+) create mode 100644 scripts/root_to_node_distance.py diff --git a/defaults/auspice_config.json b/defaults/auspice_config.json index 1bd087018..80cdbc545 100644 --- a/defaults/auspice_config.json +++ b/defaults/auspice_config.json @@ -20,6 +20,11 @@ "title": "S1 mutations", "type": "continuous" }, + { + "key": "S1_cumulative_mutations", + "title": "S1 mutations (cumulative)", + "type": "continuous" + }, { "key": "mlr_lineage_fitness", "title": "MLR lineage fitness", diff --git a/defaults/distance_maps/S1.json b/defaults/distance_maps/S1.json index be7d59bce..45c9bb4fa 100644 --- a/defaults/distance_maps/S1.json +++ b/defaults/distance_maps/S1.json @@ -676,5 +676,6 @@ "685": 1 } }, + "ignored_characters": ["-", "X", "N"], "name": "S1" } diff --git a/scripts/root_to_node_distance.py b/scripts/root_to_node_distance.py new file mode 100644 index 000000000..d9ca5e1d8 --- /dev/null +++ b/scripts/root_to_node_distance.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +""" +Sum per-branch sequence distances along the path from the root to each node. + +`augur distance --compare-to root` compares every node directly against the root +sequence, so recurrent and reverting substitutions at the same site do not +accumulate -- the metric saturates as sites are hit repeatedly. This script +instead walks the tree and, for every edge, adds the per-branch distance +computed with augur's own distance-map logic, so each substitution event along +the root-to-node path is counted. The result is a molecular-clock-like tally of +cumulative amino-acid changes rather than a net difference from the reference. +""" +import argparse +from collections import defaultdict + +from Bio import Phylo +from augur.distance import read_distance_map, get_distance_between_nodes +from augur.reconstruct_sequences import load_alignments +from augur.io import write_json + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Sum per-branch sequence distances from the root to each node.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument("--tree", required=True, help="Newick tree with named internal nodes") + parser.add_argument("--alignment", nargs="+", required=True, + help="per-gene FASTA alignment(s) including internal-node sequences") + parser.add_argument("--gene-names", nargs="+", required=True, + help="gene name for each alignment, paired positionally with --alignment") + parser.add_argument("--map", required=True, + help="distance map JSON (sites, weights, ignored_characters)") + parser.add_argument("--attribute-name", required=True, + help="name to store the cumulative distance under in the node-data JSON") + parser.add_argument("--output", required=True, help="output node-data JSON") + args = parser.parse_args() + + tree = Phylo.read(args.tree, "newick") + distance_map = read_distance_map(args.map) + + # Flatten per-gene alignments to {node_name: {gene: sequence}}, exactly as augur distance does. + alignments = load_alignments(args.alignment, args.gene_names) + sequences_by_node_and_gene = defaultdict(dict) + for gene, alignment in alignments.items(): + for record in alignment: + sequences_by_node_and_gene[record.name][gene] = str(record.seq) + + # Walk from the root, accumulating the per-branch distance on each edge. + # Iterative preorder guarantees a parent is assigned before its children and, + # unlike recursion, is safe for the deep trees ncov can produce. + cumulative = {tree.root.name: 0} + node_data = {tree.root.name: {args.attribute_name: 0}} + for node in tree.find_clades(order="preorder"): + for child in node.clades: + branch_distance = get_distance_between_nodes( + sequences_by_node_and_gene[node.name], + sequences_by_node_and_gene[child.name], + distance_map, + ) + cumulative[child.name] = cumulative[node.name] + branch_distance + node_data[child.name] = {args.attribute_name: cumulative[child.name]} + + write_json({"nodes": node_data}, args.output) diff --git a/workflow/snakemake_rules/main_workflow.smk b/workflow/snakemake_rules/main_workflow.smk index 56947163d..7941c3482 100644 --- a/workflow/snakemake_rules/main_workflow.smk +++ b/workflow/snakemake_rules/main_workflow.smk @@ -909,6 +909,31 @@ rule distances: --output {output} """ +rule cumulative_distances: + input: + tree = rules.refine.output.tree, + alignments = "results/{build_name}/translations/aligned.gene.S_withInternalNodes.fasta", + distance_maps = ["defaults/distance_maps/S1.json"] + params: + genes = 'S', + attribute_names = 'S1_cumulative_mutations' + output: + node_data = "results/{build_name}/cumulative_distances.json" + log: + "logs/cumulative_distances_{build_name}.txt" + conda: + config["conda_environment"] + shell: + r""" + python3 scripts/root_to_node_distance.py \ + --tree {input.tree} \ + --alignment {input.alignments} \ + --gene-names {params.genes} \ + --map {input.distance_maps} \ + --attribute-name {params.attribute_names} \ + --output {output.node_data} 2>&1 | tee {log} + """ + rule traits: message: """ @@ -1197,6 +1222,7 @@ def _get_node_data_by_wildcards(wildcards): rules.traits.output.node_data, rules.mlr_lineage_fitness.output.node_data, rules.distances.output.node_data, + rules.cumulative_distances.output.node_data, rules.calculate_epiweeks.output.node_data, ] From eef6e886337d184623223237e474fbaa2a2ed99c Mon Sep 17 00:00:00 2001 From: Trevor Bedford Date: Sun, 2 Aug 2026 11:54:54 +0200 Subject: [PATCH 2/2] Make S1_mutations the cumulative root-to-node count Now that the cumulative metric can be compared head-to-head against the old net-vs-root count, promote it to be the canonical S1_mutations coloring. Remove the vs-root `rule distances` (augur distance --compare-to root); the existing `rule cumulative_distances` now emits the attribute S1_mutations directly. Drop the temporary "S1 mutations (cumulative)" coloring and the stale rules.distances node-data wiring in export. Co-Authored-By: Claude Opus 4.8 --- defaults/auspice_config.json | 5 ---- workflow/snakemake_rules/main_workflow.smk | 28 +--------------------- 2 files changed, 1 insertion(+), 32 deletions(-) diff --git a/defaults/auspice_config.json b/defaults/auspice_config.json index 80cdbc545..1bd087018 100644 --- a/defaults/auspice_config.json +++ b/defaults/auspice_config.json @@ -20,11 +20,6 @@ "title": "S1 mutations", "type": "continuous" }, - { - "key": "S1_cumulative_mutations", - "title": "S1 mutations (cumulative)", - "type": "continuous" - }, { "key": "mlr_lineage_fitness", "title": "MLR lineage fitness", diff --git a/workflow/snakemake_rules/main_workflow.smk b/workflow/snakemake_rules/main_workflow.smk index 7941c3482..fb7451368 100644 --- a/workflow/snakemake_rules/main_workflow.smk +++ b/workflow/snakemake_rules/main_workflow.smk @@ -884,31 +884,6 @@ rule translate: --output {output.node_data} 2>&1 | tee {log} """ -rule distances: - input: - tree = rules.refine.output.tree, - alignments = "results/{build_name}/translations/aligned.gene.S_withInternalNodes.fasta", - distance_maps = ["defaults/distance_maps/S1.json"] - params: - genes = 'S', - comparisons = ['root'], - attribute_names = ['S1_mutations'] - output: - node_data = "results/{build_name}/distances.json" - conda: - config["conda_environment"] - shell: - r""" - augur distance \ - --tree {input.tree} \ - --alignment {input.alignments} \ - --gene-names {params.genes} \ - --compare-to {params.comparisons} \ - --attribute-name {params.attribute_names} \ - --map {input.distance_maps} \ - --output {output} - """ - rule cumulative_distances: input: tree = rules.refine.output.tree, @@ -916,7 +891,7 @@ rule cumulative_distances: distance_maps = ["defaults/distance_maps/S1.json"] params: genes = 'S', - attribute_names = 'S1_cumulative_mutations' + attribute_names = 'S1_mutations' output: node_data = "results/{build_name}/cumulative_distances.json" log: @@ -1221,7 +1196,6 @@ def _get_node_data_by_wildcards(wildcards): rules.recency.output.node_data, rules.traits.output.node_data, rules.mlr_lineage_fitness.output.node_data, - rules.distances.output.node_data, rules.cumulative_distances.output.node_data, rules.calculate_epiweeks.output.node_data, ]