tooluniverse-phylogenetics

# Phylogenetics and Sequence Analysis

Comprehensive phylogenetics and sequence analysis using PhyKIT, Biopython, and DendroPy. Designed for bioinformatics questions about multiple sequence alignments, phylogenetic trees, parsimony, molecular evolution, and comparative genomics.

**IMPORTANT**: This skill handles complex phylogenetic workflows. Most implementation details have been moved to `references/` for progressive disclosure. This document focuses on high-level decision-making and workflow orchestration.

---

## When to Use This Skill

Apply when users:
- Have FASTA alignment files and ask about parsimony informative sites, gaps, or alignment quality
- Have Newick tree files and ask about treeness, tree length, evolutionary rate, or DVMC
- Ask about treeness/RCV, RCV, or relative composition variability
- Need to compare phylogenetic metrics between groups (fungi vs animals, etc.)
- Ask about PhyKIT functions (treeness, rcv, dvmc, evo_rate, parsimony_informative, tree_length)
- Have gene family data with paired alignments and trees
- Need Mann-Whitney U tests or other statistical comparisons of phylogenetic metrics
- Ask about bootstrap support, branch lengths, or tree topology
- Need to build trees (NJ, UPGMA, parsimony) from alignments
- Ask about Robinson-Foulds distance or tree comparison

**BixBench Coverage**: 33 questions across 8 projects (bix-4, bix-11, bix-12, bix-25, bix-35, bix-38, bix-45, bix-60)

**NOT for** (use other skills instead):
- Multiple sequence alignment generation → Use external tools (MUSCLE, MAFFT, ClustalW)
- Maximum Likelihood tree construction → Use IQ-TREE, RAxML, or PhyML
- Bayesian phylogenetics → Use MrBayes or BEAST
- Ancestral state reconstruction → Use separate tools

---

## Core Principles

1. **Data-first approach** - Discover and validate all input files (alignments, trees) before any analysis
2. **PhyKIT-compatible** - Use PhyKIT functions for treeness, RCV, DVMC, parsimony, evolutionary rate (matches BixBench expected outputs)
3. **Format-flexible** - Support FASTA, PHYLIP, Nexus, Newick, and auto-detect formats
4. **Batch processing** - Process hundreds of gene alignments/trees in a single analysis
5. **Statistical rigor** - Mann-Whitney U, medians, percentiles, standard deviations with scipy.stats
6. **Precision awareness** - Match rounding to 4 decimal places (PhyKIT default) or as requested
7. **Group comparison** - Compare metrics between taxa groups (e.g., fungi vs animals)
8. **Question-driven** - Parse exactly what is asked and return the specific number/statistic

---

## Required Python Packages

```python
# Core (MUST be installed)
import numpy as np
import pandas as pd
from scipy import stats
from Bio import AlignIO, Phylo, SeqIO
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

# PhyKIT (primary computation engine)
from phykit.services.tree.treeness import Treeness
from phykit.services.tree.total_tree_length import TotalTreeLength
from phykit.services.tree.evolutionary_rate import EvolutionaryRate
from phykit.services.tree.dvmc import DVMC
from phykit.services.tree.treeness_over_rcv import TreenessOverRCV
from phykit.services.alignment.parsimony_informative_sites import ParsimonyInformative
from phykit.services.alignment.rcv import RelativeCompositionVariability

# DendroPy (for advanced tree operations)
import dendropy

# ToolUniverse (for sequence retrieval)
from tooluniverse import ToolUniverse
```

**Installation**:
```bash
pip install phykit dendropy biopython pandas numpy scipy
```

---

## High-Level Workflow Decision Tree

