Flu-UShER Pipeline

A Snakemake pipeline for building phylogenetic trees of influenza virus sequences using UShER. The pipeline processes influenza sequences by segment and subtype, creating separate trees for HA and NA segments by subtype (e.g., H1, H3, N1, N2) while combining all subtypes for other segments (PB2, PB1, PA, NP, MP, NS).

Directory Structure

flu-usher/
├── Snakefile                # Main pipeline file
├── config.yaml              # Configuration file
├── data/                    # Input data directories
│   ├── H1N1/                # Example: H1N1 sequences (all segments)
│   ├── H3N2/                # Example: H3N2 sequences (all segments)
│   ├── H5N1/                # Example: H5N1 sequences (all segments)
│   └── H7N9/                # Example: H7N9 sequences (all segments)
│       ├── sequences.fasta  # One or more FASTA files containing sequences
│       └── metadata.xls     # One or more metadata files from GISAID
├── logs/                    # Log files (created by the pipeline)
├── results/                 # Output results (organized by segment/subtype)
│   ├── HA/                  # HA segment results by subtype
│   │   ├── H1/              # H1 subtype tree and files
│   │   ├── H3/              # H3 subtype tree and files
│   │   ├── H5/              # H5 subtype tree and files
│   │   ├── H7/              # H7 subtype tree and files
│   │   └── H9/              # H9 subtype tree and files
│   ├── NA/                  # NA segment results by subtype
│   │   ├── N1/              # N1 subtype tree and files
│   │   ├── N2/              # N2 subtype tree and files
│   │   └── N9/              # N9 subtype tree and files
│   └── PB2/                 # Other segment results (PB1, PA, NP, MP, NS)
│       └── all/             # All subtypes combined
├── scripts/
│   ├── parse_gisaid_data.py                        # Parse GISAID data by segment
│   ├── download_ref_seq.py                         # Download reference sequences
│   ├── curate_and_extract_coding_seqs.py           # Curate alignments and extract coding sequences
│   ├── randomize_alignment.py                      # Randomize alignment order
│   ├── trim_dag.py                                 # Trim suboptimal trees from DAG
│   ├── convert_DAG_protobuf_to_newick_samples.py   # Sample tree from DAG
│   ├── create_root_samples_file.py                 # Create samples file for root extraction
│   ├── extract_root_sequence.py                    # Infer root sequence from tree
│   ├── simplified_host_classifier.py               # Classify hosts into groups (used by augment_metadata.py)
│   ├── augment_metadata.py                         # Add host, geographic, and temporal group columns to metadata
│   ├── create_samples_file.py                      # Create samples file for subtree extraction (generic column filter)
│   ├── prepare_host_annotation.py                  # Build 2-col CSV (isolate_id, host_group) for PastML
│   └── prepare_subtype_annotation.py               # Build 2-col CSV (isolate_id, subtype) for PastML, normalizing "A / H5N1" → "H5N1"
└── notebooks/               # Jupyter notebooks for development and analysis

Usage

1. Set up your environment

   Create and activate the conda environment:

   conda env create -f environment.yml
   conda activate flu-usher
   

   Clone larch, then build and install it from source into the environment, where $CONDA_PREFIX is the path to the active conda environment (e.g., /home/hhaddox/miniforge3/envs/flu-usher/):

   git clone https://github.com/matsengrp/larch.git
   cd larch
   git submodule update --init --recursive
   mkdir build
   cd build
   cmake -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX ..
   make -j16
   make install
   

   Other dependencies are managed per-step via conda environments in envs/, which Snakemake invokes automatically.

2. Configure the pipeline

   Edit config.yaml to:

   • Specify input directories containing your GISAID data
   • List HA subtypes to analyze (e.g., H1, H3, H5, H7, H9)
   • List NA subtypes to analyze (e.g., N1, N2, N9)
   • Set reference accession numbers for each segment-subtype combination
   • Adjust filtering thresholds for sequence curation (max_frac_gaps, max_frac_ambig)
   • Set number of randomizations for tree building (n_randomizations, default: 10)
   • Set desired number of threads
   • Specify geographic groups to extract (geographic_groups_to_extract)
   • Optionally specify rerooting nodes for final trees (reroot)

3. Prepare your GISAID data

   The pipeline expects two types of files with input data downloaded from GISAID:

   • One or more FASTA files containing nucleotide sequences
     • FASTA sequence ID format: EPI|SEGMENT|NAME|EPI_ISL|SUBTYPE
   • One or more XLS files with metadata for the above sequences

   The pipeline will automatically:

   • Aggregate all FASTA and metadata files from each input directory
   • Parse sequences by segment
   • Group HA and NA segments by subtype (defined by the subtype column in the metadata file)
   • Combine all subtypes for internal segments (PB2, PB1, PA, NP, MP, NS)

