---
name: busco-phylogeny
description: Generate phylogenies from genome assemblies using BUSCO/compleasm-based single-copy orthologs with scheduler-aware workflow generation
---

# BUSCO-based Phylogenomics Workflow Generator

This skill provides phylogenomics expertise for generating comprehensive, scheduler-aware workflows for phylogenetic inference from genome assemblies using single-copy orthologs.

## Purpose

This skill helps users generate phylogenies from genome assemblies by:
1. Handling mixed input (local files and NCBI accessions)
2. Creating scheduler-specific scripts (SLURM, PBS, cloud, local)
3. Setting up complete workflows from raw genomes to final trees
4. Providing quality control and recommendations
5. Supporting flexible software management (bioconda, Docker, custom)

## Available Resources

The skill provides access to these bundled resources:

### Scripts (`scripts/`)
- **`query_ncbi_assemblies.py`** - Query NCBI for available genome assemblies by taxon name (new!)
- **`download_ncbi_genomes.py`** - Download genomes from NCBI using BioProjects or Assembly accessions
- **`rename_genomes.py`** - Rename genome files with meaningful sample names (important!)
- **`generate_qc_report.sh`** - Generate quality control reports from compleasm results
- **`extract_orthologs.sh`** - Extract and reorganize single-copy orthologs
- **`run_aliscore.sh`** - Wrapper for Aliscore to identify randomly similar sequences (RSS)
- **`run_alicut.sh`** - Wrapper for ALICUT to remove RSS positions from alignments
- **`run_aliscore_alicut_batch.sh`** - Batch process all alignments through Aliscore + ALICUT
- **`convert_fasconcat_to_partition.py`** - Convert FASconCAT output to IQ-TREE partition format
- **`predownloaded_aliscore_alicut/`** - Pre-tested Aliscore and ALICUT Perl scripts

### Templates (`templates/`)
- **`slurm/`** - SLURM job scheduler templates
- **`pbs/`** - PBS/Torque job scheduler templates
- **`local/`** - Local machine templates (with GNU parallel)
- **`README.md`** - Complete template documentation

### References (`references/`)
- **`REFERENCE.md`** - Detailed technical reference including:
  - Sample naming best practices
  - BUSCO lineage datasets (complete list)
  - Resource recommendations (memory, CPUs, walltime)
  - Detailed step-by-step implementation guides
  - Quality control guidelines
  - Aliscore/ALICUT detailed guide
  - Tool citations and download links
  - Software installation guide
  - Common issues and troubleshooting

## Workflow Overview

The complete phylogenomics pipeline follows this sequence:

**Input Preparation** → **Ortholog Identification** → **Quality Control** → **Ortholog Extraction** → **Alignment** → **Trimming** → **Concatenation** → **Phylogenetic Inference**

## Initial User Questions

When a user requests phylogeny generation, gather the following information systematically:

### Step 1: Detect Computing Environment

Before asking questions, attempt to detect the local computing environment:

```bash
# Check for job schedulers
command -v sbatch >/dev/null 2>&1  # SLURM
command -v qsub >/dev/null 2>&1    # PBS/Torque
command -v parallel >/dev/null 2>&1  # GNU parallel
```

Report findings to the user, then confirm: **"I detected [X] on this machine. Will you be running the scripts here or on a different system?"**

### Required Information

Ask these questions to gather essential workflow parameters:

1. **Computing Environment**
   - Where will these scripts run? (SLURM cluster, PBS/Torque cluster, Cloud computing, Local machine)

2. **Input Data**
   - Local genome files, NCBI accessions, or both?
   - If NCBI: Do you already have Assembly accessions (GCA_*/GCF_*) or BioProject accessions (PRJNA*/PRJEB*/PRJDA*)?
   - If user doesn't have accessions: Offer to help find assemblies using `query_ncbi_assemblies.py` (see "STEP 0A: Query NCBI for Assemblies" below)
   - If local files: What are the file paths?

