bio-phylo-bayesian-inference

What this skill does

## Version Compatibility

Reference examples tested with: MrBayes 3.2.7+, BEAST2 2.7+, Tracer 1.7+, RevBayes 1.2+, PhyloBayes MPI 1.9+, RWTY (R package)

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `mb --version`, `beast -version`, `rb --version`, `pb --version`
- Python: `pip show biopython` then `help(module.function)` to check signatures
- R: `packageVersion('rwty')` then `?function_name` to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Bayesian Phylogenetic Inference

**"Run a Bayesian phylogenetic analysis"** -> Infer posterior distribution of trees and parameters via MCMC sampling, producing posterior probability support and enabling formal model comparison.
- CLI: `mb` (MrBayes), `beast` (BEAST2), `rb` (RevBayes), `pb`/`bpcomp` (PhyloBayes)
- Python: BioPython `Bio.Phylo` for parsing output trees; `arviz`/`pandas` for trace diagnostics

## When to Use Bayesian vs ML

| Factor | ML (IQ-TREE/RAxML-NG) | MrBayes | BEAST2 | RevBayes | PhyloBayes |
|--------|------------------------|---------|--------|----------|------------|
| Primary use | Topology, hypothesis testing | Posterior probabilities on topology | Divergence times, demographics, tip-dating | Non-standard models, full graphical model control | Deep phylogenies with compositional heterogeneity |
| Speed | Fast | Moderate | Slow | Moderate | Very slow |
| Ease of use | Easy | Easy | Moderate (BEAUti GUI) | Steep learning curve | Moderate |
| Model flexibility | Good | Good | Extensive (packages) | Maximum (scripted) | CAT/CAT-GTR |
| Default choice | Yes, for topology | When posterior probabilities needed | When time calibration needed (see divergence-dating) | When no existing model fits | When LBA from model misspecification suspected |
| Multi-run | Automatic | 2 runs x 4 chains default | Must run independently | Manual | Must run independently |

**Default recommendation:** Start with ML (modern-tree-inference). Move to Bayesian when posterior probabilities, divergence times, or site-heterogeneous models are specifically needed.

## MrBayes

### Basic Analysis

**Goal:** Obtain a Bayesian phylogeny with posterior probability support values.

**Approach:** Load a Nexus alignment into MrBayes, set the substitution model and priors, run two independent MCMC analyses with Metropolis-coupled chains (MC3), then summarize parameters and trees after discarding burn-in.

```
execute alignment.nex
lset nst=6 rates=invgamma
prset brlenspr=unconstrained:exp(10.0)
mcmc ngen=1000000 samplefreq=100 nchains=4 nruns=2
sump
sumt
```

Key settings:
- `nst=6 rates=invgamma`: GTR+I+G model. Use `nst=2` for HKY, `nst=1` for JC/F81
- `nchains=4 nruns=2`: Two independent runs each with 4 Metropolis-coupled chains (1 cold + 3 heated). This is the default and should not be reduced
- `samplefreq=100`: Sample every 100 generations. Adjust so total samples = ngen/samplefreq is at least 10,000
- `sump`: Summarizes parameter traces, reports ESS and PSRF
- `sumt`: Summarizes trees, builds majority-rule consensus with posterior probabilities

### MC3 (Metropolis-Coupled MCMC)

MrBayes runs Metropolis coupling by default: one cold chain explores the posterior while heated chains explore broader landscape and swap states with the cold chain. This improves mixing for multimodal posteriors. The acceptance rate for swaps between adjacent chains should be 20-70%; adjust `temp=0.2` (default) if swaps are too rare or too frequent.

### Mixed Models and Partitions

```
charset gene1 = 1-300
charset gene2 = 301-600
partition by_gene = 2: gene1, gene2
set partition = by_gene
lset applyto=(1) nst=6 rates=invgamma
lset applyto=(2) nst=2 rates=gamma
unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all)
prset applyto=(all) ratepr=variable
```

