Python实战：从OTU表到肠道微生物组Alpha/Beta多样性分析

很多人看到 “hack your gut microbiome” 这个标题，第一反应可能是“肠道菌群还能像写代码一样被修改？”其实，从开发者的角度看，肠道微生物组更像一个高维计数矩阵：每个样本对应若干细菌分类群的测序读数，我们要做的就是数据清洗、标准化、降维、统计检验和可视化。这篇文章就用 Python 完整演示一套从 OTU 表到 Alpha/Beta 多样性分析的流程，适合想入门生物信息学的开发者和数据分析师，也适合给已有测序数据但不知道怎么下手的同学做参考。

1. 把肠道微生物组当成一个数据分析问题

1.1 肠道微生物组是什么

肠道微生物组是指生活在人体肠道内的细菌、真菌、古菌和病毒等微生物的总称。这些微生物参与食物消化、维生素合成、免疫调节，甚至通过肠-脑轴影响神经系统状态。现代微生物组研究最常用的手段是 16S rRNA 基因测序：细菌的 16S rRNA 基因既有高度保守的区域，也有可变区域，通过扩增并测序可变区，可以大致区分不同细菌种类。

测序结论并不是“你的肠道里有 3 亿个双歧杆菌”这种直观结果，而是生成一张矩阵：每一行是一个样本，每一列是一个细菌分类群，交叉点是该分类群在这份样本中检测到的序列数目。下游数据分析的大部分工作，都是围绕这张“OTU 表”或“ASV 表”展开。

1.2 “hack” 在这里指什么

这里的 hack 不是指攻击网站或绕过安全机制，而是指一种工程化、数据驱动的改造与理解方式。传统微生物组研究依赖湿实验，周期长、成本高；而计算分析可以快速完成群落结构比较、多样性评估、差异物种筛选等任务。只要你掌握了 Python 数据处理的基本功，就能完成入门级的微生物组分析。

对开发者来说，这还是一个很友好的领域：

• 输入数据是标准的表格，处理起来和电商订单表没什么本质区别。
• 分析思路可以拆成函数和脚本，适合工程化管理。
• 绝大多数工具是开源的，公共数据集也很多。

所以，“能不能 hack 你的肠道微生物组”这个问题的答案是：能，但这里说的 hack 其实是科学的统计分析与可视化。

1.3 16S 测序产出的数据结构

完整的 16S 分析流程一般包括：

1. 原始测序数据（FASTQ）
2. 质量控制与拼接
3. OTU/ASV 聚类
4. 物种注释
5. 生成 OTU 丰度表

OTU 全称是 operational taxonomic unit，操作分类单元，通常按 97% 序列相似度聚类。ASV 则更精确，是按序列变异划分的特征。从数据结构上讲，OTU 表和 ASV 表完全一致，都可以看作样本×特征的计数矩阵。

本教程为了聚焦数据分析部分，不讲解上游序列比对，而是直接从 OTU 表开始。我们会用一份模拟数据来演示完整流程，但你完全可以把自己的测序数据整理成同样格式，套用下面的脚本。

2. 环境准备：搭建可复现的分析环境

2.1 安装 Python 与依赖

建议使用 Python 3.9 或更高版本。本教程核心依赖只有四个库：

• pandas：表格处理
• numpy：数值计算
• scipy：统计检验
• matplotlib：绘图

如果你还没有安装，可以用 pip 一次性安装完成。

pip install pandas numpy scipy matplotlib

如果你更习惯 conda 管理环境，也可以创建一个独立环境，避免和项目环境冲突。

conda create -n microbiome python=3.10 -y conda activate microbiome pip install pandas numpy scipy matplotlib

国内网络环境下，pip 下载慢时可以配置镜像源，例如清华 PyPI 镜像：

pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple

2.2 项目目录结构

建议把脚本和数据分开存放，这样后续扩展更容易。以下是一个推荐结构：

microbiome-hack/ ├── data/ │ └── microbiome_data.csv ├── output/ ├── scripts/ │ ├── generate_data.py │ ├── alpha_diversity.py │ ├── beta_diversity.py │ ├── composition.py │ └── diff_abundance.py └── README.md

所有脚本中都会使用相对路径读取data/下的文件，并把结果写入output/。

2.3 生成模拟 OTU 数据

