Model-Based Inference: Maximum Likelihood

Neighbor-joining gave you a tree from a distance matrix, and parsimony gave you the tree needing the fewest changes. Neither one asks how probable your actual alignment is under an explicit model of how sequences evolve. Maximum likelihood does exactly that: it scores every candidate tree by the probability of observing your data, then searches for the tree that scores highest. This lesson builds a working picture of that score and runs two standard tools, IQ-TREE and RAxML-NG, on the same alignment you used earlier.

What likelihood means

The likelihood of a tree is the probability of the alignment given that tree and a chosen substitution model, written L = P(alignment | topology, branch lengths, model). It is computed column by column and multiplied across all sites, so the raw number is astronomically small. We therefore always report the natural log of it, the log-likelihood or lnL, which is a negative number. Less negative is better: a tree with lnL of -8423 fits the data better than one with lnL of -8600. Note that the score depends jointly on the topology, the branch lengths, and the model, so all three are optimized together, not separately.

Model and search

A substitution model is a rate matrix describing how one nucleotide changes into another, for example GTR (general time-reversible) with a gamma term (+G) that lets sites evolve at different rates. Picking the wrong model biases branch lengths and support, so IQ-TREE's ModelFinder tests many models and keeps the one with the best Bayesian Information Criterion. The tree itself cannot be found by brute force because the number of possible topologies explodes with taxa, so both tools use a heuristic search: they start from a fast parsimony or BIONJ tree, apply local rearrangements (NNI and SPR moves), keep any move that raises the lnL, and repeat from several starting trees to avoid getting stuck in a local optimum.

# ModelFinder Plus (MFP) picks the model, then infers the ML tree
# -B 1000 runs 1000 ultrafast bootstrap replicates; -T AUTO picks thread count
iqtree2 -s aln.fasta -m MFP -B 1000 -T AUTO --prefix ml_iqtree

IQ-TREE multicore version 2.3.6 for MacOS ARM 64-bit

ModelFinder will test up to 484 DNA models (sample size: 1200) ...
Akaike Information Criterion:           GTR+F+G4
Bayesian Information Criterion:         GTR+F+G4
Best-fit model: GTR+F+G4 chosen according to BIC

Optimizing model parameters (epsilon = 0.100)
Initial log-likelihood: -8451.902
BEST SCORE FOUND : -8423.157

Analysis results written to:
  IQ-TREE report:                ml_iqtree.iqtree
  Maximum-likelihood tree:       ml_iqtree.treefile
  Likelihood distances:          ml_iqtree.mldist
  Ultrafast bootstrap approximation results written to:
    Split support values:        ml_iqtree.splits.nex
    Consensus tree:              ml_iqtree.contree
  Screen log file:               ml_iqtree.log

Open ml_iqtree.iqtree for the human-readable report: it names the best-fit model, prints the final lnL, and lists estimated model parameters. The ML tree is in ml_iqtree.treefile, a Newick string whose branch lengths are in expected substitutions per site, not raw counts. The numbers at internal nodes are ultrafast bootstrap supports on a 0 to 100 scale, and a value of 95 or higher is conventionally read as strong support for that clade. Load the treefile into FigTree or iTOL to view it.

# Cross-check with RAxML-NG using a fixed model
# --all does an ML search plus 100 Felsenstein bootstraps in one run
raxml-ng --all --msa aln.fasta --model GTR+G --prefix ml_raxml --seed 42 --bs-trees 100

RAxML-NG v1.2.2 released on 30.04.2024

Analysis options: ML tree search + bootstrapping
Starting trees: 10 random + 10 parsimony

[00:00:12] ML tree search completed, best tree logLH: -8421.984
Final LogLikelihood: -8421.984

Optimized model parameters saved to: ml_raxml.raxml.bestModel
Best ML tree saved to:               ml_raxml.raxml.bestTree
Best ML tree with Felsenstein bootstrap (FBP) support values saved to:
                                     ml_raxml.raxml.support

The two ML runs land on almost the same topology, and their final lnL values (-8423.2 versus -8422.0) are close because both fit essentially the same GTR + gamma model. The roughly one-unit gap comes from the two tools running independent heuristic searches and estimating base frequencies slightly differently (IQ-TREE's +F uses empirical frequencies counted from the alignment), not from one model fitting fundamentally better; raw lnL from two separate searches is not a clean model comparison. Comparing against the earlier results, the strongly supported clades match your NJ and parsimony trees, which is reassuring. Where they differ, trust the ML trees: parsimony often leaves short internal branches ambiguous, and uncorrected NJ distances do not account for the multiple hidden substitutions the substitution model corrects for, nor for the among-site rate variation the gamma term captures, so ML branch lengths and node resolution are more reliable on divergent sequences.