```
START: User question about phylogenetic data
│
├─ Q1: What type of analysis is needed?
│  │
│  ├─ ALIGNMENT ANALYSIS (FASTA/PHYLIP files)
│  │  ├─ Parsimony informative sites → phykit_parsimony_informative()
│  │  ├─ RCV score → phykit_rcv()
│  │  ├─ Gap percentage → alignment_gap_percentage()
│  │  ├─ GC content → alignment_statistics()
│  │  └─ See: references/sequence_alignment.md
│  │
│  ├─ TREE ANALYSIS (Newick files)
│  │  ├─ Treeness → phykit_treeness()
│  │  ├─ Tree length → phykit_tree_length()
│  │  ├─ Evolutionary rate → phykit_evolutionary_rate()
│  │  ├─ DVMC → phykit_dvmc()
│  │  ├─ Bootstrap support → extract_bootstrap_support()
│  │  └─ See: references/tree_building.md
│  │
│  ├─ COMBINED ANALYSIS (alignment + tree)
│  │  └─ Treeness/RCV → phykit_treeness_over_rcv()
│  │
│  ├─ TREE CONSTRUCTION (build from alignment)
│  │  ├─ Neighbor-Joining → build_nj_tree()
│  │  ├─ UPGMA → build_upgma_tree()
│  │  ├─ Parsimony → build_parsimony_tree()
│  │  └─ See: references/tree_building.md
│  │
│  ├─ GROUP COMPARISON (fungi vs animals, etc.)
│  │  ├─ Batch compute metrics per group
│  │  ├─ Mann-Whitney U test
│  │  ├─ Summary statistics (median, mean, percentiles)
│  │  └─ See: references/parsimony_analysis.md
│  │
│  └─ TREE COMPARISON
│     ├─ Robinson-Foulds distance → robinson_foulds_distance()
│     └─ Bootstrap consensus → bootstrap_analysis()
│
├─ Q2: What data format is available?
│  ├─ FASTA (.fa, .fasta, .faa, .fna)
│  ├─ PHYLIP (.phy, .phylip) - Use phylip-relaxed for long names
│  ├─ Nexus (.nex, .nexus)
│  ├─ Newick (.nwk, .newick, .tre, .tree)
│  └─ Auto-detect with load_alignment() or load_tree()
│
└─ Q3: Is this a batch analysis?
   ├─ Single gene → Run metric function once
   ├─ Multiple genes → Use batch_compute_metric()
   └─ Group comparison → Use discover_gene_files() + compare_groups()
```

---

## Quick Reference: Common Metrics

| Metric | Function | Input | Description |
|--------|----------|-------|-------------|
| **Treeness** | `phykit_treeness(tree_file)` | Newick | Internal branch length / Total branch length |
| **RCV** | `phykit_rcv(aln_file)` | FASTA/PHYLIP | Relative Composition Variability |
| **Treeness/RCV** | `phykit_treeness_over_rcv(tree, aln)` | Both | Treeness divided by RCV |
| **Tree Length** | `phykit_tree_length(tree_file)` | Newick | Sum of all branch lengths |
| **Evolutionary Rate** | `phykit_evolutionary_rate(tree_file)` | Newick | Total branch length / num terminals |
| **DVMC** | `phykit_dvmc(tree_file)` | Newick | Degree of Violation of Molecular Clock |
| **Parsimony Sites** | `phykit_parsimony_informative(aln_file)` | FASTA/PHYLIP | Sites with ≥2 chars appearing ≥2 times |
| **Gap Percentage** | `alignment_gap_percentage(aln_file)` | FASTA/PHYLIP | Percentage of gap characters |

See `scripts/tree_statistics.py` for implementation.

---

## Common Analysis Patterns (BixBench)

### Pattern 1: Single Metric Across Groups

**Question**: "What is the median DVMC for fungi vs animals?"

**Workflow**:
```python
# 1. Discover files
fungi_genes = discover_gene_files("data/fungi")
animal_genes = discover_gene_files("data/animals")

# 2. Compute metric
fungi_dvmc = batch_dvmc(fungi_genes)
animal_dvmc = batch_dvmc(animal_genes)

# 3. Compare
fungi_values = list(fungi_dvmc.values())
animal_values = list(animal_dvmc.values())

print(f"Fungi median DVMC: {np.median(fungi_values):.4f}")
print(f"Animal median DVMC: {np.median(animal_values):.4f}")
```

**See**: `references/parsimony_analysis.md` for full implementation

### Pattern 2: Statistical Comparison

**Question**: "What is the Mann-Whitney U statistic comparing treeness between groups?"

**Workflow**:
```python
from scipy import stats

# Compute treeness for both groups
group1_treeness = batch_treeness(group1_genes)
group2_treeness = batch_treeness(group2_genes)

# Mann-Whitney U test (two-sided)
u_stat, p_value = stats.mannwhitneyu(
    list(group1_treeness.values()),
    list(group2_treeness.values()),
    alternative='two-sided'
)

print(f"U statistic: {u_stat:.0f}")
print(f"P-value: {p_value:.4e}")
```

### Pattern 3: Filtering + Metric

