Bio Applied Phylodynamics

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

Phylogeography

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