{
  "skill": "alterlab-phylogenetics",
  "evals": [
    {
      "id": "full-ml-tree-from-fasta",
      "prompt": "I have a FASTA of about 150 unaligned protein sequences from a gene family. Align them, trim the bad columns, infer a maximum-likelihood tree with automatic model selection and 1000 bootstraps, then root it and give me a figure.",
      "expected_output": "Invokes alterlab-phylogenetics: runs MAFFT (linsi/einsi for ~150 seqs) to align, TrimAl to trim, IQ-TREE 2 with -m TEST for model selection and -B 1000 ultrafast bootstraps, then midpoint-roots and visualizes with ETE3. Full align -> trim -> ML tree -> visualize pipeline.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "IQ-TREE" },
        { "type": "behavior", "value": "Sequences MAFFT alignment, optional TrimAl, then IQ-TREE 2 with model selection (-m TEST) and -B 1000 bootstraps." }
      ]
    },
    {
      "id": "fasttree-large-dataset",
      "prompt": "I've got an alignment of around 8000 nucleotide sequences and IQ-TREE is taking forever. What's the fast way to get an approximate ML tree out of this?",
      "expected_output": "Invokes alterlab-phylogenetics: recommends FastTree for the large dataset (FastTree -nt -gtr for nucleotides), noting it is 10-100x faster than IQ-TREE for thousands of sequences, as the appropriate fast approximate-ML choice.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "FastTree" },
        { "type": "behavior", "value": "Recommends FastTree (with -nt -gtr) for the large alignment as the fast approximate-ML alternative to IQ-TREE." }
      ]
    },
    {
      "id": "mafft-method-selection",
      "prompt": "I'm aligning around 60 closely related sequences and accuracy matters more than speed. Which MAFFT mode should I use and how do I run it?",
      "expected_output": "Invokes alterlab-phylogenetics: recommends an accurate MAFFT mode for the small set (linsi via --localpair --maxiterate 1000, or einsi via --genafpair) since accuracy is preferred for <200 sequences, and shows the run command.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "MAFFT" },
        { "type": "behavior", "value": "Selects an accurate MAFFT mode (linsi/einsi) for the small sequence set per the method-selection guidance." }
      ]
    },
    {
      "id": "viral-molecular-clock-recombination",
      "prompt": "I'm doing phylodynamics on a set of dated SARS-CoV-2 genomes and want to estimate divergence dates. Anything I should check before building the tree, and how do I add the temporal calibration?",
      "expected_output": "Invokes alterlab-phylogenetics: warns to check for recombination first (RDP4/GARD) in viral sequences, then builds the tree and adds molecular-clock / temporal calibration via IQ-TREE 2 --date with dates.txt (and --date-CI). Phylodynamics/molecular-clock workflow.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Advises a recombination check (RDP4/GARD) before tree-building and uses IQ-TREE --date for temporal/molecular-clock calibration." }
      ]
    },
    {
      "id": "near-miss-etetoolkit",
      "prompt": "I already have a finished Newick tree. I just want to prune it to a subset of taxa keeping branch lengths, detect duplication versus speciation events on the gene tree, and compute the Robinson-Foulds distance to a second tree.",
      "expected_output": "Should NOT trigger alterlab-phylogenetics. Manipulating an existing tree (prune, duplication/speciation event detection via get_descendant_evol_events, robinson_foulds comparison) is alterlab-etetoolkit territory. phylogenetics is the upstream align-and-infer pipeline and defers pure tree manipulation/comparison to the ETE Toolkit skill.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "etetoolkit" }
      ]
    },
    {
      "id": "near-miss-biopython",
      "prompt": "I just need to parse a multi-record FASTA file, reverse-complement each sequence, translate it to protein, and write the results back out to a new FASTA.",
      "expected_output": "Should NOT trigger alterlab-phylogenetics. Basic sequence parsing, reverse-complement, translation, and FASTA I/O is general sequence handling for alterlab-biopython, with no alignment or tree inference involved. phylogenetics is specifically for reconstructing evolutionary trees.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "biopython" }
      ]
    }
  ]
}

# IQ-TREE Phylogenetic Inference Reference

