bio-phylo-divergence-dating

What this skill does

## Version Compatibility

Reference examples tested with: BEAST2 2.7+, MCMCTree (PAML 4.10+), TreePL 1.0+, Python 3.9+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `beast -version` then `beast -help` to confirm flags
- CLI: `mcmctree --help` or check PAML documentation for parameter names
- CLI: `treePL` to confirm installation
- Python: `pip show biopython` then `help(module.function)` to check signatures

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

# Divergence Time Estimation

**"Estimate when lineages diverged"** -> Combine phylogenetic trees with fossil calibrations and molecular clock models to infer absolute divergence times.
- CLI: BEAUti + BEAST2 for moderate data with complex models
- CLI: MCMCTree (PAML) for genome-scale data with approximate likelihood
- CLI: TreePL for very large trees with penalized likelihood

## Clock Model Selection

| Model | Assumption | When to Use |
|-------|------------|-------------|
| Strict clock | Constant rate across all lineages | Only when a molecular clock test does not reject clocklike behavior; rarely appropriate for inter-species data |
| UCLN (uncorrelated lognormal) | Branch rates drawn independently from lognormal | **Default choice**; rates vary across lineages without parent-child correlation |
| Uncorrelated exponential | Rates from exponential distribution | More extreme rate variation; when UCLN coefficient of variation is very high |
| Autocorrelated (ACLN) | Rates correlated with parent branch | When evolutionary rates change gradually (e.g., body size evolution driving rate changes) |
| Random local clock | Few distinct rate regimes | When a small number of rate shifts are suspected rather than continuous variation |

### Clock Model Diagnostics

Under UCLN in BEAST2, the posterior of `ucld.stdev` measures clocklikeness:

| ucld.stdev | Interpretation |
|------------|----------------|
| ~0 | Data is clocklike; strict clock may suffice |
| 0.1-0.5 | Moderate rate variation; UCLN appropriate |
| >1.0 | Extreme rate variation; consider model adequacy or data partitioning |

In MCMCTree, the equivalent diagnostic is `sigma2` under clock=2 (independent rates). Near-zero values suggest a strict clock; very large values indicate substantial rate heterogeneity.

## Fossil Calibration Strategies

### Node Calibration (Traditional)

Prior distributions constrain the age of internal nodes based on fossil evidence:

| Distribution | When to Use | Parameterization |
|--------------|-------------|------------------|
| Lognormal | Most common; fossil sets minimum, soft maximum | Offset = fossil minimum age, mean and stdev control tail |
| Exponential | Strong belief divergence is close to fossil age | Rate parameter controls decay from minimum |
| Uniform | Only minimum and maximum known | Dangerous with narrow bounds; truncation distorts posteriors |
| Normal | Independently known divergence age (e.g., biogeographic event) | Mean = best estimate, stdev = uncertainty |

### Critical Calibration Pitfalls

These are the most common sources of error in divergence dating, and each can silently produce incorrect results:

1. **Prior interaction**: Multiple calibration priors interact through the tree prior (birth-death or Yule). The effective prior on a node may differ substantially from the specified prior. ALWAYS run the analysis sampling from the prior only (no data) first and compare effective priors to specified priors. If they differ, adjust calibration distributions.

2. **Truncation effects**: When a calibration density extends beyond the root age, the density gets truncated and renormalized, distorting the effective prior. This is especially problematic for uniform and lognormal priors near the root.

3. **Secondary calibrations**: Using divergence time estimates from a previous study as calibrations compounds uncertainty. The original confidence intervals become point-like constraints in the new analysis. When secondary calibrations are unavoidable, use substantially wider priors than the original confidence intervals.

4. **Fossil placement**: The single largest source of calibration error is incorrect phylogenetic placement of fossils. A fossil assigned to the wrong node can bias all downstream estimates. Verify fossil placement against current morphological and phylogenetic evidence.

5. **Minimum vs maximum confusion**: Fossils provide minimum age constraints (the lineage is at least as old as the fossil). Maximum ages require additional geological or biogeographic arguments and are harder to justify.

## Tip-Dating and Fossilized Birth-Death (FBD)

### FBD vs Node Calibration

| Factor | Node Calibration | Fossilized Birth-Death |
|--------|-----------------|----------------------|
| Fossil usage | Few well-placed fossils | All available fossils as tips |
| Prior distributions | Ad hoc (lognormal, uniform) | Coherent birth-death framework |
| Morphological data | Not required | Beneficial but optional |
| Extinct taxa | Excluded | Included as tips |
| Setup complexity | Lower | Higher (morphological matrix, FBD parameters) |
| Software | BEAST2, MCMCTree | BEAST2 (SA package) |

**When to prefer FBD**: Rich fossil record available, morphological data matrix exists, wanting to include extinct taxa, or concerned about ad hoc prior distributions.

**When to prefer node calibration**: Few well-placed fossils, no morphological data, simpler setup needed, or working with MCMCTree.

## Software Decision Framework

| Scenario | Recommended Tool | Rationale |
|----------|-----------------|-----------|
| Genome-scale data (>100 loci) | MCMCTree | Approximate likelihood is 100-1000x faster than exact |
| Moderate data, complex models | BEAST2 | Richest model ecosystem, FBD support |
| Very large tree (>1000 taxa) | TreePL | Penalized likelihood scales to 10,000+ taxa |
| Fossils as tips needed | BEAST2 + SA package | Only tool with FBD tip-dating |
| Mixture models + dating | IQ2MC (IQ-TREE 3 + MCMCTree) | Combines mixture model fit with Bayesian dating |
| Quick exploratory dating | TreePL | Fast; correlates well with Bayesian estimates |

## BEAST2 Dating Workflow

**Goal:** Estimate divergence times with calibrated relaxed clock in BEAST2.

**Approach:** Prepare XML via BEAUti, run BEAST2, verify convergence, summarize results.

### Step-by-Step Workflow

1. **Prepare alignment** in BEAUti: import FASTA/NEXUS, set data type
2. **Set substitution model**: use bModelTest for model averaging (recommended) or fix a model
3. **Set clock model**: UCLN is the default choice; set initial clock rate based on prior knowledge
4. **Set tree prior**: Birth-Death for extant-only data, FBD when including fossils as tips
5. **Add fossil calibrations**: specify prior distributions with justified parameters for each calibrated node
6. **Run from prior first**: sample from the prior only (uncheck "use data" or set `sampleFromPrior="true"`) to verify effective priors match expectations
7. **Run full analysis**: at least 2 independent runs with different random seeds
8. **Check convergence**: Tracer: all ESS values >= 200; compare posterior traces between runs
9. **Check topological convergence**: use RWTY or similar
10. **Combine and summarize**: LogCombiner (remove burn-in from each run), then TreeAnnotator for MCC tree

### BEAST2 Clock Rate Prior

**Goal:** Set an informable prior on the clock rate to avoid poor mixing.

**Approach:** Estimate the expected substitution rate from prior information and set the clock rate prior accordingly.

```bash
# Rough clock rate estimation:
# rate = divergence / (2 * time)
# For mammals: ~0.01 substitutions/site/Myr for nuclear genes
# For mitochondrial: ~0.02 substitutions/site/Myr
# Set lognormal prior with M = ln(expected_rate), S = 0.5-1.0
```

### Key BEAST2 Packages

| Package | Purpose |
|---------|---------|
| SA (s