3. **Taxonomic Scope & Dataset Details**
   - What taxonomic group? (determines BUSCO lineage dataset)
   - How many taxa/genomes will be analyzed?
   - What is the approximate phylogenetic breadth? (species-level, genus-level, family-level, order-level, etc.)
   - See `references/REFERENCE.md` for complete lineage list

4. **Environment Management**
   - Use unified conda environment (default, recommended), or separate environments per tool?

5. **Resource Constraints**
   - How many CPU cores/threads to use in total? (Ask user to specify, do not auto-detect)
   - Available memory (RAM) per node/machine?
   - Maximum walltime for jobs?
   - See `references/REFERENCE.md` for resource recommendations

6. **Parallelization Strategy**

   Ask the user how they want to handle parallel processing:

   - **For job schedulers (SLURM/PBS)**:
     - Use array jobs for parallel steps? (Recommended: Yes)
     - Which steps to parallelize? (Steps 2, 5, 6, 8C recommended)

   - **For local machines**:
     - Use GNU parallel for parallel steps? (requires `parallel` installed)
     - How many concurrent jobs?

   - **For all systems**:
     - Optimize for maximum throughput or simplicity?

7. **Scheduler-Specific Configuration** (if using SLURM or PBS)
   - Account/Username for compute time charges
   - Partition/Queue to submit jobs to
   - Email notifications? (address and when: START, END, FAIL, ALL)
   - Job dependencies? (Recommended: Yes for linear workflow)
   - Output log directory? (Default: `logs/`)

8. **Alignment Trimming Preference**
   - Aliscore/ALICUT (traditional, thorough), trimAl (fast), BMGE (entropy-based), or ClipKit (modern)?

9. **Substitution Model Selection** (for IQ-TREE phylogenetic inference)

   **Context needed**: Taxonomic breadth, number of taxa, evolutionary rates

   **Action**: Fetch IQ-TREE model documentation and suggest appropriate amino acid substitution models based on dataset characteristics.

   Use the substitution model recommendation system (see "Substitution Model Recommendation" section below).

10. **Educational Goals**
   - Are you learning bioinformatics and would you like comprehensive explanations of each workflow step?
   - If yes: After completing each major workflow stage, offer to explain what the step accomplishes, why certain choices were made, and what best practices are being followed.
   - Store this preference to use throughout the workflow.

---

## Recommended Directory Structure

Organize analyses with dedicated folders for each pipeline step:

```
project_name/
├── logs/                          # All log files
├── 00_genomes/                    # Input genome assemblies
├── 01_busco_results/              # BUSCO/compleasm outputs
├── 02_qc/                         # Quality control reports
├── 03_extracted_orthologs/        # Extracted single-copy orthologs
├── 04_alignments/                 # Multiple sequence alignments
├── 05_trimmed/                    # Trimmed alignments
├── 06_concatenation/              # Supermatrix and partition files
├── 07_partition_search/           # Partition model selection
├── 08_concatenated_tree/          # Concatenated ML tree
├── 09_gene_trees/                 # Individual gene trees
├── 10_species_tree/               # ASTRAL species tree
└── scripts/                       # All analysis scripts
```

**Benefits**: Easy debugging, clear workflow progression, reproducibility, prevents root directory clutter.

---

## Template System

This skill uses a template-based system to reduce token usage and improve maintainability. Script templates are stored in the `templates/` directory and organized by computing environment.

### How to Use Templates

When generating scripts for users:

1. **Read the appropriate template** for their computing environment:
   ```
   Read("templates/slurm/02_compleasm_first.job")
   ```

2. **Replace placeholders** with user-specific values:
   - `TOTAL_THREADS` → e.g., `64`
   - `THREADS_PER_JOB` → e.g., `16`
   - `NUM_GENOMES` → e.g., `20`
   - `NUM_LOCI` → e.g., `2795`
   - `LINEAGE` → e.g., `insecta_odb10`
   - `MODEL_SET` → e.g., `LG,WAG,JTT,Q.pfam`

3. **Present the customized script** to the user with setup instructions

### Available Templates

