drug-discovery-prompts/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-phylogenetics/references/tree_building.md

---
title: "Tree Building and Analysis"
task: ""
lineage_type: import
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-phylogenetics/references/tree_building.md
upstream_sha: e2520a96
imported_at: 2026-06-26
prompt_class: prompt
upstream_changes: accepted
author: upstream
validated: false
---

# Tree Building and Analysis

Detailed guide for phylogenetic tree construction and analysis using Biopython, PhyKIT, and DendroPy.

---

## Tree File Loading

### Auto-detect Tree Format

```python
from Bio import Phylo
import io

def load_tree(filepath):
    """Load tree with auto-format detection.

    Supports: Newick, Nexus
    """
    with open(filepath, 'r') as f:
        content = f.read().strip()

    # Try Newick first (most common)
    try:
        tree = Phylo.read(io.StringIO(content), 'newick')
        return tree, 'newick'
    except Exception:
        pass

    # Try Nexus
    try:
        tree = Phylo.read(io.StringIO(content), 'nexus')
        return tree, 'nexus'
    except Exception:
        pass

    raise ValueError(f"Cannot parse tree: {filepath}")
```

### Usage

```python
# Load any tree format
tree, format_detected = load_tree("gene1.nwk")
print(f"Format: {format_detected}")
print(f"Terminals: {tree.count_terminals()}")
```

---

## PhyKIT Tree Metrics

### Treeness

```python
from phykit.services.tree.treeness import Treeness
from types import SimpleNamespace

def phykit_treeness(tree_file):
    """Calculate treeness (internal branch length / total branch length).

    Treeness measures the proportion of tree length on internal branches.
    Higher values (closer to 1) indicate stronger phylogenetic signal.

    Returns: float (0-1)
    """
    t = Treeness(SimpleNamespace(tree=tree_file))
    tree = t.read_tree_file()
    return t.calculate_treeness(tree)
```

### Tree Length

```python
from phykit.services.tree.total_tree_length import TotalTreeLength

def phykit_tree_length(tree_file):
    """Calculate total tree length (sum of all branch lengths).

    Returns: float
    """
    tl = TotalTreeLength(SimpleNamespace(tree=tree_file))
    tree = tl.read_tree_file()
    return tl.calculate_total_tree_length(tree)
```

### Evolutionary Rate

```python
from phykit.services.tree.evolutionary_rate import EvolutionaryRate

def phykit_evolutionary_rate(tree_file):
    """Calculate evolutionary rate (total branch length / number of terminals).

    Represents average substitution rate per lineage.

    Returns: float
    """
    er = EvolutionaryRate(SimpleNamespace(tree=tree_file))
    tree = er.read_tree_file()
    total_bl = tree.total_branch_length()
    num_terminals = tree.count_terminals()
    return total_bl / num_terminals
```

### Degree of Violation of Molecular Clock (DVMC)

```python
from phykit.services.tree.dvmc import DVMC

def phykit_dvmc(tree_file):
    """Calculate Degree of Violation of Molecular Clock.

    DVMC is the standard deviation of root-to-tip distances.
    Lower values indicate tree is more clock-like.

    Returns: float
    """
    d = DVMC(SimpleNamespace(tree=tree_file))
    tree = d.read_tree_file()
    return d.determine_dvmc(tree)
```

### Treeness/RCV

```python
def phykit_treeness_over_rcv(tree_file, aln_file):
    """Calculate treeness/RCV ratio.

    Combines phylogenetic signal (treeness) with compositional uniformity (RCV).
    Higher values indicate better phylogenetic quality.

    Returns: (treeness_over_rcv, treeness, rcv)
    """
    treeness = phykit_treeness(tree_file)
    rcv = phykit_rcv(aln_file)

    if rcv == 0:
        return float('inf'), treeness, rcv

    return treeness / rcv, treeness, rcv
```

---

## Distance-Based Tree Construction

