scikit-bio

scikit-bio

概述

scikit-bio 是一个全面的 Python 生物数据分析库，涵盖序列操作、比对、系统发育、微生物生态学和多变量统计。它提供专用数据结构（DNA、RNA、Protein、DistanceMatrix、TreeNode、TabularMSA），并可与更广泛的 Python 科学计算生态系统集成。

适用场景

• 计算微生物组数据的 α/β 多样性并运行 PERMANOVA
• 根据距离矩阵构建系统发育树（NJ、UPGMA）
• 对群落组成数据进行 PCoA 排序
• 读写生物学格式（FASTA、FASTQ、Newick、BIOM）
• 两两序列比对（Smith-Waterman、Needleman-Wunsch）
• 计算用于系统发育 β 多样性分析的 UniFrac 距离
• 对生态距离矩阵进行统计检验（ANOSIM、Mantel）
• 使用 QIIME 2 工件和微生物组流水线
• 对于高通量 NGS 比对/变异检测，请改用 STAR/BWA
• 对于蛋白质结构预测，请改用 AlphaFold/ESMFold

前置条件

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

快速入门

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']}")

核心 API

1. 序列操作

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. 序列比对

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. 系统发育树

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. 多样性分析

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. 排序分析

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. 统计检验

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. 文件 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. 距离矩阵

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}")

核心概念

数据结构

• DNA/RNA/Protein：经过验证的生物序列对象，包含元数据（序列级、逐碱基位置级、区间区域）
• TabularMSA：以表格形式表示的多序列比对；支持共识、保守性和位置熵分析
• TreeNode：支持遍历、剪枝、重新定根和距离计算的系统发育树
• DistanceMatrix：按 ID 索引的对称距离矩阵；可与排序分析和统计功能集成
• Table：用于微生物组计数数据的 BIOM 特征表（样本 × 特征）；可与 pandas、polars、AnnData 互操作

系统发育与非系统发育多样性

指标类型	非系统发育	系统发育（需要树）
α 多样性	Shannon、Simpson、observed_otus、Pielou's evenness	Faith's PD
β 多样性	Bray-Curtis、Jaccard	Unweighted UniFrac、Weighted UniFrac

计数必须为整数

多样性函数要求整数丰度计数，而非相对频率。如果你有比例值，请乘回原始计数：

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

常见工作流

工作流 1：微生物组多样性分析

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}")

工作流 2：系统发育树构建

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}")

工作流 3：序列处理流水线

1. 使用生成器读取 FASTQ 序列（通过 skbio.io.read 并设置 constructor=DNA）
2. 根据 positional_metadata['quality'] 中的质量分数进行筛选
3. 使用 find_with_regex() 和序列切片修剪接头
4. 使用 find_with_regex('ATG.{3}') 查找基序
5. 通过 dna.transcribe().translate() 进行翻译
6. 使用 skbio.io.write() 将输出写入 FASTA

关键参数

参数	函数	默认值	作用
metric	alpha/beta_diversity	必需	多样性指标名称（例如，'shannon'、'braycurtis'）
permutations	permanova/anosim/mantel	999	置换次数；次数越高，p 值越精确
tree	alpha/beta_diversity	None	系统发育指标（Faith PD、UniFrac）所必需
otu_ids	alpha/beta_diversity	None	将计数映射到树末端的特征 ID
method	mantel	'pearson'	相关方法：'pearson' 或 'spearman'
constructor	io.read	None	用于解析的序列类（DNA、RNA、Protein）
format	io.read/write	Auto	文件格式（fasta、fastq、newick 等）
genetic_code	translate	1	NCBI 遗传密码表（1=标准，11=细菌）
k	kmer_frequencies	必需	用于频率计算的 k-mer 长度

最佳实践

1. 大文件使用生成器 — skbio.io.read() 会返回生成器；处理数百万条序列时，避免使用 list()
2. 进行系统发育多样性分析时使用投影 CRS — 始终提供匹配的 tree 和 otu_ids；使用 tree.shear(feature_ids) 将树裁剪至特征集合
3. 将 PERMANOVA 与 PERMDISP 搭配使用 — PERMANOVA 对离散度差异较敏感；运行 PERMDISP 检查组间同质性
4. 使用 999 次或更多置换，以获得适合发表的 p 值；探索性分析使用 99 次即可
5. 大型特征表使用 HDF5 BIOM 格式，而非 JSON（I/O 更快，文件更小）
6. 反模式 — 在多样性函数中使用相对丰度：多样性函数要求整数计数，而非比例值。必要时转换回整数计数
7. 反模式 — k-mer 分析中的 k 值过小：k < 3 时判别能力有限；序列比较使用 k=5-7

常见操作示例

操作示例：计算多样性前进行稀释

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)

操作示例：部分计算 Beta 多样性（大型数据集）

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")

操作示例：BIOM 表操作

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)

故障排除

问题	原因	解决方案
ValueError: Ids must be unique	DistanceMatrix/序列中存在重复 ID	对 ID 去重：ids = list(dict.fromkeys(ids))
ValueError: Counts must be integers	将相对丰度传入多样性函数	使用整数计数；转换回整数：(proportions * 1000).astype(int)
大型 FASTA 文件导致内存错误	一次性加载了所有序列	使用生成器：for seq in skbio.io.read(...)
树末端与 OTU ID 不匹配	系统发育多样性计算失败	裁剪树：tree = tree.shear(feature_ids)
PERMANOVA p=0.001，但 PCoA 中各组重叠	离散度差异显著	运行 permdisp() 进行检查；PERMANOVA 检验位置和离散度
翻译结果错误	细菌使用了默认遗传密码（标准遗传密码）	对细菌/古菌序列设置 genetic_code=11
NJ 树构建缓慢	n 较大时复杂度为 O(n³)	对 >1000 个分类单元使用 GME 或 BME：from skbio.tree import gme

随附资源

• references/extended_api.md — 扩展 API 参考，涵盖高级比对参数（SSW、间隙罚分、CIGAR）、树构建算法（NJ、UPGMA、GME/BME）、BIOM 表操作、蛋白质嵌入，以及与 QIIME 2、pandas 和 scikit-learn 的集成模式

未从原文迁移的内容：原文中的单个 api_reference.md（749 行）已精简为 references/extended_api.md，重点介绍正文未涵盖的功能。省略的内容：与核心 API 重复的基础示例、冗长的故障排除内容（已由上方的故障排除表涵盖）。

参考资料

• scikit-bio 文档：https://scikit.bio/docs/latest/
• GitHub：https://github.com/scikit-bio/scikit-bio
• QIIME 2 论坛：https://forum.qiime2.org

相关技能

• scanpy-scrna-seq — 单细胞 RNA-seq 分析（数据领域不同，但共享排序概念）
• statistical-analysis — 用于补充生态统计的一般统计检验
• matplotlib-scientific-plotting — PCoA 图和多样性比较的可视化