scikit-bio

scikit-bio

Overview

scikit-bio is a comprehensive Python library for biological data analysis, spanning sequence manipulation, alignment, phylogenetics, microbial ecology, and multivariate statistics. It provides specialized data structures (DNA, RNA, Protein, DistanceMatrix, TreeNode, TabularMSA) that integrate with the broader Python scientific stack.

When to Use

• Calculating alpha/beta diversity and running PERMANOVA on microbiome data
• Building phylogenetic trees from distance matrices (NJ, UPGMA)
• Performing PCoA ordination on community composition data
• Reading/writing biological formats (FASTA, FASTQ, Newick, BIOM)
• Pairwise sequence alignment (Smith-Waterman, Needleman-Wunsch)
• Computing UniFrac distances for phylogenetic beta diversity
• Statistical testing on ecological distance matrices (ANOSIM, Mantel)
• Working with QIIME 2 artifacts and microbiome pipelines
• For high-throughput NGS alignment/variant calling, use STAR/BWA instead
• For protein structure prediction, use AlphaFold/ESMFold instead

Prerequisites

pip install scikit-bio
# Optional: pip install biom-format  — HDF5 BIOM table support
# Optional: pip install matplotlib seaborn  — visualization

Quick Start

import skbio
from skbio.diversity import alpha_diversity, beta_diversity
from skbio.stats.distance import permanova
from skbio.stats.ordination import pcoa
import numpy as np

# Sample OTU counts (samples × features)
counts = np.array([[10, 20, 30], [15, 25, 5], [5, 10, 40], [20, 5, 15]])
sample_ids = ['S1', 'S2', 'S3', 'S4']
grouping = ['control', 'control', 'treatment', 'treatment']

# Alpha diversity
shannon = alpha_diversity('shannon', counts, ids=sample_ids)
print(f"Shannon diversity: {shannon.values}")  # [1.09, 1.04, 0.94, 1.03]

# Beta diversity → PCoA → PERMANOVA
bc_dm = beta_diversity('braycurtis', counts, ids=sample_ids)
pcoa_results = pcoa(bc_dm)
results = permanova(bc_dm, grouping, permutations=999)
print(f"PERMANOVA p-value: {results['p-value']}")

Core API

1. Sequence Manipulation

from skbio import DNA, RNA, Protein

# Create and manipulate sequences
dna = DNA('ATCGATCGATCG', metadata={'id': 'gene1', 'description': 'test'})
rc = dna.reverse_complement()
rna = dna.transcribe()
protein = rna.translate()
print(f"DNA: {dna}, RC: {rc}, Protein: {protein}")

# Motif finding and k-mer analysis
motif_positions = dna.find_with_regex('ATG.{3}')
kmer_freqs = dna.kmer_frequencies(k=3)
print(f"3-mer frequencies: {dict(list(kmer_freqs.items())[:3])}")

# Sequence properties
print(f"Has degenerates: {dna.has_degenerates()}")
print(f"GC content: {dna.gc_content():.2f}")
degapped = dna.degap()  # Remove gap characters

# Metadata: sequence-level, positional, interval
dna = DNA('ATCGATCG', metadata={'id': 'seq1'},
          positional_metadata={'quality': [30, 35, 40, 38, 32, 36, 34, 33]})
dna.interval_metadata.add([(0, 4)], metadata={'type': 'promoter'})
print(f"Quality scores: {list(dna.positional_metadata['quality'])}")

2. Sequence Alignment

from skbio import DNA
from skbio.alignment import local_pairwise_align_ssw, TabularMSA

# Pairwise local alignment (Smith-Waterman via SSW)
seq1 = DNA('ACTCGATCGATCGATCGATCG')
seq2 = DNA('ATCGATCGATCGATCGATCGA')
alignment, score, start_end = local_pairwise_align_ssw(seq1, seq2)
print(f"Score: {score}, Positions: {start_end}")

# Multiple sequence alignment from file
msa = TabularMSA.read('alignment.fasta', constructor=DNA)
consensus = msa.consensus()
conservation = msa.conservation()
print(f"Consensus: {consensus[:20]}, Conservation: {conservation[:5]}")

3. Phylogenetic Trees

from skbio import TreeNode, DistanceMatrix
from skbio.tree import nj, upgma

# Build tree from distance matrix
data = [[0, 5, 9, 9], [5, 0, 10, 10], [9, 10, 0, 8], [9, 10, 8, 0]]
dm = DistanceMatrix(data, ids=['A', 'B', 'C', 'D'])
tree = nj(dm)
print(tree.ascii_art())

# Tree operations
subtree = tree.shear(['A', 'B', 'C'])  # Prune to subset
tips = [node.name for node in tree.tips()]
lca = tree.lowest_common_ancestor(['A', 'B'])
print(f"Tips: {tips}, LCA children: {len(list(lca.children))}")

