Amino Acid Recoding

In molecular systematics and phylogenomics, site-heterogeneous phylogenetic models and physicochemical amino-acid recoding are the primary statistical frameworks developed to overcome long-branch attraction (LBA) and compositional bias across deep evolutionary divergences—such as the root of the eukaryotic tree of life (Asgard archaea to LECA), the basal metazoan split (Ctenophora vs. Porifera), the fungal placement of Microsporidia, and the monophyly of Ecdysozoa.

When standard site-homogeneous Maximum Likelihood (ML) or Bayesian models (such as WAG+Γ, JTT+Γ, or LG+Γ) evaluate a multiple sequence alignment, they apply a single 20 × 20 empirical amino-acid exchangeability matrix (R) and a single stationary equilibrium frequency vector (π = [π₁, π₂, ..., π₂₀]) to every alignment column i, allowing columns to vary only in their scalar evolutionary rate r_i drawn from a discrete Gamma (+Γ) or FreeRate (+R) distribution. While computationally fast, this assumption violates basic protein biophysics: in any folded three-dimensional protein domain, steric packing, hydrogen bonding, and electrostatic charge constraints restrict each individual alignment site i to a tiny effective amino-acid alphabet (K_eff ≈ 2 to 5 biochemically compatible residues) rather than all 20 amino acids. Because the probability that two rapidly evolving lineages independently mutate to the same amino acid (convergent homoplasy) scales as 1 / (K_eff − 1) rather than 1 / 19, site-homogeneous matrices severely underestimate homoplasy along long branches.

To model site-specific biochemical constraints directly within the Markov substitution rate matrix Q^(k), site-heterogeneous profile mixture models assign alignment columns to distinct equilibrium amino-acid profile classes π^(k) = (π₁^(k), ..., π₂₀^(k)):

1. The Instantaneous Rate Matrix Q^(k) for Profile Class k

Within profile class k, the instantaneous substitution rate Q_ab^(k) from amino acid a to amino acid b (a ≠ b) is factored into a global symmetric exchangeability matrix R_ab (either estimated directly from the alignment as in GTR, or fixed to the empirical LG matrix) multiplied by the class-specific target amino-acid equilibrium frequency π_b^(k):

Q_ab^(k) = (1 / Z_k) · R_ab · π_b^(k) (for a ≠ b), with Q_aa^(k) = − ∑_b≠a Q_ab^(k)

where Z_k is a normalization factor ensuring that branch lengths t are measured in expected substitutions per site across the mixture. For example, a buried hydrophobic core class k_hydro has π^(k) concentrated almost exclusively on {Ile, Val, Leu, Met} (K_eff ≈ 3.2), so the transition probability matrix P^(k)(t) = exp(Q^(k) · t) correctly predicts that two long branches will frequently converge on Ile or Val by chance alone—preventing the likelihood engine from misinterpreting convergent Ile–Ile matches as synapomorphies.

2. Effective Number of Amino Acids (K_eff) from Shannon Entropy

To quantify the biochemical restriction of any alignment column i or profile class π^(k), molecular evolutionists compute the Shannon site entropy H_i (in bits or nats) and exponentiate it to obtain the effective amino-acid alphabet size K_eff (the Hill number of order 1):

H_i = − ∑_a=1..20 π_a⁽ⁱ⁾ ln( π_a⁽ⁱ⁾ ) ⟹ K_eff(i) = exp( H_i ) = ∏_a=1..20 ( π_a⁽ⁱ⁾ )^−π_a(i) ∈ [1.0, 20.0]

• Invariant Catalytic Residue (H_i = 0): K_eff = 1.0 (zero homoplasy, zero evolutionary change).
• Constrained Solenoid/Propeller Repeat Core (WD40, Kelch, TPR, Armadillo, C4 Zinc Finger): K_eff ≈ 3.1 to 4.2 (high double-hit convergence probability 1 / (K_eff − 1) ≈ 31% to 48%).
• Globular Hydrophobic Core vs. Polar Surface: K_eff ≈ 4.5 to 8.5.
• Unconstrained Uniform Distribution (π_a = 1/20): K_eff = 20.0 (theoretical random ceiling assumed by unweighted Jukes–Cantor/Poisson models).

Even when across-site heterogeneity (K_eff) is modeled via CAT or LG+C60, a second orthogonal source of systematic error—across-branch compositional heterogeneity (non-stationarity)—arises when different lineages evolve under divergent genomic GC/AT mutational biases or thermophilic adaptation:

1. The FYMINK vs. GARP Codon-Driven Amino-Acid Shift: Because genetic codons have fixed nucleotide compositions at their first and second positions, genomes with extreme AT-bias (such as intracellular Microsporidia, rhabditid nematodes, endosymbiotic Buchnera and Rickettsia, Plasmodium falciparum, and mitochondrial genomes) systematically accumulate amino acids encoded by AU-rich codons (Phe [UUY], Tyr [UAY], Met [AUG], Ile [AUH], Asn [AAY], Lys [AAR] — the FYMINK group). Conversely, genomes with high GC-bias (such as Streptomyces, halophilic archaea, choanoflagellates, or avian GC-biased gene conversion regions) systematically accumulate amino acids encoded by GC-rich codons (Gly [GGN], Ala [GCN], Arg [CGN], Pro [CCN] — the GARP group).
2. The Thermophilic IVYWREL Enrichment: Independently of GC content, hyperthermophilic archaea and bacteria growing at 80 °C to 105 °C enrich their proteomes in seven hydrophobic and salt-bridge-forming residues (Ile, Val, Tyr, Trp, Arg, Glu, Leu — the IVYWREL predictor, which correlates with optimal growth temperature T_opt at R > 0.93) while purging thermolabile Asn, Gln, Ser, and Cys.
3. Relative Composition Frequency Variability (RCFV): Introduced by Mark Wilkinson, Peter Foster, Thorsten Struck, and colleagues (Zhong et al., 2011), RCFV measures the total absolute deviation of taxon-specific amino-acid (or nucleotide) frequencies μ_ij from the alignment-wide mean frequency μ̄_i across n taxa and m character states (m = 20 for proteins, 4 for DNA):

RCFV = (1 / n) ∑_j=1..n ∑_i=1..m | μ_ij − μ̄_i |

Along with taxon-specific nRCFV (normalized for sequence length and missing data), elevated RCFV values (RCFV > 0.05–0.10 in protein partitions or > 0.15 in intrinsically disordered C-terminal tails) diagnose partitions where convergent FYMINK or GARP drift is actively pulling unrelated AT-rich or GC-rich lineages together on the tree.

When across-branch compositional drift (RCFV) is strong, how can a phylogeneticist strip away superficial intra-class substitutions (Ile ↔ Val, Lys ↔ Arg, Ala ↔ Gly ↔ Ser) while preserving deep structural substitutions that occurred hundreds of millions of years ago?

By recoding the 20-letter amino-acid alphabet into a reduced 6-letter or 4-letter physicochemical alphabet prior to running Maximum Likelihood or Bayesian inference (CAT-GTR on Dayhoff-6 states):