### Reversible-Jump Model Selection

MrBayes can average over all 203 time-reversible substitution models during MCMC:

```
lset nst=mixed rates=invgamma
```

This integrates over model uncertainty rather than fixing a single model. Preferred when model choice is uncertain.

## BEAST2

### Workflow

1. Prepare alignment in BEAUti (generates XML configuration)
2. Set substitution model, clock model, and tree prior in BEAUti
3. Run BEAST2 from command line
4. Check convergence in Tracer
5. Summarize trees with TreeAnnotator
6. Visualize in FigTree or equivalent

```bash
beast -threads 4 -seed 12345 analysis.xml
treeannotator -burnin 10 -heights median analysis.trees consensus.tree
```

### bModelTest for Bayesian Model Averaging

The bModelTest package averages over all 203 time-reversible nucleotide substitution models during MCMC, eliminating the need for separate model selection:

1. Install via BEAUti Package Manager
2. In BEAUti: Site Model tab -> select "bModelTest"
3. The analysis samples across models weighted by their posterior probability

This is analogous to MrBayes `nst=mixed` but implemented as a BEAST2 package.

### Running Multiple Independent Analyses

BEAST2 does NOT run multiple chains by default (unlike MrBayes). Always run at least two independent analyses with different seeds:

```bash
beast -threads 4 -seed 12345 analysis.xml
beast -threads 4 -seed 67890 analysis.xml
```

Then compare traces in Tracer to confirm both runs converge to the same posterior.

## MCMC Convergence Diagnostics

This is the most critical aspect of any Bayesian phylogenetic analysis. An unconverged analysis produces meaningless results regardless of the model or data.

### Effective Sample Size (ESS)

ESS measures the number of effectively independent samples after accounting for autocorrelation.

| ESS | Interpretation |
|-----|---------------|
| < 100 | Insufficient; results unreliable |
| 100-199 | Marginal; increase run length |
| >= 200 | Adequate (Tracer default threshold) |
| >= 625 | Conservative target for precise credible intervals |

**Critical rule:** If ESS < 200 for any parameter, run the chain longer. Do NOT simply increase thinning (samplefreq). Thinning discards information and does not improve ESS; it is only justified to reduce file size.

### Trace Plot Interpretation

A well-converged trace plot resembles a "hairy caterpillar," with rapid oscillation around a stable mean and no visible trend. Warning signs:
- Upward/downward trends: chain has not reached stationarity
- Steps or plateaus: chain stuck in local optima
- Different levels between runs: runs exploring different regions of parameter space
- Periodic oscillations: poor mixing, possibly due to correlated parameters

### PSRF (Potential Scale Reduction Factor)

MrBayes reports PSRF (Gelman-Rubin diagnostic) comparing variance within and between runs:

| PSRF | Interpretation |
|------|---------------|
| < 1.01 | Good convergence |
| 1.01-1.05 | Acceptable but monitor |
| > 1.05 | Not converged; run longer or check model |

### Topological Convergence (RWTY)

Convergence of continuous parameters (likelihood, branch lengths) does NOT guarantee topological convergence. The tree topology may still be poorly sampled even when ESS for likelihood is high.

**Goal:** Assess whether MCMC chains have adequately sampled tree space.

**Approach:** Use the RWTY R package to compute topological ESS and visualize treespace using multidimensional scaling of Robinson-Foulds distances between sampled trees.

```r
library(rwty)
trees_run1 <- load.trees('run1.t', type='nexus')
trees_run2 <- load.trees('run2.t', type='nexus')
rwty_output <- analyze.rwty(list(run1=trees_run1, run2=trees_run2), burnin=25)
makeplot.topology(rwty_output)
makeplot.treespace(rwty_output)
```

What to look for:
- Topological ESS should be >= 200 (same threshold as continuous ESS)
- Treespace plot: samples from independent runs should overlap in a single cluster
- If runs occupy different regions of treespace, they have not converged on topology
- Split frequency standard deviation