From 395966b467f229c3a77964d0f828a100b44af88a Mon Sep 17 00:00:00 2001
From: Victor Lin <vlin@fredhutch.org>
Date: Fri, 17 Jul 2026 13:05:00 -0700
Subject: [PATCH] Use augur subsample

Augur subsample supports more complex subsampling patterns that are not
supported by a single call to augur filter.

This is a breaking change and the old configuration will no longer work.
---
 CHANGELOG.md                               |  1 +
 phylogenetic/defaults/config.yaml          | 49 +++++++++++++++-------
 phylogenetic/rules/annotate_phylogeny.smk  |  4 +-
 phylogenetic/rules/config.smk              | 19 ++++++++-
 phylogenetic/rules/construct_phylogeny.smk |  4 +-
 phylogenetic/rules/export.smk              |  6 +--
 phylogenetic/rules/prepare_sequences.smk   | 29 ++++++-------
 7 files changed, 72 insertions(+), 40 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 06a75d7..6ee02de 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ the outputs of the workflows.
 
 ## 2026
 
+* 20 July 2026: `filter` is replaced with `subsample`. See updated default config for the expected schema. **This is a breaking change**.
 * 20 July 2026: Command-line parameters for `filter`, `refine`, `ancestral`, `traits`, and `tip_frequencies` are now defined as a single argument string under `<rule>.<build>` in config. This reverts the change from 14 July 2026 and provides a consistent and flexible pattern for passing **non-file** arguments to Augur **This is a breaking change**.
 * 14 July 2026: The refine option for clock filtering has moved from `refine.<build>: "--clock-filter-iqd <N>"` to `refine.<build>.clock_filter_iqd: <N>`. **This is a breaking change**.
 * 14 July 2026: The following refine options can be configured using `refine.<build>.<option>`: `coalescent`, `date_inference`, `timetree`, `date_confidence`, `clock_filter_iqd`, `divergence_units`.
diff --git a/phylogenetic/defaults/config.yaml b/phylogenetic/defaults/config.yaml
index 1a25445..eb155eb 100644
--- a/phylogenetic/defaults/config.yaml
+++ b/phylogenetic/defaults/config.yaml
@@ -28,21 +28,40 @@ reference: "reference.gb"
 # argument strings, but instead using their dedicated YAML config parameters.
 # They can be customized per build by referencing '{build}' in the path.
 
-# https://docs.nextstrain.org/projects/augur/en/stable/usage/cli/filter.html
-filter:
-  exclude: "{build}/exclude.txt"
-  include: "{build}/include.txt"
-  north-america: >-
-    --min-length 8000
-    --group-by country year month MuV_genotype division
-    --subsample-max-sequences 4000
-    --min-date 2006
-    --query "region=='North America' & (MuV_genotype=='G')"
-  global: >-
-    --min-length 8000
-    --group-by country year month MuV_genotype division
-    --subsample-max-sequences 4000
-    --min-date 1950
+# Subsample is different. Instead of an argument string, it takes a YAML
+# configuration which must follow a specific schema.
+#
+# https://docs.nextstrain.org/projects/augur/en/stable/usage/cli/subsample.html#configuration
+subsample:
+  north-america:
+    samples:
+      sample:
+        exclude: "north-america/exclude.txt"
+        include: "north-america/include.txt"
+        min_length: 8000
+        group_by:
+        - country
+        - year
+        - month
+        - MuV_genotype
+        - division
+        max_sequences: 4000
+        min_date: 2006
+        query: "region=='North America' & (MuV_genotype=='G')"
+  global:
+    samples:
+      sample:
+        exclude: "global/exclude.txt"
+        include: "global/include.txt"
+        min_length: 8000
+        group_by:
+        - country
+        - year
+        - month
+        - MuV_genotype
+        - division
+        max_sequences: 4000
+        min_date: 1950
 
 # https://docs.nextstrain.org/projects/augur/en/stable/usage/cli/refine.html
 refine:
diff --git a/phylogenetic/rules/annotate_phylogeny.smk b/phylogenetic/rules/annotate_phylogeny.smk
index d56ada2..a16c05f 100644
--- a/phylogenetic/rules/annotate_phylogeny.smk
+++ b/phylogenetic/rules/annotate_phylogeny.smk
@@ -3,7 +3,7 @@ This part of the workflow creates additonal annotations for the phylogenetic tre
 
 REQUIRED INPUTS:
 
-    metadata   = results/{build}/filtered.tsv
+    metadata   = results/{build}/subsampled.tsv
     alignment  = results/{build}/aligned.fasta
     tree       = results/{build}/tree.nwk
 
@@ -91,7 +91,7 @@ rule traits:
     """
     input:
         tree = "results/{build}/tree.nwk",
-        metadata = "results/{build}/filtered.tsv",
+        metadata = "results/{build}/subsampled.tsv",
     output:
         node_data = "results/{build}/traits.json",
     log:
diff --git a/phylogenetic/rules/config.smk b/phylogenetic/rules/config.smk
index 08d8483..799e638 100644
--- a/phylogenetic/rules/config.smk
+++ b/phylogenetic/rules/config.smk
@@ -5,6 +5,10 @@ OUTPUTS:
 
     results/run_config.yaml
 """