Key templates by workflow step:
- **Step 0 (setup)**: Environment setup script in `references/REFERENCE.md`
- **Step 2 (compleasm)**: `02_compleasm_first`, `02_compleasm_parallel`
- **Step 8A (partition search)**: `08a_partition_search`
- **Step 8C (gene trees)**: `08c_gene_trees_array`, `08c_gene_trees_parallel`, `08c_gene_trees_serial`

See `templates/README.md` for complete template documentation.

---

## Substitution Model Recommendation

When asked about substitution model selection (Question 9), use this systematic approach:

### Step 1: Fetch IQ-TREE Documentation

Use WebFetch to retrieve current model information:
```
WebFetch(url="https://iqtree.github.io/doc/Substitution-Models",
         prompt="Extract all amino acid substitution models with descriptions and usage guidelines")
```

### Step 2: Analyze Dataset Characteristics

Consider these factors from user responses:
- **Taxonomic Scope**: Species/genus (shallow) vs. family/order (moderate) vs. class/phylum+ (deep)
- **Number of Taxa**: <20 (small), 20-50 (medium), >50 (large)
- **Evolutionary Rates**: Fast-evolving, moderate, or slow-evolving
- **Sequence Type**: Nuclear proteins, mitochondrial, or chloroplast

### Step 3: Recommend Models

Provide 3-5 appropriate models based on dataset characteristics. For detailed model recommendation matrices and taxonomically-targeted models, see `references/REFERENCE.md` section "Substitution Model Recommendation".

**General recommendations**:
- **Nuclear proteins (most common)**: LG, WAG, JTT, Q.pfam
- **Mitochondrial**: mtREV, mtZOA, mtMAM, mtART, mtVer, mtInv
- **Chloroplast**: cpREV
- **Taxonomically-targeted**: Q.bird, Q.mammal, Q.insect, Q.plant, Q.yeast (when applicable)

### Step 4: Present Recommendations

Format recommendations with justifications and explain how models will be used in IQ-TREE steps 8A and 8C.

### Step 5: Store Model Set

Store the final comma-separated model list (e.g., "LG,WAG,JTT,Q.pfam") for use in Step 8 template placeholders.

---

## Workflow Implementation

Once required information is gathered, guide the user through these steps. For each step, use templates where available and refer to `references/REFERENCE.md` for detailed implementation.

### STEP 0: Environment Setup

**ALWAYS start by generating a setup script** for the user's environment.

Use the unified conda environment setup script from `references/REFERENCE.md` (Section: "Software Installation Guide"). This creates a single conda environment with all necessary tools:
- compleasm, MAFFT, trimming tools (trimAl, ClipKit, BMGE)
- IQ-TREE, ASTRAL, Perl with BioPerl, GNU parallel
- Downloads and installs Aliscore/ALICUT Perl scripts

**Key points**:
- Users choose between mamba (faster) or conda
- Users choose between predownloaded Aliscore/ALICUT scripts (tested) or latest from GitHub
- All subsequent steps use `conda activate phylo` (the unified environment)

See `references/REFERENCE.md` for the complete setup script template.

---

### STEP 0A: Query NCBI for Assemblies (Optional)

**Use this step when**: User wants to use NCBI data but doesn't have specific assembly accessions yet.

This optional preliminary step helps users discover available genome assemblies by taxon name before proceeding with the main workflow.

**IMPORTANT**: This step requires the NCBI datasets CLI tool (`datasets`) to be installed **on the local machine where Claude is running**. If Claude is running locally but preparing scripts for a remote HPC cluster, the `datasets` tool must be installed locally to query NCBI and find assemblies. The query functionality cannot be used if preparing scripts for remote execution without local access to the datasets CLI tool.

#### When to Offer This Step

Offer this step when:
- User wants to analyze genomes from NCBI
- User doesn't have specific Assembly or BioProject accessions
- User mentions a taxonomic group (e.g., "I want to build a phylogeny for beetles")
- The `datasets` CLI tool is available locally (check with `command -v datasets`)

