Tooluniverse Phylogenetics

Phylogenetics and Sequence Analysis

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• `scogs_fungi.zip` / `scogs_animals.zip` (BUSCO single-copy ortholog phylogenetics) → these contain the pre-computed alignments (*.faa.mafft.clipkit) and trees (*.faa.mafft.clipkit.treefile) from the original analysis. Use these directly with PhyKIT (see "BUSCO scogs questions" below). Re-running BUSCO → MAFFT → IQ-TREE from *.busco.zip files takes 1–6 hours AND produces slightly different numbers due to seed/version drift.
• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *tree*, *phykit*, *saturation*, *treeness*) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5–10× turn count).

---

BUSCO scogs questions (multi-species phylogenomics)

data folders with scogs_fungi.zip and/or scogs_animals.zip ship pre-computed per-ortholog alignments (and sometimes trees). The question asks for a metric per group, or a Mann-Whitney U / median / ratio comparison between groups.

PRIMARY SCRIPT — both groups in one pass (use this FIRST)

When the question compares animals vs fungi (Mann-Whitney U, ratio, fold-change, paired difference), the bundled paired-comparison script extracts both zips, computes the metric per ortholog for each group, and emits ALL of: per-group summary, two-tailed Mann-Whitney U + p-value (in both orderings since U is asymmetric), paired-ortholog median diff, paired-ortholog median ratio, group-median ratio, and lowest-non-zero ratios — in one run, no aggregation step needed:

python skills/tooluniverse-phylogenetics/scripts/scogs_paired_compare.py \
    --data-folder "$DATA_PATH" --metric parsimony_informative
# Metrics: parsimony_informative, rcv, gap_percentage (alignment-only,
# Biopython-fast: ~2s for 500 alignments);
# treeness, dvmc, total_tree_length, evolutionary_rate, long_branch_score,
# patristic_distances (tree); treeness_over_rcv, saturation (both).

Output blocks (parse in Python or grep):

# SUMMARY group=animals: n=... mean=... median=... min=... max=... p25=... p75=... lowest_nonzero=... n_nonzero=...
# SUMMARY group=fungi:   n=... mean=... median=... min=... max=... p25=... p75=... lowest_nonzero=... n_nonzero=...
# MWU animals_vs_fungi: U=... p=...
# MWU fungi_vs_animals: U=... p=...        <-- U(a,b) + U(b,a) = n_a*n_b
# PAIRED n_common=N: median_diff(animals-fungi)=...  median_diff(fungi-animals)=...
# PAIRED RATIO median(animals/fungi)=... (n=...)    <-- for each common ortholog: a_val/b_val, then median
# PAIRED RATIO median(fungi/animals)=... (n=...)
# GROUP_MEDIAN_RATIO animals/fungi=...               <-- median(group_a) / median(group_b)
# GROUP_MEDIAN_RATIO fungi/animals=...
# GROUP_MEDIAN_DIFF animals-fungi=...
# LOWEST_NONZERO animals=... fungi=...
# LOWEST_NONZERO_RATIO animals/fungi=...
# LOWEST_NONZERO_RATIO fungi/animals=...

For long_branch_score and patristic_distances (multi-value-per-tree metrics), pass --per-tree-stat mean or --per-tree-stat median to choose the per-tree summary BEFORE the cross-tree MWU. The question wording "comparing median long branch scores" means per-tree summary = median; "comparing mean long branch scores" means per-tree summary = mean. Run TWICE (once with each) if uncertain.

Single-group script (when only one group is asked about)

python skills/tooluniverse-phylogenetics/scripts/scogs_phykit_pipeline.py \
    --data-folder "$DATA_PATH" --group fungi --metric treeness --out /tmp/f.tsv
# Auto-falls-back to .faa.mafft when .faa.mafft.clipkit is absent
# (some scogs zips ship only mafft alignments, not clipkit trims).

phykit parsimony_informative is NOT a valid CLI subcommand

PhyKIT's CLI exposes parsimony-informative-site count as parsimony_informative_sites (alias pis). Calling phykit parsimony_informative <file> returns the help banner with non-zero exit and silently produces zero values. The bundled scripts translate parsimony_informative → parsimony_informative_sites automatically. The output is <n_pi>\t<n_total>\t<percent> — column THREE is the percentage that questions usually ask for.

Group-median ratio vs paired ratio (read this carefully)

When a question phrases tree-length / RCV / DVMC comparisons as "ratio of fungal to animal X across orthologs", there are TWO distinct quantities:

1. GROUP_MEDIAN_RATIO = median(values_fungi) / median(values_animals). Use ALL orthologs in each group independently. This is what group-comparison published numbers usually report (n_fungi can differ from n_animals, and "across" is a population statement, not a paired one).

2. PAIRED RATIO median = for each ortholog present in BOTH groups, compute value_fungi / value_animals, then take the median across common orthologs. Smaller denominator (intersection only) and a different number when the groups have different size.

Default to GROUP_MEDIAN_RATIO unless the question explicitly says "matched ortholog", "paired ortholog", "per-ortholog ratio", or "for each ortholog". If the answer phrasing is ambiguous, BOTH numbers are in the script's output — pick the one matching the question's "across" / "paired" / "ratio of medians" phrasing.

Total amino-acid count across single-copy orthologs — single representative, not all species

When a BUSCO single-copy ortholog dataset (single_copy_busco_sequences/) is present and the question asks "how many total amino acids are present in all single-copy ortholog sequences", count one representative sequence per ortholog, not the sum across all species/copies.

Each <ortholog_id>.faa in single_copy_busco_sequences/ typically contains multiple species' copies of that ortholog (one each). Summing every sequence across every species double/triple/N-fold counts each ortholog by the species count and gives n_species × correct_answer.

❌ WRONG: for f in *.faa: sum(len(rec.seq) for rec in SeqIO.parse(f, 'fasta')) then sum across files (multi-species sum)

✅ RIGHT: for f in *.faa: first_rec = next(SeqIO.parse(f, 'fasta')); total += len(first_rec.seq) (one representative per ortholog)

If your answer is n_species × GT (e.g. 32228 when GT looks like 13809 = 32228/2.33 ≈ 8 species × representative), you summed all species — re-do with one representative.

Lowest-non-zero ratios

For metrics that can legitimately equal 0 for highly conserved or very short alignments (parsimony informative %, RCV on near-identical seqs), "lowest" in a question typically means "lowest non-zero". The paired script emits LOWEST_NONZERO_RATIO for both orderings — use that line when the raw min in a group is 0.

File-layout fallback (alignment naming)

scogs zips ship in two shapes:

• Full: <gene>.faa, <gene>.faa.mafft, <gene>.faa.mafft.clipkit,

<gene>.faa.mafft.clipkit.treefile, plus iqtree/bionj/log/mldist. Use clipkit alignment + treefile for tree-paired metrics.

• Alignment-only: just <gene>.faa + <gene>.faa.mafft. No

trees, no clipkit. Used for parsimony, RCV, gap-percentage questions. Use the .faa.mafft (NOT raw .faa) — the published metric was computed on the MAFFT-aligned file.

Both bundled scripts auto-detect the layout and use the best available alignment per ortholog. Do NOT re-run MAFFT or ClipKit yourself; the shipped files are canonical.

Anti-pattern: running phykit on the raw *.busco.zip extracted ortholog FASTAs and aligning/tree-building yourself. The pre-computed files in scogs_*.zip are the canonical inputs.

---

PhyKIT, Biopython, and DendroPy for alignment/tree analysis, evolutionary metrics, and comparative genomics.

LOOK UP, DON'T GUESS

When uncertain about any scientific fact, SEARCH databases first.

---

When to Use

FASTA/PHYLIP/Nexus/Newick files; treeness, RCV, DVMC, evolutionary rate, parsimony sites, tree length, bootstrap; group comparisons (Mann-Whitney U); tree construction (NJ/UPGMA/parsimony); Robinson-Foulds distance.

De novo alignment / tree from your own sequences: to align raw sequences (not pre-computed files), call EBI_msa_align (Clustal Omega / MUSCLE / MAFFT / Kalign / T-Coffee via EMBL-EBI), then pass its data.aligned_fasta string as the aligned_sequences argument of EBI_build_phylogenetic_tree (note the arg name differs from the output key) for a neighbour-joining or UPGMA tree (Newick). Feed that Newick / alignment straight into the PhyKIT metrics below.

Still NOT for: maximum-likelihood trees (IQ-TREE/RAxML) or Bayesian inference (MrBayes/BEAST) — EBI_build_phylogenetic_tree only does distance-based NJ/UPGMA. For publication ML/Bayesian phylogenies, run dedicated tooling; use the pre-computed scogs_* trees when available.