> Commands below call `iqtree`. bioconda's `iqtree` package is IQ-TREE 3 (3.1.x) and installs
> `iqtree` and `iqtree3`; the older 2.x package installed `iqtree2`. The options here are the
> same on both.

## Basic Command Syntax

```bash
iqtree -s alignment.fasta --prefix output -B 1000 -T AUTO --redo   # -m MFP is the default
```

## Key Parameters

| Flag | Description | Default |
|------|-------------|---------|
| `-s` | Input alignment file | Required |
| `--prefix` | Output file prefix | alignment name |
| `-m` | Substitution model (or TEST) | GTR+G |
| `-B` | Ultrafast bootstrap replicates | Off |
| `-b` | Standard bootstrap replicates (slow) | Off |
| `-T` | Number of threads (or AUTO) | 1 |
| `-o` | Outgroup taxa name(s) | None (unrooted) |
| `--redo` | Overwrite existing results | Off |
| `-alrt` | SH-aLRT test replicates | Off |

## Model Selection

```bash
# Full model testing (automatically selects best model)
iqtree -s alignment.fasta -m TEST --prefix test_run -B 1000 -T 4

# Specify model explicitly
iqtree -s alignment.fasta -m GTR+G4 --prefix gtr_run -B 1000

# Protein sequences
iqtree -s protein.fasta -m TEST --prefix prot_tree -B 1000

# Codon-based analysis
iqtree -s codon.fasta -m GY --prefix codon_tree -B 1000
```

## Bootstrapping Methods

### Ultrafast Bootstrap (UFBoot, recommended)
```bash
iqtree -s alignment.fasta -B 1000  # 1000 replicates
# Values ≥95 are reliable
# ~10× faster than standard bootstrap
```

### Standard Bootstrap
```bash
iqtree -s alignment.fasta -b 100  # 100 replicates (very slow)
```

### SH-aLRT Test (fast alternative)
```bash
iqtree -s alignment.fasta -alrt 1000 -B 1000  # Both SH-aLRT and UFBoot
# SH-aLRT ≥80 AND UFBoot ≥95 = well-supported branch
```

## Branch Support Interpretation

| Bootstrap Value | Interpretation |
|----------------|----------------|
| ≥ 95 | Well-supported (strongly supported) |
| 70–94 | Moderately supported |
| 50–69 | Weakly supported |
| < 50 | Unreliable (not supported) |

## Output Files

| File | Description |
|------|-------------|
| `{prefix}.treefile` | Best ML tree in Newick format |
| `{prefix}.iqtree` | Full analysis report |
| `{prefix}.log` | Computation log |
| `{prefix}.contree` | Consensus tree from bootstrap |
| `{prefix}.splits.nex` | Network splits |
| `{prefix}.bionj` | BioNJ starting tree |
| `{prefix}.model.gz` | Saved model parameters |

## Advanced Analyses

### Molecular Clock (Dating)

```bash
# Temporal analysis with sampling dates
iqtree -s alignment.fasta -m GTR+G \
        --date dates.tsv \           # Tab-separated: taxon_name  YYYY-MM-DD
        --clock-test \               # Test for clock-like evolution
        --date-CI 95 \              # 95% CI for node dates
        --prefix dated_tree
```

### Concordance Factors

```bash
# Gene concordance factor (gCF) - requires multiple gene alignments
iqtree --gcf gene_trees.nwk \
        --tree main_tree.treefile \
        --cf-verbose \
        --prefix cf_analysis
```

### Ancestral Sequence Reconstruction

```bash
iqtree -s alignment.fasta -m LG+G4 \
        -asr \                      # Marginal ancestral state reconstruction
        --prefix anc_tree
# Output: {prefix}.state (ancestral sequences per node)
```

### Partition Model (Multi-Gene)

```bash
# Create partition file (partitions.txt):
# DNA, gene1 = 1-500
# DNA, gene2 = 501-1000

iqtree -s concat_alignment.fasta \
        -p partitions.txt \
        -m TEST \
        -B 1000 \
        --prefix partition_tree
```

## IQ-TREE Log Parsing

