Phylogenetics and Sequence Analysis

Four traps that produce a confidently wrong number

Each of these was observed producing a wrong answer while the correct guidance was already present further down this file. Check them before you answer.

1. PhyKIT prints more than one column, and for saturation the two conventions disagree — state which you used. phykit saturation prints saturation <TAB> |saturation-1|. Its own --help is explicit: "The first value is the saturation value and the second column is the absolute value of saturation minus 1." But several published analyses (and some reference answers derived from them) report the second column as "the saturation value". The two always sum to 1.0000, which is the tell that you may be looking at the wrong one — on the fungal scogs set the medians are 0.39 (col 1) and 0.61 (col 2).

   So: follow phykit and use column 1 unless the question or source defines saturation the other way, and say in your answer which column you read. Do not silently pick the one that looks closer to an expected number.

   treeness_over_rcv has no such ambiguity: it gives ratio <TAB> treeness <TAB> RCV and the ratio is first.

2. "Gap percentage" means the fraction of alignment COLUMNS containing at least one gap, not the fraction of residues that are gaps. On the fungal scogs set the residue definition maxes out at 0.556, so a ">70% gaps" filter selects nothing and the question looks unanswerable; by columns, three orthologs qualify (max 0.783).

3. treeness_over_rcv and rcv take the UNTRIMMED .faa.mafft, while saturation takes the trimmed .clipkit. RCV measures variability across columns, so trimming changes it: median 0.2683 untrimmed against 0.3050 trimmed, and among >70%-gap genes the maximum is 0.2572 untrimmed against 0.4174 trimmed.

4. Never loop PhyKIT per file. phykit_batch_analysis is parallel and does ~250 trees in about 35 seconds; a shell loop takes ~9 minutes and runs out of turns mid-way, producing no answer at all. It also selects the right column for every function, which removes trap 1 entirely.

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.

Question phrasing	Count
"total amino acids in all single-copy ortholog sequences"	Sum of ONE sequence per ortholog (either the FIRST entry per file or the median-length entry)
"total amino acids across N species' single-copy orthologs"	Sum across species explicitly (multi-species sum)
"average length of single-copy orthologs"	Mean per-ortholog length (one per ortholog)

❌ 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.
• 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.

Which alignment goes with which metric (this changes the answer)

The tree is always the ClipKit-derived .faa.mafft.clipkit.treefile. The alignment argument depends on the metric:

metric	alignment to pass
treeness, dvmc, total_tree_length, long_branch_score	tree only — no alignment
saturation	.faa.mafft.clipkit (trimmed)
treeness_over_rcv / rcv	.faa.mafft (untrimmed)
parsimony-informative sites, gap percentage	.faa.mafft (untrimmed)

RCV measures compositional variability across the alignment's columns, so trimming changes it materially — and treeness_over_rcv divides by RCV, so the trimmed alignment shifts the ratio for every gene. Verified on the fungal scogs set (249 orthologs, canonical shipped files):

median treeness/RCV   untrimmed .faa.mafft = 0.2683    trimmed .clipkit = 0.3050
max treeness/RCV      (over the 3 genes with >70% gapped columns:
                       1260807at2759 0.0861, 1567796at2759 0.1866, 939345at2759 0.2572)
                      untrimmed .faa.mafft = 0.2572    trimmed .clipkit = 0.4174

Plain treeness needs no alignment and is unaffected — it reproduces exactly (median 0.0501 on the same 249 files), which is how the alignment choice was isolated as the cause rather than the tree set or the tool.

phykit_batch_analysis takes the two independently, so pass them explicitly:

tu run phykit_batch_analysis '{"operation":"batch","function":"treeness_over_rcv",
  "directory":"<dir>","extension":".faa.mafft",
  "tree_directory":"<dir>","tree_extension":".faa.mafft.clipkit.treefile"}'

Gap percentage in these questions means the fraction of alignment columns containing at least one gap, not the fraction of all residues that are gaps. The two differ by an order of magnitude: with the residue definition no fungal ortholog exceeds 70% gaps, so a ">70% gaps" filter silently selects nothing.

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

Metric	Input	Description
Treeness	Newick	Internal / total branch length
RCV	FASTA/PHYLIP	Relative Composition Variability
Treeness/RCV	Both	Signal quality ratio
Tree Length	Newick	Sum of all branch lengths
Evolutionary Rate	Newick	Total length / num terminals
DVMC	Newick	Degree of Violation of Molecular Clock
Parsimony Sites	FASTA/PHYLIP	Sites with >=2 chars appearing >=2 times

