bio-applied-phylodynamics

name	bio-applied-phylodynamics
description	Build time-scaled phylogenies with TreeTime/Augur, validate clock via root-to-tip regression, interpret BEAST2 skyline plots and phylogeography. Use when dating an outbreak, estimating TMRCA, R0, or Ne(t).
tool_type	python
primary_tool	TreeTime

Viral Phylodynamics and Molecular Epidemiology

When to Use

• Dating a viral outbreak or estimating time to most recent common ancestor (TMRCA) from an alignment + sampling dates
• Validating whether a sequence set has clock-like signal before running BEAST2 or augur refine
• Building a Nextstrain (Augur/Auspice) pipeline from FASTA + metadata to a time-scaled tree
• Interpreting a Bayesian Skyline Plot (Ne(t)) or discrete/continuous phylogeography output from BEAST2
• Estimating R0 or growth rate from a coalescent tree or genomic surveillance data

Version Compatibility

TreeTime ≥0.11, Nextstrain Augur ≥24 (Python ≥3.10), BEAST2 ≥2.7 (+Tracer ≥1.7.2 for convergence), IQ-TREE2 ≥2.3 with LSD2, scipy ≥1.11, Biopython ≥1.83.

Prerequisites

• pip install treetime biopython scipy pandas matplotlib and conda install -c bioconda nextstrain-augur (or the nextstrain CLI)
• A ML tree (see bio-phylogenetics-modern-tree-inference) and a MAFFT/MUSCLE alignment as inputs
• Understand basic tree I/O (bio-phylogenetics-tree-io) — Newick parsing, tip/node metadata

Root-to-Tip Regression (Clock Signal Validation)

Goal: confirm the sequence set evolves clock-like before committing to a Bayesian timetree analysis. Approach: regress root-to-tip genetic distance against sampling date; R² > 0.70 and a positive slope indicate usable clock signal. Always inspect the residuals — a decent R² can still hide recombination or selection.

import numpy as np
from scipy import stats

def clock_signal(sampling_dates, root_to_tip, genome_len=29903):
    """Validate molecular clock signal via root-to-tip regression.

    Args:
        sampling_dates: array of decimal-year collection dates.
        root_to_tip: array of genetic distances from the tree root (subs/site).
        genome_len: genome length used to convert rate to expected SNPs.

    Returns:
        dict with slope (subs/site/year), r_squared, p-value, and an
        estimated SNPs-per-2-weeks figure for sanity-checking against
        known virus rates (e.g. ~1-2 SNPs/2wk for SARS-CoV-2).
    """
    slope, intercept, r, p, se = stats.linregress(sampling_dates, root_to_tip)
    return {
        "rate_subs_per_site_per_year": slope,
        "r_squared": r ** 2,
        "p_value": p,
        "snps_per_2weeks": slope * (14 / 365) * genome_len,
        "clock_signal_ok": (r ** 2 > 0.70) and (slope > 0),
    }

# sampling_dates: decimal years per tip; root_to_tip: from ete3/Bio.Phylo tree traversal
result = clock_signal(np.array([2020.1, 2020.4, 2020.9, 2021.3]),
                       np.array([0.0001, 0.0004, 0.0009, 0.0013]))
print(result)

Time-Scaled Phylogenies (Nextstrain Augur / TreeTime)

Goal: convert a divergence tree into a dated timetree with TMRCA estimates for downstream Auspice visualization or transmission-cluster analysis. Approach: filter/subsample to avoid geographic and temporal oversampling, align, build an ML tree, then time-scale with augur refine (a TreeTime wrapper).

# 1. Filter and subsample to reduce sampling bias
augur filter --sequences seqs.fasta --metadata metadata.tsv \
    --min-date 2020-01-01 --subsample-max-sequences 300 \
    --output filtered.fasta

# 2. Align to reference
augur align --sequences filtered.fasta \
    --reference-sequence reference.gb --output aligned.fasta