```python
def parse_iqtree_log(log_file: str) -> dict:
    """Extract key results from IQ-TREE log file."""
    results = {}
    with open(log_file) as f:
        for line in f:
            if "Best-fit model" in line:
                results["best_model"] = line.split(":")[1].strip()
            elif "Log-likelihood of the tree:" in line:
                results["log_likelihood"] = float(line.split(":")[1].strip())
            elif "Number of free parameters" in line:
                results["free_params"] = int(line.split(":")[1].strip())
            elif "Akaike information criterion" in line:
                results["AIC"] = float(line.split(":")[1].strip())
            elif "Bayesian information criterion" in line:
                results["BIC"] = float(line.split(":")[1].strip())
            elif "Total CPU time used" in line:
                results["cpu_time"] = line.split(":")[1].strip()
    return results

# Example:
# results = parse_iqtree_log("output.log")
# print(f"Best model: {results['best_model']}")
# print(f"Log-likelihood: {results['log_likelihood']:.2f}")
```

## Common Issues and Solutions

| Issue | Likely Cause | Solution |
|-------|-------------|---------|
| All bootstrap values = 0 | Too few taxa | Need ≥4 taxa for bootstrap |
| Very long branches | Alignment artifacts | Re-trim alignment; check for outliers |
| Memory error | Too many sequences | Use FastTree; or reduce `-T` to 1 |
| Poor model fit | Wrong alphabet | Check nucleotide vs. protein specification |
| Identical sequences | Duplicate sequences | Remove duplicates before alignment |

## MAFFT Alignment Guide

```bash
# Accurate (< 200 sequences)
mafft --localpair --maxiterate 1000 input.fasta > aligned.fasta

# Medium (200-1000 sequences)
mafft --auto input.fasta > aligned.fasta

# Fast (> 1000 sequences)
mafft --fftns input.fasta > aligned.fasta

# Very large (> 10000 sequences)
mafft --retree 1 input.fasta > aligned.fasta

# Using multiple threads
mafft --thread 8 --auto input.fasta > aligned.fasta
```

"""
Phylogenetic Analysis Pipeline
===============================
Complete workflow: MAFFT alignment → IQ-TREE 3 tree → ETE3 visualization.

Requirements (CLI binaries are NOT on PyPI):
    conda install -c bioconda mafft iqtree fasttree   # or Homebrew: brew install mafft fasttree; brew install brewsci/bio/iqtree
    uv pip install "ete3==3.1.3"                       # visualization only (optional)

Usage:
    python phylogenetic_analysis.py sequences.fasta --type nt --threads 4
    python phylogenetic_analysis.py proteins.fasta --type aa --fasttree
"""

import argparse
import os
import subprocess
import sys
from pathlib import Path
from shutil import which

# bioconda's `iqtree` package is IQ-TREE 3 and ships `iqtree` + `iqtree3`; the 2.x
# package shipped `iqtree2`. Resolve whichever is on PATH instead of hardcoding one.
IQTREE_CANDIDATES = ("iqtree", "iqtree3", "iqtree2")


def find_iqtree():
    """Return the first IQ-TREE executable on PATH, or None."""
    for name in IQTREE_CANDIDATES:
        if which(name) is not None:
            return name
    return None


def check_dependencies(use_fasttree: bool = False):
    """Check that the required CLI tools are on PATH (binaries are not pip-installable)."""
    tools = {"mafft": "conda install -c bioconda mafft  (or: brew install mafft)"}
    if use_fasttree:
        tools["FastTree"] = "conda install -c bioconda fasttree  (or: brew install fasttree)"

    missing = [f"  {tool}: {hint}" for tool, hint in tools.items() if which(tool) is None]
    if not use_fasttree and find_iqtree() is None:
        missing.append(
            "  iqtree (or iqtree3/iqtree2): "
            "conda install -c bioconda iqtree  (or: brew install brewsci/bio/iqtree)"
        )
    if missing:
        print("Missing dependencies:")
        for m in missing:
            print(m)
        sys.exit(1)
    print("All dependencies found.")


def count_sequences(fasta_file: str) -> int:
    """Count sequences in a FASTA file."""
    with open(fasta_file) as f:
        return sum(1 for line in f if line.startswith('>'))