---

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

Pattern	Method
"median X"	np.median(values)
"maximum X"	np.max(values)
"difference in median"	abs(np.median(a) - np.median(b))
"Mann-Whitney U"	stats.mannwhitneyu(a, b)[0]
"fold-change"	np.median(a) / np.median(b)

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

---

Interpretation

Metric	Good	Acceptable	Poor
Treeness	>0.8	0.5-0.8	<0.5
RCV	<0.2	0.2-0.5	>0.5
Treeness/RCV	>2.0	1.0-2.0	<1.0
Bootstrap	>95%	70-95%	<70%
Parsimony sites	>30%	10-30%	<10%

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.

The batch tool is parallel: ~250 trees finish in about 35 seconds. A per-tree shell loop takes ~9 minutes for the same work and is the single most common way these questions end with no answer at all — the run hits its turn or time budget mid-loop and reports "I'll report when it finishes" instead of a number. If you find yourself writing for f in *.treefile, stop and call the batch tool.

Supported function values include treeness, saturation, dvmc, long_branch_score, total_tree_length, parsimony_informative, treeness_over_rcv (alias toverr). dvmc and long_branch_score are covered — you do not need to loop for those.

Two-group comparisons (Mann-Whitney U, differences of medians). Questions comparing fungi against animals need one batch call per group, then the test on the two value lists — not a per-tree loop over both groups:

tu run phykit_batch_analysis '{"operation":"batch","function":"dvmc","directory":"<fungi>","extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"dvmc","directory":"<animals>","extension":".treefile"}'
# then scipy.stats.mannwhitneyu(fungi_values, animal_values)

Ask for values in the result when you need the full list for a test; the batch tool returns them for sets up to 50 and summary statistics always. For larger sets, compute the statistic from the per-group summaries the tool returns rather than re-deriving every value by hand.

PhyKIT column conventions — take the right one

Several PhyKIT subcommands print more than one number per file, and the value the question wants is usually not the first:

subcommand	prints	the value asked for
saturation	saturation <TAB> |saturation-1|	column 1 per phykit's docs; some sources report col 2 — say which you used
treeness_over_rcv	treeness/RCV <TAB> treeness <TAB> RCV	column 1, the ratio
parsimony_informative_sites	n_pi <TAB> n_total <TAB> %PIS	column 3 for a percentage

Taking saturation's first column gives exactly 1 - answer: a fungal set whose saturation is 0.6146 reports 0.3854 instead, and the two sum to 1.0000, which is the tell. phykit_batch_analysis already selects the right column for each function — another reason to call it rather than run the CLI yourself.

Commit the value you computed

Two failures in this benchmark came from computing the right number and then answering a different one:

• a tree-length ratio computed as 2.1775, then answered as 1.9 after re-reading "paired orthologs";
• an average treeness that listed 19 among the alternatives, then committed 10.

When a question is ambiguous, compute the reading you judge most literal, state the alternative in one clause, and answer with the value you actually computed. Do not replace a computed result with a re-derived one at the last step — if two readings are both defensible, give the computed number first and name the other, rather than silently switching.

PhyKIT column-position cheat sheet (parse output carefully)

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

Command	Output columns	Column to take
phykit saturation	saturation_value <TAB> abs(saturation-1)	col 1 is the "saturation value" (1 = no saturation; closer to 1 = less saturated). col 2 = |saturation - 1| (distance from no-saturation; higher = MORE saturated, less signal retained). Use col 1 for "saturation value" questions; col 2 for "distance from saturation"
phykit toverr (a.k.a. treeness_over_rcv)	treeness/RCV <TAB> treeness <TAB> RCV	col 1 (treeness/RCV ratio)
phykit long_branch_score -v (verbose)	taxon <TAB> score per line	aggregate scores per tree (mean)
phykit long_branch_score (no -v)	mean <TAB> median <TAB> 25%ile <TAB> 75%ile <TAB> min <TAB> max <TAB> std <TAB> var <TAB> n	col 1 (mean) for "mean LB score"
phykit patristic_distances (no -v)	summary stats line (same shape as LB)	col 1 (mean) for "mean patristic distance"

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