Build interactive Nextstrain phylogenetic trees of influenza HA from different subtypes and annotate by MHC II usage.
Specifically, these trees are designed to show a broad range of sequences from different species for each subtype. The dropdowns at the left enable the trees to be colored by:
- The directly measured titers of pseudovirus with tested HAs on sialic-acid versus MHC-II versus no-sialic-acid cells. These are the properties with names like log10__titers.
- The identity of amino acids at each protein site (color by Genotype, see links below to understand site numbering scheme).
- The predicted tufted duck MHC-II usage as determined by summing effects of all mutations as measured in H5 or H7 HA pseudovirus deep mutational scanning, quantified as the total effect of all mutations of the max magnitude effect of the mutation with the single largest impact. Note that these predictions assume mutations are non-epistatic and have the same effects across genetic backgrounds, which will not always be true.
Interactive Nextstrain trees are at https://nextstrain.org/community/jbloomlab/flu-ha-mhcii-usage-trees, specifically:
- H1 HAs:
- H2 HAs:
- H3 HAs:
- H5 HAs:
- H7 HAs:
Trees are built for different HA subtypes. The configuration for these trees is specified under trees in config.yaml. The trees are built subsampling sequences from Genbank, and you can also specify specific accessions to include under manual_adds and accessions_to_include, and accessions to exclude under exclude.
Molecular clock outliers are excluded according to clock-filter-iqd; the value of bad_dates_in_keep_accessions_action specifies what to do if strains specified for inclusion fail this filter.
The trees are annotated by the tufted duck MHC II usage predicted under an additive model from the measurements in data/mutation_effect_data; these predictions are applied to all sequences.
The trees are also annotated by direct measurements of MHC II usage (as relative titers for different strains) for specific strains in the data/trees/{subtype}/strain_annotations.tsv files.; the value of missing_strain_annotations_action specifies what to do if strains with specified titers are not found.
See ./data/README.md for more explanation of the input data.
The final trees are placed in ./auspice where they can be viewed as Nextstrain Community Builds at the links above.
Structure of repo:
- config.yaml: configuration
- environment.yaml: conda environment
- Snakefile: Snakemake file with pipeline
- ./scripts/: scripts used by pipeline
- ./data/: input data (see data/README.md for details on the subdirectory structure, mutation effect data, and the
build_protein_sites_and_mutation_effects.pyscript that generatesprotein_sites.tsvandmutation_effects.tsvfor each tree) ./results/: results created by pipeline- ./auspice/: final JSON trees for visualization via https://auspice.us/ or Nextstrain community builds.
To run the pipeline, build and activate the conda environment in environment.yaml, and then run:
snakemake --software-deployment-method conda -j <n_cpus>
Download all Influenza A virus sequences and metadata from GenBank using ncbi-datasets-cli.
The download includes genomic sequences, CDS nucleotide sequences, gene annotations, and BioSample metadata as a single zip archive.
The taxon ID for the downloaded sequences is specified in config.yaml as datasets_taxid.
The datestamp file in results/ncbi_dataset/ncbi_dataset_download_date.txt records the download date.
Extract genome metadata, genomic FASTA headers, and per-gene annotation from the NCBI zip file using dataformat.
This processes the large zip once so downstream per-subtype rules can work from smaller TSV files.
For each tree defined in config.yaml, extract full-length hemagglutinin (HA) coding sequences and associated metadata.
HA accessions are identified by segment 4 in the genome metadata or gene-name of "HA", "ha", or "HA1" in the annotation report (some older NCBI records lack the segment field but have the gene annotation).
Subtypes are parsed from genomic FASTA headers and filtered by the subtype regex in the config (e.g., H5N\d{1,2} for all H5Nx) and the CDS length range specified in cds_length_range YAML pointed to in config (either bound can be null for no limit).
For accessions with missing collection dates, the year is parsed from the strain name as a fallback (e.g., "1930" from "A/swine/Iowa/15/1930"); sequences with no date even after this fallback are dropped, as are sequences with ambiguous nucleotides or that are not valid coding sequences.
Outputs per-tree FASTA and metadata TSV files in results/trees/{tree}/.
Download the NCBI taxonomy dump files for use by taxonkit.
These files are used to look up host taxonomy (class, order) from NCBI taxonomy IDs.
NCBI's structured host / host_tax_id fields are empty for a small fraction of records. For each such row, this rule parses the 2nd /-delimited token from the strain name (e.g. "chicken" in A/chicken/Scotland/1959) and, if that token is mapped to a tax ID in data/strain_token_host_map.tsv, fills in host and host_tax_id. Tokens that are locations, environmental samples, processed products, or ambiguous are explicitly listed in the map with a blank tax ID so they remain empty-host. Unknown tokens (not in the map) are logged as warnings but also left empty — the default is conservative.
For each tree, annotate the metadata with host taxonomy information using taxonkit to look up lineage from NCBI taxonomy IDs.
Renames host to host_specific (the NCBI Taxonomy scientific name), and adds host_general (classifying each host as "human", "avian", "swine", "bovine", "equine", "carnivore", or "other mammal") and host_order (the taxonomic order).
Outputs the annotated metadata TSV in results/trees/{tree}/.
Subsample the HA sequences for each tree using augur subsample.
The subsampling strategy (grouping, max sequences, date filtering, include/exclude lists) is configured per tree in config.yaml under augur_subsample.
Outputs subsampled FASTA and metadata TSV files in results/trees/{tree}/.
For each tree, collapse host orders that appear at frequency below collapse_low_freq_host_order (default 0.005) in the subsampled data to "other".
This reduces the number of colors needed when visualizing by host order in auspice.
Outputs modified subsampled metadata in results/trees/{tree}/.
Align the subsampled HA CDS sequences for each tree using augur align (which wraps mafft).
A per-tree reference sequence (configured as reference_sequence in config.yaml) guides the alignment; the reference is removed from the output.
Outputs the aligned FASTA in results/trees/{tree}/.
Infer a maximum-likelihood phylogenetic tree from the alignment using augur tree (which wraps IQ-TREE).
Uses a fixed seed for reproducibility and collapses zero-length branches.
Outputs a raw Newick tree in results/trees/{tree}/ (before any temporal refinement).
Refine the raw tree using augur refine (which wraps TreeTime) to build a time-resolved phylogeny.
Additional augur refine flags are configured per tree in config.yaml under augur_refine; this can include a molecular clock filter (clock-filter-iqd) to remove tips that deviate too far from the root-to-tip regression, and keep-ids exempts specific accessions from clock-filter pruning so that manually requested strains are retained even if they are clock outliers.
Outputs a refined Newick tree, a node-data JSON with branch lengths and inferred dates, and the combined stdout+stderr from augur refine (for downstream parsing) in results/trees/{tree}/.
Parse the augur refine output to extract any strains whose dates were reset because they were flagged as molecular clock outliers (kept in the tree via keep-ids but with date constraints removed).
Outputs a TSV with accession, metadata date, and inferred date for each outlier strain, or an empty table if none.
The script includes safety checks to error (rather than silently report zero outliers) if the augur refine output format has changed.
Infer ancestral nucleotide sequences using augur ancestral.
Mutations are reported relative to the same reference sequence used for alignment.
Outputs a node-data JSON with nucleotide mutations per branch in results/trees/{tree}/.
Translate nucleotide mutations to amino acid mutations using augur translate.
Uses a GFF3 annotation file to define gene regions (SigPep, HA1, HA2) for translation.
Outputs a node-data JSON with amino acid mutations per branch in results/trees/{tree}/.
Score each tree node based on the cumulative effect of its amino acid mutations on measured phenotypes.
Per-mutation effects are provided in a TSV file (configured as mutation_effects in config.yaml).
The script validates that proteins in the effects file match the genes in the tree's amino acid mutations.
For each phenotype (configured under phenotypes), computes the total effect, the max magnitude effect, and the identity of the max magnitude mutation.
Outputs a node-data JSON with phenotype scores in results/trees/{tree}/.
Generate an auspice config JSON that defines continuous viridis color scales for each phenotype's total effect and max magnitude effect scores, and a categorical coloring for the max magnitude mutation.
The color scale is configured per tree as phenotype_color_scale in config.yaml.
Outputs an auspice config JSON in results/trees/{tree}/.
For each tree, read the per-strain annotations TSV (configured as strain_annotations in config.yaml) and map each non-accession column onto the corresponding tree tip.
If the TSV is empty, nothing is annotated.
Numeric columns are colored continuously using the shared phenotype_color_scale from config.yaml; non-numeric columns become categorical colorings and are added to the Auspice filters.
Behavior when an accession in the TSV is not a tip in the final refined tree is controlled by missing_strain_annotations_action in config.yaml (ignore warns and drops those rows; error fails the pipeline).
Outputs a node-data JSON and an auspice config JSON (both consumed by augur export) in results/trees/{tree}/.
For each tree, format the shared description markdown file by replacing {tree} with the tree name (e.g., "H5").
This produces a tree-specific description file in results/trees/{tree}/ that is passed to augur export.
Export interactive auspice v2 JSONs using augur export.
Each tree uses an auspice config file (configured as auspice_config in config.yaml) that defines colorings, filters, display defaults, and metadata; a generated phenotype auspice config with color scales for mutation effect scores; a generated strain annotations auspice config with colorings (and filters) for columns from strain_annotations.tsv; a markdown description file for the tree sidebar; and a title for the tree visualization.
The output files are placed in auspice/ with names like {auspice_prefix}_{tree}.json (the prefix is set in config.yaml, typically matching the repo name for Nextstrain community builds).
For each entry in pruned_trees in config.yaml, post-process the source tree's exported Auspice JSON down to only the tips whose rows in strain_annotations.tsv match a column/value filter (e.g. keep only H5 tips with has_titers == "yes").
Internal nodes left with a single descendant after pruning are collapsed, with their branch mutations merged onto the surviving child so no mutation information is lost.
An optional title can override the Auspice title on the pruned JSON.
Outputs a pruned Auspice JSON in auspice/ alongside the full trees.
For each tree, verify that every accession listed in data/trees/{tree}/accessions_to_include.txt (comments and blank lines are ignored) and every accession in data/trees/{tree}/manual_add_metadata.tsv appears as a leaf in the exported Auspice JSON.
Writes results/trees/{tree}/include_validation.tsv with one row per expected accession (accession, source, in_tree) and errors the pipeline if any are missing.
These TSVs are final rule all targets.