# Bio Applied Phylodynamics

> 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).

- Skill: `pavel-kravchenko/bio-applied-phylodynamics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-phylodynamics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-phylodynamics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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).

```bash
# 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`.

```python
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