**Question**: "What is the treeness/RCV for alignments with <5% gaps?"

**Workflow**:
```python
# 1. Filter by gap percentage
valid_genes = []
for entry in gene_files:
    if 'aln_file' in entry:
        gap_pct = alignment_gap_percentage(entry['aln_file'])
        if gap_pct < 5.0:
            valid_genes.append(entry)

# 2. Compute metric on filtered set
results = batch_treeness_over_rcv(valid_genes)

# 3. Report
values = [r[0] for r in results.values()]  # treeness/rcv ratio
print(f"Median treeness/RCV: {np.median(values):.4f}")
```

### Pattern 4: Specific Gene Lookup

**Question**: "What is the evolutionary rate for gene X?"

**Workflow**:
```python
# Find gene file
gene_files = discover_gene_files("data/")
gene_entry = [g for g in gene_files if g['gene_id'] == 'X'][0]

# Compute metric
evo_rate = phykit_evolutionary_rate(gene_entry['tree_file'])

print(f"Evolutionary rate for gene X: {evo_rate:.4f}")
```

---

## Choosing Methods: When to Use What

### Alignment Methods

**When building alignments** (use external tools, not this skill):

| Method | Speed | Accuracy | Use Case |
|--------|-------|----------|----------|
| **ClustalW** | Slow | Medium | Small datasets (<100 sequences), educational |
| **MUSCLE** | Fast | High | Medium datasets (100-1000 sequences) |
| **MAFFT** | Very Fast | Very High | **Recommended** - Large datasets (>1000 sequences) |

**For this skill**: Work with pre-aligned sequences. Use `load_alignment()` to read any format.

### Tree Building Methods

**When to use which tree method:**

| Method | Speed | Accuracy | Use Case |
|--------|-------|----------|----------|
| **Neighbor-Joining** | Fast | Medium | Quick trees, large datasets, exploratory |
| **UPGMA** | Fast | Low | Assumes molecular clock, special cases only |
| **Maximum Parsimony** | Medium | Medium | Small datasets, discrete characters |
| **Maximum Likelihood** | Slow | High | **Use external tools** (IQ-TREE, RAxML) for production |

**Implementation in this skill**:
```python
# Fast distance-based trees
tree = build_nj_tree("alignment.fa")  # Neighbor-Joining
tree = build_upgma_tree("alignment.fa")  # UPGMA

# Parsimony (for small alignments)
tree = build_parsimony_tree("alignment.fa")
```

**For production ML trees**: Use IQ-TREE or RAxML externally, then analyze with this skill.

See `references/tree_building.md` for detailed implementations.

---

## Batch Processing

### Discovering Gene Files

```python
# Auto-discover paired alignment + tree files
gene_files = discover_gene_files("data/")

# Result: list of dicts with 'gene_id', 'aln_file', 'tree_file'
# [
#   {'gene_id': 'gene1', 'aln_file': 'gene1.fa', 'tree_file': 'gene1.nwk'},
#   {'gene_id': 'gene2', 'aln_file': 'gene2.fa', 'tree_file': 'gene2.nwk'},
#   ...
# ]
```

### Computing Metrics in Batch

```python
# Tree metrics
treeness_results = batch_treeness(gene_files)
tree_length_results = batch_tree_length(gene_files)
dvmc_results = batch_dvmc(gene_files)
evo_rate_results = batch_evolutionary_rate(gene_files)

# Alignment metrics
rcv_results = batch_rcv(gene_files)
pi_results = batch_parsimony_informative(gene_files)
gap_results = batch_gap_percentage(gene_files)

# Combined metrics
treeness_rcv_results = batch_treeness_over_rcv(gene_files)

# All return dict: {gene_id: value}
```

### Statistical Analysis

```python
# Summary statistics
stats = summary_stats(list(treeness_results.values()))
# Returns: {'mean': ..., 'median': ..., 'std': ..., 'min': ..., 'max': ...}

# Group comparison
comparison = compare_groups(
    list(fungi_treeness.values()),
    list(animal_treeness.values()),
    group1_name="Fungi",
    group2_name="Animals"
)
# Returns: {'u_statistic': ..., 'p_value': ..., 'Fungi': {...}, 'Animals': {...}}
```

See `references/parsimony_analysis.md` for full workflow.

---

## Answer Extraction for BixBench