在开始正式分析前，先生成一份模拟数据。下面脚本会生成 30 个样本、15 个细菌分类群，高纤维组和高脂组各 15 个样本，并把两组样本的群落结构设置得差异明显，方便后续观察分析方法的效果。

# scripts/generate_data.py """ 生成一份模拟的肠道微生物组 OTU 计数表。 每组 15 个样本，共 30 个样本，15 个细菌分类群（用 g1-g15 表示）。 高纤维组与高脂组的群落结构被设置为明显不同，方便后续分析看出差异。 """ from pathlib import Path import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parents[1] DATA_DIR = ROOT / 'data' DATA_DIR.mkdir(exist_ok=True) np.random.seed(42) N_SAMPLES_PER_GROUP = 15 N_TAXA = 15 groups = ['fiber'] * N_SAMPLES_PER_GROUP + ['fat'] * N_SAMPLES_PER_GROUP # 基础平均丰度：近似均匀分布，总深度约 3000 条序列 base = np.random.dirichlet(np.ones(N_TAXA)) * 3000 otu = np.zeros((len(groups), N_TAXA)) for i, g in enumerate(groups): effects = np.ones(N_TAXA) * 0.5 if g == 'fiber': effects[0:4] = 3.2 # 模拟短链脂肪酸产生菌占优 effects[8:10] = 0.01 # 模拟低丰度分类群 else: effects[10:13] = 4.0 # 模拟高脂饮食相关分类群占优 effects[0:2] = 0.01 mu = base * effects otu[i, :] = np.random.poisson(mu) df = pd.DataFrame( otu, index=[f's{i+1:02d}' for i in range(len(groups))], columns=[f'g{j+1}' for j in range(N_TAXA)], ) df.insert(0, 'group', groups) df.to_csv(DATA_DIR / 'microbiome_data.csv', index_label='sample') print(f'模拟数据已写入 {DATA_DIR / "microbiome_data.csv"}') print(df.shape)

运行脚本：

cd microbiome-hack python scripts/generate_data.py

生成的microbiome_data.csv第一列是分组信息，后续列是各分类群的原始计数。后续所有分析脚本都基于它进行。

3. 分析基础：从 OTU 表到多样性指数

在写正式分析脚本之前，需要先理解几个核心概念。这些概念会直接影响代码的实现方式。

3.1 OTU 计数矩阵与数据标准化

OTU 表里存的是测序得到的序列条数，数值受两个因素影响：一是样本中真实菌群丰度，二是测序深度。两个样本即使菌群结构完全相同，只要测序深度不同，计数就可能差很多。因此，做多样性分析前通常需要标准化。

最朴素的方法是相对丰度标准化，即每个样本内各分类群计数除以该样本总计数：

relative_abundance = otu.div(otu.sum(axis=1), axis=0)

这一操作把每行总和变成 1，方便比较组成比例。另一种传统做法是抽平，把每个样本随机抽到相同总深度，但在现代流程中，相对丰度结合合适的统计模型已经足够入门使用。

3.2 Alpha 多样性

Alpha 多样性描述的是一个样本内部物种的丰富度和均匀度，常见指标有 Shannon、Simpson、Chao1。

Shannon 指数公式：

H = -Σ p_i * ln(p_i)

其中 p_i 是分类群 i 的相对丰度。Shannon 指数越高，说明群落越多样。

Simpson 指数公式：

D = 1 - Σ p_i^2

它反映随机抽取两个个体属于不同物种的概率，越接近 1 说明多样性越高。

Chao1 估计的是物种总数：

Chao1 = S_obs + n1^2 / (2 * n2)

其中 S_obs 是观测到的分类群数量，n1 是出现次数为 1 的分类群数，n2 是出现次数为 2 的分类群数。Chao1 本质上是把隐藏在样本里却可能没被检测到的物种估算进来。

3.3 Beta 多样性

Beta 多样性比较的是不同样本之间的群落差异。最常用的距离是 Bray-Curtis 距离：

BC = Σ |a_i - b_i| / Σ (a_i + b_i)

其中 a_i 和 b_i 分别是两个样本中分类群 i 的丰度。Bray-Curtis 距离取值在 0 到 1 之间，0 表示完全相同，1 表示完全不相同。

3.4 PCoA 降维原理