+def main():
+    write_subsample_config()
+    write_config("results/run_config.yaml")
+
 
 def conditional(option, argument):
     """Used for config-defined arguments whose presence necessitates a command-line option
@@ -25,4 +29,17 @@ def conditional(option, argument):
     raise WorkflowError(f"Workflow function conditional() received an argument value of unexpected type: {type(argument).__name__}")
 
 
-write_config("results/run_config.yaml")
+def write_subsample_config():
+    for build in config["builds"]:
+        if "custom_subsample" in config:
+            section = ["custom_subsample", build]
+        else:
+            section = ["subsample", build]
+        write_config(f"results/{build}/subsample_config.yaml", section=section)
+
+
+try:
+    main()
+except InvalidConfigError as e:
+    print(f"ERROR: {e}", file=sys.stderr)
+    exit(1)
diff --git a/phylogenetic/rules/construct_phylogeny.smk b/phylogenetic/rules/construct_phylogeny.smk
index fcbe9ac..972b307 100644
--- a/phylogenetic/rules/construct_phylogeny.smk
+++ b/phylogenetic/rules/construct_phylogeny.smk
@@ -3,7 +3,7 @@ This part of the workflow constructs the phylogenetic tree.
 
 REQUIRED INPUTS:
 
-    metadata   = results/{build}/filtered.tsv
+    metadata   = results/{build}/subsampled.tsv
     alignment  = results/{build}/aligned.fasta
 
 OUTPUTS:
@@ -45,7 +45,7 @@ rule refine:
     input:
         tree = "results/{build}/tree_raw.nwk",
         alignment = "results/{build}/aligned.fasta",
-        metadata = "results/{build}/filtered.tsv"
+        metadata = "results/{build}/subsampled.tsv"
     output:
         tree = "results/{build}/tree.nwk",
         node_data = "results/{build}/branch_lengths.json",
diff --git a/phylogenetic/rules/export.smk b/phylogenetic/rules/export.smk
index b62af25..9af6b1e 100644
--- a/phylogenetic/rules/export.smk
+++ b/phylogenetic/rules/export.smk
@@ -4,7 +4,7 @@ export a Nextstrain dataset.
 
 REQUIRED INPUTS:
 
-    metadata        = results/{build}/filtered.tsv
+    metadata        = results/{build}/subsampled.tsv
     tree            = results/{build}/tree.nwk
     branch_lengths  = results/{build}/branch_lengths.json
     node_data       = results/{build}/*.json
@@ -53,7 +53,7 @@ rule export:
     """Exporting data files for for auspice"""
     input:
         tree = "results/{build}/tree.nwk",
-        metadata = "results/{build}/filtered.tsv",
+        metadata = "results/{build}/subsampled.tsv",
         branch_lengths = "results/{build}/branch_lengths.json",
         traits = "results/{build}/traits.json",
         nt_muts = "results/{build}/nt_muts.json",
@@ -93,7 +93,7 @@ rule tip_frequencies:
     """
     input:
         tree = "results/{build}/tree.nwk",
-        metadata = "results/{build}/filtered.tsv",
+        metadata = "results/{build}/subsampled.tsv",
     output:
         tip_freq = "auspice/mumps_{build}_tip-frequencies.json"
     log:
diff --git a/phylogenetic/rules/prepare_sequences.smk b/phylogenetic/rules/prepare_sequences.smk
index c97ace2..6e7dc53 100644
--- a/phylogenetic/rules/prepare_sequences.smk
+++ b/phylogenetic/rules/prepare_sequences.smk
@@ -9,7 +9,7 @@ REQUIRED INPUTS:
 
 OUTPUTS:
 
-    metadata  = results/{build}/filtered.tsv
+    metadata  = results/{build}/subsampled.tsv
     alignment = results/{build}/aligned.fasta
 
 This part of the workflow usually includes the following steps:
@@ -21,39 +21,34 @@ This part of the workflow usually includes the following steps:
 
 See Augur's usage docs for these commands for more details.
 """
+from augur.subsample import get_referenced_files
 
-rule filter:
-    """
-    Filtering sequences
-    """
+rule subsample:
     input:
+        config = "results/{build}/subsample_config.yaml",
         sequences = "results/sequences.fasta",
         metadata = "results/metadata.tsv",
-        exclude = resolve_config_path(config["filter"]["exclude"]),
-        include = resolve_config_path(config["filter"]["include"]),
+        referenced_files = lambda w: get_referenced_files(f"results/{w.build}/subsample_config.yaml"),
     output:
-        sequences = "results/{build}/filtered.fasta",
-        metadata = "results/{build}/filtered.tsv",
+        sequences = "results/{build}/subsampled.fasta",
+        metadata = "results/{build}/subsampled.tsv",
     log:
-        "logs/{build}/filtered.txt",
+        "logs/{build}/subsample.txt",
     benchmark:
-        "benchmarks/{build}/filtered.txt",
+        "benchmarks/{build}/subsample.txt",
     params:
-        args = lambda w: config['filter'][w.build],
         strain_id = config.get("strain_id_field", "strain"),
     shell:
         r"""
         exec &> >(tee {log:q})
 
-        augur filter \
+        augur subsample \
             --sequences {input.sequences:q} \
             --metadata {input.metadata:q} \
             --metadata-id-columns {params.strain_id:q} \
-            --exclude {input.exclude:q} \
-            --include {input.include:q} \
             --output-sequences {output.sequences:q} \
             --output-metadata {output.metadata:q} \
-            {params.args}
+            --config {input.config}
         """
 
 rule align:
@@ -62,7 +57,7 @@ rule align:
       - filling gaps with N
     """
     input:
-        sequences = "results/{build}/filtered.fasta",
+        sequences = "results/{build}/subsampled.fasta",
         reference = resolve_config_path(config['reference']),
     output:
         alignment = "results/{build}/aligned.fasta",