### Neighbor-Joining (NJ)

```python
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

def build_nj_tree(aln_file, model='identity'):
    """Build Neighbor-Joining tree from alignment.

    NJ is fast and reasonably accurate for most cases.

    Args:
        aln_file: path to alignment file
        model: distance model ('identity' or 'blosum62')

    Returns: Biopython Phylo.BaseTree
    """
    alignment, _ = load_alignment(aln_file)
    calculator = DistanceCalculator(model)
    dm = calculator.get_distance(alignment)
    constructor = DistanceTreeConstructor()
    return constructor.nj(dm)
```

### UPGMA

```python
def build_upgma_tree(aln_file, model='identity'):
    """Build UPGMA tree from alignment.

    UPGMA assumes molecular clock. Only use for ultrametric data.

    Args:
        aln_file: path to alignment file
        model: distance model ('identity' or 'blosum62')

    Returns: Biopython Phylo.BaseTree
    """
    alignment, _ = load_alignment(aln_file)
    calculator = DistanceCalculator(model)
    dm = calculator.get_distance(alignment)
    constructor = DistanceTreeConstructor()
    return constructor.upgma(dm)
```

### Saving Trees

```python
# Build tree
tree = build_nj_tree("alignment.fa")

# Save as Newick
Phylo.write(tree, "tree.nwk", "newick")

# Save as Nexus
Phylo.write(tree, "tree.nex", "nexus")
```

---

## Maximum Parsimony Tree Construction

### Simple Parsimony Search

```python
from Bio.Phylo.TreeConstruction import ParsimonyScorer, NNITreeSearcher, ParsimonyTreeConstructor

def build_parsimony_tree(aln_file):
    """Build Maximum Parsimony tree from alignment.

    Uses Nearest Neighbor Interchange (NNI) search.
    Best for small alignments (<50 sequences).

    Returns: Biopython Phylo.BaseTree
    """
    alignment, _ = load_alignment(aln_file)

    # Start with NJ tree
    calculator = DistanceCalculator('identity')
    dm = calculator.get_distance(alignment)
    constructor = DistanceTreeConstructor()
    starting_tree = constructor.nj(dm)

    # Parsimony search
    scorer = ParsimonyScorer()
    searcher = NNITreeSearcher(scorer)
    pars_constructor = ParsimonyTreeConstructor(searcher, starting_tree)

    return pars_constructor.build_tree(alignment)
```

### Getting Parsimony Score

```python
from Bio.Phylo.TreeConstruction import ParsimonyScorer

def get_parsimony_score(tree, alignment):
    """Calculate parsimony score for a tree.

    Lower scores indicate fewer required substitutions.

    Returns: int (total substitutions)
    """
    scorer = ParsimonyScorer()
    return scorer.get_score(tree, alignment)
```

---

## Bootstrap Analysis

### Bootstrap Support

```python
from Bio.Phylo.Consensus import bootstrap_consensus, bootstrap_trees, majority_consensus

def bootstrap_analysis(aln_file, n_replicates=100, model='identity'):
    """Perform bootstrap analysis on alignment.

    Generates bootstrap replicates and builds consensus tree.

    Args:
        aln_file: path to alignment
        n_replicates: number of bootstrap replicates (default 100)
        model: distance model

    Returns: consensus tree with support values
    """
    alignment, _ = load_alignment(aln_file)
    calculator = DistanceCalculator(model)
    constructor = DistanceTreeConstructor(calculator, 'nj')

    # Generate bootstrap trees
    trees = list(bootstrap_trees(alignment, n_replicates, constructor))

    # Build majority consensus
    consensus = majority_consensus(trees, cutoff=0.5)

    return consensus
```

### Extracting Bootstrap Support Values