得到样本间的距离矩阵后，每个样本相当于处在高维空间中的一个点。为了可视化，需要使用 PCoA 把坐标降维到二维平面。

PCoA 的核心步骤如下：

1. 对距离矩阵取平方。
2. 利用中心化矩阵做 Gower 变换。
3. 对变换后的矩阵进行特征分解。
4. 取特征值最大的几个特征向量作为主坐标。

最终每个样本会得到一个二维坐标，画出来就是散点图。如果两种分组在 PCoA 图中明显分开，说明两组微生物群落结构存在差异。

4. 完整实操：用 Python 分析微生物组数据

下面我们编写完整分析脚本。每个脚本都可以独立运行，运行前请确保当前工作目录在项目根目录下。

4.1 数据加载与预处理

先把数据读进来，确认基本结构。

# scripts/load_data.py import pandas as pd from pathlib import Path ROOT = Path(__file__).resolve().parents[1] df = pd.read_csv(ROOT / 'data' / 'microbiome_data.csv', index_col=0) group = df['group'] otu = df.drop(columns='group') print('样本数：', otu.shape[0]) print('分类群数：', otu.shape[1]) print('分组分布：') print(group.value_counts()) print(otu.head())

预期输出中可以看到 30 个样本、15 个分类群，fiber 和 fat 组各 15 个样本。

4.2 Alpha 多样性组间比较

编写alpha_diversity.py，计算三个 Alpha 多样性指标，并做 Welch t 检验比较两组的 Shannon 指数是否存在显著差异。

# scripts/alpha_diversity.py import numpy as np import pandas as pd from scipy import stats import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from pathlib import Path ROOT = Path(__file__).resolve().parents[1] OUTPUT_DIR = ROOT / 'output' OUTPUT_DIR.mkdir(exist_ok=True) df = pd.read_csv(ROOT / 'data' / 'microbiome_data.csv', index_col=0) group = df['group'] otu = df.drop(columns='group') def shannon_index(counts): counts = np.asarray(counts, dtype=float) counts = counts[counts > 0] p = counts / counts.sum() return float(-np.sum(p * np.log(p))) def simpson_index(counts): counts = np.asarray(counts, dtype=float) p = counts / counts.sum() return float(1 - np.sum(p ** 2)) def chao1_index(counts): counts = np.asarray(counts, dtype=int) observed = int(np.sum(counts > 0)) f1 = int(np.sum(counts == 1)) f2 = int(np.sum(counts == 2)) if f1 == 0: return float(observed) if f2 == 0: f2 = 1 return observed + (f1 * (f1 - 1)) / (2.0 * (f2 + 1)) alpha = pd.DataFrame({ 'shannon': otu.apply(shannon_index, axis=1), 'simpson': otu.apply(simpson_index, axis=1), 'chao1': otu.apply(chao1_index, axis=1), }, index=otu.index) alpha['group'] = group.values print('Alpha 多样性前 5 行：') print(alpha.head()) fiber = alpha[alpha['group'] == 'fiber']['shannon'] fat = alpha[alpha['group'] == 'fat']['shannon'] t_stat, p_value = stats.ttest_ind(fiber, fat, equal_var=False) print(f'Welch t 检验 p 值: {p_value:.4g}') fig, ax = plt.subplots(figsize=(6, 5)) alpha.boxplot(column='shannon', by='group', ax=ax) ax.set_ylabel('Shannon index') ax.set_title('Alpha diversity between groups') fig.suptitle('') plt.savefig(OUTPUT_DIR / 'alpha_shannon_boxplot.png', dpi=150, bbox_inches='tight') print('箱线图已保存到 output/alpha_shannon_boxplot.png')

运行结果中，如果 p 值小于 0.05，说明两组 Alpha 多样性存在统计学差异。这里需要记住，p 值只说明差异是否显著，并不说明差异有多大，所以还要结合箱线图看实际分布。

4.3 Beta 多样性与 PCoA 可视化

这部分会计算 Bray-Curtis 距离矩阵，然后进行 PCoA 降维。