If the datasets tool is not available locally, inform the user they can:
1. Install it via conda: `conda install -c conda-forge ncbi-datasets-cli`
2. Download it manually from https://www.ncbi.nlm.nih.gov/datasets/docs/v2/download-and-install/
3. Use the NCBI website to manually search for assemblies

#### Workflow

1. **Check for datasets CLI tool availability**:
   ```bash
   command -v datasets >/dev/null 2>&1 && echo "datasets CLI found" || echo "datasets CLI not found"
   ```

2. **Ask for focal taxon**: Request the taxonomic group of interest
   - Examples: "Coleoptera", "Drosophila", "Apis mellifera"
   - Can be at any taxonomic level (order, family, genus, species)

3. **Suggest quality filtering criteria**: Before querying, ask the user about quality preferences to help select the best assemblies for phylogenomics:

   **Assembly Quality Considerations:**
   - **Assembly Level**: Chromosome-level assemblies are preferred over Scaffold over Contig
     - Chromosome: Most complete, best for phylogenomics
     - Scaffold: Medium quality, may have gaps between contigs
     - Contig: Lower quality, more fragmented

   - **Contig N50**: Higher N50 indicates better assembly contiguity
     - Excellent: N50 > 10 Mbp
     - Good: N50 > 1 Mbp
     - Poor: N50 < 100 kbp

   - **Annotation Status**: Annotated assemblies have gene predictions and BUSCO scores
     - Useful for assessing completeness
     - BUSCO > 95% is excellent for phylogenomics

   - **RefSeq vs GenBank**:
     - RefSeq (GCF_*): Curated, higher quality, recommended
     - GenBank (GCA_*): Submitted assemblies, variable quality

   **Suggested Questions for Users:**
   - "Do you prefer chromosome-level assemblies only, or are scaffold-level acceptable?"
   - "Should I filter for high-quality assemblies (e.g., contig N50 > 1 Mbp)?"
   - "Do you want to see RefSeq assemblies only (higher quality)?"
   - "Would you like to see quality metrics (N50, BUSCO scores) to help you decide?"

4. **Query NCBI using the script**: Use `scripts/query_ncbi_assemblies.py` with appropriate filters based on user preferences:

   ```bash
   # Basic query with quality metrics
   python scripts/query_ncbi_assemblies.py --taxon "Coleoptera" --show-quality

   # Filter for chromosome-level assemblies
   python scripts/query_ncbi_assemblies.py --taxon "Drosophila" --assembly-level Chromosome --show-quality

   # Filter for high-quality assemblies (N50 > 1 Mbp)
   python scripts/query_ncbi_assemblies.py --taxon "Apis" --min-contig-n50 1000000 --show-quality

   # RefSeq annotated assemblies only (highest quality)
   python scripts/query_ncbi_assemblies.py --taxon "Felidae" --refseq-only --annotated --show-quality

   # Combined filters: Chromosome-level, high N50, annotated
   python scripts/query_ncbi_assemblies.py --taxon "Coleoptera" \
     --assembly-level Chromosome \
     --min-contig-n50 1000000 \
     --annotated \
     --show-quality \
     --save assembly_accessions.txt
   ```

   **Available Filtering Options:**
   - `--assembly-level {Chromosome,Scaffold,Contig}` - Filter by assembly level
   - `--min-contig-n50 N` - Minimum contig N50 in base pairs (e.g., 1000000 for 1 Mbp)
   - `--annotated` - Only return assemblies with gene annotations
   - `--refseq-only` - Only return RefSeq (GCF_*) assemblies
   - `--show-quality` - Display quality metrics (N50, BUSCO scores) in output
   - `--max-results N` - Maximum number of results (default: 20)
   - `--save FILE` - Save accessions to file for later download

5. **Present results to user**: The script displays:
   - Assembly accession (GCA_* or GCF_*)
   - Organism name
   - Assembly level (Chromosome, Scaffold, Contig)
   - Assembly name
   - **If --show-quality:** Contig N50, BUSCO completeness %, Annotation status
   - **Quality summary:** Count of assemblies by level