# 3. Build ML divergence tree (IQ-TREE under the hood)
augur tree --alignment aligned.fasta --output tree_raw.nwk

# 4. Time-scale with TreeTime (coalescent prior, relaxed by default)
augur refine --tree tree_raw.nwk --alignment aligned.fasta \
    --metadata metadata.tsv --timetree --coalescent opt \
    --output-tree timetree.nwk --output-node-data branch_lengths.json

# 5. Export for Auspice
augur export v2 --tree timetree.nwk --node-data branch_lengths.json \
    --output auspice.json

Use a strict clock only for short time spans within one host species; switch to a relaxed (UCLN) clock in BEAST2 (--coalescent opt in Augur approximates this) whenever the dataset spans multiple years or hosts, since substitution rates vary across lineages.

Effective Population Size and R0

Goal: turn a BEAST2 Bayesian Skyline Plot or coalescent growth rate into an epidemiological signal. Approach: Ne(t) tracks coalescent rate, not census population size; a rising Ne means exponential growth, a plateau means endemic equilibrium, a sharp decline means an intervention or immunity effect. R0 can be approximated from the exponential-growth-phase coalescent rate r and the serial interval T.

import numpy as np

def r0_from_growth_rate(growth_rate_per_year, serial_interval_days):
    """Estimate R0 from an exponential-phase coalescent/molecular-clock growth rate.

    R0 = 1 + r * T_generation (valid only in the early exponential phase;
    breaks down once susceptible depletion or interventions occur).
    """
    serial_interval_years = serial_interval_days / 365
    return 1 + growth_rate_per_year * serial_interval_years

r0 = r0_from_growth_rate(growth_rate_per_year=45, serial_interval_days=5.2)
print(f"R0 estimate: {r0:.2f}")  # sanity-check against 2-4 for most respiratory viruses

Ne Trend	Epidemiological Meaning
Rising Ne	Exponential epidemic growth
Plateau	Endemic equilibrium
Sharp decline	Intervention, seasonal end, or population immunity

Phylogeography

Method	Tool	Use Case
Discrete (DTA)	BEAST2	Location states per node, directional migration rates
Continuous	BEAST2	Lat/lon Brownian motion, diffusion coefficients
Parsimony/ML	Nextstrain/Augur	Fast, visual — animated Auspice maps

Key metrics: migration rate matrix (transitions/year), Bayes factor (evidence for a route), geographic diffusion coefficient (km²/year, continuous model only).

Pitfalls

• R² > 0.70 is necessary but not sufficient: strong root-to-tip regression can still arise from recombination or selection — always inspect residuals, not just R².
• Strict clock on long time scales: rates vary across lineages; use a relaxed clock (UCLN) in BEAST2, or --coalescent opt in Augur, for multi-year or cross-host datasets.
• Ne ≠ actual viral population size: BSP estimates the effective population size from coalescent timing; heavy sampling density inflates apparent Ne independent of true prevalence.
• Sampling bias distorts phylogeography: oversampling one location fabricates spurious migration routes toward it — always subsample with augur filter first.
• BEAST2 MCMC convergence: require ESS > 200 for every parameter in Tracer before trusting a posterior; insufficient chain length is the most common failure mode.
• Recombination breaks the clock assumption: HIV and coronaviruses recombine — screen with RDP4 or ClonalFrameML before fitting molecular clock models.
• R0 = 1 + r·T only holds in the early exponential phase: once susceptible depletion or interventions kick in, use birth-death models (BDEI/BDSIR) instead.

See Also

• bio-phylogenetics-modern-tree-inference — build the ML divergence tree that feeds augur refine/TreeTime
• bio-phylogenetics-tree-io — Newick/Nexus parsing for tree objects used in clock regressions
• bio-epidemiological-genomics-transmission-inference — cluster and pair transmission events from timetrees
• bio-workflows-outbreak-pipeline — end-to-end outbreak genomics pipeline that wraps this analysis