Bio Core Phylogenetics

Phylogenetics: Reconstructing Evolutionary History

When to Use

• Building a tree (NJ or UPGMA) from an aligned FASTA/MSA of DNA or protein sequences.
• Computing and comparing distance models (p-distance, Jukes-Cantor, Kimura 2-parameter) for a set of sequences.
• Reading, writing, or navigating Newick-format trees (Bio.Phylo), e.g. finding a common ancestor or pairwise tree distance.
• Interpreting or generating bootstrap support values for clades.
• Deciding whether a quick distance-based tree suffices or a full ML/Bayesian pipeline (IQ-TREE, MrBayes) is warranted.

Version Compatibility

BioPython ≥1.81, NumPy ≥1.24, Python ≥3.9. IQ-TREE2 ≥2.2 / RAxML-NG ≥1.2 for the external ML reference commands below.

Prerequisites

pip install biopython numpy

Prior concepts/skills: a multiple sequence alignment already computed (see bio-alignment-msa-parsing), basic sequence I/O (bio-sequence-io-read-sequences).

Distance Models

Goal: turn a set of aligned sequences into a corrected pairwise distance matrix. Approach: compute raw p-distance, then apply the JC69 or K2P correction, which inflates the distance to account for unseen multiple hits at the same site.

import numpy as np

def p_distance(seq1, seq2):
    """Proportion of differing sites, ignoring gap columns."""
    pairs = [(a, b) for a, b in zip(seq1, seq2) if a != '-' and b != '-']
    diffs = sum(1 for a, b in pairs if a != b)
    return diffs / len(pairs) if pairs else 0.0

def jukes_cantor(seq1, seq2):
    """JC69-corrected distance; undefined (returns inf) once p >= 0.75 (saturation)."""
    p = p_distance(seq1, seq2)
    if p >= 0.75:
        return float('inf')
    return -0.75 * np.log(1 - 4 * p / 3)

def kimura_k2p(seq1, seq2):
    """Kimura 2-parameter distance: separates transitions (S) from transversions (V)."""
    purines, pyrimidines = set('AG'), set('CT')
    pairs = [(a, b) for a, b in zip(seq1, seq2) if a != '-' and b != '-']
    n = len(pairs)
    S = sum(1 for a, b in pairs if a != b and
            ((a in purines and b in purines) or (a in pyrimidines and b in pyrimidines))) / n
    V = sum(1 for a, b in pairs if a != b and
            not ((a in purines and b in purines) or (a in pyrimidines and b in pyrimidines))) / n
    t1, t2 = 1 - 2 * S - V, 1 - 2 * V
    if t1 <= 0 or t2 <= 0:
        return float('inf')
    return -0.5 * np.log(t1) - 0.25 * np.log(t2)

def build_distance_matrix(sequences, method=jukes_cantor):
    """Pairwise distance matrix (NumPy array) for a list of aligned sequence strings."""
    n = len(sequences)
    matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            d = method(sequences[i], sequences[j])
            matrix[i, j] = matrix[j, i] = d
    return matrix

Tree Construction with BioPython

Goal: go from an alignment file to a tree object you can query and draw. Approach: load the alignment with AlignIO, get a DistanceMatrix via DistanceCalculator, then build with DistanceTreeConstructor — prefer .nj() (Neighbor-Joining, no clock assumption) over .upgma() (assumes a molecular clock) unless you specifically need an ultrametric tree.

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

def build_tree(fasta_alignment_path, method="nj", model="identity"):
    """
    Build a phylogenetic tree from an aligned FASTA file.
    model: 'identity' for DNA, 'blosum62'/'pam250' for protein.
    method: 'nj' (Neighbor-Joining, preferred) or 'upgma' (assumes molecular clock).
    """
    aln = AlignIO.read(fasta_alignment_path, "fasta")
    calc = DistanceCalculator(model)
    dm = calc.get_distance(aln)

    constructor = DistanceTreeConstructor()
    tree = constructor.nj(dm) if method == "nj" else constructor.upgma(dm)
    return tree

