scikit-bio

scikit-bio - 生物信息学与生态学

scikit-bio 提供严谨生物学分析所需的数据结构和统计方法。它擅长计算 α/β 多样性、进行排序分析（PCoA），以及处理复杂的系统发育树。

适用场景

• 分析微生物组数据（分类组成、群落结构）。
• 计算生态多样性指标（Shannon、Simpson、UniFrac）。
• 进行排序分析以实现可视化（PCoA、DCA）。
• 进行高层次的序列操作（带元数据的 DNA、RNA、Protein）。
• 读取和写入系统发育树（Newick 格式）。
• 分析成对及多序列比对。
• 对群落差异进行统计检验（PERMANOVA、ANOSIM）。

参考文档

官方文档：http://scikit-bio.org/
教程：http://scikit-bio.org/docs/latest/
搜索模式：skbio.sequence、skbio.stats.distance、skbio.diversity、skbio.stats.ordination

核心原则

生物序列语法

scikit-bio 不使用原始字符串，而是使用带类型的对象（DNA、RNA、Protein）。这些对象了解自己的字母表，能够处理质量分数（Phred），并支持生物学操作（转录、翻译）。

距离矩阵

微生物组研究通常归结为比较样本。DistanceMatrix 对象是核心，便于进行索引、子集选取和统计检验（例如 PERMANOVA）。

多样性指标

提供生态学中使用的数百种多样性指标的标准化实现，确保不同研究结果具有可重复性。

快速参考

安装

pip install scikit-bio

标准导入

import skbio
import numpy as np
import pandas as pd
from skbio import DNA, RNA, Protein, Sequence
from skbio.stats.distance import DistanceMatrix
from skbio.diversity import alpha, beta
from skbio.stats.ordination import pcoa

基本模式 - 序列操作

from skbio import DNA

# 1. Create a sequence with metadata
seq = DNA("ACC--GTT", metadata={'id': 'sample1', 'desc': 'gene A'})

# 2. Biological operations
rc = seq.reverse_complement()
degapped = seq.degap()

# 3. Validation
print(f"Is valid? {seq.is_valid()}") # Checks against DNA alphabet

关键规则

✅ 应该做

• 使用带类型的序列 - 始终使用 DNA、RNA 或 Protein，而非通用的 Sequence，以便使用特定字母表的方法。
• 设置元数据 - 使用 metadata 和 positional_metadata（用于质量分数）属性，使数据保持自包含。
• 检查字母表 - 从不可信来源导入数据时，使用 .is_valid() 检查是否含有非 IUPAC 字符。
• 投影距离矩阵 - 使用 pcoa() 可视化高维群落数据。
• 为 UniFrac 指定参考 - 计算 UniFrac 多样性时，确保系统发育树包含丰度表中的所有分类单元。
• 验证树 - 使用 skbio.tree.TreeNode 操作树，因为它提供了可靠的遍历方法。

❌ 不应该做

• 混用序列类型 - 不要尝试将 DNA 与 RNA 对象进行比对。
• 忽略质量分数 - 如果使用 FASTQ 数据，请将 Phred 分数存储在 positional_metadata 中。
• 手动计算距离 - 不要使用原始 NumPy 计算生物学距离；请使用 skbio.diversity.beta 获取 UniFrac 或 Bray-Curtis 距离。
• 忘记匹配 ID - 进行排序分析或 PERMANOVA 时，确保距离矩阵中的 ID 与元数据映射中的 ID 完全匹配。

反模式（绝对不要）

from skbio import DNA

# ❌ BAD: Manual reverse complement with string logic
rc_str = seq_str[::-1].replace('A', 't').replace('T', 'a')... # Fragile

# ✅ GOOD: Built-in validated method
rc_seq = DNA(seq_str).reverse_complement()

# ❌ BAD: Using raw lists for community analysis
# data = [[1, 0, 5], [2, 1, 0]]
# dist = manual_bray_curtis(data)

# ✅ GOOD: Using DistanceMatrix
from skbio.stats.distance import DistanceMatrix
dm = DistanceMatrix(matrix_data, ids=['S1', 'S2', 'S3'])

# ❌ BAD: Stripping metadata to run scikit-learn
# vals = seq.values # Metadata lost!

# ✅ GOOD: Process within skbio or use standard IO
seq.write('output.fasta')

序列分析（skbio.sequence）

高级序列操作

from skbio import DNA, Protein

# Sequence with quality scores
seq = DNA("ACGT", positional_metadata={'quality': [30, 35, 40, 20]})

# Slicing preserves metadata
sub = seq[1:3] 
print(sub.positional_metadata) # {'quality': [35, 40]}

# K-mer frequencies
freqs = seq.kmer_frequencies(k=2)

# Translation (DNA -> Protein)
protein = DNA("ATGCGA").translate()

多样性分析（skbio.diversity）

α 多样性（样本内）

from skbio.diversity import alpha

counts = [10, 0, 5, 2, 20] # Species abundances
otus = ['OTU1', 'OTU2', 'OTU3', 'OTU4', 'OTU5']

