Bio Phylo Distance Calculations

Version Compatibility

Reference examples tested with: BioPython 1.83+. Model-corrected alternatives: ape 5.8+, phangorn, FastME, scikit-bio.

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

• Python: pip show biopython then help(module.function) to check signatures
• R: packageVersion('ape') then ?dist.dna to verify model strings
• CLI: fastme --version then fastme --help to confirm flags

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

Bio.Phylo DistanceCalculator offers only identity / substitution-matrix distances and has NO model-based JC/K80/TN93; for model-corrected DNA distances use ape dist.dna(model=...) (introspect its model list and gamma= support if the version differs).

Distance Calculations -- Distance Erases the Information ML Keeps; the Correction Is the Model

"Build a tree from my alignment" -> Estimate model-corrected pairwise distances, then cluster them into a tree, knowing the correction is where the biology lives and the matrix discards what ML uses.

• Python: DistanceCalculator(...).get_distance(aln) then DistanceTreeConstructor().nj(dm) (Bio.Phylo, identity-only)
• R: dist.dna(x, model='TN93') then fastme.bal(d) / nj(d) (ape, model-corrected)

Scope: computing evolutionary distance matrices and building/bootstrapping distance trees (NJ, BIONJ, FastME, UPGMA), and testing saturation before trusting them. ML inference, model selection, and where NJ seeds an ML search -> modern-tree-inference. Alignment quality and trimming, which gate every distance -> alignment/alignment-io. Reading/writing/converting the trees these methods emit -> tree-io. Rooting and pruning -> tree-manipulation.

The Single Most Important Modern Insight

A distance-based phylogeny is built in two lossy steps and both losses are silent. First, a distance is a model-corrected estimate of expected substitutions per site, not the observed proportion of differing sites; second, the matrix collapses each sequence pair to one scalar and discards the site-pattern information that maximum likelihood evaluates. Distance methods are therefore fast but doubly limited -- by how good the correction is, and by what a matrix of scalars can encode -- and they are structurally incapable of using the per-site signal that makes ML the gold standard for hard problems. Three load-bearing facts:

1. The correction is the whole biological content. The raw p-distance (proportion of differing sites) systematically UNDERCOUNTS true divergence because every site that mutated twice (back-substitution), in parallel on both lineages, or to a third state (multiple hit) is missed, and the gap grows without bound toward saturation. The algorithm downstream (NJ, FastME) is just arithmetic on the numbers handed to it; feed it p-distances on divergent data and it returns a confidently wrong tree, fast.
2. The matrix throws away what ML keeps. Two completely different alignments can produce the same distance matrix and therefore the same tree. ML never collapses the data, so it cannot be fooled this way; once the matrix exists, which sites drove the divergence, among-site rate heterogeneity, and competing site-pattern signal are gone.
3. Consistency is conditional, not a guarantee. NJ is statistically consistent -- it returns the true tree as sequence length grows -- but only if the input distances are additive and correctly estimated. A consistent algorithm fed biased (saturated, misspecified) distances converges confidently on the WRONG tree. "I used NJ so it is consistent" is not a defense; the correction is exactly what most users get wrong.

Distance Corrections

Every nucleotide correction inverts the probability that a site differs under a model, so as p approaches the equilibrium difference (0.75 for equal base frequencies under JC) the correction diverges to infinity and its variance explodes -- the mathematical face of saturation.

LogDet / paralinear (Lockhart et al. 1994; Lake 1994) is the non-obvious one: JC/K80/TN93 all assume stationarity (constant base composition across the tree), and when composition drifts (thermophiles, AT-rich insect mtDNA, GC-rich chloroplasts) every stationary correction groups taxa by base composition rather than ancestry -- compositional attraction, a sibling of long-branch attraction that fools NJ and stationary-model ML alike. LogDet is computed from the determinant of the 4x4 pairwise divergence matrix and is consistent under the general Markov model, at the cost of needing long sequences and not easily taking gamma. If a topology flips between JC and LogDet, suspect compositional heterogeneity.

Algorithm Decision

