003 Name Skill 7fabca72

---
name: bio-ecological-genomics-species-delimitation
description: Delimits species boundaries from molecular data using distance-based (ASAP), tree-based (bPTP, GMYC), and coalescent (BPP) methods. Compares multiple delimitation results with delimtools. Use when delineating putative species from DNA barcoding data, resolving cryptic species complexes, or validating taxonomic assignments. Emphasizes multi-method consensus following integrative taxonomy best practice.
tool_type: mixed
primary_tool: ASAP
---

# Species Delimitation

Delimits putative species from molecular data using complementary distance-based, tree-based, and coalescent methods.

## Overview of Methods

| Method | Input | Approach | Strengths | Limitations |
|--------|-------|----------|-----------|-------------|
| ASAP | Aligned sequences | Distance-based partitioning | Fast, automatic, ranks partitions | Single locus only |
| bPTP | Rooted phylogeny | Bayesian branching rates | Bayesian support values | Over-splits with structure |
| GMYC | Ultrametric tree | Coalescent/speciation transition | Well-established theory | Requires ultrametric tree |
| BPP | Multi-locus alignment | Full coalescent model | Multi-locus, statistically rigorous | Computationally intensive |

## ASAP (Assemble Species by Automatic Partitioning)

Distance-based method that finds the optimal partition of sequences into putative species:

```bash
# ASAP web interface: https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html
# Upload aligned FASTA, select distance model

# Command-line ASAP (if available)
ASAP -i aligned_sequences.fasta -d K2P -o asap_results/
# -d K2P: Kimura 2-parameter distance (standard for COI barcoding)
# -d p: uncorrected p-distance (for closely related taxa)
# -d JC: Jukes-Cantor (simplest substitution model)
```

### Interpreting ASAP Results

```
# ASAP ranks partitions by asap-score (lower = better)
# Reports: partition rank, number of groups, asap-score, p-value, threshold distance
# Always examine top 5-10 partitions, not just the best one
# Large gaps between consecutive asap-scores indicate robust partitioning
```

### Python ASAP-Style Analysis

```python
from Bio import AlignIO, Phylo
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform
import numpy as np

alignment = AlignIO.read('aligned_sequences.fasta', 'fasta')

calculator = DistanceCalculator('identity')
dm = calculator.get_distance(alignment)

names = [r.id for r in alignment]
n = len(names)
dist_array = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        dist_array[i][j] = dm[names[i], names[j]]

condensed = squareform(dist_array)
Z = linkage(condensed, method='average')

# Scan thresholds to identify barcode gap
# COI barcode gap typically at 2-3% for animals
thresholds = np.arange(0.01, 0.10, 0.005)
for t in thresholds:
    clusters = fcluster(Z, t=t, criterion='distance')
    n_clusters = len(set(clusters))
    print(f'Threshold {t:.3f}: {n_clusters} groups')
```

## bPTP (Bayesian Poisson Tree Processes)

Tree-based method that models speciation as a branching rate shift:

```bash
# iTaxoTools PTP (Python 3 compatible)
pip install PTP-pyqt5

# Run bPTP from command line
python -m PTP.PTP -t rooted_tree.nwk -o bptp_results \
    --type bayesian --ngen 100000 --burnin 0.1 --seed 42

# --type bayesian: Bayesian MCMC (preferred over ML for support values)
# --ngen 100000: MCMC generations; increase for large trees
# --burnin 0.1: discard first 10% of samples
```

### Interpreting bPTP Output

```
# Support values for each species partition:
#   > 0.95: strong support
#   0.80-0.95: moderate support
#   < 0.80: weak support, boundary uncertain
# bPTP can over-split when populations have strong geographic structure
# Compare against distance-based methods to identify over-splitting
```

## GMYC (Generalized Mixed Yule Coalescent)

Models the transition between interspecific (Yule) and intraspecific (coalescent) branching rates on an ultrametric tree:

```r
library(splits)
library(ape)

# Step 1: Prepare ultrametric tree
# GMYC requires an ultrametric (time-calibrated) tree
tree <- read.tree('rooted_tree.nwk')

# chronos() for quick ultrametric conversion (penalized likelihood)
# lambda=1: smoothing parameter; higher = more clock-like
ultrametric_tree <- chronos(tree, lambda = 1)
class(ultrametric_tree) <- 'phylo'

# Verify ultrametric property
is.ultrametric(ultrametric_tree)

# Step 2: Run single-threshold GMYC
gmyc_result <- gmyc(ultrametric_tree, method = 'single')

# Summary: number of entities, confidence interval, likelihood ratio test
summary(gmyc_result)

# Number of putative species (ML entities)
cat('GMYC species:', gmyc_result$entity[1], '\n')

# Likelihood ratio test p-value
# p < 0.05: significant speciation-coalescent transition detected
cat('LR test p-value:', gmyc_result$p.value[1], '\n')

# Step 3: Extract species assignments
species_list <- spec.list(gmyc_result)
cat('Species assignments:\n')
print(species_list)
```

