Phylogenetics and Sequence Analysis

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

# 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:

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

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:

# 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:

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:

# 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:

# 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):

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

Tree Building Methods

When to use which tree method:

Implementation in this skill:

# 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

# 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

# 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

# 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

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

ToolUniverse Integration

Sequence Retrieval

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.