# Tree comparison
tree2 = upgma(dm)
rf_dist = tree.compare_rfd(tree2)
cophenetic_dm = tree.cophenetic_matrix()
print(f"Robinson-Foulds distance: {rf_dist}")

4. Diversity Analysis

from skbio.diversity import alpha_diversity, beta_diversity
import numpy as np

counts = np.array([[10, 20, 30, 0], [15, 25, 5, 10], [5, 10, 40, 2]])
sample_ids = ['S1', 'S2', 'S3']

# Alpha diversity (multiple metrics)
for metric in ['shannon', 'simpson', 'observed_otus', 'pielou_e']:
    alpha = alpha_diversity(metric, counts, ids=sample_ids)
    print(f"{metric}: {alpha.values.round(3)}")

# Beta diversity
bc_dm = beta_diversity('braycurtis', counts, ids=sample_ids)
jaccard_dm = beta_diversity('jaccard', counts, ids=sample_ids)
print(f"Bray-Curtis S1-S2: {bc_dm['S1', 'S2']:.3f}")

# Phylogenetic diversity (requires tree + OTU IDs)
from skbio.diversity import alpha_diversity, beta_diversity

faith_pd = alpha_diversity('faith_pd', counts, ids=sample_ids,
                           tree=tree, otu_ids=feature_ids)
unifrac_dm = beta_diversity('unweighted_unifrac', counts,
                            ids=sample_ids, tree=tree, otu_ids=feature_ids)
w_unifrac_dm = beta_diversity('weighted_unifrac', counts,
                              ids=sample_ids, tree=tree, otu_ids=feature_ids)
print(f"Faith PD: {faith_pd.values}")

5. Ordination

from skbio.stats.ordination import pcoa, cca

# PCoA from distance matrix
pcoa_results = pcoa(bc_dm)
pc1 = pcoa_results.samples['PC1']
pc2 = pcoa_results.samples['PC2']
prop = pcoa_results.proportion_explained
print(f"PC1 explains {prop.iloc[0]:.1%}, PC2 explains {prop.iloc[1]:.1%}")

# CCA with environmental variables (constrained ordination)
# species_matrix: samples × species counts
# env_matrix: samples × environmental variables
cca_results = cca(species_matrix, env_matrix)
biplot_scores = cca_results.biplot_scores
print(f"Biplot scores shape: {biplot_scores.shape}")

# Save/load ordination results
pcoa_results.write('pcoa_results.txt')
loaded = skbio.OrdinationResults.read('pcoa_results.txt')

6. Statistical Testing

from skbio.stats.distance import permanova, anosim, permdisp, mantel

grouping = ['control', 'control', 'treatment']

# PERMANOVA — test group differences
perm_results = permanova(bc_dm, grouping, permutations=999)
print(f"PERMANOVA: F={perm_results['test statistic']:.3f}, "
      f"p={perm_results['p-value']:.4f}")

# ANOSIM — alternative group test
anosim_results = anosim(bc_dm, grouping, permutations=999)
print(f"ANOSIM: R={anosim_results['test statistic']:.3f}, "
      f"p={anosim_results['p-value']:.4f}")

# PERMDISP — test homogeneity of dispersions
permdisp_results = permdisp(bc_dm, grouping, permutations=999)

# Mantel test — correlation between distance matrices
r, p_value, n = mantel(bc_dm, jaccard_dm, method='pearson', permutations=999)
print(f"Mantel: r={r:.3f}, p={p_value:.4f}")

7. File I/O

import skbio

# Read various formats
seq = skbio.DNA.read('input.fasta', format='fasta')
tree = skbio.TreeNode.read('tree.nwk')
dm = skbio.DistanceMatrix.read('distances.txt')

# Memory-efficient reading (generator for large files)
for seq in skbio.io.read('large.fasta', format='fasta', constructor=skbio.DNA):
    print(f"{seq.metadata['id']}: {len(seq)} bp")

# Write and convert formats
seq.write('output.fasta', format='fasta')
seqs = list(skbio.io.read('input.fastq', format='fastq', constructor=skbio.DNA))
skbio.io.write(seqs, format='fasta', into='output.fasta')

8. Distance Matrices

from skbio import DistanceMatrix
import numpy as np

# Create from array
data = np.array([[0, 1, 2], [1, 0, 3], [2, 3, 0]])
dm = DistanceMatrix(data, ids=['A', 'B', 'C'])

# Access and slice
print(f"A-B distance: {dm['A', 'B']}")
subset = dm.filter(['A', 'C'])
condensed = dm.condensed_form()  # scipy-compatible
df = dm.to_data_frame()
print(f"Shape: {df.shape}")