• Dayhoff-6 (6-State) Recoding (Derived from Margaret Dayhoff's PAM250 Log-Odds Matrix, 1978): Groups the 20 amino acids into six biochemically conservative exchangeability classes within which rapid homoplastic and GC/AT-driven substitutions occur, so that only inter-group transversions (changes in residue volume, charge, or hydrophobicity) are scored on the phylogenetic tree:
  1. 0 — Small / Neutral / Turn-Forming (A, G, P, S, T): Alanine, Glycine, Proline, Serine, Threonine.
  2. 1 — Acidic and Polar Amide (D, E, N, Q): Aspartate, Glutamate, Asparagine, Glutamine.
  3. 2 — Basic / Positively Charged (H, K, R): Histidine, Lysine, Arginine (masks thermophilic/GC-driven Lys ↔ Arg swaps!).
  4. 3 — Branched-Chain Aliphatic / Hydrophobic (I, L, M, V): Isoleucine, Leucine, Methionine, Valine (masks AT-driven Ile ↔ Val ↔ Leu swaps!).
  5. 4 — Aromatic (F, W, Y): Phenylalanine, Tryptophan, Tyrosine (masks AT-driven Phe ↔ Tyr swaps!).
  6. 5 — Sulfhydryl / Disulfide (C): Cysteine (kept in its own singleton class due to the structural permanence of disulfide bridges and metal coordination).
• Susko–Roger 4-State (SR4) Recoding (Susko & Roger, 2007): Optimizes bin assignments mathematically by minimizing within-bin substitution rate variance across empirical alignments into four states: {A, G, N, P, S, T}, {C, H, W, Y}, {D, E, K, Q, R}, and {F, I, L, M, V} (enabling direct use of standard 4-state nucleotide GTR+Γ and non-stationary GTR models in IQ-TREE and PhyloBayes).
• Kosiol 6-State (KGB6) and Sub-Domain Partitioned Recoding: Because recoding an entire alignment into 6 states discards ~40%–50% of total sequence information, applying Dayhoff-6 recoding indiscriminately to slowly evolving catalytic cores can reduce bootstrap support on short internal branches. Modern sub-domain partitioning protocols resolve this trade-off by keeping conserved globular catalytic cores and synapomorphic Short Linear Motifs (SLiMs) in the full 20-state LG+C60-PMSF model while selectively stripping or Dayhoff-6 recoding high-RCFV disordered tails and low-K_eff tandem repeat domains.

Why is PMSF (Posterior Mean Site Frequency) 20- to 35-fold faster than standard LG+C60 in IQ-TREE?

In a standard 60-component empirical profile mixture model (LG+C60+F+Γ4) with 4 Gamma rate categories, computing the likelihood of a single candidate tree topology during tree search requires running Felsenstein's pruning algorithm 60 × 4 = 240 separate times at every alignment column (and storing 240 conditional likelihood vectors at every internal node of the tree), which consumes massive RAM and CPU time. In 2018, Huai-Chun Wang, Bui Quang Minh, Edward Susko, and Andrew Roger introduced the PMSF approximation: first, IQ-TREE fits the full LG+C60+F+Γ mixture onto a fast initial guide tree (BioNJ or LG+F+Γ tree) and uses Bayes' theorem to compute the posterior probability p(k | X_i) that alignment column i belongs to each of the 60 profile classes k. It then averages the 60 profiles weighted by those posterior probabilities into a single custom 20-element equilibrium frequency vector π̂⁽ⁱ⁾ = ∑_k=1..60 p(k | X_i) π^(k) unique to site i. During all subsequent NNI/SPR tree searches and 1,000-replicate bootstrapping, each site i only needs 1 equilibrium vector (4 Gamma categories instead of 240)—cutting CPU time and memory by ~20- to 35-fold while matching or outperforming full LG+C60 and Bayesian CAT-GTR in suppressing long-branch attraction!

How do Posterior Predictive Simulations test whether a phylogenetic model is adequate for a dataset?

Standard model-selection criteria such as the Bayesian Information Criterion (BIC) or Akaike Information Criterion (AIC) only tell you which model is the best relative performer among a list of candidates—they do not tell you whether the winning model actually fits the real biological data (absolute model adequacy). Posterior Predictive Simulations (implemented in PhyloBayes and RevBayes) test absolute adequacy by taking parameter draws from the fitted model (e.g., LG+Γ vs. CAT-GTR), simulating synthetic sequence alignments of the exact same size along the inferred tree, and comparing summary statistics of the real alignment against the null distribution of simulated alignments:

1. Across-Site Biochemical Diversity Test (Mean Biochemical Diversity / Number of Distinct Amino Acids per Site): Site-homogeneous LG+Γ consistently predicts far too many distinct amino acids per column (e.g., 7.8 amino acids/site simulated vs. 4.2 observed in the real alignment, Z > 10 SD failure), whereas CAT-GTR and LG+C60-PMSF accurately match the real biochemical diversity.
2. Across-Branch Compositional Homogeneity Test (Maximum Pairwise Compositional Chi-Square / RCFV): Compares the observed RCFV across taxa to the simulated stationary null; if Z > 3.0 SD, the alignment violates compositional stationarity and requires Dayhoff-6 recoding, high-RCFV sub-domain stripping, or branch-heterogeneous (NDCH / COaLA) models.

1. Lartillot, N., & Philippe, H. (2004). "A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process." Molecular Biology and Evolution, 21(6), 1095–1109. https://doi.org/10.1093/molbev/msh112
2. Susko, E., & Roger, A. J. (2007). "On reduced amino acid alphabets for phylogenetic inference." Molecular Biology and Evolution, 24(9), 2139–2150. https://doi.org/10.1093/molbev/msm144
3. Quang, L. S., Gascuel, O., & Lartillot, N. (2008). "Empirical profile mixture models for phylogenetic reconstruction." Bioinformatics, 24(20), 2317–2323. https://doi.org/10.1093/bioinformatics/btn445
4. Wang, H. C., Minh, B. Q., Susko, E., & Roger, A. J. (2018). "Modeling site heterogeneity with posterior mean site frequency profiles accelerates accurate phylogenomic estimation." Systematic Biology, 67(2), 216–235. https://doi.org/10.1093/sysbio/syx068