4. Run the pipeline

   snakemake --cores <number_of_cores> --use-conda
   

5. Output

   The pipeline organizes outputs by segment and subtype:

   For HA and NA segments (e.g., results/HA/H5/ or results/NA/N1/):

   • raw_sequences.fasta.xz: Parsed sequences for this segment-subtype
   • reference/: Reference sequence and Nextclade dataset files
   • msa.fasta.xz: Multiple sequence alignment from Nextclade
   • msa.tsv.xz: Nextclade alignment metadata
   • curated_msa.fasta.xz: Quality-filtered alignment (e.g., gaps < 3%, ambiguities < 0%)
   • unaligned_coding_seqs/: Directory containing per-gene unaligned coding sequences extracted from curated sequences
   • curated_reference.fasta/txt/gff/gtf: Reference files matching the curated alignment
   • curated_root.fasta: Root sequence (reference or inferred from rerooting)
   • randomized_{n}/: Directory for each randomized alignment (n = 0, 1, 2, ...)
     • msa.fasta.xz: Randomized alignment
     • msa.vcf: Variant call format file
     • preopt_tree.pb.gz: Initial parsimony tree
     • opt_tree.pb.gz: Optimized tree
     • dag.pb: DAG representation of the tree
   • larch_merged_dag.pb: Merged DAG from all randomizations
   • trimmed_dag.pb: Trimmed DAG with suboptimal trees removed
   • sampled_tree.nh: Newick format tree sampled from trimmed DAG
   • sampled_tree.pb.gz: MAT protobuf of sampled tree
   • final_tree.pb.gz: Final tree (rerooted if specified in config)
   • final_tree.jsonl.gz: Interactive Taxonium visualization file
   • geographic_trees/: Geographic subtree visualizations
     • {geo_group}_samples.txt: Sample list for each geographic group
     • {geo_group}_tree.pb.gz: Extracted subtree for each geographic group
     • {geo_group}_tree.jsonl.gz: Taxonium visualization for each geographic group
   • host_ancestral/: Per-node host inference (PastML / DOWNPASS)
     • combined_ancestral_states.tab: Tab-separated file of inferred host_group for every node (leaves + internals); the node column joins to internal node IDs in final_tree.pb.gz. Ambiguous internals may appear on multiple rows (one per equally-parsimonious state).
     • host_tree.html: Interactive PastML visualization of the ancestral reconstruction.
     • named.tree_final_tree.nwk: PastML's labeled-tree output (identical names to final_tree.nwk here).
   • subtype_ancestral/: Per-node subtype inference (PastML / DOWNPASS)
     • combined_ancestral_states.tab: Tab-separated file of inferred subtype (H*N* form, normalized from the raw GISAID subtype column) for every node (leaves + internals); same join + ambiguity semantics as host_ancestral/. On HA per-subtype trees the H letter is fixed across all tips and only the inferred N partner varies (analogously for NA trees); on internal-segment trees neither letter is constrained, so the full H*N* can vary along the tree.
     • subtype_tree.html: Interactive PastML visualization of the ancestral reconstruction.

   For the other segments (e.g., results/PB2/all/ or results/NP/all/):

   • Same outputs as above, but combining all influenza subtypes

   Global outputs:

   • results/combined_metadata.csv: Aggregated metadata from all input files
   • results/combined_metadata_augmented.csv: Metadata with host_group, geographic_group, and temporal_group columns added
   • results/host_annotation.csv: 2-column (isolate_id, host_group) annotation table for PastML
   • results/subtype_annotation.csv: 2-column (isolate_id, subtype) annotation table for PastML (subtype normalized to H*N* from the raw subtype column)
   • results/input_data_md5sums.txt: Provenance manifest with md5 hashes for every input FASTA / XLS file under the configured input_dirs, one file per line in standard md5sum format (<hash> <path>)
   • results/notebooks/: Executed analysis notebooks
     • analyze_metadata.html: Metadata analysis report
     • analyze_alignments.html: Alignment statistics report

Pipeline Steps

1. Parse GISAID data (parse_gisaid_data.py):

   • Aggregates sequences from multiple input directories
   • Splits by segment and subtype
   • Creates unified metadata file

2. Download reference (download_ref_seq.py):

   • Fetches appropriate reference sequence for each segment-subtype
   • Creates Nextclade dataset configuration

3. Align sequences (Nextclade):

   • Performs codon-aware alignment to reference
   • Includes reference in the output alignment
   • Generates alignment metadata TSV

4. Curate alignment and extract coding sequences (curate_and_extract_coding_seqs.py):

   • Filters sequences by quality (gaps, ambiguities, terminal gaps)
   • Extracts only coding regions based on GFF annotations
   • Sanitizes sequence IDs for UShER compatibility
   • Filters duplicate sequences
   • Validates CDS for all genes (sequences must pass ALL validations)
   • Extracts per-gene unaligned coding sequences from curated alignment
   • Preserves original nucleotides removed during alignment