tree = build_tree("alignment.fasta", method="nj")
Phylo.draw_ascii(tree)
tree.distance("Human", "Chimp")            # pairwise distance along tree branches
tree.common_ancestor("Human", "Chimp")      # MRCA node
tree.count_terminals(), tree.total_branch_length()

Newick I/O and Navigation

Goal: read/write trees in the standard interchange format and query their structure. Approach: Bio.Phylo.read/.write handle Newick, NEXUS, and phyloXML via the format argument; trees are navigated as nested Clade objects.

import io
from Bio import Phylo

newick = "((Human:0.01,Chimp:0.012):0.02,(Mouse:0.25,Rat:0.23):0.3,Zebrafish:0.6);"
tree = Phylo.read(io.StringIO(newick), "newick")

tree.get_terminals()          # leaf (tip) Clade objects
tree.get_nonterminals()       # internal Clade objects
tree.distance("Human", "Chimp")
tree.common_ancestor("Human", "Chimp")

out = io.StringIO()
Phylo.write(tree, out, "newick")   # or format="nexus" / "phyloxml"

Bootstrap Support

Goal: quantify how well-supported each clade is given sampling variance in the alignment. Approach: resample alignment columns with replacement N times, rebuild a tree each time, and count how often each clade from the original tree recurs. Bio.Phylo.Consensus automates this on top of DistanceTreeConstructor.

import random
from Bio import AlignIO
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo.Consensus import bootstrap_trees, get_support

def bootstrap_tree(fasta_alignment_path, n_replicates=100, model="identity", method="nj"):
    """
    Build an NJ/UPGMA tree annotated with bootstrap support (0-100) on each internal node.
    """
    aln = AlignIO.read(fasta_alignment_path, "fasta")
    calc = DistanceCalculator(model)
    constructor = DistanceTreeConstructor(calc, method)

    target_tree = constructor.build_tree(aln)
    replicate_trees = list(bootstrap_trees(aln, n_replicates, constructor))
    support_tree = get_support(target_tree, replicate_trees)  # confidence in [0, 1] per clade
    return support_tree

random.seed(0)
support_tree = bootstrap_tree("alignment.fasta", n_replicates=100)
for clade in support_tree.get_nonterminals():
    if clade.confidence is not None:
        print(f"clade with {clade.count_terminals()} tips: support={clade.confidence * 100:.0f}%")

External ML/Bayesian Tools

Run with automatic model selection and 1000 ultrafast bootstraps:

iqtree2 -s alignment.fasta -m TEST -B 1000 -T AUTO

Pitfalls

• Tip order does not imply closeness: Only shared branching points (nodes) indicate relationships. Rotating branches around any internal node produces an equivalent tree. Always look at which node is shared, not which tips are adjacent.
• UPGMA assumes a molecular clock: Produces rooted, ultrametric trees (all tips equidistant from root). This assumption is frequently violated. Neighbor-Joining (NJ) does not assume a clock — prefer NJ for most datasets.
• Branch lengths have units: In sequence-based trees, branch lengths = expected substitutions per site. They are not time units unless the tree is explicitly time-calibrated.
• Bootstrap values are NOT probabilities: Bootstrap 95 means 95% of replicates recovered that clade — not 95% probability of being correct. Thresholds: ≥70 for NJ/parsimony, ≥80 for ML.
• Model selection matters for ML: Default to IQ-TREE's built-in model selection (-m TEST) rather than hardcoding JC69. GTR+Γ+I is common for nucleotide data; LG or WAG for protein.
• p-distance underestimates true distance: Multiple hits at the same site are invisible. JC69 correction becomes undefined at p ≥ 0.75 (saturation).

See Also

• bio-phylogenetics-tree-io — reading/writing additional tree formats (NEXUS, phyloXML)
• bio-phylogenetics-tree-manipulation — rerooting, pruning, and editing existing trees
• bio-phylogenetics-modern-tree-inference — ML/Bayesian inference with IQ-TREE, ete3
• bio-alignment-msa-parsing — producing the aligned input this skill consumes