From 991f1b02c954a9692af8cb0004798883232cfa97 Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Wed, 5 Aug 2026 17:24:47 -0700 Subject: [PATCH 1/4] Add separate config section for filter_for_pre_subsample_alignment Similar to "Add separate frequencies config" (0b221851), the filter_for_pre_subsample_alignment rule shouldn't rely on config from another rule. --- CHANGELOG.md | 1 + config.schema.yaml | 27 +++++++++++++++++++++++++++ config/configfile.yaml | 19 +++++++++++++++++++ workflow/snakemake_rules/core.smk | 6 +++--- 4 files changed, 50 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 98b7b83..0fa7b8e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,4 +4,5 @@ We use this CHANGELOG to document breaking changes, new features, bug fixes, and ## 2026 +* TBD: Phylogenetic workflow configuration for initial quality filtering has moved to its own section, `filter_for_pre_subsample_alignment`. * 11 August 2026: Phylogenetic workflow configuration is now validated against a strict schema. The workflow will error if your configuration has extraneous entries that were previously ignored. diff --git a/config.schema.yaml b/config.schema.yaml index 774060b..0b1e5b2 100644 --- a/config.schema.yaml +++ b/config.schema.yaml @@ -194,6 +194,33 @@ properties: type: string accession: type: string + filter_for_pre_subsample_alignment: + type: object + additionalProperties: false + properties: + group_by: + type: string + min_length: + <<: *per_build_map + patternProperties: + "^.*$": + type: integer + min_coverage: + <<: *per_build_map + patternProperties: + "^.*$": + type: number + resolutions: + <<: *per_resolution_map + patternProperties: + "^.*$": + type: object + additionalProperties: false + properties: + min_date: + type: string + background_min_date: + type: string f_dms_data: type: string f_dms_antibodies: diff --git a/config/configfile.yaml b/config/configfile.yaml index 02f67c3..48467a2 100644 --- a/config/configfile.yaml +++ b/config/configfile.yaml @@ -107,6 +107,25 @@ nextclade_attributes: reference_name: "hRSV/B/Australia/VIC-RCH056/2019" accession: "EPI_ISL_1653999" +filter_for_pre_subsample_alignment: + min_length: + genome: 10000 + G: 600 + F: 1200 + F-antibody-escape: 1200 + min_coverage: + genome: 0.3 + G: 0.3 + F: 0.3 + F-antibody-escape: 0.75 + resolutions: + all-time: + min_date: "1975-01-01" + 6y: + min_date: 6Y + 3y: + min_date: 3Y + # configuration specific to the F deep mutational scanning antibody escape data f_dms_data: dms-data/all_antibodies.csv f_dms_antibodies: # columns in `f_dms_data` with per-mutation escape diff --git a/workflow/snakemake_rules/core.smk b/workflow/snakemake_rules/core.smk index 96f2fe9..750f117 100644 --- a/workflow/snakemake_rules/core.smk +++ b/workflow/snakemake_rules/core.smk @@ -235,10 +235,10 @@ rule filter_for_pre_subsample_alignment: benchmark: "benchmarks/filter_for_pre_subsample_alignment_{a_or_b}_{build_name}_{resolution}.txt" params: - min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter"]["min_coverage"][w.build_name]}', - min_length=lambda w: config["filter"]["min_length"][w.build_name], + min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter_for_pre_subsample_alignment"]["min_coverage"][w.build_name]}', + min_length=lambda w: config["filter_for_pre_subsample_alignment"]["min_length"][w.build_name], strain_id=config["strain_id_field"], - min_date=lambda w: config["filter"]["resolutions"][w.resolution]["min_date"], + min_date=lambda w: config["filter_for_pre_subsample_alignment"]["resolutions"][w.resolution]["min_date"], shell: r""" exec &> >(tee {log:q}) From ae0fce0c9beadc870c6fe18cd74ce63172f94420 Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Wed, 5 Aug 2026 17:08:53 -0700 Subject: [PATCH 2/4] Use augur subsample The previous subsampling implementation was fixed to a two-sample recent+background split with some hardcoded parameters. Replacing it with augur subsample allows for more flexible configuration. In Snakemake, implementation is mostly copied from pathogen repos that have switched over to augur subsample. One notable difference is that the combine_samples rule must stay to handle output from the enrich_antibody_escape rule. In the config YAML, the subsampling configuration is much more verbose as a byproduct of increased flexibility. It was generated using a script, which I'll add in another commit since it makes additional changes. --- CHANGELOG.md | 3 +- config.schema.yaml | 66 +--- config/configfile.yaml | 572 ++++++++++++++++++++++++++-- workflow/snakemake_rules/config.smk | 16 + workflow/snakemake_rules/core.smk | 105 +---- 5 files changed, 586 insertions(+), 176 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0fa7b8e..0386aec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,5 +4,6 @@ We use this CHANGELOG to document breaking changes, new features, bug fixes, and ## 2026 -* TBD: Phylogenetic workflow configuration for initial quality filtering has moved to its own section, `filter_for_pre_subsample_alignment`. +* TBD: The `filter` section in phylogenetic workflow configuration has been replaced by `subsample`/`custom_subsample` for subsampling, and `filter_for_pre_subsample_alignment` for initial quality filtering. **This is a breaking change**. + * NOTE: The workflow does not yet support proximal samples. * 11 August 2026: Phylogenetic workflow configuration is now validated against a strict schema. The workflow will error if your configuration has extraneous entries that were previously ignored. diff --git a/config.schema.yaml b/config.schema.yaml index 0b1e5b2..94859fe 100644 --- a/config.schema.yaml +++ b/config.schema.yaml @@ -5,6 +5,13 @@ description: >- configuration file. $defs: + per_full_build_map: &per_full_build_map + type: object + additionalProperties: false + propertyNames: + title: Full build name + description: Full build name (e.g. 'a/genome/all-time') + per_subtype_map: &per_subtype_map type: object additionalProperties: false @@ -79,52 +86,19 @@ properties: type: string display_strain_field: type: string - filter: - type: object - additionalProperties: false - properties: - group_by: - type: string - min_coverage: - <<: *per_build_map - patternProperties: - "^.*$": - type: number - min_length: - <<: *per_build_map - patternProperties: - "^.*$": - type: integer - resolutions: - <<: *per_resolution_map - patternProperties: - "^.*$": - type: object - additionalProperties: false - properties: - min_date: - type: string - background_min_date: - type: string - subsample_max_sequences: - <<: *per_build_map - patternProperties: - "^.*$": - type: integer - exclude_where: - type: object - additionalProperties: false - properties: - recent: - type: array - items: - type: string - background: - type: array - items: - type: string - missing_data_threshold: - type: integer + subsample: &subsample_config + <<: *per_full_build_map + description: >- + Subsampling configuration. When using --configfile, it is recommended to + use 'custom_subsample' instead to ignore default subsampling configuration. + patternProperties: + "^.*$": + $ref: "https://nextstrain.org/schemas/augur/subsample-config-unaligned/v1" + custom_subsample: + <<: *subsample_config + description: >- + Custom subsampling configuration. When using --configfile, this is + recommended over 'subsample' to ignore default subsampling configuration. files: type: object additionalProperties: false diff --git a/config/configfile.yaml b/config/configfile.yaml index 48467a2..5939779 100644 --- a/config/configfile.yaml +++ b/config/configfile.yaml @@ -29,41 +29,543 @@ strain_id_field: "accession" display_strain_field: "strain" -filter: - group_by: "year country" - min_coverage: - genome: 0.3 - G: 0.3 - F: 0.3 - F-antibody-escape: 0.75 - - min_length: - genome: 10000 - G: 600 - F: 1200 - F-antibody-escape: 1200 - resolutions: - all-time: - min_date: "1975-01-01" - 6y: - min_date: 6Y - background_min_date: 12Y - 3y: - min_date: 3Y - background_min_date: 12Y - - subsample_max_sequences: - genome: 3000 - G: 3000 - F: 3000 - F-antibody-escape: 2000 - - exclude_where: - recent: ["qc.overallStatus=bad"] - background: ["qc.overallStatus=bad", "qc.overallStatus=mediocre"] - - missing_data_threshold: 1000 - +subsample: + a/genome/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + a/genome/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/genome/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/G/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + a/G/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/G/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/F/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + a/F/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/F/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/F-antibody-escape/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: '1975-01-01' + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + a/F-antibody-escape/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: 6Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 200 + min_date: 12Y + max_date: 6Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 & clade.str.startswith("A.D", na=False) + a/F-antibody-escape/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: 3Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + background: + include: config/include_a.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 200 + min_date: 12Y + max_date: 3Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 & clade.str.startswith("A.D", na=False) + b/genome/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + b/genome/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/genome/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 10000 + query: genome_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/G/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + b/G/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/G/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 600 + query: G_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/F/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: '1975-01-01' + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + b/F/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 6Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 6Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/F/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 3000 + min_date: 3Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 300 + min_date: 12Y + max_date: 3Y + min_length: 1200 + query: F_coverage>0.3 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/F-antibody-escape/all-time: + samples: + sample: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: '1975-01-01' + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + b/F-antibody-escape/6y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: 6Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 200 + min_date: 12Y + max_date: 6Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 & clade.str.startswith("B.D", na=False) + b/F-antibody-escape/3y: + samples: + recent: + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + group_by: + - year + - country + max_sequences: 2000 + min_date: 3Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 + background: + include: config/include_b.txt + exclude: config/outliers_ppx.txt + exclude_where: + - qc.overallStatus=bad + - qc.overallStatus=mediocre + group_by: + - year + - country + max_sequences: 200 + min_date: 12Y + max_date: 3Y + min_length: 1200 + query: F_coverage>0.75 & missing_data<1000 & clade.str.startswith("B.D", na=False) files: auspice_config: "config/auspice_config.json" auspice_config_additional_colorings: "config/auspice_config_additional_colorings.json" diff --git a/workflow/snakemake_rules/config.smk b/workflow/snakemake_rules/config.smk index ab284ff..1b537b6 100644 --- a/workflow/snakemake_rules/config.smk +++ b/workflow/snakemake_rules/config.smk @@ -14,6 +14,7 @@ def main(): "results/run_config.yaml", Path(workflow.basedir) / "config.schema.yaml" ) + write_subsample_config() # TODO: move this to nextstrain/shared @@ -39,6 +40,21 @@ def dump_and_validate(dump_path, schema_path): except ValidateError as e: raise InvalidConfigError(str(e)) from e + +def write_subsample_config(): + # TODO: Support custom build names in the workflow and infer from + # config["builds"]. + for a_or_b in ["a", "b"]: + for build_name in ["genome", "G", "F", "F-antibody-escape"]: + for resolution in ["all-time", "6y", "3y"]: + build = f"{a_or_b}/{build_name}/{resolution}" + 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: diff --git a/workflow/snakemake_rules/core.smk b/workflow/snakemake_rules/core.smk index 750f117..b528426 100644 --- a/workflow/snakemake_rules/core.smk +++ b/workflow/snakemake_rules/core.smk @@ -3,7 +3,7 @@ This part of the workflow expects input files sequences = "data/sequences.fasta" metadata = "data/metadata.tsv" """ - +from augur.subsample import get_referenced_files rule index_sequences: @@ -56,121 +56,38 @@ rule newreference: """ -rule filter_recent: - """ - filtering sequences - """ +rule subsample: input: sequences="results/{a_or_b}/sequences.fasta", metadata="results/{a_or_b}/metadata.tsv", sequence_index=rules.index_sequences.output, - exclude=config["exclude"], + config="results/{a_or_b}/{build_name}/{resolution}/subsample_config.yaml", + referenced_files=lambda w: get_referenced_files(f"results/{w.a_or_b}/{w.build_name}/{w.resolution}/subsample_config.yaml"), output: - sequences=build_dir - + "/{a_or_b}/{build_name}/{resolution}/filtered_recent.fasta", + sequences=build_dir + "/{a_or_b}/{build_name}/{resolution}/subsampled.fasta", log: - "logs/filter_recent_{a_or_b}_{build_name}_{resolution}.txt" + "logs/subsample_{a_or_b}_{build_name}_{resolution}.txt", benchmark: - "benchmarks/filter_recent_{a_or_b}_{build_name}_{resolution}.txt" + "benchmarks/subsample_{a_or_b}_{build_name}_{resolution}.txt", params: - group_by=config["filter"]["group_by"], - min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter"]["min_coverage"][w.build_name]}', - min_length=lambda w: config["filter"]["min_length"][w.build_name], - subsample_max_sequences=lambda w: config["filter"][ - "subsample_max_sequences" - ][w.build_name], strain_id=config["strain_id_field"], - min_date=lambda w: config["filter"]["resolutions"][w.resolution]["min_date"], - exclude_where=config["filter"]["exclude_where"]["recent"], - missing_data_threshold=config["filter"]["missing_data_threshold"], shell: r""" exec &> >(tee {log:q}) - augur filter \ + augur subsample \ --sequences {input.sequences} \ --sequence-index {input.sequence_index} \ --metadata {input.metadata} \ --metadata-id-columns {params.strain_id} \ - --exclude {input.exclude} \ - --exclude-where {params.exclude_where:q} \ - --min-date {params.min_date} \ - --min-length {params.min_length} \ - --output {output.sequences} \ - --group-by {params.group_by} \ - --subsample-max-sequences {params.subsample_max_sequences} \ - --query '({params.min_coverage}) & missing_data<{params.missing_data_threshold}' - """ - - -rule filter_background: - """ - filtering sequences - """ - input: - sequences="results/{a_or_b}/sequences.fasta", - metadata="results/{a_or_b}/metadata.tsv", - sequence_index=rules.index_sequences.output, - include="config/include_{a_or_b}.txt", - exclude=config["exclude"], - output: - sequences=build_dir - + "/{a_or_b}/{build_name}/{resolution}/filtered_background.fasta", - metadata=build_dir - + "/{a_or_b}/{build_name}/{resolution}/filtered_background_metadata.tsv", - log: - "logs/filter_background_{a_or_b}_{build_name}_{resolution}.txt" - benchmark: - "benchmarks/filter_background_{a_or_b}_{build_name}_{resolution}.txt" - params: - group_by=config["filter"]["group_by"], - min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter"]["min_coverage"][w.build_name]}', - min_length=lambda w: config["filter"]["min_length"][w.build_name], - subsample_max_sequences=lambda w: int( - config["filter"]["subsample_max_sequences"][w.build_name], - ) - // 10, - strain_id=config["strain_id_field"], - max_date=lambda w: config["filter"]["resolutions"][w.resolution]["min_date"], - min_date=lambda w: config["filter"]["resolutions"][w.resolution][ - "background_min_date" - ], - exclude_where=config["filter"]["exclude_where"]["background"], - missing_data_threshold=config["filter"]["missing_data_threshold"], - clade_prefix=lambda w: f"{w.a_or_b.upper()}.D", - shell: - r""" - exec &> >(tee {log:q}) - - augur filter \ - --sequences {input.sequences} \ - --sequence-index {input.sequence_index} \ - --metadata {input.metadata} \ - --metadata-id-columns {params.strain_id} \ - --include {input.include} \ - --exclude {input.exclude} \ - --exclude-where {params.exclude_where:q} \ - --min-date {params.min_date} \ - --max-date {params.max_date} \ - --min-length {params.min_length} \ - --output-sequences {output.sequences} \ - --output-metadata {output.metadata} \ - --group-by {params.group_by} \ - --subsample-max-sequences {params.subsample_max_sequences} \ - --query '({params.min_coverage}) & missing_data<{params.missing_data_threshold} & clade.str.startswith("{params.clade_prefix}", na=False)' + --config {input.config} \ + --output-sequences {output.sequences} """ rule combine_samples: input: subsamples=lambda w: ( - ( - [ - rules.filter_recent.output.sequences, - rules.filter_background.output.sequences, - ] - if "background_min_date" in config["filter"]["resolutions"][w.resolution] - else [rules.filter_recent.output.sequences] - ) + [rules.subsample.output.sequences] # potentially add sequences sampled to include maximum escape sequences + ( [ From c33cf9220b5a80a7546f5e09f1a798b176852ab4 Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Wed, 5 Aug 2026 17:22:55 -0700 Subject: [PATCH 3/4] Generate default config YAML The script restores some of the logic that was originally in the hardcoded filter implementation, allowing easier bulk edits to the subsample config. One downside is that the generated file is less readable with a strict YAML style and no comments. Comments have been moved to the script, but ideally they'd live in a schema which is used to generate user-facing docs. --- .gitattributes | 5 + README.md | 9 +- config/configfile.yaml | 135 +++++++-------- scripts/generate_default_config.py | 259 +++++++++++++++++++++++++++++ 4 files changed, 334 insertions(+), 74 deletions(-) create mode 100644 scripts/generate_default_config.py diff --git a/.gitattributes b/.gitattributes index 46610df..7c5d1ec 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,3 +1,8 @@ # Allow Git to decide if file is text or binary # Always use LF line endings even on Windows. * text=auto eol=lf + +# This is a large generated file that, while text, it is not useful to +# routinely show the diff of. A diff can be forced as needed, e.g. with `git +# diff --text`. +/config/configfile.yaml -diff diff --git a/README.md b/README.md index fe171cc..13bc190 100644 --- a/README.md +++ b/README.md @@ -46,7 +46,14 @@ Once you've run the build, you can view the results with: ## Configuration -The default configuration is in [`config/configfile.yaml`](./config/configfile.yaml). +The default configuration is generated by +[scripts/generate_default_config.py](./scripts/generate_default_config.py). It +can be run with Nextstrain CLI: + +```sh +nextstrain shell . -c 'python scripts/generate_default_config.py' +``` + The workflow is contained in the [Snakefile](Snakefile) with included [rules](workflow/snakemake_rules/). Each rule specifies its file inputs and outputs and pulls its parameters from the config. There is little redirection and each diff --git a/config/configfile.yaml b/config/configfile.yaml index 5939779..38d15c0 100644 --- a/config/configfile.yaml +++ b/config/configfile.yaml @@ -1,34 +1,31 @@ -conda_environment: "workflow/envs/nextstrain.yaml" - -genesforglycosylation: ["G", "F"] - -# if the build name has a "-" in it, then the pipeline assumes -# the gene is the part before the first "-", so "F-antibody-escape" -# is parsed as the gene "F". -builds_to_run: ["genome", "G", "F", "F-antibody-escape"] - -resolutions_to_run: ["all-time", "6y", "3y"] - -subtypes: ['a', 'b'] - -# Both files must have a {a_or_b} expandable field to be replaced by "a" or "b" -# depending on if they are specified in the `subtypes` param above +# [DO NOT EDIT] This file was generated by scripts/generate_default_config.py. +conda_environment: workflow/envs/nextstrain.yaml +genesforglycosylation: +- G +- F +builds_to_run: +- genome +- G +- F +- F-antibody-escape +resolutions_to_run: +- all-time +- 6y +- 3y +subtypes: +- a +- b inputs: - - name: ppx_open - metadata: "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata.tsv.gz" - sequences: "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences.fasta.xz" - - name: ppx_restricted - metadata: "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata_restricted.tsv.gz" - sequences: "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences_restricted.fasta.xz" - -exclude: "config/outliers_ppx.txt" - -description: "config/description.md" - -strain_id_field: "accession" -display_strain_field: "strain" - - +- name: ppx_open + metadata: https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata.tsv.gz + sequences: https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences.fasta.xz +- name: ppx_restricted + metadata: https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata_restricted.tsv.gz + sequences: https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences_restricted.fasta.xz +exclude: config/outliers_ppx.txt +description: config/description.md +strain_id_field: accession +display_strain_field: strain subsample: a/genome/all-time: samples: @@ -567,48 +564,41 @@ subsample: min_length: 1200 query: F_coverage>0.75 & missing_data<1000 & clade.str.startswith("B.D", na=False) files: - auspice_config: "config/auspice_config.json" - auspice_config_additional_colorings: "config/auspice_config_additional_colorings.json" - auspice_config_f_antibody_escape: "config/auspice_config_dms-defaults.json" - auspice_config_non-genome_builds: "config/auspice_config_non-genome.json" - + auspice_config: config/auspice_config.json + auspice_config_additional_colorings: config/auspice_config_additional_colorings.json + auspice_config_f_antibody_escape: config/auspice_config_dms-defaults.json + auspice_config_non-genome_builds: config/auspice_config_non-genome.json refine: - coalescent: "opt" - date_inference: "marginal" + coalescent: opt + date_inference: marginal clock_filter_iqd: 4 - divergence_units: "mutations-per-site" - + divergence_units: mutations-per-site ancestral: - inference: "joint" - + inference: joint cds: - F: "F" - G: "G" - genome: "F" - F-antibody-escape: "F" - + F: F + G: G + genome: F + F-antibody-escape: F traits: - columns: "country region" - + columns: country region frequencies: resolutions: all-time: - min_date: "1975-01-01" + min_date: '1975-01-01' 6y: min_date: 6Y 3y: min_date: 3Y - nextclade_attributes: a: - name: "RSV-A NextClade using real-time tree" - reference_name: "hRSV/A/England/397/2017" - accession: "EPI_ISL_412866" + name: RSV-A NextClade using real-time tree + reference_name: hRSV/A/England/397/2017 + accession: EPI_ISL_412866 b: - name: "RSV-B NextClade using real-time tree" - reference_name: "hRSV/B/Australia/VIC-RCH056/2019" - accession: "EPI_ISL_1653999" - + name: RSV-B NextClade using real-time tree + reference_name: hRSV/B/Australia/VIC-RCH056/2019 + accession: EPI_ISL_1653999 filter_for_pre_subsample_alignment: min_length: genome: 10000 @@ -622,24 +612,23 @@ filter_for_pre_subsample_alignment: F-antibody-escape: 0.75 resolutions: all-time: - min_date: "1975-01-01" + min_date: '1975-01-01' 6y: min_date: 6Y 3y: min_date: 3Y - -# configuration specific to the F deep mutational scanning antibody escape data f_dms_data: dms-data/all_antibodies.csv -f_dms_antibodies: # columns in `f_dms_data` with per-mutation escape - - Clesrovimab-Fab - - Clesrovimab-IgG - - Nirsevimab-Fab - - Nirsevimab-IgG -dms_only_positive_escape: true # for DMS escape values, set any values < 0 to 0 -enrich_antibody_escape: # additional filtering for antibody escape sequences added to tree - F-antibody-escape: # for this build, enrich by these criteria - nseqs_per_antibody_scoretype: 500 # add this many sequences for each antibody for total and max escape - group_by: [country, year] # group by these variables - max_identical_f_prot_muts: 2 # for each group, no more than this many w identical F protein mutations - max_identical_max_escape_mut: 6 # for each group, no more than this many w same top F escape mutation - +f_dms_antibodies: +- Clesrovimab-Fab +- Clesrovimab-IgG +- Nirsevimab-Fab +- Nirsevimab-IgG +dms_only_positive_escape: true +enrich_antibody_escape: + F-antibody-escape: + nseqs_per_antibody_scoretype: 500 + group_by: + - country + - year + max_identical_f_prot_muts: 2 + max_identical_max_escape_mut: 6 diff --git a/scripts/generate_default_config.py b/scripts/generate_default_config.py new file mode 100644 index 0000000..2576c52 --- /dev/null +++ b/scripts/generate_default_config.py @@ -0,0 +1,259 @@ +"""Generate the default configfile for the RSV workflow.""" + +from pathlib import Path +import yaml + + +class NoAliasDumper(yaml.SafeDumper): + def ignore_aliases(self, data): + return True + + +def main(): + config = generate_config() + + path = Path(__file__).resolve().parent.parent / "config" / "configfile.yaml" + + with open(path, "w") as f: + print("# [DO NOT EDIT] This file was generated by scripts/generate_default_config.py.", file=f) + yaml.dump(config, f, Dumper=NoAliasDumper, sort_keys=False, width=1000) + + +def generate_config(): + return { + "conda_environment": "workflow/envs/nextstrain.yaml", + + "genesforglycosylation": ["G", "F"], + # if the build name has a "-" in it, then the pipeline assumes + # the gene is the part before the first "-", so "F-antibody-escape" + # is parsed as the gene "F". + "builds_to_run": ["genome", "G", "F", "F-antibody-escape"], + + "resolutions_to_run": ["all-time", "6y", "3y"], + + "subtypes": ["a", "b"], + + # Both files must have a {a_or_b} expandable field to be replaced by "a" or "b" + # depending on if they are specified in the `subtypes` param above + "inputs": [ + { + "name": "ppx_open", + "metadata": "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata.tsv.gz", + "sequences": "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences.fasta.xz", + }, + { + "name": "ppx_restricted", + "metadata": "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/metadata_restricted.tsv.gz", + "sequences": "https://data.nextstrain.org/files/workflows/rsv/{a_or_b}/sequences_restricted.fasta.xz", + }, + ], + + "exclude": "config/outliers_ppx.txt", + + "description": "config/description.md", + + "strain_id_field": "accession", + "display_strain_field": "strain", + + "subsample": generate_subsample_config(), + + "files": { + "auspice_config": "config/auspice_config.json", + "auspice_config_additional_colorings": "config/auspice_config_additional_colorings.json", + "auspice_config_f_antibody_escape": "config/auspice_config_dms-defaults.json", + "auspice_config_non-genome_builds": "config/auspice_config_non-genome.json", + }, + + "refine": { + "coalescent": "opt", + "date_inference": "marginal", + "clock_filter_iqd": 4, + "divergence_units": "mutations-per-site", + }, + + "ancestral": { + "inference": "joint", + }, + + "cds": { + "F": "F", + "G": "G", + "genome": "F", + "F-antibody-escape": "F", + }, + + "traits": { + "columns": "country region", + }, + + "frequencies": { + "resolutions": { + "all-time": {"min_date": "1975-01-01"}, + "6y": {"min_date": "6Y"}, + "3y": {"min_date": "3Y"}, + }, + }, + + "nextclade_attributes": { + "a": { + "name": "RSV-A NextClade using real-time tree", + "reference_name": "hRSV/A/England/397/2017", + "accession": "EPI_ISL_412866", + }, + "b": { + "name": "RSV-B NextClade using real-time tree", + "reference_name": "hRSV/B/Australia/VIC-RCH056/2019", + "accession": "EPI_ISL_1653999", + }, + }, + + "filter_for_pre_subsample_alignment": { + "min_length": { + "genome": 10000, + "G": 600, + "F": 1200, + "F-antibody-escape": 1200, + }, + "min_coverage": { + "genome": 0.3, + "G": 0.3, + "F": 0.3, + "F-antibody-escape": 0.75, + }, + "resolutions": { + "all-time": {"min_date": "1975-01-01"}, + "6y": {"min_date": "6Y"}, + "3y": {"min_date": "3Y"}, + }, + }, + + # configuration specific to the F deep mutational scanning antibody escape data + "f_dms_data": "dms-data/all_antibodies.csv", + "f_dms_antibodies": [ # columns in `f_dms_data` with per-mutation escape + "Clesrovimab-Fab", + "Clesrovimab-IgG", + "Nirsevimab-Fab", + "Nirsevimab-IgG", + ], + "dms_only_positive_escape": True, # for DMS escape values, set any values < 0 to 0 + "enrich_antibody_escape": { # additional filtering for antibody escape sequences added to tree + "F-antibody-escape": { # for this build, enrich by these criteria + "nseqs_per_antibody_scoretype": 500, # add this many sequences for each antibody for total and max escape + "group_by": ["country", "year"], # group by these variables + "max_identical_f_prot_muts": 2, # for each group, no more than this many w identical F protein mutations + "max_identical_max_escape_mut": 6, # for each group, no more than this many w same top F escape mutation + }, + }, + } + + +def generate_subsample_config(): + """Generate the expanded subsample configuration for all builds.""" + subsample = {} + + subtypes = ["a", "b"] + + builds = { + "genome": { + "coverage_col": "genome_coverage", + "min_coverage": 0.3, + "min_length": 10000, + "recent_max_seqs": 3000, + }, + "G": { + "coverage_col": "G_coverage", + "min_coverage": 0.3, + "min_length": 600, + "recent_max_seqs": 3000, + }, + "F": { + "coverage_col": "F_coverage", + "min_coverage": 0.3, + "min_length": 1200, + "recent_max_seqs": 3000, + }, + "F-antibody-escape": { + "coverage_col": "F_coverage", + "min_coverage": 0.75, + "min_length": 1200, + "recent_max_seqs": 2000, + }, + } + + group_by = ["year", "country"] + exclude_file = "config/outliers_ppx.txt" + recent_exclude = ["qc.overallStatus=bad"] + background_exclude = ["qc.overallStatus=bad", "qc.overallStatus=mediocre"] + + for a_or_b in subtypes: + for build_name, build_info in builds.items(): + recent_query = f"{build_info['coverage_col']}>{build_info['min_coverage']} & missing_data<1000" + background_query = f'{build_info["coverage_col"]}>{build_info["min_coverage"]} & missing_data<1000 & clade.str.startswith("{a_or_b.upper()}.D", na=False)' + + resolutions = { + "all-time": { + "sample": { + "exclude": exclude_file, + "exclude_where": recent_exclude, + "group_by": group_by, + "max_sequences": build_info["recent_max_seqs"], + "min_date": "1975-01-01", + "min_length": build_info["min_length"], + "query": recent_query, + }, + }, + "6y": { + "recent": { + "exclude": exclude_file, + "exclude_where": recent_exclude, + "group_by": group_by, + "max_sequences": build_info["recent_max_seqs"], + "min_date": "6Y", + "min_length": build_info["min_length"], + "query": recent_query, + }, + "background": { + "include": f"config/include_{a_or_b}.txt", + "exclude": exclude_file, + "exclude_where": background_exclude, + "group_by": group_by, + "max_sequences": int(build_info["recent_max_seqs"] / 10), + "min_date": "12Y", + "max_date": "6Y", + "min_length": build_info["min_length"], + "query": background_query, + }, + }, + "3y": { + "recent": { + "exclude": exclude_file, + "exclude_where": recent_exclude, + "group_by": group_by, + "max_sequences": build_info["recent_max_seqs"], + "min_date": "3Y", + "min_length": build_info["min_length"], + "query": recent_query, + }, + "background": { + "include": f"config/include_{a_or_b}.txt", + "exclude": exclude_file, + "exclude_where": background_exclude, + "group_by": group_by, + "max_sequences": int(build_info["recent_max_seqs"] / 10), + "min_date": "12Y", + "max_date": "3Y", + "min_length": build_info["min_length"], + "query": background_query, + }, + }, + } + + for resolution, samples in resolutions.items(): + key = f"{a_or_b}/{build_name}/{resolution}" + subsample[key] = {"samples": samples} + + return subsample + + +if __name__ == "__main__": + main() From d96e18a335c1c811795f56440768f0c04e8a452f Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Mon, 10 Aug 2026 15:18:13 -0700 Subject: [PATCH 4/4] Rename rule to filter_for_f_antibody_escape The new name makes it more obvious that this rule is only used when build_name=F-antibody-escape. --- CHANGELOG.md | 2 +- config.schema.yaml | 2 +- config/configfile.yaml | 2 +- scripts/generate_default_config.py | 2 +- workflow/snakemake_rules/core.smk | 14 +++++++------- 5 files changed, 11 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0386aec..253c007 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,6 @@ We use this CHANGELOG to document breaking changes, new features, bug fixes, and ## 2026 -* TBD: The `filter` section in phylogenetic workflow configuration has been replaced by `subsample`/`custom_subsample` for subsampling, and `filter_for_pre_subsample_alignment` for initial quality filtering. **This is a breaking change**. +* TBD: The `filter` section in phylogenetic workflow configuration has been replaced by `subsample`/`custom_subsample` for subsampling, and `filter_for_f_antibody_escape` for initial quality filtering. **This is a breaking change**. * NOTE: The workflow does not yet support proximal samples. * 11 August 2026: Phylogenetic workflow configuration is now validated against a strict schema. The workflow will error if your configuration has extraneous entries that were previously ignored. diff --git a/config.schema.yaml b/config.schema.yaml index 94859fe..ae26a49 100644 --- a/config.schema.yaml +++ b/config.schema.yaml @@ -168,7 +168,7 @@ properties: type: string accession: type: string - filter_for_pre_subsample_alignment: + filter_for_f_antibody_escape: type: object additionalProperties: false properties: diff --git a/config/configfile.yaml b/config/configfile.yaml index 38d15c0..f9430e4 100644 --- a/config/configfile.yaml +++ b/config/configfile.yaml @@ -599,7 +599,7 @@ nextclade_attributes: name: RSV-B NextClade using real-time tree reference_name: hRSV/B/Australia/VIC-RCH056/2019 accession: EPI_ISL_1653999 -filter_for_pre_subsample_alignment: +filter_for_f_antibody_escape: min_length: genome: 10000 G: 600 diff --git a/scripts/generate_default_config.py b/scripts/generate_default_config.py index 2576c52..5d28fb4 100644 --- a/scripts/generate_default_config.py +++ b/scripts/generate_default_config.py @@ -107,7 +107,7 @@ def generate_config(): }, }, - "filter_for_pre_subsample_alignment": { + "filter_for_f_antibody_escape": { "min_length": { "genome": 10000, "G": 600, diff --git a/workflow/snakemake_rules/core.smk b/workflow/snakemake_rules/core.smk index b528426..e4a9af4 100644 --- a/workflow/snakemake_rules/core.smk +++ b/workflow/snakemake_rules/core.smk @@ -137,7 +137,7 @@ rule get_nextclade_dataset: """ -rule filter_for_pre_subsample_alignment: +rule filter_for_f_antibody_escape: """ Do the quality filtering applied to each sequence set before subsampling """ @@ -148,14 +148,14 @@ rule filter_for_pre_subsample_alignment: output: sequences=build_dir + "/{a_or_b}/{build_name}/{resolution}/pre_subsample/filtered_for_alignment.fasta", log: - "logs/filter_for_pre_subsample_alignment_{a_or_b}_{build_name}_{resolution}.txt" + "logs/filter_for_f_antibody_escape_{a_or_b}_{build_name}_{resolution}.txt" benchmark: - "benchmarks/filter_for_pre_subsample_alignment_{a_or_b}_{build_name}_{resolution}.txt" + "benchmarks/filter_for_f_antibody_escape_{a_or_b}_{build_name}_{resolution}.txt" params: - min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter_for_pre_subsample_alignment"]["min_coverage"][w.build_name]}', - min_length=lambda w: config["filter_for_pre_subsample_alignment"]["min_length"][w.build_name], + min_coverage=lambda w: f'{w.build_name.split("-")[0]}_coverage>{config["filter_for_f_antibody_escape"]["min_coverage"][w.build_name]}', + min_length=lambda w: config["filter_for_f_antibody_escape"]["min_length"][w.build_name], strain_id=config["strain_id_field"], - min_date=lambda w: config["filter_for_pre_subsample_alignment"]["resolutions"][w.resolution]["min_date"], + min_date=lambda w: config["filter_for_f_antibody_escape"]["resolutions"][w.resolution]["min_date"], shell: r""" exec &> >(tee {log:q}) @@ -178,7 +178,7 @@ rule align_pre_subsample_sequences: Aligning all pre-subsampled quality-filtered sequences """ input: - sequences=rules.filter_for_pre_subsample_alignment.output.sequences, + sequences=rules.filter_for_f_antibody_escape.output.sequences, dataset=rules.get_nextclade_dataset.output.dataset, output: alignment=build_dir + "/{a_or_b}/{build_name}/{resolution}/pre_subsample/sequences.aligned.fasta",