```python
def extract_bootstrap_support(tree_file):
    """Extract bootstrap support values from internal nodes.

    Returns: dict with support statistics
    """
    tree, _ = load_tree(tree_file)

    supports = []
    for clade in tree.get_nonterminals():
        if clade.confidence is not None:
            supports.append(clade.confidence)

    if not supports:
        return {'supports': [], 'mean': None, 'median': None}

    return {
        'supports': supports,
        'mean': float(np.mean(supports)),
        'median': float(np.median(supports)),
        'min': float(np.min(supports)),
        'max': float(np.max(supports)),
        'n_nodes': len(supports),
        'above_70': sum(1 for s in supports if s >= 70),
        'above_90': sum(1 for s in supports if s >= 90),
    }
```

---

## Branch Length Analysis

### Branch Statistics

```python
def tree_branch_stats(tree_file):
    """Compute branch length statistics from a tree.

    Returns: dict with branch statistics
    """
    tree, _ = load_tree(tree_file)

    internal_lengths = []
    terminal_lengths = []

    for clade in tree.find_clades():
        if clade.branch_length is not None:
            if clade.is_terminal():
                terminal_lengths.append(clade.branch_length)
            else:
                internal_lengths.append(clade.branch_length)

    all_lengths = internal_lengths + terminal_lengths

    return {
        'total_length': sum(all_lengths),
        'n_internal': len(internal_lengths),
        'n_terminal': len(terminal_lengths),
        'internal_sum': sum(internal_lengths),
        'terminal_sum': sum(terminal_lengths),
        'mean_branch': np.mean(all_lengths) if all_lengths else 0,
        'max_branch': max(all_lengths) if all_lengths else 0,
        'min_branch': min(all_lengths) if all_lengths else 0,
    }
```

### Example

```python
stats = tree_branch_stats("tree.nwk")

print(f"Total length: {stats['total_length']:.4f}")
print(f"Internal branches: {stats['n_internal']}")
print(f"Terminal branches: {stats['n_terminal']}")
print(f"Mean branch length: {stats['mean_branch']:.4f}")
```

---

## Tree Comparison

### Robinson-Foulds Distance

```python
import dendropy
from dendropy.calculate import treecompare

def robinson_foulds_distance(tree_file1, tree_file2):
    """Calculate Robinson-Foulds distance between two trees.

    RF distance measures topological difference between trees.
    0 = identical topology, higher = more different.

    Returns: int (number of differing splits)
    """
    tree1 = dendropy.Tree.get(path=tree_file1, schema="newick")
    tree2 = dendropy.Tree.get(
        path=tree_file2,
        schema="newick",
        taxon_namespace=tree1.taxon_namespace  # Must share namespace
    )

    return treecompare.symmetric_difference(tree1, tree2)
```

### Weighted Robinson-Foulds

```python
def weighted_robinson_foulds(tree_file1, tree_file2):
    """Calculate weighted RF distance (includes branch lengths)."""
    tree1 = dendropy.Tree.get(path=tree_file1, schema="newick")
    tree2 = dendropy.Tree.get(
        path=tree_file2,
        schema="newick",
        taxon_namespace=tree1.taxon_namespace
    )

    return treecompare.weighted_robinson_foulds_distance(tree1, tree2)
```

---

## Batch Processing

### Tree Metric Batching

```python
def batch_treeness(gene_files):
    """Compute treeness for all genes."""
    return batch_compute_metric(gene_files, phykit_treeness, requires='tree')

def batch_tree_length(gene_files):
    """Compute tree length for all genes."""
    return batch_compute_metric(gene_files, phykit_tree_length, requires='tree')

def batch_evolutionary_rate(gene_files):
    """Compute evolutionary rate for all genes."""
    return batch_compute_metric(gene_files, phykit_evolutionary_rate, requires='tree')

def batch_dvmc(gene_files):
    """Compute DVMC for all genes."""
    return batch_compute_metric(gene_files, phykit_dvmc, requires='tree')

def batch_treeness_over_rcv(gene_files):
    """Compute treeness/RCV for all genes (requires both tree and alignment)."""
    results = {}
    for entry in gene_files:
        if 'tree_file' in entry and 'aln_file' in entry:
            try:
                results[entry['gene_id']] = phykit_treeness_over_rcv(
                    entry['tree_file'], entry['aln_file']
                )
            except Exception:
                pass
    return results
```