# scripts/beta_diversity.py import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from pathlib import Path ROOT = Path(__file__).resolve().parents[1] OUTPUT_DIR = ROOT / 'output' OUTPUT_DIR.mkdir(exist_ok=True) df = pd.read_csv(ROOT / 'data' / 'microbiome_data.csv', index_col=0) group = df['group'].values otu = df.drop(columns='group').values def bray_curtis(a, b): a = np.asarray(a, dtype=float) b = np.asarray(b, dtype=float) num = np.abs(a - b).sum() den = a.sum() + b.sum() return num / den if den > 0 else 0.0 n = otu.shape[0] dm = np.zeros((n, n)) for i in range(n): for j in range(i + 1, n): d = bray_curtis(otu[i], otu[j]) dm[i, j] = d dm[j, i] = d print('Bray-Curtis 距离矩阵前 5 行：') print(pd.DataFrame(dm, index=df.index, columns=df.index).iloc[:5, :5]) def pcoa(distance_matrix): n = distance_matrix.shape[0] A = distance_matrix ** 2 J = np.eye(n) - np.ones((n, n)) / n B = -0.5 * J @ A @ J eigvals, eigvecs = np.linalg.eigh(B) idx = np.argsort(eigvals)[::-1] eigvals = eigvals[idx] eigvecs = eigvecs[:, idx] valid = eigvals > 1e-10 eigvals = eigvals[valid] eigvecs = eigvecs[:, valid] coords = eigvecs * np.sqrt(eigvals)[None, :] explained = eigvals / eigvals.sum() return coords, explained coords, explained = pcoa(dm) print(f'前两个主坐标解释方差比例: {explained[:2] * 100:.2f}%') colors = {'fiber': '#2b8cbe', 'fat': '#de2d26'} plt.figure(figsize=(7, 6)) for label in np.unique(group): mask = group == label plt.scatter( coords[mask, 0], coords[mask, 1], label=label, color=colors[label], alpha=0.8, s=60 ) plt.axhline(0, color='gray', linewidth=0.8, linestyle='--') plt.axvline(0, color='gray', linewidth=0.8, linestyle='--') plt.xlabel(f'PCo1 ({explained[0]*100:.1f}%)') plt.ylabel(f'PCo2 ({explained[1]*100:.1f}%)') plt.title('Beta diversity (Bray-Curtis) PCoA') plt.legend() plt.savefig(OUTPUT_DIR / 'pcoa.png', dpi=150, bbox_inches='tight') print('PCoA 图已保存到 output/pcoa.png')

如果两组在图中明显分离，说明肠道微生物群落组成存在明显差异。

4.4 群落组成可视化

堆叠柱状图是最直观展示样本中分类群相对丰度的方式。

# scripts/composition.py import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from pathlib import Path ROOT = Path(__file__).resolve().parents[1] OUTPUT_DIR = ROOT / 'output' OUTPUT_DIR.mkdir(exist_ok=True) df = pd.read_csv(ROOT / 'data' / 'microbiome_data.csv', index_col=0) group = df['group'] otu = df.drop(columns='group') rel = otu.div(otu.sum(axis=1), axis=0) sample_labels = rel.index.tolist() n = len(rel) bottom = np.zeros(n) plt.figure(figsize=(12, 5)) for taxon in rel.columns: values = rel[taxon].values plt.bar(np.arange(n), values, bottom=bottom, label=taxon, linewidth=0) bottom += values plt.xticks(np.arange(n), sample_labels, rotation=90) plt.ylabel('Relative abundance') plt.xlabel('Sample') plt.title('Microbial community composition') plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left', fontsize=8) plt.tight_layout() plt.savefig(OUTPUT_DIR / 'composition_stacked_bar.png', dpi=150) print('堆叠柱状图已保存到 output/composition_stacked_bar.png')

从图中能直观看到，fiber 组和 fat 组的优势分类群明显不同，这正是生成数据时预设的效果。

4.5 差异分类群初筛

最后，用 Mann-Whitney U 检验对每个分类群做组间比较，并做简单的 Bonferroni 校正。

# scripts/diff_abundance.py import numpy as np import pandas as pd from scipy.stats import mannwhitneyu from pathlib import Path ROOT = Path(__file__).resolve().parents[1] df = pd.read_csv(ROOT / 'data' / 'microbiome_data.csv', index_col=0) group = df['group'] otu = df.drop(columns='group') fiber_mask = group == 'fiber' rows = [] for taxon in otu.columns: fiber_vals