### Multiple-Threshold GMYC

```r
# Allows different coalescent-to-speciation transitions across the tree
# Better for datasets with variable population sizes or divergence times
gmyc_multi <- gmyc(ultrametric_tree, method = 'multiple')
summary(gmyc_multi)
```

### GMYC Visualization

```r
# Color-coded tree by species partition
pdf('gmyc_species_tree.pdf', width = 10, height = 12)
plot(gmyc_result, cex = 0.6)
dev.off()

# Support surface: likelihood across threshold positions
pdf('gmyc_support.pdf', width = 8, height = 5)
plot(gmyc_result, type = 'support')
dev.off()
```

## BPP (Bayesian Phylogenetics and Phylogeography)

Full multispecies coalescent model for rigorous species delimitation:

```
# BPP control file for species delimitation (A10 analysis)
# A10: joint species delimitation and species tree estimation

seed = 12345

seqfile = alignment.phy
Imapfile = imap.txt      # maps individuals to putative species
outfile = bpp_results.txt
mcmcfile = bpp_mcmc.txt

speciesdelimitation = 1   # 1 = delimitation mode
speciestree = 1           # 1 = estimate species tree
speciesmodelprior = 1     # 1 = uniform prior on all rooted trees

# species&tree: species names, sample sizes, starting tree
species&tree = 4  speciesA speciesB speciesC speciesD
                  10 8 12 6
                  ((speciesA, speciesB), (speciesC, speciesD));

# Prior on theta (ancestral population size)
# thetaprior = alpha beta: inverse-gamma(alpha, beta)
# alpha=3, beta=0.004: mean=0.002 (typical for moderate diversity)
thetaprior = 3 0.004

# Prior on tau (divergence time)
# tauprior = alpha beta: inverse-gamma(alpha, beta)
# alpha=3, beta=0.002: mean=0.001 (relatively recent divergence)
tauprior = 3 0.002

# MCMC settings
nsample = 100000
sampfreq = 2
burnin = 50000
```

### Interpreting BPP Results

```
# Posterior probability for each delimitation model
# PP > 0.95: strong support for that delimitation
# PP 0.50-0.95: moderate support
# PP < 0.50: uncertain, consider lumping species
# Always run multiple independent chains to check convergence
```

## Multi-Method Comparison

Compare delimitation results across methods to identify consensus species:

```r
library(ape)

# Partition vectors from each method (species assignment per individual)
asap_partition <- c(1, 1, 1, 2, 2, 3, 3, 3, 4, 4)
bptp_partition <- c(1, 1, 1, 2, 2, 3, 3, 4, 5, 5)
gmyc_partition <- c(1, 1, 1, 2, 2, 3, 3, 3, 4, 4)

# Compare with adjusted Rand index
# ARI = 1: perfect agreement; ARI = 0: random; ARI < 0: less than random
library(fossil)
cat('ASAP vs bPTP ARI:', adj.rand.index(asap_partition, bptp_partition), '\n')
cat('ASAP vs GMYC ARI:', adj.rand.index(asap_partition, gmyc_partition), '\n')
cat('bPTP vs GMYC ARI:', adj.rand.index(bptp_partition, gmyc_partition), '\n')

# Consensus: individuals assigned to the same species by >= 2 methods
comparison <- data.frame(
    individual = paste0('ind_', 1:10),
    ASAP = asap_partition,
    bPTP = bptp_partition,
    GMYC = gmyc_partition
)

# Count number of species per method
cat('\nSpecies counts:\n')
cat('  ASAP:', length(unique(asap_partition)), '\n')
cat('  bPTP:', length(unique(bptp_partition)), '\n')
cat('  GMYC:', length(unique(gmyc_partition)), '\n')
```

## Ultrametric Tree Preparation

Most tree-based methods require ultrametric input:

```r
library(ape)
library(phytools)

tree <- read.tree('ml_tree.nwk')

# Method 1: chronos (penalized likelihood, fastest)
# model='correlated': rates correlated across branches (relaxed clock)
ultra_tree <- chronos(tree, model = 'correlated')

# Method 2: BEAST2 (Bayesian, most rigorous but slow)
# Produces ultrametric tree with credibility intervals
# Run externally, import the MCC tree

# Verify
is.ultrametric(ultra_tree)

# Fix near-ultrametric trees with small rounding errors
ultra_tree <- force.ultrametric(ultra_tree, method = 'extend')
```

## Related Skills

- edna-metabarcoding - Generate barcode sequences from eDNA
- conservation-genetics - Population-level genetic assessment
- phylogenetics/tree-io - Tree input/output for delimitation methods
- phylogenetics/modern-tree-inference - Phylogenetic tree construction
- database-access/entrez-fetch - Retrieve barcode sequences from GenBank