---

Required Packages

import numpy as np, pandas as pd
from scipy import stats
from Bio import AlignIO, Phylo, SeqIO
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
import dendropy

---

Workflow Decision Tree

ALIGNMENT ANALYSIS (FASTA/PHYLIP):
  Parsimony sites → phykit_parsimony_informative()
  RCV → phykit_rcv()
  Gap % → alignment_gap_percentage()

TREE ANALYSIS (Newick):
  Treeness → phykit_treeness()
  Tree length → phykit_tree_length()
  Evolutionary rate → phykit_evolutionary_rate()
  DVMC → phykit_dvmc()
  Bootstrap → extract_bootstrap_support()

COMBINED: Treeness/RCV → phykit_treeness_over_rcv(tree, aln)

TREE CONSTRUCTION: NJ → build_nj_tree(); UPGMA → build_upgma_tree(); Parsimony → build_parsimony_tree()

GROUP COMPARISON: batch metrics → Mann-Whitney U → summary stats

TREE COMPARISON: Robinson-Foulds → robinson_foulds_distance()

---

Quick Reference

---

Common Patterns

Single Metric Across Groups

fungi_dvmc = batch_dvmc(discover_gene_files("data/fungi"))
animal_dvmc = batch_dvmc(discover_gene_files("data/animals"))
print(f"Fungi median: {np.median(list(fungi_dvmc.values())):.4f}")

Statistical Comparison

u_stat, p_value = stats.mannwhitneyu(list(g1.values()), list(g2.values()), alternative='two-sided')

Filtering + Metric

Filter by gap percentage < 5%, then compute treeness/RCV on filtered set.

Batch Processing

gene_files = discover_gene_files("data/")  # → [{gene_id, aln_file, tree_file}]
treeness_results = batch_treeness(gene_files)  # → {gene_id: value}

---

Answer Extraction

Rounding: PhyKIT default 4 decimals. U stats = integer. Question wording overrides.

---

Interpretation

Completeness Checklist

All files identified; group structure detected; correct PhyKIT function; ALL genes processed (not sample); correct test; 4-decimal rounding; specific statistic (median/max/U/p); Mann-Whitney alternative='two-sided'.

---

Analysis conventions

MANDATORY: Use phykit_batch_analysis tool for batch computations

For ANY question asking for statistics across multiple trees/alignments (median treeness, mean saturation, DVMC percentage, gap percentage, long branch scores), use the ToolUniverse tool:

tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"./trees","extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"saturation","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"gap_percentage","directory":"./alignments","extension":".fa"}'

Do NOT run phykit manually in a loop — the tool handles all files and returns correct summary statistics.

PhyKIT column-position cheat sheet (parse output carefully)

When parsing PhyKIT stdout for batch metrics, the column you want depends on the metric:

Rule of thumb: phykit toverr and saturation produce multi-column lines per alignment. Don't grep the value that "looks like the answer" — count columns from the header in phykit <metric> --help. If your batch median is wildly off the published number (e.g., median treeness/RCV ≈ 0.20 when expected ≈ 0.26), you almost certainly picked the wrong column.

Preferred: don't parse phykit output by hand — call the phykit_batch_analysis tool, which already returns the correct column for each metric. Supported function values are treeness, saturation, dvmc, long_branch_score, total_tree_length, parsimony_informative:

tu run phykit_batch_analysis '{"operation":"batch","function":"saturation","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'

For treeness_over_rcv (toverr / treeness/RCV ratio) the tool has no matching function; use the bundled scogs_*.py scripts below, which compute it directly.

Sanity targets for biological scogs trees: median saturation ~0.4–0.7, median treeness/RCV ~0.2–0.4, median treeness ~0.05–0.15. Values an order of magnitude off these mean wrong column.

Bundled script: BUSCO target_orthologs intersection