def run_mafft(input_fasta: str, output_fasta: str, n_threads: int = 4,
               method: str = "auto") -> str:
    """Run MAFFT multiple sequence alignment."""
    n_seqs = count_sequences(input_fasta)
    print(f"MAFFT: Aligning {n_seqs} sequences...")

    # Auto-select method based on dataset size
    if method == "auto":
        if n_seqs <= 200:
            cmd = ["mafft", "--localpair", "--maxiterate", "1000",
                   "--thread", str(n_threads), "--inputorder", input_fasta]
        elif n_seqs <= 1000:
            cmd = ["mafft", "--auto", "--thread", str(n_threads),
                   "--inputorder", input_fasta]
        else:
            cmd = ["mafft", "--fftns", "--thread", str(n_threads),
                   "--inputorder", input_fasta]
    else:
        cmd = ["mafft", f"--{method}", "--thread", str(n_threads),
               "--inputorder", input_fasta]

    with open(output_fasta, 'w') as out:
        result = subprocess.run(cmd, stdout=out, stderr=subprocess.PIPE, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"MAFFT failed:\n{result.stderr[:500]}")

    print(f"  Alignment complete → {output_fasta}")
    return output_fasta


def run_iqtree(aligned_fasta: str, prefix: str, seq_type: str = "nt",
                bootstrap: int = 1000, n_threads: int = 4,
                outgroup: str = None) -> str:
    """Run IQ-TREE phylogenetic inference (IQ-TREE 2 or 3, whichever is installed)."""
    exe = find_iqtree()
    if exe is None:
        raise RuntimeError(
            "No IQ-TREE executable found on PATH "
            f"(looked for {', '.join(IQTREE_CANDIDATES)}). "
            "Install with: conda install -c bioconda iqtree"
        )
    print(f"{exe}: Building maximum likelihood tree...")

    cmd = [
        exe,
        "-s", aligned_fasta,
        "--prefix", prefix,
        "-m", "MFP",            # ModelFinder Plus (IQ-TREE's default model search)
        "-B", str(bootstrap),   # Ultrafast bootstrap
        "-T", str(n_threads),
        "--redo",
        "-alrt", "1000",        # SH-aLRT test
    ]

    if outgroup:
        cmd += ["-o", outgroup]

    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"IQ-TREE failed:\n{result.stderr[:500]}")

    tree_file = f"{prefix}.treefile"

    # Extract best model from log
    log_file = f"{prefix}.log"
    if os.path.exists(log_file):
        with open(log_file) as f:
            for line in f:
                if "Best-fit model" in line:
                    print(f"  {line.strip()}")

    print(f"  Tree saved → {tree_file}")
    return tree_file


def run_fasttree(aligned_fasta: str, output_tree: str, seq_type: str = "nt") -> str:
    """Run FastTree (faster alternative for large datasets)."""
    print("FastTree: Building approximate ML tree (faster)...")

    if seq_type == "nt":
        cmd = ["FastTree", "-nt", "-gtr", "-gamma", aligned_fasta]
    else:
        cmd = ["FastTree", "-lg", "-gamma", aligned_fasta]

    with open(output_tree, 'w') as out:
        result = subprocess.run(cmd, stdout=out, stderr=subprocess.PIPE, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"FastTree failed:\n{result.stderr[:500]}")

    print(f"  Tree saved → {output_tree}")
    return output_tree


def visualize_tree(tree_file: str, output_png: str, outgroup: str = None) -> None:
    """Visualize the phylogenetic tree with ETE3."""
    try:
        from ete3 import Tree, TreeStyle
    except ImportError:
        print("ETE3 not installed. Skipping visualization.")
        print('  Install: uv pip install "ete3==3.1.3"')
        return

    t = Tree(tree_file)

    # Root the tree
    if outgroup and outgroup in [leaf.name for leaf in t.get_leaves()]:
        t.set_outgroup(outgroup)
        print(f"  Rooted at outgroup: {outgroup}")
    else:
        # Midpoint rooting
        t.set_outgroup(t.get_midpoint_outgroup())
        print("  Applied midpoint rooting")

    # Style
    ts = TreeStyle()
    ts.show_leaf_name = True
    ts.show_branch_support = True
    ts.mode = "r"  # rectangular

    try:
        t.render(output_png, tree_style=ts, w=800, units="px")
        print(f"  Visualization saved → {output_png}")
    except Exception as e:
        print(f"  Visualization failed (display issue?): {e}")
        # Save tree in Newick format as fallback
        rooted_nwk = output_png.replace(".png", "_rooted.nwk")
        t.write(format=1, outfile=rooted_nwk)
        print(f"  Rooted tree saved → {rooted_nwk}")