Key Concepts

Data Structures

• DNA/RNA/Protein: Validated biological sequence objects with metadata (sequence-level, positional per-base, interval regions)
• TabularMSA: Multiple sequence alignment as a table; supports consensus, conservation, position entropy
• TreeNode: Phylogenetic tree with traversal, pruning, rerooting, distance calculations
• DistanceMatrix: Symmetric distance matrix with ID-based indexing; integrates with ordination and stats
• Table: BIOM feature table (samples × features) for microbiome count data; interoperates with pandas, polars, AnnData

Phylogenetic vs Non-Phylogenetic Diversity

Metric Type	Non-Phylogenetic	Phylogenetic (requires tree)
Alpha diversity	Shannon, Simpson, observed_otus, Pielou's evenness	Faith's PD
Beta diversity	Bray-Curtis, Jaccard	Unweighted UniFrac, Weighted UniFrac

Counts Must Be Integers

Diversity functions require integer abundance counts, not relative frequencies. If you have proportions, multiply back:

# Wrong: relative abundance
# counts = np.array([0.1, 0.5, 0.4])
# Correct: integer counts
counts = np.array([10, 50, 40])

Common Workflows

Workflow 1: Microbiome Diversity Analysis

import skbio
from skbio import TreeNode
from skbio.diversity import alpha_diversity, beta_diversity
from skbio.stats.ordination import pcoa
from skbio.stats.distance import permanova
import numpy as np

# 1. Load data
counts = np.loadtxt('otu_table.tsv', delimiter='\t', dtype=int, skiprows=1)
sample_ids = ['S1', 'S2', 'S3', 'S4', 'S5', 'S6']
feature_ids = ['OTU1', 'OTU2', 'OTU3', 'OTU4']
tree = TreeNode.read('phylogeny.nwk')
grouping = ['control', 'control', 'control', 'treatment', 'treatment', 'treatment']

# 2. Alpha diversity
shannon = alpha_diversity('shannon', counts, ids=sample_ids)
faith = alpha_diversity('faith_pd', counts, ids=sample_ids,
                        tree=tree, otu_ids=feature_ids)
print(f"Mean Shannon - Control: {shannon[:3].mean():.2f}, "
      f"Treatment: {shannon[3:].mean():.2f}")

# 3. Beta diversity + PCoA
unifrac_dm = beta_diversity('unweighted_unifrac', counts,
                            ids=sample_ids, tree=tree, otu_ids=feature_ids)
pcoa_results = pcoa(unifrac_dm)
print(f"PC1: {pcoa_results.proportion_explained.iloc[0]:.1%} variance")

# 4. Statistical testing
results = permanova(unifrac_dm, grouping, permutations=999)
print(f"PERMANOVA p={results['p-value']:.4f}")

Workflow 2: Phylogenetic Tree Construction

from skbio import DNA, DistanceMatrix
from skbio.tree import nj
from skbio.alignment import local_pairwise_align_ssw
import numpy as np

# 1. Read sequences
seqs = list(skbio.io.read('sequences.fasta', format='fasta', constructor=DNA))
n = len(seqs)
print(f"Loaded {n} sequences")

# 2. Compute pairwise distances (Hamming-like via alignment scores)
dist_data = np.zeros((n, n))
for i in range(n):
    for j in range(i+1, n):
        _, score, _ = local_pairwise_align_ssw(seqs[i], seqs[j])
        max_len = max(len(seqs[i]), len(seqs[j]))
        dist_data[i, j] = dist_data[j, i] = 1 - (score / max_len)

# 3. Build tree
ids = [s.metadata.get('id', f'seq{i}') for i, s in enumerate(seqs)]
dm = DistanceMatrix(dist_data, ids=ids)
tree = nj(dm)
print(tree.ascii_art())

# 4. Analyze tree
tree.write('output.nwk', format='newick')
cophenetic = tree.cophenetic_matrix()
print(f"Cophenetic matrix shape: {cophenetic.shape}")

Workflow 3: Sequence Processing Pipeline

1. Read FASTQ sequences using generator (skbio.io.read with constructor=DNA)
2. Filter by quality scores from positional_metadata['quality']
3. Trim adapters using find_with_regex() and sequence slicing
4. Find motifs with find_with_regex('ATG.{3}')
5. Translate via dna.transcribe().translate()
6. Write output FASTA with skbio.io.write()

Key Parameters

