tooluniverse-phylogenetics

What this skill does

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

```bash
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)

```bash
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.
  Use clipkit alignment + treefile for tree-paired metrics.
- **Alignment-only**: just `<gene>.faa` + `<gene>.faa.mafft`. No
  trees, no clipkit. Used for parsim