| Question Pattern | Extraction Method |
|-----------------|-------------------|
| "What is the median X?" | `np.median(values)` |
| "What is the maximum X?" | `np.max(values)` |
| "What is the difference between median X for A vs B?" | `abs(np.median(a) - np.median(b))` |
| "What percentage of X have Y above Z?" | `sum(v > Z for v in values) / len(values) * 100` |
| "What is the Mann-Whitney U statistic?" | `stats.mannwhitneyu(a, b)[0]` |
| "What is the p-value?" | `stats.mannwhitneyu(a, b)[1]` |
| "What is the X value for gene Y?" | `results[gene_id]` |
| "What is the fold-change in median X?" | `np.median(a) / np.median(b)` |
| "multiplied by 1000" | `round(value * 1000)` |

### Rounding Rules

- **PhyKIT default**: 4 decimal places
- **Percentages**: Match question format (e.g., "35%" → integer, "3.5%" → 1 decimal)
- **P-values**: Scientific notation for very small values
- **U statistics**: Integer (no decimals)
- **Always check question wording**: "rounded to 3 decimal places" overrides defaults

---

## BixBench Question Coverage

| Project | Questions | Metrics |
|---------|-----------|---------|
| **bix-4** | 7 | DVMC analysis (fungi vs animals) |
| **bix-11** | 6 | Treeness analysis (median, percentages, Mann-Whitney U) |
| **bix-12** | 5 | Parsimony informative sites (counts, percentages, ratios) |
| **bix-25** | 2 | Treeness/RCV with gap filtering |
| **bix-35** | 4 | Evolutionary rate (specific genes, comparisons) |
| **bix-38** | 5 | Tree length (fold-change, variance, paired ratios) |
| **bix-45** | 4 | RCV (Mann-Whitney U, medians, paired differences) |
| **bix-60** | 1 | Average treeness across multiple trees |

---

## ToolUniverse Integration

### Sequence Retrieval

```python
from tooluniverse import ToolUniverse

tu = ToolUniverse()
tu.load_tools()

# Get sequences from NCBI
result = tu.tools.NCBI_get_sequence(accession="NP_000546")

# Get gene tree from Ensembl
tree_result = tu.tools.EnsemblCompara_get_gene_tree(gene="ENSG00000141510")

# Get species tree from OpenTree
tree_result = tu.tools.OpenTree_get_induced_subtree(ott_ids="770315,770319")
```

---

## File Structure

```
tooluniverse-phylogenetics/
├── SKILL.md                           # This file (workflow orchestration)
├── QUICK_START.md                     # Quick reference
├── test_phylogenetics.py             # Comprehensive test suite
├── references/
│   ├── sequence_alignment.md         # Alignment analysis details
│   ├── tree_building.md              # Tree construction methods
│   ├── parsimony_analysis.md         # Statistical comparison workflows
│   └── troubleshooting.md            # Common issues and solutions
└── scripts/
    ├── format_alignment.py           # Alignment format conversion
    └── tree_statistics.py            # Core metric implementations
```

---

## Completeness Checklist

Before returning your answer, verify:

- [ ] Identified all input files (alignments and/or trees)
- [ ] Detected group structure (fungi/animals/etc.) if applicable
- [ ] Used correct PhyKIT function for the requested metric
- [ ] Processed ALL genes in each group (not just a sample)
- [ ] Applied correct statistical test if comparison requested
- [ ] Used correct rounding (4 decimals default, or as specified)
- [ ] Returned the specific statistic asked for (median, max, U stat, p-value, etc.)
- [ ] For percentage questions, confirmed whether answer is integer or decimal
- [ ] For "difference" questions, confirmed direction (A - B vs abs difference)
- [ ] For Mann-Whitney U, used `alternative='two-sided'` (default in scipy)

---

## Next Steps

- For detailed alignment analysis workflows → See `references/sequence_alignment.md`
- For tree construction methods → See `references/tree_building.md`
- For statistical comparison examples → See `references/parsimony_analysis.md`
- For common errors and solutions → See `references/troubleshooting.md`
- For script implementations → See `scripts/tree_statistics.py`

---

## Support

For issues with:
- **PhyKIT functions**: Check PhyKIT documentation at https://jlsteenwyk.com/PhyKIT/
- **Biopython tree/alignment parsing**: See https://biopython.org/wiki/Phylo
- **DendroPy operations**: See https://dendropy.org/
- **ToolUniverse integration**: Check ToolUniverse documentation

## License

Same as ToolUniverse framework license.