# Bio Core Phylogenetics

> Build NJ/UPGMA trees from aligned FASTA with BioPython Bio.Phylo, score p-distance/JC69/K2P models, and parse Newick with bootstrap support. Use for tree building from an MSA, Newick I/O, or bootstrap-value interpretation.

- Skill: `pavel-kravchenko/bio-core-phylogenetics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-phylogenetics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-phylogenetics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

```bash
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

| Model | Formula | Assumption |
|---|---|---|
| p-distance | diffs / compared_sites | No correction |
| JC69 | −3/4 · ln(1 − 4p/3) | Equal rates, equal base freq |
| K2P | −½·ln(1−2S−V) − ¼·ln(1−2V) | Distinguishes transitions (S) from transversions (V) |

**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.

```python
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.

```python
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.

```python
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`.

```python
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

| Tool | Algorithm | Use case |
|---|---|---|
| BioPython `DistanceTreeConstructor` | NJ, UPGMA | Quick trees from a distance matrix |
| IQ-TREE2 | ML + model selection | Publication-quality ML trees |
| RAxML-NG | ML | Large datasets, fast |
| MrBayes | Bayesian | Posterior probability support |
| FastTree | Approximate ML | Very large alignments |

Run with automatic model selection and 1000 ultrafast bootstraps:

```bash
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