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 d4d7144..f72125a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,9 @@ # CHANGELOG We use this CHANGELOG to document breaking changes, new features, bug fixes, and config value changes that may affect both the usage of the workflows and the outputs of the workflows. + +## 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. +* TBD: 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 new file mode 100644 index 0000000..80ea2a4 --- /dev/null +++ b/config.schema.yaml @@ -0,0 +1,227 @@ +$schema: "https://json-schema.org/draft/2020-12/schema" +title: RSV Phylogenetic Workflow Configuration +description: >- + This is the schema for the Nextstrain rsv phylogenetic workflow's + 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 + propertyNames: + title: Subtype + description: Subtype ('a' or 'b') + + per_build_map: &per_build_map + type: object + additionalProperties: false + propertyNames: + title: Build name + description: Build name ('genome', 'G', 'F', or 'F-antibody-escape') + + per_resolution_map: &per_resolution_map + type: object + additionalProperties: false + propertyNames: + title: Resolution name + description: Resolution name ('all-time', '6y', or '3y') + + input_item: + type: object + additionalProperties: false + required: + - name + anyOf: + - required: [metadata] + - required: [sequences] + properties: + name: + type: string + metadata: + type: string + sequences: + type: string + +type: object +additionalProperties: false +properties: + conda_environment: + type: string + genesforglycosylation: + type: array + items: + type: string + builds_to_run: + type: array + items: + type: string + resolutions_to_run: + type: array + items: + type: string + subtypes: + type: array + items: + type: string + inputs: + type: array + items: + $ref: "#/$defs/input_item" + additional_inputs: + type: array + items: + $ref: "#/$defs/input_item" + exclude: + type: string + description: + type: string + strain_id_field: + type: string + display_strain_field: + type: string + 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/v1#/$defs/schemaForUnalignedSequences" + 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 + properties: + auspice_config: + type: string + auspice_config_additional_colorings: + type: string + auspice_config_f_antibody_escape: + type: string + auspice_config_non-genome_builds: + type: string + refine: + type: object + additionalProperties: false + properties: + coalescent: + type: string + date_inference: + type: string + clock_filter_iqd: + type: number + divergence_units: + type: string + ancestral: + type: object + additionalProperties: false + properties: + inference: + type: string + cds: + <<: *per_build_map + patternProperties: + "^.*$": + type: string + traits: + type: object + additionalProperties: false + properties: + columns: + type: [string, array] + items: + type: string + frequencies: + type: object + additionalProperties: false + properties: + resolutions: + <<: *per_resolution_map + patternProperties: + "^.*$": + type: object + additionalProperties: false + properties: + min_date: + type: string + nextclade_attributes: + <<: *per_subtype_map + patternProperties: + "^.*$": + type: object + additionalProperties: false + properties: + name: + type: string + reference_name: + 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: + type: array + items: + type: string + dms_only_positive_escape: + type: boolean + enrich_antibody_escape: + <<: *per_build_map + patternProperties: + "^.*$": + type: object + additionalProperties: false + properties: + nseqs_per_antibody_scoretype: + type: integer + group_by: + type: array + items: + type: string + max_identical_f_prot_muts: + type: integer + max_identical_max_escape_mut: + type: integer + custom_rules: + type: array + description: Custom Snakemake rule files to include. If used, this will disable config schema validation. + items: + type: string diff --git a/config/configfile.yaml b/config/configfile.yaml index 3aa617b..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..dabdf24 --- /dev/null +++ b/scripts/generate_default_config.py @@ -0,0 +1,272 @@ +"""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.""" + subtypes = ["a", "b"] + + builds = { + "genome": { + "gene_coverage": "genome_coverage", + "min_coverage": 0.3, + "min_length": 10000, + "recent_max_seqs": 3000, + "bg_max_seqs": 300, + }, + "G": { + "gene_coverage": "G_coverage", + "min_coverage": 0.3, + "min_length": 600, + "recent_max_seqs": 3000, + "bg_max_seqs": 300, + }, + "F": { + "gene_coverage": "F_coverage", + "min_coverage": 0.3, + "min_length": 1200, + "recent_max_seqs": 3000, + "bg_max_seqs": 300, + }, + "F-antibody-escape": { + "gene_coverage": "F_coverage", + "min_coverage": 0.75, + "min_length": 1200, + "recent_max_seqs": 2000, + "bg_max_seqs": 200, + }, + } + + subsample = {} + for a_or_b in subtypes: + clade_prefix = f"{a_or_b.upper()}.D" + include_file = f"config/include_{a_or_b}.txt" + + for build_name, build_info in builds.items(): + coverage_var = build_info["gene_coverage"] + min_cov = build_info["min_coverage"] + min_len = build_info["min_length"] + recent_seqs = build_info["recent_max_seqs"] + bg_seqs = build_info["bg_max_seqs"] + + recent_query = f"{coverage_var}>{min_cov} & missing_data<1000" + bg_query = f"{coverage_var}>{min_cov} & missing_data<1000 & clade.str.startswith(\"{clade_prefix}\", na=False)" + + resolutions = { + "all-time": { + "sample": { + "exclude": "config/outliers_ppx.txt", + "exclude_where": ["qc.overallStatus=bad"], + "group_by": ["year", "country"], + "max_sequences": recent_seqs, + "min_date": "1975-01-01", + "min_length": min_len, + "query": recent_query, + }, + }, + "6y": { + "recent": { + "exclude": "config/outliers_ppx.txt", + "exclude_where": ["qc.overallStatus=bad"], + "group_by": ["year", "country"], + "max_sequences": recent_seqs, + "min_date": "6Y", + "min_length": min_len, + "query": recent_query, + }, + "background": { + "include": include_file, + "exclude": "config/outliers_ppx.txt", + "exclude_where": [ + "qc.overallStatus=bad", + "qc.overallStatus=mediocre", + ], + "group_by": ["year", "country"], + "max_sequences": bg_seqs, + "min_date": "12Y", + "max_date": "6Y", + "min_length": min_len, + "query": bg_query, + }, + }, + "3y": { + "recent": { + "exclude": "config/outliers_ppx.txt", + "exclude_where": ["qc.overallStatus=bad"], + "group_by": ["year", "country"], + "max_sequences": recent_seqs, + "min_date": "3Y", + "min_length": min_len, + "query": recent_query, + }, + "background": { + "include": include_file, + "exclude": "config/outliers_ppx.txt", + "exclude_where": [ + "qc.overallStatus=bad", + "qc.overallStatus=mediocre", + ], + "group_by": ["year", "country"], + "max_sequences": bg_seqs, + "min_date": "12Y", + "max_date": "3Y", + "min_length": min_len, + "query": bg_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 199b4a9..1b537b6 100644 --- a/workflow/snakemake_rules/config.smk +++ b/workflow/snakemake_rules/config.smk @@ -5,4 +5,58 @@ OUTPUTS: results/run_config.yaml """ -write_config("results/run_config.yaml") +import sys +from augur.validate import load_json_schema_locally, validate_json, ValidateError +from pathlib import Path + +def main(): + dump_and_validate( + "results/run_config.yaml", + Path(workflow.basedir) / "config.schema.yaml" + ) + write_subsample_config() + + +# TODO: move this to nextstrain/shared +# Copied from measles with some cleanups +# +def dump_and_validate(dump_path, schema_path): + """ + Write Snakemake's 'config' variable to a file, then validate it against the + schema. Do both in the same function so that the validation output can + easily reference the path of the dumped config for inspection. + """ + global config + + write_config(dump_path) + + if "custom_rules" in config: + print("WARNING: Skipping config schema validation because custom rules are defined.", file=sys.stderr) + return + + try: + validator = load_json_schema_locally(schema_path) + validate_json(config, validator, dump_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: + print(f"ERROR: {e}", file=sys.stderr) + exit(1) 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",