5. Randomize alignment (randomize_alignment.py):

   • Creates multiple randomized versions of the alignment (keeping reference at top)
   • Each randomization uses a different seed for variation in tree building

6. Create VCF (faToVcf):

   • Converts each randomized FASTA alignment to variant format
   • Includes reference and handles ambiguous nucleotides

7. Build initial tree (usher-sampled):

   • Creates parsimony-based phylogenetic tree for each randomization
   • Uses empty starting tree and builds incrementally

8. Optimize tree (matOptimize):

   • Refines tree topology to minimize parsimony score for each randomization
   • Multiple optimization rounds to improve tree quality

9. Convert to DAG (larch-usher):

   • Converts each optimized tree to a DAG (Directed Acyclic Graph) representation
   • Allows for representing multiple equally parsimonious tree topologies

10. Merge DAGs (larch-dagutil):

    • Combines DAGs from all randomizations into a single merged DAG
    • Trims redundant structures during merge

11. Trim DAG (trim_dag.py):

    • Removes suboptimal trees from the merged DAG
    • Retains only the most parsimonious tree topologies

12. Sample tree from DAG (convert_DAG_protobuf_to_newick_samples.py):

    • Samples a representative tree from the trimmed DAG
    • Outputs in Newick format

13. Create MAT protobuf (matOptimize):

    • Converts the sampled Newick tree to MAT protobuf format
    • Required for downstream matUtils operations

14. Reroot tree (matUtils extract):

    • Reroots the tree at a specified node if configured
    • Otherwise creates a symlink to the sampled tree

15. Create root sequence (extract_root_sequence.py or symlink):

    • If rerooted: Infers root sequence from tree mutations
    • If not rerooted: Uses reference sequence as root

16. Augment metadata (augment_metadata.py):

    • Adds host_group column (human, avian, swine, bovine, other)
    • Adds geographic_group column (north_america, europe, asia, other)
    • Adds temporal_group column (early, late, unknown based on global median date)

17. Extract geographic subtrees (matUtils extract):

    • Creates separate subtrees for each configured geographic region (filter on the augmented geographic_group column)
    • Each subtree includes the root sequence plus matching samples

18. Infer per-node host states (PastML, prepare_host_annotation.py):

    • Builds a 2-column annotation table from combined_metadata_augmented.csv mapping isolate_id → host_group
    • Runs PastML with the DOWNPASS parsimony method on final_tree.nwk to reconstruct host_group for every internal node
    • Inputs: final_tree.nwk, combined_metadata_augmented.csv
    • Outputs: host_ancestral/combined_ancestral_states.tab (per-node host_group; node IDs match final_tree.pb.gz), host_ancestral/host_tree.html (interactive viz), host_ancestral/named.tree_final_tree.nwk
    • Ambiguous internals appear on multiple rows of combined_ancestral_states.tab (one per equally-parsimonious state)

19. Infer per-node subtype states (PastML, prepare_subtype_annotation.py):

    • Builds a 2-column annotation table from combined_metadata_augmented.csv mapping isolate_id → subtype, normalizing the raw GISAID subtype ("A / H5N1") to H*N* form ("H5N1") inside the prep script (unknown for non-matches)
    • Runs PastML with DOWNPASS on final_tree.nwk to reconstruct subtype for every internal node
    • Inputs: final_tree.nwk, combined_metadata_augmented.csv
    • Outputs: subtype_ancestral/combined_ancestral_states.tab, subtype_ancestral/subtype_tree.html
    • On HA per-subtype trees the H letter is fixed across all tips and only the inferred N partner varies (analogously for NA trees); on internal-segment trees (PB2/all, PB1/all, PA/all, NP/all, MP/all, NS/all) neither letter is constrained, so the full H*N* can vary along the tree and ancestral subtype tracks reassortment events

20. Create visualizations (usher_to_taxonium):

    • Converts final tree and all subtrees to Taxonium format
    • Incorporates metadata (including host, geographic, and temporal groups) for interactive exploration

21. Execute analysis notebooks (jupyter nbconvert):

    • Runs analysis notebooks after all pipeline outputs are complete
    • Generates HTML reports in results/notebooks/
    • Includes metadata analysis and alignment statistics

22. Record input data md5 sums (md5sum):

    • Walks every .fasta and .xls file under the configured input_dirs and writes their md5 hashes to results/input_data_md5sums.txt as a provenance manifest
    • Files are listed as explicit Snakemake inputs, so the manifest is regenerated whenever any input data file changes

Requirements

The main conda environment (environment.yml) provides Snakemake and the build dependencies needed to compile larch. Per-step dependencies are managed via separate conda environments in envs/, which Snakemake invokes automatically using the --use-conda flag.