Phylogenetics & Evolution

Lesson 7 of 13 · 11 min

Substitution Models: JC to GTR

A phylogenetic tree is only as trustworthy as the substitution model used to score it. Such a model gives the probability that one nucleotide turns into another over a branch, and every maximum-likelihood or Bayesian method needs one to convert an alignment into branch lengths and a topology. Picking a model that fits your data is not cosmetic tuning; it changes the tree you infer.

Why Raw Distances Fail

If two sequences differ at 30 percent of their sites, the true number of substitutions per site is higher than 0.30, because some positions were hit more than once, reverted to the original base, or changed in parallel. These hidden and multiple substitutions make the raw p-distance saturate and underestimate real divergence, worst of all between distant taxa. A substitution model corrects for the unseen events by estimating expected substitutions per site instead of counting visible mismatches.

JC69 To GTR

The classic DNA models form a nested ladder of rising realism. JC69 assumes equal base frequencies and one rate for every change; K80 adds a separate transition/transversion ratio called kappa; HKY combines that ratio with unequal, estimated base frequencies; and GTR gives a distinct exchange rate to all six reversible substitution pairs on top of unequal frequencies. Exchange rates are always relative, so one is fixed as a reference, which is why GTR contributes five free rate parameters rather than six. Extra parameters fit better but spend degrees of freedom, and that trade-off is exactly what model selection must referee.

text
Model  IQ-TREE  Base frequencies      Substitution rates                Free params
JC69   JC       equal (0.25 each)     1 (all changes equal)             0
K80    K2P      equal (0.25 each)     2 (transition vs transversion)    1
HKY85  HKY      unequal (estimated)   2 (transition vs transversion)    4
GTR    GTR      unequal (estimated)   6 (all reversible pairs)          8

Rate-variation add-ons (attach to any model above):
+G4    among-site gamma variation, 4 rate categories (shape alpha)   +1 param
+I     proportion of invariable, never-changing sites                +1 param

The two add-ons in the table model the fact that sites evolve at different speeds. +G4 replaces a single rate with a discrete gamma distribution split into four rate categories, tuned by a shape parameter alpha, where a small alpha means strong rate variation across sites. +I sets aside a proportion of invariable sites that never change, and although +I and +G4 can be combined, their parameters are partly correlated and can be hard to estimate together.

Automated Model Selection

bash
# ModelFinder Plus: score candidate models, then build the ML tree
iqtree2 -s example.phy -m MFP -B 1000 -T AUTO --prefix primates

# The ranked model table is written to primates.iqtree
grep -A 12 "sorted by BIC" primates.iqtree
ModelFinder will test up to 286 DNA models (sample size: 898 nucleotides) ...
 ...
Akaike Information Criterion:           GTR+F+I+G4
Corrected Akaike Information Criterion:  GTR+F+I+G4
Bayesian Information Criterion:          GTR+F+G4
Best-fit model according to BIC: GTR+F+G4

List of models sorted by BIC scores:

Model                  LogL         AIC      w-AIC        AICc     w-AICc         BIC      w-BIC
GTR+F+G4           -6183.253   12426.506 +   0.413   12428.651 +   0.430   12570.506 +   0.879
GTR+F+I+G4         -6181.900   12425.800 +   0.587   12428.091 +   0.569   12574.600 +   0.114
HKY+F+G4           -6201.744   12455.488 -   0.000   12457.100 -   0.000   12580.288 +   0.007
GTR+F+I            -6190.402   12440.804 -   0.000   12442.949 -   0.000   12584.804 -   0.001
K2P+G4             -6233.187   12512.374 -   0.000   12513.637 -   0.000   12622.774 -   0.000
JC+G4              -6258.921   12561.842 -   0.000   12562.999 -   0.000   12667.442 -   0.000

The command runs ModelFinder Plus (-m MFP), which scores hundreds of models before building the tree, with 1000 ultrafast bootstrap replicates (-B 1000) and automatic threading (-T AUTO). In the report the +F suffix marks models whose base frequencies are counted straight from the alignment, which is why HKY and GTR appear as HKY+F and GTR+F while the equal-frequency JC and K2P do not. Lower AIC, AICc, and BIC all mean a better fit, and IQ-TREE selects by BIC by default; notice that AIC and AICc here prefer the richer GTR+F+I+G4, while BIC's heavier penalty for the extra invariable-sites parameter picks the simpler GTR+F+G4. Trust BIC for a conservative choice, rerun with -m MF for selection only, or add -merit AICc to switch the ranking criterion.

Try it yourself: Grab a small single-gene alignment (the example.phy shipped with IQ-TREE, or any FASTA of one gene from 8 to 15 species) and run iqtree2 -s your_alignment -m MFP -B 1000 -T AUTO. Open the resulting .iqtree file, find the line 'Best-fit model according to BIC', and check whether including +I changed the winner versus AIC. Then force the simplest model with -m JC and compare its log-likelihood to confirm the fitted model really scores better.