def tree_summary(tree_file: str) -> dict:
    """Print summary statistics for the tree."""
    try:
        from ete3 import Tree
        t = Tree(tree_file)
        t.set_outgroup(t.get_midpoint_outgroup())

        leaves = t.get_leaves()
        branch_lengths = [n.dist for n in t.traverse() if n.dist > 0]

        stats = {
            "n_taxa": len(leaves),
            "total_branch_length": sum(branch_lengths),
            "mean_branch_length": sum(branch_lengths) / len(branch_lengths) if branch_lengths else 0,
            "max_branch_length": max(branch_lengths) if branch_lengths else 0,
        }

        print("\nTree Summary:")
        for k, v in stats.items():
            if isinstance(v, float):
                print(f"  {k}: {v:.6f}")
            else:
                print(f"  {k}: {v}")

        return stats
    except Exception as e:
        print(f"Could not compute tree stats: {e}")
        return {}


def main():
    parser = argparse.ArgumentParser(description="Phylogenetic analysis pipeline")
    parser.add_argument("input", help="Input FASTA file (unaligned)")
    parser.add_argument("--type", choices=["nt", "aa"], default="nt",
                        help="Sequence type: nt (nucleotide) or aa (amino acid)")
    parser.add_argument("--threads", type=int, default=4, help="Number of threads")
    parser.add_argument("--bootstrap", type=int, default=1000,
                        help="Bootstrap replicates for IQ-TREE")
    parser.add_argument("--fasttree", action="store_true",
                        help="Use FastTree instead of IQ-TREE (faster, less accurate)")
    parser.add_argument("--outgroup", help="Outgroup taxon name for rooting")
    parser.add_argument("--mafft-method", default="auto",
                        choices=["auto", "linsi", "einsi", "fftnsi", "fftns"],
                        help="MAFFT alignment method")
    parser.add_argument("--output-dir", default="phylo_results",
                        help="Output directory")

    args = parser.parse_args()

    check_dependencies(use_fasttree=args.fasttree)

    # Setup
    os.makedirs(args.output_dir, exist_ok=True)
    prefix = os.path.join(args.output_dir, Path(args.input).stem)

    print("=" * 60)
    print("Phylogenetic Analysis Pipeline")
    print("=" * 60)
    print(f"Input: {args.input}")
    print(f"Sequence type: {args.type}")
    print(f"Output dir: {args.output_dir}")

    # Step 1: Multiple Sequence Alignment
    print("\n[Step 1/3] Multiple Sequence Alignment (MAFFT)")
    aligned = run_mafft(
        args.input,
        f"{prefix}_aligned.fasta",
        n_threads=args.threads,
        method=args.mafft_method
    )

    # Step 2: Tree Inference
    print("\n[Step 2/3] Tree Inference")
    if args.fasttree:
        tree_file = run_fasttree(aligned, f"{prefix}.tree", seq_type=args.type)
    else:
        tree_file = run_iqtree(
            aligned, prefix,
            seq_type=args.type,
            bootstrap=args.bootstrap,
            n_threads=args.threads,
            outgroup=args.outgroup
        )

    # Step 3: Visualization
    print("\n[Step 3/3] Visualization (ETE3)")
    visualize_tree(tree_file, f"{prefix}_tree.png", outgroup=args.outgroup)
    tree_summary(tree_file)

    print("\n" + "=" * 60)
    print("Analysis complete!")
    print(f"Key outputs:")
    print(f"  Aligned sequences: {aligned}")
    print(f"  Tree file: {tree_file}")
    print(f"  Visualization: {prefix}_tree.png")


if __name__ == "__main__":
    main()