When the data folder has *.busco.zip files + target_orthologs.txt, use the bundled script — do NOT enumerate single_copy_busco_sequences/*.faa across all zips manually:

python skills/tooluniverse-phylogenetics/scripts/busco_target_orthologs.py \
  --data-folder /path/to/data

The default run prints FIVE summary lines covering every common interpretation of "total amino acids":

# SUMMARY: n_targets=K, n_intersected=N (single-copy in ALL S species), intersected_total_aa=A, sum_all_aa=B
# SUMMARY group=all: intersected n=N total_aa=A, sum_all total_aa=B
# SUMMARY group=animals: sum_all total_aa=X        <-- per-group sum (animal species only)
# SUMMARY group=fungi:   sum_all total_aa=Y        <-- per-group sum (fungal species only)

Picking the right SUMMARY line (read carefully)

Match the question phrasing to the summary line:

Default rule when the data folder contains BOTH animal AND fungal busco zips: published "total amino acids" answers almost always refer to ONE group (the analysis group), NOT the cross-group union. Use group=animals: sum_all or group=fungi: sum_all. Do NOT pick the union number (sum_all_aa) unless the question explicitly says "across all 8 species" or "fungi and animals combined".

The script emits the per-group sums BEFORE the union sum on stdout for this exact reason — read the output line by line and stop at the group=animals / group=fungi line that matches the analysis group implied by the question.

Single-copy orthologs across species — comparison set + intersection

Two-step rule when counting across BUSCO single_copy_busco_sequences/ data:

1. Find the comparison set first. If a target_orthologs.txt (or similar named subset list) exists in the data folder, that file IS the comparison set — restrict to those ortholog IDs only. Do not enumerate every BUSCO single-copy file across species. Do not assume "all" means the whole BUSCO output when a target list is provided.

2. Then apply the intersection rule. "Single-copy ortholog" across species means single-copy in EVERY species in the comparison set. If an ortholog is missing from one species' single_copy_busco_sequences/, exclude it from the count entirely — do not partially count the species that do have it.

Sanity check: if any species shows a much smaller per-ortholog count than others (e.g., one species at ~600 aa while others are 4000+ aa for the same ortholog set), the missing-from-some orthologs are inflating the per-ortholog average — drop them first.

Worked example. data folder has 8 species (4 animal, 4 fungal) *.busco.zip + target_orthologs.txt listing 10 ortholog IDs:

• Wrong: enumerate all single_copy_busco_sequences/*.faa across all 8 species → ≈80 files → sum AA → answer 32228 (treats every per-species copy independently).
• Right: for each of the 10 target IDs, check it appears as single_copy in all 8 species → keep only intersected IDs (often 5/10 — some target IDs are multi-copy/missing in one species) → for kept IDs, sum AA across the 8 species → 13809.
• "5 trees" semantics: when a question says "5 trees" but you find 10 treefiles, the GT used the intersected subset (orthologs single-copy in all species) — not all 10. Re-derive the intersection before averaging.

Process the FULL set, not a sample (batch metrics)

When a question asks for a median/percentile/mean across orthologs, your batch must include EVERY ortholog in the relevant comparison set:

• scogs_fungi.zip ships ~255 fungal alignments+trees; scogs_animals.zip ships ~241. Median computed from a 10-file sample is NOT the published answer.
• For phykit_batch_analysis, always point at the extracted scogs directory containing all per-ortholog files, not a hand-picked subset.
• If your computed RCV/treeness/DVMC median diverges from a sanity-check target by >10%, count files first — you likely processed a subset.

Filter THEN compute (don't compute then filter)

Questions of the form "max X in genes with >70% gaps" require the filter to be applied before the max:

# 1. Compute gap% per alignment
# 2. Keep only alignments with gap% > 70
# 3. Compute treeness/RCV ON THE FILTERED SET
# 4. Take max

Computing the metric across all genes and then taking max returns the global max, which is wrong.

Animals vs fungi — long branch score aggregation

PhyKIT's long_branch_score -v outputs per-taxon LB scores (one row per leaf in the tree). For per-tree summaries: 1. Per-tree: run phykit long_branch_score -v <tree> → list of per-taxon scores. 2. Per-tree summary: collapse to ONE number per tree using either the mean or the median of those per-taxon scores. 3. Per-group summary: aggregate per-tree numbers (median/mean/MWU U + p-value).

Match the per-tree summary to the question phrasing:

The bundled scogs_paired_compare.py --metric long_branch_score --per-tree-stat {mean,median} does steps 1+2 for both groups in one pass and emits the cross-group MWU U + p-value directly.

Common error: averaging the four animal species and four fungal species directly without going through the per-tree step — this conflates species LB and ortholog LB and yields the wrong delta.

Treeness/RCV: use the right input file

phykit toverr (a.k.a. treeness_over_rcv) takes BOTH alignment and tree. Use the trimmed alignment (*.faa.mafft.clipkit) paired with its treefile (*.faa.mafft.clipkit.treefile), not the raw .faa.mafft. ClipKit-trimmed alignments are what produced the canonical tree, so the RCV must be computed on the same trimmed alignment for the ratio to match published numbers.

Parsimony informative sites

• Exclude gap-only columns before counting — a column that is all gaps is not informative.
• A site is parsimony informative when ≥2 different non-gap characters each appear in ≥2 taxa.
• Use Biopython AlignIO or the AMAS tool to iterate columns and count.

Treeness (RCV ratio)

Treeness = sum of internal branch lengths / total tree length. Internal branches are those that do not lead to a leaf (tip).

PhyKIT usage

PhyKIT (pip install phykit) provides command-line functions for tree and alignment statistics. Common functions:

• phykit treeness <tree_file> — outputs treeness (RCV) value
• phykit saturation <alignment_file> -t <tree_file> — outputs saturation value
• phykit dvmc <tree_file> — degree of violation of the molecular clock
• phykit long_branch_score <tree_file> — long-branch score (LBS)
• phykit alignment_length <alignment_file> — alignment length
• phykit parsimony_informative <alignment_file> — count parsimony informative sites

When running PhyKIT on multiple gene trees/alignments, use the bundled batch script:

# Treeness across all trees
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir scogs_fungi --function treeness --ext .treefile --stat median

# Saturation with paired alignment+tree
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir alignments --function saturation --tree-dir trees \
  --ext .fa --tree-ext .treefile --stat median

# Long branch score (mean per tree, then median across trees)
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function long_branch_score --ext .treefile \
  --per-tree-stat mean --stat median

# DVMC
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function dvmc --ext .treefile --stat all

# Gap percentage across all alignments
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir alignments --function gap_percentage --ext .fa

# Evolutionary rate (median across trees)
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function evolutionary_rate --ext .treefile --stat median

# Mean patristic distance per tree, then mean across trees
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function patristic_distances --ext .treefile --stat mean

Preferred: use the `phykit_batch_analysis` ToolUniverse tool instead of running PhyKIT manually:

# Via CLI
tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"/path/to/trees","extension":".treefile"}'

# Via SDK
tu.run_one_function({"name": "phykit_batch_analysis", "arguments": {"operation": "batch", "function": "saturation", "directory": "/path/to/alignments", "extension": ".fa", "tree_directory": "/path/to/trees"}})

# Gap percentage
tu run phykit_batch_analysis '{"operation":"gap_percentage","directory":"/path/to/alignments","extension":".fa"}'

Key rules: 1. Process ALL files — don't stop at a subset. The tool handles this automatically 2. Gap percentage: total gaps / total positions across all alignments (not per-file average) 3. Long branch score: each tree produces per-taxon scores → summarize per tree (mean) → then summarize across trees (median). Use "per_tree_stat":"mean" 4. Fungi vs animal comparisons: match genes by ortholog ID (filename stem), not by file order. Run the tool on each organism's directory separately, then compare medians

References

references/sequence_alignment.md, references/tree_building.md, references/parsimony_analysis.md, scripts/tree_statistics.py

• PhyKIT: https://jlsteenwyk.com/PhyKIT/
• Biopython Phylo: https://biopython.org/wiki/Phylo
• DendroPy: https://dendropy.org/

Quick Start: Phylogenetics Skill

Fast reference for common phylogenetic analysis tasks.

---

scogs animal-vs-fungi questions — START HERE

When the data folder has scogs_animals.zip and/or scogs_fungi.zip, use the bundled paired-comparison script. It computes the metric per ortholog for BOTH groups, plus Mann-Whitney U + p-value, paired ratios, group-median ratios, and lowest-non-zero ratios — all in one shot. Avoids the timeouts that come from re-running per-group batches.

# parsimony_informative %, RCV, gap_percentage: ~2s for 500 alignments
# (Biopython, no subprocess fork-per-file)
python skills/tooluniverse-phylogenetics/scripts/scogs_paired_compare.py \
    --data-folder "$DATA_PATH" --metric parsimony_informative

# Long branch score / patristic distances: choose per-tree summary
python skills/tooluniverse-phylogenetics/scripts/scogs_paired_compare.py \
    --data-folder "$DATA_PATH" --metric long_branch_score --per-tree-stat mean
python skills/tooluniverse-phylogenetics/scripts/scogs_paired_compare.py \
    --data-folder "$DATA_PATH" --metric long_branch_score --per-tree-stat median

The output emits both MWU animals_vs_fungi AND MWU fungi_vs_animals because U is asymmetric. Grep the line whose ordering matches the question wording, e.g. "comparing animals and fungi" → use MWU animals_vs_fungi, "comparing fungi vs animals" → use the other.

For "ratio" questions: GROUP_MEDIAN_RATIO uses all orthologs (median / median), PAIRED_RATIO uses common-ortholog ratios then median. The script prints both — read the question for "matched" / "paired" / "per-ortholog" cues; default to GROUP_MEDIAN_RATIO otherwise.

---

Installation

pip install phykit dendropy biopython pandas numpy scipy

---

Common Tasks

1. Compute Single Tree Metric

from scripts.tree_statistics import phykit_treeness, phykit_tree_length

# Treeness
treeness = phykit_treeness("tree.nwk")
print(f"Treeness: {treeness:.4f}")

# Tree length
length = phykit_tree_length("tree.nwk")
print(f"Tree length: {length:.4f}")

2. Compute Single Alignment Metric

from scripts.tree_statistics import phykit_parsimony_informative, phykit_rcv

# Parsimony informative sites
pi_count, aln_len, pi_pct = phykit_parsimony_informative("alignment.fa")
print(f"PI sites: {pi_count} / {aln_len} ({pi_pct:.2f}%)")

# RCV
rcv = phykit_rcv("alignment.fa")
print(f"RCV: {rcv:.4f}")

3. Batch Analysis

from scripts.format_alignment import discover_gene_files
from scripts.tree_statistics import batch_treeness, batch_dvmc
import numpy as np

# Discover files
gene_files = discover_gene_files("data/")
print(f"Found {len(gene_files)} genes")

# Compute metric for all genes
treeness_results = batch_treeness(gene_files)
dvmc_results = batch_dvmc(gene_files)

# Get median
print(f"Median treeness: {np.median(list(treeness_results.values())):.4f}")
print(f"Median DVMC: {np.median(list(dvmc_results.values())):.4f}")

4. Compare Two Groups

from scipy import stats

# Get metrics for both groups
fungi_genes = discover_gene_files("data/", group_name="fungi")
animal_genes = discover_gene_files("data/", group_name="animals")

fungi_treeness = batch_treeness(fungi_genes)
animal_treeness = batch_treeness(animal_genes)

# Mann-Whitney U test
u_stat, p_value = stats.mannwhitneyu(
    list(fungi_treeness.values()),
    list(animal_treeness.values())
)

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

5. Filter and Analyze

from scripts.format_alignment import filter_by_gap_threshold
from scripts.tree_statistics import batch_treeness_over_rcv

# Filter by gap percentage
valid_genes = filter_by_gap_threshold(gene_files, max_gap_pct=5.0)
print(f"Genes with <5% gaps: {len(valid_genes)}")

# Compute treeness/RCV for valid genes
results = batch_treeness_over_rcv(valid_genes)

# Extract ratios
ratios = [r[0] for r in results.values()]
print(f"Median treeness/RCV: {np.median(ratios):.4f}")

6. Convert Alignment Format

# Single file
python scripts/format_alignment.py convert input.phy output.fa --format fasta

# Batch convert
python scripts/format_alignment.py batch-convert input_dir/ output_dir/ --format phylip-relaxed

---

Quick Reference: Metrics

---

Directory Structure

data/
├── fungi/
│   ├── gene1.fa
│   ├── gene1.nwk
│   ├── gene2.fa
│   └── gene2.nwk
└── animals/
    ├── gene1.fa
    ├── gene1.nwk
    ├── gene2.fa
    └── gene2.nwk

---

Common Workflows

Workflow 1: DVMC Comparison

# Question: "What is the median DVMC for fungi?"

fungi_genes = discover_gene_files("data/", group_name="fungi")
fungi_dvmc = batch_dvmc(fungi_genes)
median_dvmc = np.median(list(fungi_dvmc.values()))

print(f"Median DVMC: {median_dvmc:.4f}")

Workflow 2: Parsimony Sites Ratio

# Question: "What is the ratio of minimum PI sites (fungi / animals)?"

fungi_pi = batch_parsimony_informative(fungi_genes)
animal_pi = batch_parsimony_informative(animal_genes)

fungi_counts = [r[0] for r in fungi_pi.values()]
animal_counts = [r[0] for r in animal_pi.values()]

ratio = min(fungi_counts) / min(animal_counts)

print(f"Ratio: {ratio:.4f}")

Workflow 3: Percentage Above Threshold

# Question: "What percentage of trees have treeness > 0.45?"

treeness_results = batch_treeness(gene_files)
values = list(treeness_results.values())

above = sum(1 for v in values if v > 0.45)
percentage = (above / len(values)) * 100

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

---

Command-Line Tools

Discover Files

python scripts/format_alignment.py discover data/ --group fungi

Filter Alignments

python scripts/format_alignment.py filter data/ --max-gap 5.0 --min-seqs 4

Compute Tree Statistics

python scripts/tree_statistics.py tree.nwk alignment.fa

Remove Gappy Sequences

python scripts/format_alignment.py remove-gappy-seqs input.fa output.fa --max-gap 50.0

Remove Gappy Columns

python scripts/format_alignment.py remove-gappy-cols input.fa output.fa --max-gap 50.0

---

Troubleshooting

"Cannot parse alignment/tree"

• Check file format (FASTA must be >name\nseq, Newick must end with ;)
• See references/troubleshooting.md

No files found

• Check directory structure
• Run: python scripts/format_alignment.py discover data/

Different results from expected

• Check rounding (PhyKIT default: 4 decimals)
• Check PhyKIT version: pip show phykit

---

Next Steps

• Detailed workflows: See SKILL.md
• Alignment analysis: See references/sequence_alignment.md
• Tree construction: See references/tree_building.md
• Statistical comparisons: See references/parsimony_analysis.md
• Troubleshooting: See references/troubleshooting.md

---

Import Statements

# Core imports for most tasks
import numpy as np
import pandas as pd
from scipy import stats

# Load functions from scripts
from scripts.format_alignment import discover_gene_files
from scripts.tree_statistics import (
    phykit_treeness,
    phykit_tree_length,
    phykit_evolutionary_rate,
    phykit_dvmc,
    phykit_rcv,
    phykit_parsimony_informative,
    phykit_treeness_over_rcv,
    batch_treeness,
    batch_dvmc,
    batch_rcv,
    summary_stats,
    compare_groups
)

Parsimony Analysis and Statistical Comparisons

Comprehensive guide for batch phylogenetic analysis and statistical comparisons between groups.

---

File Discovery

Auto-discover Gene Files

import glob
import os

def discover_gene_files(data_dir, group_name=None):
    """Discover paired alignment and tree files for a group.

    Handles common naming patterns:
    - gene1.fa + gene1.nwk
    - gene1.aligned.fa + gene1.treefile
    - gene1_alignment.fa + gene1_tree.nwk

    Args:
        data_dir: root directory to search
        group_name: optional subdirectory (e.g., "fungi", "animals")

    Returns: list of dicts with 'gene_id', 'aln_file', 'tree_file'
    """
    search_dir = os.path.join(data_dir, group_name) if group_name else data_dir

    # Find alignments
    aln_files = {}
    aln_extensions = ['*.fa', '*.fasta', '*.faa', '*.fna', '*.phy', '*.phylip', '*.nex']
    for ext in aln_extensions:
        for f in glob.glob(os.path.join(search_dir, '**', ext), recursive=True):
            gene_id = os.path.splitext(os.path.basename(f))[0]
            # Remove common suffixes
            for suffix in ['.aligned', '.aln', '.msa', '_aligned', '_alignment']:
                gene_id = gene_id.replace(suffix, '')
            aln_files[gene_id] = f

    # Find trees
    tree_files = {}
    tree_extensions = ['*.nwk', '*.newick', '*.tre', '*.tree', '*.treefile']
    for ext in tree_extensions:
        for f in glob.glob(os.path.join(search_dir, '**', ext), recursive=True):
            gene_id = os.path.splitext(os.path.basename(f))[0]
            for suffix in ['.treefile', '_tree', '.rooted', '_rooted']:
                gene_id = gene_id.replace(suffix, '')
            tree_files[gene_id] = f

    # Match pairs
    results = []
    all_gene_ids = set(aln_files.keys()) | set(tree_files.keys())
    for gene_id in sorted(all_gene_ids):
        entry = {'gene_id': gene_id}
        if gene_id in aln_files:
            entry['aln_file'] = aln_files[gene_id]
        if gene_id in tree_files:
            entry['tree_file'] = tree_files[gene_id]
        results.append(entry)

    return results

Usage Examples

# Find all genes in directory
all_genes = discover_gene_files("data/")

# Find genes for specific group
fungi_genes = discover_gene_files("data/", group_name="fungi")
animal_genes = discover_gene_files("data/", group_name="animals")

# Check what was found
print(f"Total genes: {len(all_genes)}")
print(f"With alignments: {sum(1 for g in all_genes if 'aln_file' in g)}")
print(f"With trees: {sum(1 for g in all_genes if 'tree_file' in g)}")
print(f"With both: {sum(1 for g in all_genes if 'aln_file' in g and 'tree_file' in g)}")

---

Batch Metric Computation

Generic Batch Function

def batch_compute_metric(gene_files, metric_func, requires='tree'):
    """Compute a metric across all genes in a group.

    Args:
        gene_files: list from discover_gene_files()
        metric_func: function that takes file path(s) and returns a number
        requires: 'tree', 'alignment', or 'both'

    Returns: dict mapping gene_id -> metric_value (skipping failures)
    """
    results = {}
    for entry in gene_files:
        gene_id = entry['gene_id']
        try:
            if requires == 'tree' and 'tree_file' in entry:
                results[gene_id] = metric_func(entry['tree_file'])
            elif requires == 'alignment' and 'aln_file' in entry:
                results[gene_id] = metric_func(entry['aln_file'])
            elif requires == 'both' and 'tree_file' in entry and 'aln_file' in entry:
                results[gene_id] = metric_func(entry['tree_file'], entry['aln_file'])
        except Exception as e:
            # Skip genes that fail (common with malformed files)
            pass
    return results

Specialized Batch Functions

# Tree metrics
def batch_treeness(gene_files):
    return batch_compute_metric(gene_files, phykit_treeness, requires='tree')

def batch_tree_length(gene_files):
    return batch_compute_metric(gene_files, phykit_tree_length, requires='tree')

def batch_evolutionary_rate(gene_files):
    return batch_compute_metric(gene_files, phykit_evolutionary_rate, requires='tree')

def batch_dvmc(gene_files):
    return batch_compute_metric(gene_files, phykit_dvmc, requires='tree')

# Alignment metrics
def batch_rcv(gene_files):
    return batch_compute_metric(gene_files, phykit_rcv, requires='alignment')

def batch_gap_percentage(gene_files):
    return batch_compute_metric(gene_files, alignment_gap_percentage, requires='alignment')

def batch_parsimony_informative(gene_files):
    """Returns dict of gene_id -> (count, aln_len, percentage)"""
    results = {}
    for entry in gene_files:
        if 'aln_file' in entry:
            try:
                results[entry['gene_id']] = phykit_parsimony_informative(entry['aln_file'])
            except Exception:
                pass
    return results

# Combined metrics
def batch_treeness_over_rcv(gene_files):
    """Returns dict of gene_id -> (treeness_over_rcv, treeness, rcv)"""
    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

---

Summary Statistics

Basic Statistics

import numpy as np

def summary_stats(values):
    """Compute standard summary statistics for a list of values.

    Returns: dict with count, mean, median, std, var, min, max, quartiles
    """
    arr = np.array(values)
    return {
        'count': len(arr),
        'mean': float(np.mean(arr)),
        'median': float(np.median(arr)),
        'std': float(np.std(arr, ddof=1)),
        'var': float(np.var(arr, ddof=1)),
        'min': float(np.min(arr)),
        'max': float(np.max(arr)),
        'q25': float(np.percentile(arr, 25)),
        'q75': float(np.percentile(arr, 75)),
    }

Usage

# Compute metric
dvmc_results = batch_dvmc(gene_files)
dvmc_values = list(dvmc_results.values())

# Get statistics
stats = summary_stats(dvmc_values)

print(f"Count: {stats['count']}")
print(f"Mean: {stats['mean']:.4f}")
print(f"Median: {stats['median']:.4f}")
print(f"Std Dev: {stats['std']:.4f}")
print(f"Min: {stats['min']:.4f}")
print(f"Max: {stats['max']:.4f}")
print(f"Q1: {stats['q25']:.4f}")
print(f"Q3: {stats['q75']:.4f}")

---

Group Comparisons

Mann-Whitney U Test

from scipy import stats

def compare_groups(group1_values, group2_values, group1_name="Group1", group2_name="Group2"):
    """Compare two groups using Mann-Whitney U test.

    Mann-Whitney U is a non-parametric test for comparing distributions.
    Does not assume normality.

    Args:
        group1_values: list of values for group 1
        group2_values: list of values for group 2
        group1_name: label for group 1
        group2_name: label for group 2

    Returns: dict with U statistic, p-value, and summary stats for each group
    """
    arr1 = np.array(group1_values)
    arr2 = np.array(group2_values)

    # Mann-Whitney U test (two-sided, default)
    u_stat, p_value = stats.mannwhitneyu(arr1, arr2, alternative='two-sided')

    return {
        'u_statistic': float(u_stat),
        'p_value': float(p_value),
        group1_name: summary_stats(arr1.tolist()),
        group2_name: summary_stats(arr2.tolist()),
        'median_difference': float(np.median(arr1) - np.median(arr2)),
    }

Usage

# Compute metrics for both groups
fungi_treeness = batch_treeness(fungi_genes)
animal_treeness = batch_treeness(animal_genes)

# Compare
comparison = compare_groups(
    list(fungi_treeness.values()),
    list(animal_treeness.values()),
    group1_name="Fungi",
    group2_name="Animals"
)

print(f"Mann-Whitney U: {comparison['u_statistic']:.0f}")
print(f"P-value: {comparison['p_value']:.4e}")
print(f"Fungi median: {comparison['Fungi']['median']:.4f}")
print(f"Animals median: {comparison['Animals']['median']:.4f}")
print(f"Difference: {comparison['median_difference']:.4f}")

---

Paired Comparisons

Paired Gene Analysis

def paired_comparison(group1_dict, group2_dict):
    """Compare matched pairs (same gene IDs in both groups).

    Useful when comparing orthologs between species.

    Args:
        group1_dict: dict of gene_id -> value for group 1
        group2_dict: dict of gene_id -> value for group 2

    Returns: dict with paired differences and statistics
    """
    common_genes = set(group1_dict.keys()) & set(group2_dict.keys())

    diffs = []
    ratios = []
    for gene in sorted(common_genes):
        v1 = group1_dict[gene]
        v2 = group2_dict[gene]
        diffs.append(v1 - v2)
        if v2 != 0:
            ratios.append(v1 / v2)

    result = {
        'n_pairs': len(common_genes),
        'differences': summary_stats(diffs),
        'median_difference': float(np.median(diffs)),
    }

    if ratios:
        result['ratios'] = summary_stats(ratios)
        result['median_ratio'] = float(np.median(ratios))

    return result

Usage

# Compute tree length for both groups
fungi_lengths = batch_tree_length(fungi_genes)
animal_lengths = batch_tree_length(animal_genes)

# Paired comparison
paired = paired_comparison(fungi_lengths, animal_lengths)

print(f"Paired genes: {paired['n_pairs']}")
print(f"Median difference: {paired['median_difference']:.4f}")
print(f"Median ratio: {paired['median_ratio']:.4f}")

---

Complete Workflow Examples

Example 1: DVMC Comparison

# Question: "What is the median DVMC for fungi, and the Mann-Whitney U statistic comparing fungi vs animals?"

# 1. Discover files
fungi_genes = discover_gene_files("data/", group_name="fungi")
animal_genes = discover_gene_files("data/", group_name="animals")

print(f"Fungi genes: {len(fungi_genes)}")
print(f"Animal genes: {len(animal_genes)}")

# 2. Compute DVMC for both groups
fungi_dvmc = batch_dvmc(fungi_genes)
animal_dvmc = batch_dvmc(animal_genes)

print(f"Computed DVMC for {len(fungi_dvmc)} fungi genes")
print(f"Computed DVMC for {len(animal_dvmc)} animal genes")

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

# 4. Compute median for fungi
fungi_median = np.median(fungi_values)
print(f"\nFungi median DVMC: {fungi_median:.4f}")

# 5. Mann-Whitney U test
u_stat, p_value = stats.mannwhitneyu(fungi_values, animal_values, alternative='two-sided')
print(f"Mann-Whitney U statistic: {u_stat:.0f}")
print(f"P-value: {p_value:.4e}")

# Answer: fungi_median (rounded to 4 decimals), u_stat (integer)

Example 2: Treeness with Filtering

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

# 1. Discover files
gene_files = discover_gene_files("data/")

# 2. Filter by gap percentage
valid_genes = []
for entry in gene_files:
    if 'aln_file' in entry and 'tree_file' in entry:
        try:
            gap_pct = alignment_gap_percentage(entry['aln_file'])
            if gap_pct < 5.0:
                valid_genes.append(entry)
        except Exception:
            pass

print(f"Genes with <5% gaps: {len(valid_genes)}/{len(gene_files)}")

# 3. Compute treeness/RCV
results = batch_treeness_over_rcv(valid_genes)

# 4. Extract ratios
ratios = [result[0] for result in results.values()]  # First element is ratio

# 5. Compute median
median_ratio = np.median(ratios)
print(f"Median treeness/RCV: {median_ratio:.4f}")

# Answer: median_ratio (rounded to 4 decimals)

Example 3: Parsimony Sites Comparison

# Question: "What is the ratio of minimum PI sites (fungi / animals)?"

# 1. Compute PI sites
fungi_pi = batch_parsimony_informative(fungi_genes)
animal_pi = batch_parsimony_informative(animal_genes)

# 2. Extract counts (first element of tuple)
fungi_counts = [result[0] for result in fungi_pi.values()]
animal_counts = [result[0] for result in animal_pi.values()]

# 3. Find minimums
fungi_min = min(fungi_counts)
animal_min = min(animal_counts)

# 4. Compute ratio
ratio = fungi_min / animal_min

print(f"Fungi min PI sites: {fungi_min}")
print(f"Animals min PI sites: {animal_min}")
print(f"Ratio: {ratio:.4f}")

# Answer: ratio (rounded to 4 decimals)

Example 4: Percentage Above Threshold

# Question: "What percentage of animal genes have treeness > 0.45?"

# 1. Compute treeness
animal_treeness = batch_treeness(animal_genes)
treeness_values = list(animal_treeness.values())

# 2. Count above threshold
above_threshold = sum(1 for v in treeness_values if v > 0.45)

# 3. Compute percentage
percentage = (above_threshold / len(treeness_values)) * 100

print(f"Total genes: {len(treeness_values)}")
print(f"Above 0.45: {above_threshold}")
print(f"Percentage: {percentage:.2f}%")

# Answer: percentage (rounded to 2 decimals if question asks for "XX.XX%")

Example 5: Specific Gene Lookup

# Question: "What is the evolutionary rate for gene ENSG00000141510?"

# 1. Discover files
gene_files = discover_gene_files("data/")

# 2. Find specific gene
target_gene = [g for g in gene_files if g['gene_id'] == 'ENSG00000141510'][0]

# 3. Compute evolutionary rate
evo_rate = phykit_evolutionary_rate(target_gene['tree_file'])

print(f"Evolutionary rate for ENSG00000141510: {evo_rate:.4f}")

# Answer: evo_rate (rounded to 4 decimals)

---

Answer Extraction Patterns

Common Patterns

# Median
answer = np.median(values)

# Mean
answer = np.mean(values)

# Maximum
answer = np.max(values)

# Minimum
answer = np.min(values)

# Standard deviation
answer = np.std(values, ddof=1)

# Variance
answer = np.var(values, ddof=1)

# Percentile (e.g., 75th)
answer = np.percentile(values, 75)

# Count above threshold
answer = sum(1 for v in values if v > threshold)

# Percentage above threshold
answer = (sum(1 for v in values if v > threshold) / len(values)) * 100

# Difference in medians
answer = np.median(group1) - np.median(group2)

# Absolute difference
answer = abs(np.median(group1) - np.median(group2))

# Fold change (ratio)
answer = np.median(group1) / np.median(group2)

# Mann-Whitney U statistic
answer = stats.mannwhitneyu(group1, group2)[0]

# Mann-Whitney p-value
answer = stats.mannwhitneyu(group1, group2)[1]

Rounding Guide

# PhyKIT default: 4 decimals
answer = round(value, 4)

# Percentages: match question format
# "35%" -> integer
answer = round(percentage)

# "3.5%" -> 1 decimal
answer = round(percentage, 1)

# P-values: scientific notation for small values
if p_value < 0.0001:
    answer_str = f"{p_value:.2e}"
else:
    answer_str = f"{p_value:.4f}"

# U statistics: integer
answer = int(round(u_stat))

# Special: "multiplied by 1000"
answer = round(value * 1000)

---

Troubleshooting

No Files Found

Issue: discover_gene_files() returns empty list.

Solution: Check directory structure and file extensions:

import os
print(os.listdir("data/"))
print(os.listdir("data/fungi/"))

Mismatched Gene Counts

Issue: Different numbers of genes between groups.

Solution: This is normal. Use paired comparison only for matched genes:

paired = paired_comparison(fungi_results, animal_results)
print(f"Matched genes: {paired['n_pairs']}")

Failed Metric Computation

Issue: Some genes fail during batch processing.

Solution: This is expected (malformed files, missing data). Failures are silently skipped:

# Check success rate
total_genes = len(gene_files)
successful = len(results)
print(f"Success rate: {successful}/{total_genes} ({successful/total_genes*100:.1f}%)")

---

See Also

• sequence_alignment.md - Alignment-specific metrics
• tree_building.md - Tree-specific metrics
• troubleshooting.md - Common errors and solutions
• scripts/tree_statistics.py - Implementation code

Sequence Alignment Analysis

Detailed guide for analyzing multiple sequence alignments using Biopython and PhyKIT.

---

File Loading and Format Detection

Auto-detect Alignment Format

from Bio import AlignIO

def load_alignment(filepath):
    """Load alignment with auto-format detection.

    Supports: FASTA, PHYLIP, PHYLIP-relaxed, Nexus, Clustal, Stockholm
    """
    # Try phylip-relaxed BEFORE phylip to avoid misparse of long names
    formats_to_try = ['fasta', 'phylip-relaxed', 'phylip', 'nexus', 'clustal', 'stockholm']

    for fmt in formats_to_try:
        try:
            alignment = AlignIO.read(filepath, fmt)
            return alignment, fmt
        except Exception:
            continue

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

Usage

# Load any alignment format
alignment, format_detected = load_alignment("gene1.fa")
print(f"Format: {format_detected}")
print(f"Sequences: {len(alignment)}")
print(f"Length: {alignment.get_alignment_length()}")

---

Parsimony Informative Sites

PhyKIT Implementation

from phykit.services.alignment.parsimony_informative_sites import ParsimonyInformative
from types import SimpleNamespace

def phykit_parsimony_informative(aln_file):
    """Calculate parsimony informative sites.

    Parsimony informative site: A site with at least 2 different characters
    appearing at least 2 times each (excluding gaps).

    Returns: (pi_sites_count, alignment_length, pi_percentage)
    """
    pi = ParsimonyInformative(SimpleNamespace(alignment=aln_file))
    alignment, _, _ = pi.get_alignment_and_format()
    return pi.calculate_parsimony_informative_sites(alignment)

Example

# Compute parsimony informative sites
pi_count, aln_len, pi_pct = phykit_parsimony_informative("alignment.fa")

print(f"Parsimony informative sites: {pi_count}")
print(f"Alignment length: {aln_len}")
print(f"Percentage: {pi_pct:.4f}%")

Batch Processing

def batch_parsimony_informative(gene_files):
    """Compute PI sites for multiple genes."""
    results = {}
    for entry in gene_files:
        if 'aln_file' in entry:
            try:
                results[entry['gene_id']] = phykit_parsimony_informative(entry['aln_file'])
            except Exception:
                pass  # Skip genes that fail
    return results

# Usage
pi_results = batch_parsimony_informative(gene_files)

# Extract just counts
pi_counts = {gene: result[0] for gene, result in pi_results.items()}

# Extract percentages
pi_percentages = {gene: result[2] for gene, result in pi_results.items()}

---

Relative Composition Variability (RCV)

PhyKIT Implementation

from phykit.services.alignment.rcv import RelativeCompositionVariability

def phykit_rcv(aln_file):
    """Calculate Relative Composition Variability.

    RCV measures compositional heterogeneity across sequences.
    Lower values indicate more uniform base composition.

    Returns: float (RCV score)
    """
    rcv = RelativeCompositionVariability(SimpleNamespace(alignment=aln_file))
    return rcv.calculate_rcv()

Example

# Single alignment
rcv_score = phykit_rcv("alignment.fa")
print(f"RCV: {rcv_score:.4f}")

# Batch processing
def batch_rcv(gene_files):
    return batch_compute_metric(gene_files, phykit_rcv, requires='alignment')

rcv_results = batch_rcv(gene_files)

---

Alignment Gap Analysis

Gap Percentage

import numpy as np

def alignment_gap_percentage(aln_file):
    """Calculate overall gap percentage in alignment.

    Gap characters: '-', '.', '?'

    Returns: percentage (0-100)
    """
    alignment, fmt = load_alignment(aln_file)
    n_seqs = len(alignment)
    aln_len = alignment.get_alignment_length()
    total_chars = n_seqs * aln_len

    # Convert to numpy array
    arr = np.array([[c for c in str(rec.seq)] for rec in alignment])

    # Count gaps
    gap_count = np.sum(np.isin(arr, ['-', '.', '?']))

    return (gap_count / total_chars) * 100

Per-Sequence Gap Statistics

def per_sequence_gaps(aln_file):
    """Calculate gap percentage for each sequence."""
    alignment, _ = load_alignment(aln_file)

    results = []
    for record in alignment:
        seq_str = str(record.seq)
        gap_count = seq_str.count('-') + seq_str.count('.') + seq_str.count('?')
        gap_pct = (gap_count / len(seq_str)) * 100
        results.append({
            'seq_id': record.id,
            'gap_count': gap_count,
            'gap_percentage': gap_pct
        })

    return results

---

Comprehensive Alignment Statistics

All-in-One Function

def alignment_statistics(aln_file):
    """Comprehensive alignment statistics.

    Returns dict with:
    - n_seqs: Number of sequences
    - aln_len: Alignment length
    - gap_pct: Gap percentage
    - gc_pct: GC content
    - pi_sites: Parsimony informative sites count
    - pi_pct: Parsimony informative sites percentage
    - variable_sites: Number of variable sites
    - variable_pct: Variable sites percentage
    """
    alignment, fmt = load_alignment(aln_file)
    n_seqs = len(alignment)
    aln_len = alignment.get_alignment_length()

    # Convert to numpy array
    arr = np.array([[c.upper() for c in str(rec.seq)] for rec in alignment])

    # Gap percentage
    gap_count = np.sum(np.isin(arr, ['-', '.', '?']))
    gap_pct = (gap_count / (n_seqs * aln_len)) * 100

    # GC content (excluding gaps)
    non_gap = arr[~np.isin(arr, ['-', '.', '?', 'N'])]
    gc_count = np.sum(np.isin(non_gap, ['G', 'C']))
    gc_pct = (gc_count / len(non_gap)) * 100 if len(non_gap) > 0 else 0

    # Parsimony informative sites
    pi_sites, _, pi_pct = phykit_parsimony_informative(aln_file)

    # Variable sites
    variable_count = 0
    for i in range(aln_len):
        col = arr[:, i]
        non_gap_col = col[~np.isin(col, ['-', '.', '?'])]
        if len(np.unique(non_gap_col)) > 1:
            variable_count += 1

    return {
        'n_seqs': n_seqs,
        'aln_len': aln_len,
        'gap_pct': round(gap_pct, 4),
        'gc_pct': round(gc_pct, 4),
        'pi_sites': pi_sites,
        'pi_pct': round(pi_pct, 4),
        'variable_sites': variable_count,
        'variable_pct': round((variable_count / aln_len) * 100, 4)
    }

Example

# Get comprehensive stats
stats = alignment_statistics("alignment.fa")

print(f"Sequences: {stats['n_seqs']}")
print(f"Length: {stats['aln_len']}")
print(f"Gaps: {stats['gap_pct']:.2f}%")
print(f"GC: {stats['gc_pct']:.2f}%")
print(f"PI sites: {stats['pi_sites']} ({stats['pi_pct']:.2f}%)")
print(f"Variable sites: {stats['variable_sites']} ({stats['variable_pct']:.2f}%)")

---

Format Conversion

Convert Between Formats

def convert_alignment_format(input_file, output_file, output_format='fasta'):
    """Convert alignment to different format.

    Supported formats: fasta, phylip, phylip-relaxed, nexus, clustal
    """
    alignment, input_format = load_alignment(input_file)
    AlignIO.write(alignment, output_file, output_format)
    return output_format

Example

# Convert PHYLIP to FASTA
convert_alignment_format("alignment.phy", "alignment.fa", "fasta")

# Convert FASTA to Nexus
convert_alignment_format("alignment.fa", "alignment.nex", "nexus")

---

Filtering Alignments

By Gap Threshold

def filter_by_gap_threshold(gene_files, max_gap_pct=5.0):
    """Filter alignments by maximum gap percentage.

    Returns: list of gene entries passing threshold
    """
    valid_genes = []
    for entry in gene_files:
        if 'aln_file' in entry:
            try:
                gap_pct = alignment_gap_percentage(entry['aln_file'])
                if gap_pct <= max_gap_pct:
                    valid_genes.append(entry)
            except Exception:
                pass
    return valid_genes

By Minimum Sequences

def filter_by_min_sequences(gene_files, min_seqs=4):
    """Filter alignments by minimum number of sequences."""
    valid_genes = []
    for entry in gene_files:
        if 'aln_file' in entry:
            try:
                alignment, _ = load_alignment(entry['aln_file'])
                if len(alignment) >= min_seqs:
                    valid_genes.append(entry)
            except Exception:
                pass
    return valid_genes

Combined Filtering

# Filter by multiple criteria
valid_genes = gene_files
valid_genes = filter_by_gap_threshold(valid_genes, max_gap_pct=5.0)
valid_genes = filter_by_min_sequences(valid_genes, min_seqs=4)

print(f"Genes passing filters: {len(valid_genes)}/{len(gene_files)}")

---

Common Patterns

Pattern: Median Parsimony Sites

# Question: "What is the median number of parsimony informative sites?"

# 1. Batch compute
pi_results = batch_parsimony_informative(gene_files)

# 2. Extract counts
pi_counts = [result[0] for result in pi_results.values()]

# 3. Compute median
median_pi = np.median(pi_counts)

print(f"Median PI sites: {median_pi:.4f}")

Pattern: Percentage Above Threshold

# Question: "What percentage of alignments have >100 PI sites?"

pi_counts = [result[0] for result in pi_results.values()]
above_threshold = sum(1 for count in pi_counts if count > 100)
percentage = (above_threshold / len(pi_counts)) * 100

print(f"Percentage > 100 PI sites: {percentage:.2f}%")

Pattern: Ratio of Minimums

# Question: "What is the ratio of minimum PI sites (fungi / animals)?"

fungi_pi = batch_parsimony_informative(fungi_genes)
animal_pi = batch_parsimony_informative(animal_genes)

fungi_counts = [r[0] for r in fungi_pi.values()]
animal_counts = [r[0] for r in animal_pi.values()]

ratio = min(fungi_counts) / min(animal_counts)

print(f"Ratio of minimums: {ratio:.4f}")

---

Troubleshooting

"Cannot parse alignment" Error

Cause: Format detection failed.

Solution: Try specifying format explicitly:

alignment = AlignIO.read("file.phy", "phylip-relaxed")

Long Sequence Names in PHYLIP

Issue: PHYLIP format truncates names to 10 characters.

Solution: Use phylip-relaxed format:

alignment = AlignIO.read("file.phy", "phylip-relaxed")

Empty PI Sites

Cause: Alignment lacks informative sites (too few sequences or too conserved).

Solution: Check alignment quality:

stats = alignment_statistics("alignment.fa")
print(f"Variable sites: {stats['variable_sites']}")

---

See Also

• tree_building.md - Tree construction from alignments
• parsimony_analysis.md - Statistical comparison workflows
• scripts/format_alignment.py - Format conversion utilities

Tree Building and Analysis

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

---

Tree File Loading

Auto-detect Tree Format

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

# 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

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

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

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)

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

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)

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

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

# 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

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

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

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

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

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

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

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

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

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

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

tree = build_upgma_tree("alignment.fa")

For Small Datasets

Use Maximum Parsimony:

• Good for discrete characters
• Interpretable (minimum substitutions)
• Slow for large datasets

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:

# After running IQ-TREE externally
treeness = phykit_treeness("iqtree_output.treefile")
dvmc = phykit_dvmc("iqtree_output.treefile")

---

Common Patterns

Pattern: Median Treeness

# 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

# 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

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

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

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:

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:

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

Troubleshooting Guide

Common issues and solutions for phylogenetics skill.

---

File Loading Issues

"Cannot parse alignment" Error

Symptoms: ValueError: Cannot parse alignment: file.fa

Causes: 1. Invalid file format 2. Corrupted file 3. Wrong file extension

Solutions:

# Try specifying format explicitly
from Bio import AlignIO
alignment = AlignIO.read("file.phy", "phylip-relaxed")

# Check file contents
with open("file.fa", 'r') as f:
    print(f.read()[:200])  # First 200 characters

# Try different formats
for fmt in ['fasta', 'phylip', 'phylip-relaxed', 'nexus', 'clustal']:
    try:
        aln = AlignIO.read("file.fa", fmt)
        print(f"Success with format: {fmt}")
        break
    except Exception as e:
        print(f"Failed with {fmt}: {e}")

"Cannot parse tree" Error

Symptoms: ValueError: Cannot parse tree: file.nwk

Causes: 1. Missing semicolon (Newick must end with ;) 2. Invalid Newick syntax 3. File encoding issues

Solutions:

# Check tree syntax
with open("file.nwk", 'r') as f:
    content = f.read()
    print(content)
    if not content.strip().endswith(';'):
        print("ERROR: Newick tree must end with semicolon")

# Fix missing semicolon
with open("file.nwk", 'r') as f:
    content = f.read().strip()
if not content.endswith(';'):
    with open("file_fixed.nwk", 'w') as f:
        f.write(content + ';\n')

# Try loading
from Bio import Phylo
tree = Phylo.read("file_fixed.nwk", "newick")

---

PhyKIT Errors

"PhyKIT: AttributeError: 'NoneType' object has no attribute..."

Symptoms: PhyKIT functions crash with NoneType errors

Causes: 1. Tree has no branch lengths 2. Invalid tree structure 3. Empty alignment

Solutions:

# Check if tree has branch lengths
from Bio import Phylo
tree = Phylo.read("tree.nwk", "newick")

has_branch_lengths = all(
    clade.branch_length is not None
    for clade in tree.find_clades()
)

if not has_branch_lengths:
    print("ERROR: Tree missing branch lengths")
    print("Solution: Rebuild tree with distance method (NJ, UPGMA)")

# Rebuild with NJ
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
alignment = AlignIO.read("alignment.fa", "fasta")
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(alignment)
constructor = DistanceTreeConstructor()
tree_with_lengths = constructor.nj(dm)
Phylo.write(tree_with_lengths, "tree_fixed.nwk", "newick")

PhyKIT Returns Unexpected Values

Symptoms: Treeness > 1.0 or negative values

Causes: 1. Invalid input data 2. PhyKIT version mismatch

Solutions:

# Check PhyKIT version
import phykit
print(f"PhyKIT version: {phykit.__version__}")

# Validate tree
tree = Phylo.read("tree.nwk", "newick")
total_length = sum(c.branch_length for c in tree.find_clades() if c.branch_length)
print(f"Total tree length: {total_length}")

if total_length <= 0:
    print("ERROR: Tree has zero or negative branch lengths")

---

Data Discovery Issues

No Files Found by discover_gene_files()

Symptoms: Function returns empty list

Causes: 1. Wrong directory path 2. Unexpected file extensions 3. Subdirectory structure not as expected

Solutions:

import os
import glob

# Check directory exists
data_dir = "data/"
if not os.path.exists(data_dir):
    print(f"ERROR: Directory {data_dir} not found")

# List files manually
print("Files in directory:")
for root, dirs, files in os.walk(data_dir):
    for f in files:
        print(os.path.join(root, f))

# Check for specific extensions
aln_files = glob.glob(os.path.join(data_dir, '**', '*.fa'), recursive=True)
tree_files = glob.glob(os.path.join(data_dir, '**', '*.nwk'), recursive=True)
print(f"Found {len(aln_files)} .fa files")
print(f"Found {len(tree_files)} .nwk files")

# Adjust extensions if needed
if len(aln_files) == 0:
    # Try other extensions
    for ext in ['*.fasta', '*.faa', '*.phy']:
        files = glob.glob(os.path.join(data_dir, '**', ext), recursive=True)
        if files:
            print(f"Found {len(files)} {ext} files")

Mismatched Gene IDs

Symptoms: Alignment and tree files don't match

Causes: 1. Different naming conventions 2. Suffixes not handled by discover_gene_files()

Solutions:

# Debug gene ID extraction
import os

def debug_gene_id_extraction(filepath):
    basename = os.path.basename(filepath)
    gene_id = os.path.splitext(basename)[0]
    print(f"File: {basename} -> Gene ID: {gene_id}")

# Test on your files
debug_gene_id_extraction("gene1.aligned.fa")
debug_gene_id_extraction("gene1.nwk")

# If IDs don't match, add custom suffix handling
def discover_gene_files_custom(data_dir):
    # Add your custom suffixes
    custom_suffixes = ['.aligned', '.aln', '.msa', '_aligned',
                      '.treefile', '_tree', '.YOUR_SUFFIX']
    # ... rest of function

---

Batch Processing Issues

Some Genes Fail Silently

Symptoms: Fewer results than input genes

Causes: 1. Malformed files (expected, handled by try/except) 2. Missing required files 3. Computation errors

Solutions:

# Track failures explicitly
def batch_compute_metric_verbose(gene_files, metric_func, requires='tree'):
    results = {}
    failures = []

    for entry in gene_files:
        gene_id = entry['gene_id']
        try:
            if requires == 'tree' and 'tree_file' in entry:
                results[gene_id] = metric_func(entry['tree_file'])
            elif requires == 'alignment' and 'aln_file' in entry:
                results[gene_id] = metric_func(entry['aln_file'])
        except Exception as e:
            failures.append({
                'gene_id': gene_id,
                'error': str(e),
                'file': entry.get('tree_file') or entry.get('aln_file')
            })

    print(f"Success: {len(results)}/{len(gene_files)}")
    print(f"Failures: {len(failures)}")

    if failures:
        print("\nFirst 5 failures:")
        for fail in failures[:5]:
            print(f"  {fail['gene_id']}: {fail['error'][:50]}")

    return results, failures

# Use verbose version for debugging
results, failures = batch_compute_metric_verbose(
    gene_files, phykit_treeness, requires='tree'
)

Out of Memory Errors

Symptoms: Process crashes with large datasets

Causes: 1. Too many genes processed at once 2. Large alignment files

Solutions:

# Process in batches
def batch_compute_in_chunks(gene_files, metric_func, requires='tree', chunk_size=100):
    all_results = {}

    for i in range(0, len(gene_files), chunk_size):
        chunk = gene_files[i:i+chunk_size]
        print(f"Processing chunk {i//chunk_size + 1} ({len(chunk)} genes)...")

        chunk_results = batch_compute_metric(chunk, metric_func, requires)
        all_results.update(chunk_results)

    return all_results

# Use chunked processing
results = batch_compute_in_chunks(gene_files, phykit_treeness, chunk_size=50)

---

Statistical Analysis Issues

Mann-Whitney U Test Fails

Symptoms: ValueError: All numbers are identical or similar

Causes: 1. No variation in data (all values identical) 2. Too few samples

Solutions:

# Check for variation
import numpy as np

def safe_mannwhitneyu(group1, group2):
    """Mann-Whitney U with error handling."""
    arr1 = np.array(group1)
    arr2 = np.array(group2)

    # Check for variation
    if len(np.unique(arr1)) == 1 and len(np.unique(arr2)) == 1:
        if arr1[0] == arr2[0]:
            print("WARNING: No variation in either group")
            return None, 1.0  # p=1.0 means no difference

    # Check sample size
    if len(arr1) < 3 or len(arr2) < 3:
        print("WARNING: Sample size too small for Mann-Whitney U")
        return None, None

    try:
        return stats.mannwhitneyu(arr1, arr2, alternative='two-sided')
    except Exception as e:
        print(f"ERROR: {e}")
        return None, None

# Usage
u_stat, p_value = safe_mannwhitneyu(group1_values, group2_values)
if u_stat is not None:
    print(f"U: {u_stat:.0f}, p: {p_value:.4e}")

---

Rounding and Precision Issues

Answer Doesn't Match Expected

Symptoms: Your answer is close but not exact

Causes: 1. Wrong rounding precision 2. Different PhyKIT version 3. Float precision differences

Solutions:

# Check PhyKIT default precision
value = phykit_treeness("tree.nwk")
print(f"Raw value: {value}")
print(f"4 decimals: {round(value, 4)}")
print(f"3 decimals: {round(value, 3)}")
print(f"2 decimals: {round(value, 2)}")

# Match question rounding
# "rounded to 3 decimal places" -> round(value, 3)
# "rounded to nearest integer" -> round(value)
# "multiplied by 1000" -> round(value * 1000)

# For percentages, check format
percentage = 34.5678
print(f"Integer: {round(percentage)}")  # 35
print(f"1 decimal: {round(percentage, 1)}")  # 34.6
print(f"2 decimals: {round(percentage, 2)}")  # 34.57

---

Tree Construction Issues

NJ/UPGMA Tree Construction Fails

Symptoms: Error during tree building

Causes: 1. Alignment too short 2. Sequences too similar (zero distances) 3. Invalid alignment

Solutions:

# Check alignment quality
from Bio import AlignIO

alignment = AlignIO.read("alignment.fa", "fasta")
print(f"Sequences: {len(alignment)}")
print(f"Length: {alignment.get_alignment_length()}")

if len(alignment) < 3:
    print("ERROR: Need at least 3 sequences for tree")

if alignment.get_alignment_length() < 10:
    print("WARNING: Alignment very short")

# Check distance matrix
from Bio.Phylo.TreeConstruction import DistanceCalculator

calculator = DistanceCalculator('identity')
dm = calculator.get_distance(alignment)
print("Distance matrix:")
print(dm)

# If all distances are 0, sequences are identical

Parsimony Tree Takes Too Long

Symptoms: build_parsimony_tree() doesn't finish

Causes: 1. Too many sequences (parsimony is slow) 2. Long alignment

Solutions:

# Check size
alignment = AlignIO.read("alignment.fa", "fasta")
n_seqs = len(alignment)
aln_len = alignment.get_alignment_length()

print(f"Size: {n_seqs} sequences x {aln_len} sites")

if n_seqs > 50:
    print("WARNING: Parsimony will be very slow for >50 sequences")
    print("Recommendation: Use NJ instead, or use external tools (IQ-TREE, RAxML)")

# Use NJ as faster alternative
tree = build_nj_tree("alignment.fa")

---

Performance Issues

Slow Batch Processing

Symptoms: Processing takes very long

Causes: 1. Many genes 2. Large files 3. Inefficient I/O

Solutions:

# Add progress tracking
from tqdm import tqdm

def batch_compute_with_progress(gene_files, metric_func, requires='tree'):
    results = {}
    for entry in tqdm(gene_files, desc="Processing genes"):
        gene_id = entry['gene_id']
        try:
            if requires == 'tree' and 'tree_file' in entry:
                results[gene_id] = metric_func(entry['tree_file'])
            elif requires == 'alignment' and 'aln_file' in entry:
                results[gene_id] = metric_func(entry['aln_file'])
        except Exception:
            pass
    return results

# Use with progress bar
results = batch_compute_with_progress(gene_files, phykit_treeness)

---

Common Issues

Answer Format Mismatch

Symptoms: Answer doesn't match expected format

Causes: 1. Wrong type (float vs int) 2. Wrong precision 3. Scientific notation vs decimal

Solutions:

# Common answer formats

# Integer (U statistic, counts)
answer = int(round(value))

# Float with 4 decimals (PhyKIT default)
answer = round(value, 4)

# Percentage as integer
answer = int(round(percentage))

# Percentage with decimals
answer = round(percentage, 2)

# Scientific notation for small p-values
if p_value < 0.0001:
    answer = f"{p_value:.2e}"  # e.g., "1.23e-05"
else:
    answer = round(p_value, 4)

# Multiplied values
answer = int(round(value * 1000))

---

Getting Help

Enable Debug Output

# For PhyKIT issues
import logging
logging.basicConfig(level=logging.DEBUG)

# For Biopython issues
from Bio import BiopythonWarning
import warnings
warnings.simplefilter('always', BiopythonWarning)

Check Versions

import phykit
import Bio
import dendropy
import numpy
import scipy

print(f"PhyKIT: {phykit.__version__}")
print(f"Biopython: {Bio.__version__}")
print(f"DendroPy: {dendropy.__version__}")
print(f"NumPy: {numpy.__version__}")
print(f"SciPy: {scipy.__version__}")

Test with Simple Example

# Create minimal test files
test_alignment = ">seq1\nACGT\n>seq2\nACGT\n>seq3\nTCGA\n"
with open("test.fa", 'w') as f:
    f.write(test_alignment)

# Test alignment loading
alignment = AlignIO.read("test.fa", "fasta")
print(f"Loaded {len(alignment)} sequences")

# Test tree building
tree = build_nj_tree("test.fa")
Phylo.write(tree, "test.nwk", "newick")

# Test PhyKIT
treeness = phykit_treeness("test.nwk")
print(f"Treeness: {treeness:.4f}")

---

See Also

• sequence_alignment.md - Alignment analysis details
• tree_building.md - Tree construction methods
• parsimony_analysis.md - Statistical workflows
• PhyKIT documentation: https://jlsteenwyk.com/PhyKIT/
• Biopython Phylo tutorial: https://biopython.org/wiki/Phylo