UPGMA is the key warning: it forces every tip equidistant from the root, which molecular data essentially never satisfies because lineages evolve at different rates. A fast-evolving lineage gets pulled toward the tips and grouped by total divergence rather than true ancestry -- a clustering-by-rate artifact analogous to long-branch attraction. Flat rule: do not use UPGMA for molecular phylogeny unless a clock is independently established. NJ does not assume a clock because its Q-matrix corrects each pairwise distance for that taxon's average divergence to everyone else. NJ/BIONJ/FastME can seed ML searches: some ML programs build a BIONJ or ME starting tree in O(n^3) before likelihood NNI/SPR (though IQ-TREE2 and RAxML-NG default to parsimony starting trees), so this machinery still has a place in an all-ML pipeline.

Tool Taxonomy

Loud flag, the second flat rule: Bio.Phylo and scikit-bio do NOT correct for multiple hits -- their built-in distance is p-distance or score-based. Presenting Bio.Phylo's 'identity' distance as a "Jukes-Cantor tree" is a common and wrong shortcut. For a real correction, compute the matrix in ape/FastME (or compute it yourself) and pass it in.

Build a Model-Corrected Distance Matrix and NJ Tree

Goal: Turn an alignment into a tree, contrasting Bio.Phylo's identity-only path against ape's model-corrected path.

Approach: In Python use Bio.Phylo for the NJ/UPGMA algorithm and pure-Python I/O, but treat its distance as an uncorrected p-distance; for any divergent data move the distance step to ape dist.dna, then cluster with FastME (best) or NJ.

from Bio import AlignIO
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

aln = AlignIO.read('alignment.fasta', 'fasta')
calc = DistanceCalculator('identity')        # identity-only: this is a p-distance, NOT a JC/K80 correction
dm = calc.get_distance(aln)                   # multiple/back/parallel hits are NOT corrected here
tree = DistanceTreeConstructor().nj(dm)       # NJ algorithm is correct; the DISTANCES are the limitation
# For a model-corrected distance, build the matrix in ape (below) and pass it to skbio/Bio.Phylo, or stay in R.

library(ape)
aln <- read.dna('alignment.fasta', format = 'fasta')
d <- dist.dna(aln, model = 'TN93', gamma = 0.5)   # model-corrected; gamma applies ASRV (alpha < 1 = strong)
tree <- fastme.bal(d, nni = TRUE, spr = TRUE)     # balanced minimum evolution: the modern best distance tree
# nj(d) / bionj(d) are the faster single-pass alternatives; bionj seeds ML searches

Pre-flight: Test Substitution Saturation

Goal: Decide whether the data retain phylogenetic signal before trusting any deep distance tree.

Approach: Compute Xia's index of substitution saturation Iss and compare it to the simulation-derived critical value Iss.c; a saturation plot of transitions against a corrected distance is the visual companion.

library(ape)                    # the formal entropy-based Iss vs Iss.c test lives in DAMBE; this is the ape/base-R saturation plot
d_jc  <- dist.dna(aln, model = 'JC69')
ts_tv <- dist.dna(aln, model = 'TS')   # transitions; plot against d_jc -- a PLATEAU means saturated, signal erased
plot(d_jc, ts_tv)                       # unsaturated = roughly linear; bent-over transition curve = drop those sites
# Interpretation gate (Xia 2003): Iss < Iss.c => signal retained (usable); Iss >= Iss.c => substantially saturated, do not use.

Bootstrap a Distance Tree

Goal: Quantify how reproducible each clade is under resampling -- precision, not accuracy.

Approach: Resample alignment columns with replacement, recompute the matrix with the SAME correction, rebuild with the SAME algorithm, and summarize clade frequencies on a majority-rule consensus.

library(ape)
boot <- boot.phylo(tree, aln, function(x) fastme.bal(dist.dna(x, model = 'TN93')), B = 500)
# 100-1000 replicates standard; support measures SAMPLING stability only -- it cannot detect a bias in the distances

In Python, Bio.Phylo's bootstrap_consensus(aln, 100, DistanceTreeConstructor(calc, 'nj'), majority_consensus) does the same on the identity distance (same correction caveat).