shannon = alpha.shannon(counts)
simpson = alpha.simpson(counts)
observed_otus = alpha.observed_otus(counts)

print(f"Shannon Index: {shannon:.3f}")

β 多样性（样本间）

from skbio.diversity import beta
import numpy as np

# Abundance table (samples x taxa)
data = np.array([[10, 20, 0], 
                 [5, 15, 2], 
                 [0, 1, 30]])
ids = ['Sample1', 'Sample2', 'Sample3']

# Bray-Curtis distance
bc_dm = beta.pw_distances(data, ids=ids, metric='braycurtis')

# Weighted UniFrac (Requires a tree)
# unifrac_dm = beta.weighted_unifrac(data, otu_ids, tree)

系统发育（skbio.tree）

处理树

from skbio import TreeNode
from io import StringIO

# Load Newick tree
tree_str = "((A:0.1, B:0.2)C:0.3, D:0.4)E;"
tree = TreeNode.read(StringIO(tree_str))

# Traverse
for node in tree.tips():
    print(f"Leaf: {node.name}, dist: {node.length}")

# Find Lowest Common Ancestor
lca = tree.find_lca(['A', 'B'])
print(f"LCA of A and B is {lca.name}")

# Rooting
tree.root_at('D')

统计与排序分析

PCoA（主坐标分析）

from skbio.stats.ordination import pcoa
from skbio.stats.distance import DistanceMatrix

# Create DM
dm = DistanceMatrix([[0, 0.5, 0.8], [0.5, 0, 0.2], [0.8, 0.2, 0]], 
                    ids=['A', 'B', 'C'])

# Perform PCoA
results = pcoa(dm)

# View proportion of variance explained
print(results.proportion_explained)

# Access coordinates for plotting
coords = results.samples

群落检验（PERMANOVA）

from skbio.stats.distance import permanova

# metadata linking IDs to groups
metadata = pd.DataFrame({
    'BodySite': ['Gut', 'Gut', 'Skin', 'Skin']}, 
    index=['S1', 'S2', 'S3', 'S4'])

# Assuming dm is a DistanceMatrix of S1..S4
# results = permanova(dm, metadata, column='BodySite', permutations=999)
# print(results['p-value'])

实用工作流

1. 微生物组距离可视化

def plot_microbiome_pcoa(abundance_df, metadata_df, group_col):
    """Full workflow: Abundance -> Distance -> PCoA -> Table."""
    from skbio.diversity import beta
    from skbio.stats.ordination import pcoa
    
    # 1. Calculate Bray-Curtis distance
    dm = beta.pw_distances(abundance_df.values, ids=abundance_df.index, metric='braycurtis')
    
    # 2. PCoA
    pc = pcoa(dm)
    
    # 3. Merge with metadata for plotting
    plot_data = pc.samples[['PC1', 'PC2']].join(metadata_df)
    
    return plot_data # Ready for Seaborn/Plotly

2. 序列 QC 与过滤

def filter_low_quality_sequences(sequences, min_avg_qual=30):
    """Filters DNA sequences based on Phred scores."""
    valid_seqs = []
    for s in sequences:
        # Assuming quality is in positional_metadata
        avg_q = np.mean(s.positional_metadata['quality'])
        if avg_q >= min_avg_qual:
            valid_seqs.append(s)
    return valid_seqs

3. 系统发育距离计算

def get_tip_distances(tree):
    """Returns a distance matrix of all tips in a tree."""
    dm = tree.tip_tip_distances()
    return dm

性能优化

向量化多样性计算

为多个样本计算多样性时，将整个二维数组传递给 beta.pw_distances，而不是循环遍历样本对。

树遍历

使用 tree.preorder() 或 tree.postorder()，而不是手动递归，这样在大型（尖端数超过 10,000）树上性能会好得多。

常见陷阱与解决方案

距离矩阵中的“缺失 ID”

如果元数据索引与 DistanceMatrix ID 不匹配，PERMANOVA 和 PCoA 将无法运行。

# ✅ Solution: Ensure alignment
common_ids = set(dm.ids).intersection(metadata.index)
dm_sub = dm.filter(common_ids)
metadata_sub = metadata.loc[list(common_ids)]

UniFrac 的内存使用量

UniFrac 计算可能非常耗费内存。对于超大型数据集，可以考虑使用优化版 skbio.diversity.beta.unweighted_unifrac，或使用 Unifrac-Binaries 等外部工具。

IUPAC 模糊字符

标准 DNA 对象不允许任意字符。

# ❌ Problem: DNA("ACGN") raises error if 'N' isn't handled
# ✅ Solution: scikit-bio DNA supports IUPAC (N, R, Y, etc.) by default.
# But if you have non-standard characters:
from skbio import Sequence
custom = Sequence("ACG-X") # Generic sequence permits anything

scikit-bio 是现代微生物组和进化研究的数学核心。通过强制执行严格的生物学类型约束并提供经过验证的生态学指标，它确保你的生物学洞见建立在严谨的统计基础之上。