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/CHANGELOG.md b/CHANGELOG.md index 98b7b83..253c007 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,4 +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_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/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.schema.yaml b/config.schema.yaml index 774060b..ae26a49 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 @@ -194,6 +168,33 @@ properties: type: string accession: type: string + filter_for_f_antibody_escape: + 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..f9430e4 100644 --- a/config/configfile.yaml +++ b/config/configfile.yaml @@ -1,124 +1,634 @@ -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" - - -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 - +- 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: + 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" - 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" - -# configuration specific to the F deep mutational scanning antibody escape data + name: RSV-B NextClade using real-time tree + reference_name: hRSV/B/Australia/VIC-RCH056/2019 + accession: EPI_ISL_1653999 +filter_for_f_antibody_escape: + 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 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..5d28fb4 --- /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_f_antibody_escape": { + "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() 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 96f2fe9..e4a9af4 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 + ( [ @@ -220,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 """ @@ -231,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"]["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_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"]["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}) @@ -261,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",