---

## Decision Guide: Which Tree Method?

### For Exploratory Analysis

**Use Neighbor-Joining (NJ)**:
- Fast and reasonably accurate
- No assumptions about clock-like evolution
- Good for initial tree estimation

```python
tree = build_nj_tree("alignment.fa")
```

### For Clock-Like Data

**Use UPGMA**:
- Assumes molecular clock (equal evolutionary rates)
- Faster than NJ
- Only if data truly clock-like

```python
tree = build_upgma_tree("alignment.fa")
```

### For Small Datasets

**Use Maximum Parsimony**:
- Good for discrete characters
- Interpretable (minimum substitutions)
- Slow for large datasets

```python
tree = build_parsimony_tree("alignment.fa")
```

### For Production Analysis

**Use External Tools**:
- **IQ-TREE**: Fast Maximum Likelihood, automatic model selection
- **RAxML**: Maximum Likelihood, bootstrap support
- **PhyML**: Fast ML for moderate datasets
- **MrBayes**: Bayesian inference

Then analyze resulting trees with this skill:
```python
# After running IQ-TREE externally
treeness = phykit_treeness("iqtree_output.treefile")
dvmc = phykit_dvmc("iqtree_output.treefile")
```

---

## Common Patterns

### Pattern: Median Treeness

```python
# Question: "What is the median treeness for fungi?"

# 1. Discover and compute
fungi_genes = discover_gene_files("data/fungi")
fungi_treeness = batch_treeness(fungi_genes)

# 2. Extract values
treeness_values = list(fungi_treeness.values())

# 3. Compute median
median_treeness = np.median(treeness_values)

print(f"Median treeness: {median_treeness:.4f}")
```

### Pattern: Percentage Above Threshold

```python
# Question: "What percentage of trees have treeness > 0.5?"

treeness_values = list(batch_treeness(gene_files).values())
above_threshold = sum(1 for v in treeness_values if v > 0.5)
percentage = (above_threshold / len(treeness_values)) * 100

print(f"Percentage > 0.5: {percentage:.2f}%")
```

### Pattern: Fold Change

```python
# Question: "What is the fold change in median tree length (fungi / animals)?"

fungi_lengths = list(batch_tree_length(fungi_genes).values())
animal_lengths = list(batch_tree_length(animal_genes).values())

fold_change = np.median(fungi_lengths) / np.median(animal_lengths)

print(f"Fold change: {fold_change:.4f}")
```

---

## Troubleshooting

### "Cannot parse tree" Error

**Cause**: Format detection failed or malformed Newick.

**Solution**: Check Newick syntax:
```python
# Valid Newick must end with semicolon
# (A:0.1,B:0.2,(C:0.3,D:0.4):0.5);
```

### Missing Branch Lengths

**Issue**: Tree has no branch lengths for PhyKIT metrics.

**Solution**: PhyKIT requires branch lengths. Rebuild tree with distances:
```python
tree = build_nj_tree("alignment.fa")  # NJ always has branch lengths
```

### Bootstrap Values Not Found

**Issue**: `extract_bootstrap_support()` returns empty list.

**Solution**: Check if tree has confidence values:
```python
tree, _ = load_tree("tree.nwk")
for clade in tree.get_nonterminals():
    print(f"Confidence: {clade.confidence}")
```

If None, tree doesn't have bootstrap support. Run bootstrap analysis:
```python
consensus = bootstrap_analysis("alignment.fa", n_replicates=100)
```

---

## See Also

- `sequence_alignment.md` - Alignment analysis
- `parsimony_analysis.md` - Statistical workflows
- `scripts/tree_statistics.py` - Implementation code