Parameter	Function	Default	Effect
metric	alpha/beta_diversity	Required	Diversity metric name (e.g., 'shannon', 'braycurtis')
permutations	permanova/anosim/mantel	999	Number of permutations; higher = more precise p-values
tree	alpha/beta_diversity	None	Required for phylogenetic metrics (Faith PD, UniFrac)
otu_ids	alpha/beta_diversity	None	Feature IDs mapping counts to tree tips
method	mantel	'pearson'	Correlation method: 'pearson' or 'spearman'
constructor	io.read	None	Sequence class for parsing (DNA, RNA, Protein)
format	io.read/write	Auto	File format (fasta, fastq, newick, etc.)
genetic_code	translate	1	NCBI genetic code (1=standard, 11=bacterial)
k	kmer_frequencies	Required	k-mer length for frequency calculation

Best Practices

1. Use generators for large files — skbio.io.read() returns a generator; avoid list() on millions of sequences
2. Use projected CRS for phylogenetic diversity — always provide matching tree and otu_ids; prune tree to feature set with tree.shear(feature_ids)
3. Pair PERMANOVA with PERMDISP — PERMANOVA is sensitive to dispersion differences; run PERMDISP to check group homogeneity
4. Use 999+ permutations for publication-quality p-values; 99 only for exploratory analysis
5. Use HDF5 BIOM format over JSON for large feature tables (faster I/O, smaller files)
6. Anti-pattern — relative abundance in diversity functions: diversity functions require integer counts, not proportions. Convert back if needed
7. Anti-pattern — small k for k-mer analysis: k < 3 provides little discriminatory power; use k=5-7 for sequence comparison

Common Recipes

Recipe: Rarefaction Before Diversity

from skbio.diversity import subsample_counts, alpha_diversity
import numpy as np

# Rarefy to minimum depth
min_depth = min(counts.sum(axis=1))
rarefied = np.array([subsample_counts(row, n=min_depth) for row in counts])
print(f"Rarefied to {min_depth} counts per sample")

# Calculate diversity on rarefied data
shannon = alpha_diversity('shannon', rarefied, ids=sample_ids)

Recipe: Partial Beta Diversity (Large Datasets)

from skbio.diversity import partial_beta_diversity
import itertools

# Only compute specific pairs (saves time on large matrices)
pairs = list(itertools.combinations(sample_ids[:10], 2))
partial_dm = partial_beta_diversity('braycurtis', counts,
                                    ids=sample_ids, id_pairs=pairs)
print(f"Computed {len(pairs)} pairwise distances")

Recipe: BIOM Table Operations

from skbio import Table
import numpy as np

# Read BIOM table
table = Table.read('feature_table.biom')
print(f"Samples: {table.shape[1]}, Features: {table.shape[0]}")

# Filter low-abundance features
filtered = table.filter(lambda row, id_, md: row.sum() > 10, axis='observation')

# Convert to pandas
df = table.to_dataframe()
# Normalize to relative abundance
rel_abundance = df.div(df.sum(axis=0), axis=1)

Troubleshooting

Problem	Cause	Solution
ValueError: Ids must be unique	Duplicate IDs in DistanceMatrix/sequences	Deduplicate IDs: ids = list(dict.fromkeys(ids))
ValueError: Counts must be integers	Relative abundance passed to diversity	Use integer counts; multiply back: (proportions * 1000).astype(int)
Memory error on large FASTA	Loading all sequences at once	Use generator: for seq in skbio.io.read(...)
Tree tip / OTU ID mismatch	Phylogenetic diversity fails	Prune tree: tree = tree.shear(feature_ids)
PERMANOVA p=0.001 but groups overlap in PCoA	Significant dispersion difference	Run permdisp() to check; PERMANOVA tests location AND dispersion
Wrong translation	Default genetic code (standard) used for bacteria	Set genetic_code=11 for bacterial/archaeal sequences
Slow NJ tree construction	O(n³) for large n	Use GME or BME for >1000 taxa: from skbio.tree import gme

Bundled Resources

• references/extended_api.md — Extended API reference covering advanced alignment parameters (SSW, gap penalties, CIGAR), tree construction algorithms (NJ vs UPGMA vs GME/BME), BIOM table manipulation, protein embeddings, and integration patterns with QIIME 2, pandas, and scikit-learn

Not migrated from original: the original's single api_reference.md (749 lines) was condensed into references/extended_api.md with focus on capabilities not covered inline. Omitted content: basic examples duplicating Core API, verbose troubleshooting (covered in Troubleshooting table above).

References

• scikit-bio documentation: https://scikit.bio/docs/latest/
• GitHub: https://github.com/scikit-bio/scikit-bio
• QIIME 2 forum: https://forum.qiime2.org

Related Skills

• scanpy-scrna-seq — single-cell RNA-seq analysis (different data domain but shares ordination concepts)
• statistical-analysis — general statistical testing to complement ecological stats
• matplotlib-scientific-plotting — visualization of PCoA plots and diversity comparisons