Bio Comparative Genomics Ancestral Reconstruction

Version Compatibility

Reference examples tested with: PAML 4.10.7+, IQ-TREE 2.3.6+, GRASP 2024+ (web/CLI), FastML 3.11+, RevBayes 1.2.4+, BayesTraits V4.1+, R 4.4+, ape 5.8+, phytools 2.3+, corHMM 2.9+, geiger 2.0.11+, phangorn 2.12+, RERconverge 0.3.0+, BioPython 1.84+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('corHMM') then ?ancRECON, ?make.simmap, ?ace
• CLI: codeml (no --version; check paml -h or examine Phylip.tre example), iqtree2 --version, revbayes --version
• Python: pip show biopython; check Bio.Phylo.PAML.codeml API

If code throws AttributeError, ImportError, missing slot errors on R S4 objects, or PAML mlc parsing failures, introspect the installed package (? in R, help() in Python) and adapt the example rather than retrying. PAML output formats are stable across 4.9 -> 4.10; IQ-TREE's --ancestral flag replaced -asr in v2.0+.

Ancestral State Reconstruction

"What did this gene / trait / genome look like at an internal node?" -> Choose the reconstruction framework that matches the data class (sequence / discrete trait / continuous trait / gene content) and the inference question (point estimate vs full posterior; marginal vs joint vs scaled-conditional). The single most common mistake is reconstructing under a site-independent or trait-stationary model when the underlying biology demands a hidden-rate or epistatic model -- the resulting "ancestor" is mathematically optimal under the wrong model and is silently wrong (Beaulieu & O'Meara 2016 Syst Biol 65:583; Boyko & Beaulieu 2021 MEE 12:468).

• Sequence ASR (protein resurrection): PAML codeml RateAncestor=1; IQ-TREE2 --ancestral; GRASP (graph-based, handles indels); FastML (Bayesian)
• Discrete traits: R ape::ace(type='discrete'); corHMM::corHMM() (rate categories); phytools::make.simmap() stochastic mapping; BayesTraits MultiState/Discrete
• Continuous traits: phytools::fastAnc(); phytools::contMap(); geiger::fitContinuous(model='BM'|'OU'|'EB'); RPANDA fit_t_env()
• Ancestral gene content (presence/absence): ape::ace(type='discrete', model='ARD'); Dollo parsimony in phangorn; ALE/GeneRax for full DTL (see [[gene-tree-species-tree-reconciliation]])

Algorithmic Taxonomy

Methodology evolves; verify the latest corHMM / phytools vignettes before locking on a single approach. For continuous-trait macroevolution, consult Cooper 2016 model-adequacy reviews.

Decision Tree by Experimental Scenario

Per-Tool Failure Modes

Long-branch attraction (LBA) at deep nodes

Trigger: Tree with two long terminal branches separated by a short internal branch; mixed amino-acid compositions across taxa.

Mechanism: Site-homogeneous models (LG, WAG, JTT) assume constant amino-acid equilibrium frequencies across the tree. When real compositions differ (e.g. thermophiles vs mesophiles), models underestimate the probability of convergent substitutions at compositionally-constrained sites, producing apparent shared derived states between long branches (Szánthó 2023 Syst Biol 72:767). The reconstructed ancestor at the deep node is biased toward whichever long-branch composition the model favors.

Symptom: Bootstrap support at the contested node remains high under site-homogeneous models but collapses under CAT-GTR / CAT-PMSF. Posterior predictive checks for compositional homogeneity reject the model (Foster 2004 Syst Biol 53:485).

Fix: Move to PhyloBayes-MPI with CAT-GTR or IQ-TREE2 with CAT-PMSF (-m LG+C60+F+R then -ft <tree> for posterior mean site frequencies). For ASR specifically, use ancestral reconstruction only when the model passes compositional adequacy. Slow-fast site removal (recoded amino acids; Susko & Roger 2007 MBE 24:2139) is an alternative.

Epistasis breaking site-independent ASR

Trigger: Multiple sites in the same protein evolve under coupled constraints (compensatory pairs in RNA secondary structure; buried-residue covariance; allosteric networks).

Mechanism: ML/Bayesian ASR assumes sites are independent given the tree; the joint ancestral sequence is the product of per-site posteriors. Real proteins evolve through compensatory substitutions where a destabilizing mutation at site i is compensated by a mutation at site j (Pollock 2012 PNAS 109:E1352; Shah 2015 PNAS 112:E3226). The ML ancestral sequence can contain a never-tested combination of states.

Symptom: The reconstructed protein fails to fold or is non-functional when expressed; positions flagged ambiguous (P < 0.9) are non-random and cluster in 3D space when mapped to structure.

Fix: Use GRASP indel-aware reconstruction; design 4-8 alternative constructs varying ambiguous positions (P < 0.9), prioritizing residues that are structurally coupled to high-confidence ML states; experimentally test each construct; report the range of functional reconstructions, not the single ML sequence. Hochberg & Thornton 2017 Annu Rev Biophys 46:247 review epistasis strategies.

Root placement error cascading to deep ancestors

Trigger: Trees rooted by midpoint, outgroup with very long branch, or --prefix auto-root.

Mechanism: Marginal ASR posteriors at internal nodes depend on the root's position because the root defines the time direction of substitution. A wrong root flips state inferences for deep nodes (especially when ancestral state is asymmetric, e.g. presence -> absence is more common than reverse).

Symptom: Re-rooting the tree changes the inferred ancestral state at the deepest node by > 0.2 posterior probability; STRIDE rooting (Emms 2017 MBE 34:3267) disagrees with outgroup rooting.

Fix: Run ASR over a set of candidate roots; report robust nodes (state invariant) and root-sensitive nodes separately. For phylogenomic-scale data, use STRIDE / MAD rooting (Tria 2017 Nat Eco Evo 1:0193) or ALE-rooting (Williams 2017 PNAS 114:E4602) and document the rooting strategy.

BiSSE false-positive in state-dependent diversification

Trigger: Testing whether a discrete trait influences speciation/extinction using BiSSE (Maddison 2007 Syst Biol 56:701).

Mechanism: BiSSE attributes ALL rate variation to the focal trait. When real diversification heterogeneity is caused by a hidden character correlated with the focal trait, BiSSE reports a spurious significant association (Rabosky-Goldberg 2015 Syst Biol 64:340 -- ~40% Type-I rate at moderate trees).

Symptom: BiSSE LRT highly significant but biological mechanism unclear; HiSSE rejects BiSSE in favor of a hidden-state model with the focal trait neutral.

Fix: Run HiSSE as the required null model (Beaulieu & O'Meara 2016 Syst Biol 65:583). Report BiSSE only if HiSSE-null is rejected. For traits with deep clade structure, use FiSSE (Rabosky & Goldberg 2017 Evolution 71:1432) which is robust to model misspecification by design.

Parsimony vs ML on asymmetric rates (Felsenstein bias)

Trigger: Trait has strongly asymmetric forward vs reverse rates (e.g. gene loss > gene gain).

Mechanism: Parsimony minimizes total changes, implicitly assuming symmetric rates. ML/Bayesian methods estimate the rate matrix from the data and reconstruct accordingly. Under strong asymmetry, parsimony over-reconstructs the rarer state at ancestors (Cunningham 1999 Syst Biol 48:665).

Symptom: Parsimony and ML reconstructions disagree at > 30% of nodes; ML rates fit AIC-better with ARD (all-rates-different) than ER (equal-rates).

Fix: Always run ER vs SYM vs ARD model comparison via ape::ace(model='...') AIC; use ARD when asymmetry is supported. For gene-content evolution, Dollo parsimony (gain rare, loss common) is often the better prior than equal-rates ML.

Continuous-trait BM-only model with non-BM evolution

Trigger: Fitting phytools::fastAnc() (which assumes BM) to a trait with strong directional or stabilizing selection.

Mechanism: fastAnc returns the BM-MLE ancestral state, which is a weighted mean of descendant values with weights from the BM covariance matrix. If the trait evolved under OU (stabilizing), real ancestor values were closer to the optimum than fastAnc returns; if under EB (early burst), real ancestors were more variable than fastAnc returns.

Symptom: BM model fits with geiger::fitContinuous(model='BM') give AIC > 4 above OU or EB; phylogenetic signal Pagel's lambda < 0.5; Blomberg's K significantly < 1 (Blomberg 2003 Evolution 57:717).

Fix: Run fitContinuous with multiple models (BM/OU/EB/lambda/kappa/delta); use the best AIC model's ancestral reconstruction. For OU, use OUwie::ace(); for regime shifts, bayou::bayou.mcmc. Always report Pagel's lambda alongside ancestral estimates as a phylogenetic-signal indicator. Boettiger 2012 Evolution 66:2240 and Cooper 2016 Biol J Linn Soc 118:64 detail model-adequacy testing.

Alignment error propagating to ASR

Trigger: Using MUSCLE or default MAFFT alignment on highly diverged sequences (< 30% identity).

Mechanism: Misaligned columns place non-homologous residues into the same column. ML ASR treats those residues as states of the same character, producing impossible ancestral inferences (Vialle 2018 MBE 35:1783).

Symptom: Ambiguous regions of the alignment correspond to low-confidence ASR sites; gappy columns dominate the low-confidence set; alignment scoring (TCS, Guidance2) marks the same regions as poorly aligned.

Fix: Filter alignment with HmmCleaner (Di Franco 2019 BMC Evol Biol 19:21) or PREQUAL (Whelan et al 2018 Bioinformatics 34:3929) before ASR. Segment-level filtering outperforms block-level filtering (Gblocks, trimAl) for downstream evolutionary inference. For ASR specifically, mask ambiguous columns (treat as missing) rather than removing them, to preserve coordinates.

Quantitative Thresholds

PAML codeml Ancestral Reconstruction

Goal: Reconstruct ancestral codon or protein states at internal nodes using ML under a stationary codon/protein model.

Approach: Build a codon-aware alignment (PRANK or MACSE) -> infer rooted ML tree with the same model intended for ASR -> create codeml control file with RateAncestor=1 -> run codeml -> parse rst for per-site posteriors and ancestral sequences.

'''PAML codeml ancestral reconstruction with site-confidence parsing'''

import subprocess
import re
import os
from collections import defaultdict


def write_codeml_ctl(alignment, tree, out_dir, seqtype='codon', model='M0'):
    '''Write codeml control file for ancestral reconstruction.

    seqtype: 1=codon, 2=aa, 3=codon translated
    For protein resurrection use seqtype=2 with model=3 (empirical+gamma).
    RateAncestor=1 produces rst file with per-site posteriors.
    '''
    if seqtype == 'codon':
        ctl = f'''
      seqfile = {alignment}
     treefile = {tree}
      outfile = {out_dir}/mlc
      runmode = 0
      seqtype = 1
    CodonFreq = 2
        model = 0
      NSsites = 0
        kappa = 2
    fix_omega = 0
        omega = 0.4
 RateAncestor = 1
    cleandata = 0
   getSE = 1
'''
    else:
        ctl = f'''
      seqfile = {alignment}
     treefile = {tree}
      outfile = {out_dir}/mlc
      runmode = 0
      seqtype = 2
        model = 3
   aaRatefile = lg.dat
       alpha = 0.5
        ncatG = 4
 RateAncestor = 1
    cleandata = 0
'''
    ctl_path = os.path.join(out_dir, 'codeml.ctl')
    open(ctl_path, 'w').write(ctl)
    return ctl_path


def parse_rst_posteriors(rst_file):
    '''Parse rst file. Returns per-node, per-site (state, prob) records.
    PAML rst section headers: "Prob distribution at node X" and "List of extant and reconstructed sequences".
    '''
    text = open(rst_file).read()
    node_posteriors = defaultdict(list)

    for block in re.finditer(r'Prob distribution at node (\d+)(.*?)(?=Prob distribution at node|\Z)', text, re.DOTALL):
        node = int(block.group(1))
        for ln in block.group(2).splitlines():
            m = re.match(r'\s*(\d+)\s+\S+:\s*(\w)\((\d\.\d+)\)', ln)
            if m:
                site, state, prob = int(m.group(1)), m.group(2), float(m.group(3))
                node_posteriors[node].append({'site': site, 'state': state, 'prob': prob})
    return node_posteriors


def summarize_node_confidence(node_posteriors, ambig_cutoff=0.8):
    '''Per-node confidence summary for protein-resurrection construct design.'''
    summary = {}
    for node, sites in node_posteriors.items():
        probs = [s['prob'] for s in sites]
        n_high = sum(p >= 0.95 for p in probs)
        n_amb = sum(p < ambig_cutoff for p in probs)
        summary[node] = {
            'mean_post': sum(probs) / len(probs),
            'frac_high': n_high / len(probs),
            'n_ambiguous': n_amb,
            'ambiguous_sites': [s['site'] for s in sites if s['prob'] < ambig_cutoff]
        }
    return summary

IQ-TREE2 Marginal Reconstruction

Goal: Faster ASR with native model selection, for protein/DNA alignments where PAML is too slow.

Approach: Run IQ-TREE with -m TEST to pick model -> --ancestral produces .state table of per-site posteriors; output rooted by -o <outgroup>.

iqtree2 -s alignment.fasta -m MFP --ancestral -o outgroup_taxon -B 1000 -nt AUTO --prefix asr_iqtree
# .state columns: Node | Site | State | p_A | p_C | p_G | p_T  (DNA)
# or            : Node | Site | State | p_A p_R p_N ...        (protein)

import pandas as pd

def load_iqtree_state(state_file):
    '''Returns dict[node] -> DataFrame[site, state, max_post, all_post].'''
    df = pd.read_csv(state_file, sep='\t', comment='#')
    df.columns = [c.strip() for c in df.columns]
    state_cols = [c for c in df.columns if c.startswith('p_')]
    df['max_post'] = df[state_cols].max(axis=1)
    return df.groupby('Node')

Stochastic Mapping (phytools::make.simmap)

Goal: Sample full character histories along branches for a discrete trait; quantify ancestral state probabilities with proper uncertainty.

Approach: Fit Mk rate matrix -> sample nsim simulated maps under the posterior -> summarize state probabilities per node and transitions per branch.

library(phytools)
library(corHMM)

tree <- read.tree('species_tree.nwk')
traits <- read.csv('traits.csv', row.names = 1)
x <- setNames(traits$state, rownames(traits))

# Fit Mk model and pick ARD vs SYM vs ER by AIC
fit_er <- fitMk(tree, x, model = 'ER')
fit_sym <- fitMk(tree, x, model = 'SYM')
fit_ard <- fitMk(tree, x, model = 'ARD')
aic <- sapply(list(fit_er, fit_sym, fit_ard), AIC)
best_model <- c('ER', 'SYM', 'ARD')[which.min(aic)]

# Stochastic mapping under best model, nsim >= 1000 (raise for asymmetric rates)
smaps <- make.simmap(tree, x, model = best_model, nsim = 1000, pi = 'fitzjohn')
node_pp <- summary(smaps)$ace          # posterior probabilities per state per node

# Hidden-rate alternative for clade-rate heterogeneity
hmm_fit <- corHMM(phy = tree, data = data.frame(species = names(x), trait = x),
                   rate.cat = 2, model = 'ARD', node.states = 'marginal')
hmm_fit$states                          # marginal ancestral state matrix (rows: nodes)

pi='fitzjohn' uses the Fitzjohn 2009 Syst Biol 58:595 root prior (root-state prior weighted by the likelihood of the observed tip data given each root state), preferable to pi='estimated' which fixes the prior to the estimated equilibrium distribution and can over-fit, or pi='equal' which can bias toward the rarer state when data are asymmetric.

Continuous-Trait ASR with Model Adequacy

Goal: Reconstruct ancestral values for a continuous trait while honestly reporting which model class the data support.

Approach: Fit BM/OU/EB/lambda models -> AIC compare -> reconstruct under best model -> report Pagel's lambda and Blomberg's K as signal quality.

library(phytools); library(geiger)

trait <- setNames(traits$mass, rownames(traits))
tree_p <- multi2di(tree)               # phytools/geiger require fully bifurcating

models <- c('BM', 'OU', 'EB', 'lambda', 'kappa', 'delta')
fits <- lapply(models, function(m) fitContinuous(tree_p, trait, model = m))
aic_tbl <- data.frame(model = models, AIC = sapply(fits, function(f) f$opt$aic))
best <- aic_tbl$model[which.min(aic_tbl$AIC)]

# Phylogenetic signal
lambda_fit <- phylosig(tree_p, trait, method = 'lambda', test = TRUE)
K_fit <- phylosig(tree_p, trait, method = 'K', test = TRUE)

if (best == 'BM') {
    asr <- fastAnc(tree_p, trait, CI = TRUE)        # BM ML ancestral with 95% CI
} else if (best == 'OU') {
    # Use OUwie for proper OU ancestral with regime
    library(OUwie)
    df <- data.frame(species = names(trait), regime = rep(1, length(trait)), trait = trait)
    asr_ou <- OUwie.anc(OUwie(tree_p, df, model = 'OU1'), data = df)
} else {
    # phytools::contMap can use lambda/EB rescaled tree
    asr_contmap <- contMap(tree_p, trait, method = 'fastAnc', plot = FALSE)
}

Report Pagel's lambda alongside ancestral values: lambda > 0.7 indicates BM-like, signal supports tree-based reconstruction; lambda < 0.3 indicates ecology rather than phylogeny drives variance and ASR is poorly constrained.

GRASP Indel-Aware Sequence ASR (Protein Resurrection Workflow)

Goal: Reconstruct ancestral protein sequence INCLUDING indel states, for experimental resurrection.

Approach: GRASP (Foley 2022 PLoS Comp Biol 18:e1010633) uses a partial-order alignment graph to model indels probabilistically; outputs alternative reconstructions and explicit indel uncertainty.

grasp -aln alignment.fasta -tree species.nwk -out grasp_run --inference joint --threads 8
# Produces:
#   grasp_run/asr.fasta      ancestral sequences at every internal node
#   grasp_run/asr_posterior  per-site, per-state posterior
#   grasp_run/asr_indels     indel posterior per node per gap-block

Construct-design protocol: (1) Take ML sequence as primary construct. (2) For each site with max posterior < 0.8 and a second state with posterior > 0.2, build a single-mutant alternative. (3) For each indel block with posterior < 0.8, build a present and an absent alternative. (4) Order constructs by structural compactness (avoid surface loops first). Typical batch: 4-8 constructs per ancestral node. Hochberg & Thornton 2017 Annu Rev Biophys 46:247 review the operational pipeline.

Reconciliation: When Methods Disagree

Operational rule for publication: Reconstruct under multiple models (at minimum: ER/SYM/ARD for discrete; BM/OU/EB for continuous; site-homogeneous + CAT-PMSF or PhyloBayes for sequence ASR at deep nodes); report only nodes whose state is invariant across models OR explicitly flag model-sensitive nodes. Single-model claims should be downgraded.

Cohort Gotchas

• Polyploid species in continuous-trait analyses: body size, genome size, gene count are confounded with ploidy; assign subgenomes (see [[whole-genome-duplication]]) and treat as separate tips, or use multilabel-tree methods.
• Hybrid taxa break the bifurcating-tree assumption: ape and phytools assume strict bifurcations; for hybrids, use phylogenetic networks (phangorn, RevBayes admixture) or remove hybrid tips before ASR.
• Tip-dated trees from molecular clock require careful root prior: RevBayes / BEAST2 with calibrated tip dates produce trees in absolute time; PAML/IQ-TREE work in relative substitutions. Match the time unit when integrating across tools.
• OrthoFinder species trees from gene-tree summary (STAG): branch lengths are coalescent-units when used for ASR; convert to substitutions via concatenated alignment if downstream tools require it.

Anticipated Reviewer Pushback

Common Errors

Tool Installation Notes

# CLI
conda install -c bioconda paml iqtree
# RevBayes via source build (https://revbayes.github.io/download) or homebrew on macOS
# GRASP from https://github.com/bodenlab/GRASP (Java)
# FastML web at fastml.tau.ac.il or CLI binary

# R
install.packages(c('ape', 'phytools', 'geiger', 'corHMM', 'phangorn', 'OUwie', 'bayou', 'RERconverge'))
remotes::install_github('thej022214/hisse')

# Python
pip install biopython ete4

For protein resurrection workflows, GRASP is the modern standard; for selection-context codon ASR, PAML remains the reference; for trait macroevolution, phytools + corHMM + OUwie are the working tier; for Bayesian rigor with model uncertainty, RevBayes is required.

References

• Yang Z et al 1995 Genetics 141:1641 (marginal ASR likelihood framework)
• Pupko T et al 2000 MBE 17:890 (joint ASR efficient algorithm)
• Nielsen R 2002 Syst Biol 51:729 (stochastic mapping)
• Huelsenbeck JP et al 2003 Syst Biol 52:131 (Bayesian stochastic mapping)
• Felsenstein J 1985 Am Nat 125:1 (phylogenetic independent contrasts)
• Felsenstein J 2012 Am Nat 179:145 (threshold model)
• Pagel M 1999 Nature 401:877 (lambda phylogenetic signal)
• Blomberg SP et al 2003 Evolution 57:717 (K statistic)
• Beaulieu JM et al 2013 Syst Biol 62:725 (hidden Markov state-rate decoupling)
• Beaulieu JM & O'Meara BC 2016 Syst Biol 65:583 (HiSSE)
• Rabosky DL & Goldberg EE 2015 Syst Biol 64:340 (BiSSE Type-I rates)
• Boyko JD & Beaulieu JM 2021 MEE 12:468 (generalized HMM corHMM)
• Bollback JP 2006 BMC Bioinf 7:88 (SIMMAP)
• Pollock DD et al 2012 PNAS 109:E1352 (compensatory epistasis)
• Shah P et al 2015 PNAS 112:E3226 (contingency and entrenchment epistasis)
• Hochberg GKA & Thornton JW 2017 Annu Rev Biophys 46:247 (ASR for protein resurrection)
• Foley G et al 2022 PLoS Comp Biol 18:e1010633 (GRASP indel-aware ASR)
• Szánthó LL et al 2023 Syst Biol 72(4):767-780 (compositional LBA; CAT-PMSF) -- DOI 10.1093/sysbio/syad013
• Boettiger C et al 2012 Evolution 66:2240 (model adequacy for continuous-trait macroevolution)
• Cooper N et al 2016 Biol J Linn Soc 118:64 (cautionary note OU)
• Cunningham CW 1999 Syst Biol 48:665 (asymmetric rate parsimony bias)
• Felsenstein J 1978 Syst Zool 27:401 (long branch attraction)
• Maddison WP et al 2007 Syst Biol 56:701 (BiSSE)
• Beaulieu JM et al 2012 Evolution 66:2369 (OUwie)
• Uyeda JC & Harmon LJ 2014 Syst Biol 63:902 (bayou)
• FitzJohn RG 2009 Syst Biol 58:595 (root prior)
• Whelan S et al 2018 Bioinformatics 34:3929 (PREQUAL)
• Di Franco A et al 2019 BMC Evol Biol 19:21 (HmmCleaner)
• Emms DM & Kelly S 2017 MBE 34:3267 (STRIDE rooting)
• Tria FDK et al 2017 Nat Eco Evo 1:0193 (MAD rooting)
• Williams TA et al 2017 PNAS 114:E4602 (ALE-rooting of deep phylogenies)
• Hu Z et al 2019 MBE 36:1086 (PhyloAcc)
• Fukushima K & Pollock DD 2023 Nat Eco Evo 7:155 (CSUBST)
• Muffato M et al 2023 Nat Eco Evo 7:355 (AGORA)

Related Skills

• comparative-genomics/positive-selection - Branch- and site-level selection inference on ancestral branches
• comparative-genomics/ortholog-inference - Define orthogroups whose alignments feed ASR
• comparative-genomics/gene-tree-species-tree-reconciliation - DTL-aware ancestral gene-content inference; root inference via ALE
• comparative-genomics/whole-genome-duplication - Ks-dating provides time scale for ancestral state inference
• comparative-genomics/comparative-annotation-projection - Project ancestral CDS to descendants for validation
• phylogenetics/modern-tree-inference - Generate rooted ML/Bayesian trees as ASR scaffold
• phylogenetics/bayesian-inference - RevBayes / MrBayes priors for Bayesian ASR
• phylogenetics/divergence-dating - Time-calibrated trees as input for absolute-time ASR
• alignment/multiple-alignment - PRANK / MACSE indel-aware alignment before sequence ASR
• alignment/alignment-trimming - PREQUAL / HmmCleaner filtering before ASR