6. **Help user select assemblies**: Review results with user and recommend best choices:
   - Prioritize chromosome-level assemblies with high N50
   - Prefer RefSeq (GCF_*) over GenBank (GCA_*) when available
   - Look for BUSCO completeness > 95% when available
   - Aim for consistent quality across all taxa in the phylogeny
   - Consider phylogenetic breadth (species coverage)
   - Consider data quality (RefSeq > GenBank when available)

7. **Collect selected accessions**: Compile the list of chosen assembly accessions

8. **Proceed to STEP 1**: Use the selected accessions with `download_ncbi_genomes.py`

#### Tips for Assembly Selection

- **Assembly Level**: Chromosome-level assemblies are most complete, followed by Scaffold, then Contig
- **RefSeq vs GenBank**: RefSeq (GCF_*) assemblies undergo additional curation; GenBank (GCA_*) are submitter-provided
- **Taxonomic Sampling**: For phylogenetics, aim for representative sampling across the taxonomic group
- **Quality over Quantity**: Better to have 20 high-quality assemblies than 100 poor-quality ones

---

### STEP 1: Download NCBI Genomes (if applicable)

If user provided NCBI accessions, use `scripts/download_ncbi_genomes.py`:

**For BioProjects**:
```bash
python scripts/download_ncbi_genomes.py --bioprojects PRJNA12345 -o genomes.zip
unzip genomes.zip
```

**For Assembly Accessions**:
```bash
python scripts/download_ncbi_genomes.py --assemblies GCA_123456789.1 -o genomes.zip
unzip genomes.zip
```

**IMPORTANT**: After download, genomes must be renamed with meaningful sample names (format: `[ACCESSION]_[SPECIES_NAME]`). Sample names appear in final phylogenetic trees.

Generate a script that:
1. Finds all downloaded FASTA files in ncbi_dataset directory structure
2. Moves/renames files to main genomes directory with meaningful names
3. Includes any local genome files
4. Creates final genome_list.txt with ALL genomes (local + downloaded)

See `references/REFERENCE.md` section "Sample Naming Best Practices" for detailed guidelines.

---

### STEP 2: Ortholog Identification with compleasm

Activate the unified environment and run compleasm on all genomes to identify single-copy orthologs.

**Key considerations**:
- First genome must run alone to download lineage database
- Remaining genomes can run in parallel
- Thread allocation: Miniprot scales well up to ~16-32 threads per genome

**Threading guidelines**: See `references/REFERENCE.md` for recommended thread allocation table.

**Generate scripts using templates**:
- **SLURM**: Read templates `02_compleasm_first.job` and `02_compleasm_parallel.job`
- **PBS**: Read templates `02_compleasm_first.job` and `02_compleasm_parallel.job`
- **Local**: Read templates `02_compleasm_first.sh` and `02_compleasm_parallel.sh`

Replace placeholders: `TOTAL_THREADS`, `THREADS_PER_JOB`, `NUM_GENOMES`, `LINEAGE`

For detailed implementation examples, see `references/REFERENCE.md` section "Ortholog Identification Implementation".

---

### STEP 3: Quality Control

After compleasm completes, generate QC report using `scripts/generate_qc_report.sh`:

```bash
bash scripts/generate_qc_report.sh qc_report.csv
```

Provide interpretation:
- **>95% complete**: Excellent, retain
- **90-95% complete**: Good, retain
- **85-90% complete**: Acceptable, case-by-case
- **70-85% complete**: Questionable, consider excluding
- **<70% complete**: Poor, recommend excluding

See `references/REFERENCE.md` section "Quality Control Guidelines" for detailed assessment criteria.

---

### STEP 4: Ortholog Extraction

Use `scripts/extract_orthologs.sh` to extract single-copy orthologs:

```bash
bash scripts/extract_orthologs.sh LINEAGE_NAME
```

This generates per-locus unaligned FASTA files in `single_copy_orthologs/unaligned_aa/`.

---

### STEP 5: Alignment with MAFFT

Activate the unified environment (`conda activate phylo`) which contains MAFFT.

Create locus list, then generate alignment scripts:
```bash
cd sin