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 分析流程一般包括原始测序数据FASTQ质量控制与拼接OTU/ASV 聚类物种注释生成 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 python3.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/simple2.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_okTrue) 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[fs{i1:02d} for i in range(len(groups))], columns[fg{j1} for j in range(N_TAXA)], ) df.insert(0, group, groups) df.to_csv(DATA_DIR / microbiome_data.csv, index_labelsample) 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(axis1), axis0)这一操作把每行总和变成 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 的核心步骤如下对距离矩阵取平方。利用中心化矩阵做 Gower 变换。对变换后的矩阵进行特征分解。取特征值最大的几个特征向量作为主坐标。最终每个样本会得到一个二维坐标画出来就是散点图。如果两种分组在 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_col0) group df[group] otu df.drop(columnsgroup) 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_okTrue) df pd.read_csv(ROOT / data / microbiome_data.csv, index_col0) group df[group] otu df.drop(columnsgroup) def shannon_index(counts): counts np.asarray(counts, dtypefloat) counts counts[counts 0] p counts / counts.sum() return float(-np.sum(p * np.log(p))) def simpson_index(counts): counts np.asarray(counts, dtypefloat) p counts / counts.sum() return float(1 - np.sum(p ** 2)) def chao1_index(counts): counts np.asarray(counts, dtypeint) 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, axis1), simpson: otu.apply(simpson_index, axis1), chao1: otu.apply(chao1_index, axis1), }, indexotu.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_varFalse) print(fWelch t 检验 p 值: {p_value:.4g}) fig, ax plt.subplots(figsize(6, 5)) alpha.boxplot(columnshannon, bygroup, axax) ax.set_ylabel(Shannon index) ax.set_title(Alpha diversity between groups) fig.suptitle() plt.savefig(OUTPUT_DIR / alpha_shannon_boxplot.png, dpi150, bbox_inchestight) 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_okTrue) df pd.read_csv(ROOT / data / microbiome_data.csv, index_col0) group df[group].values otu df.drop(columnsgroup).values def bray_curtis(a, b): a np.asarray(a, dtypefloat) b np.asarray(b, dtypefloat) 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, indexdf.index, columnsdf.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], labellabel, colorcolors[label], alpha0.8, s60 ) plt.axhline(0, colorgray, linewidth0.8, linestyle--) plt.axvline(0, colorgray, linewidth0.8, linestyle--) plt.xlabel(fPCo1 ({explained[0]*100:.1f}%)) plt.ylabel(fPCo2 ({explained[1]*100:.1f}%)) plt.title(Beta diversity (Bray-Curtis) PCoA) plt.legend() plt.savefig(OUTPUT_DIR / pcoa.png, dpi150, bbox_inchestight) 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_okTrue) df pd.read_csv(ROOT / data / microbiome_data.csv, index_col0) group df[group] otu df.drop(columnsgroup) rel otu.div(otu.sum(axis1), axis0) 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, bottombottom, labeltaxon, linewidth0) bottom values plt.xticks(np.arange(n), sample_labels, rotation90) plt.ylabel(Relative abundance) plt.xlabel(Sample) plt.title(Microbial community composition) plt.legend(bbox_to_anchor(1.02, 1), locupper left, fontsize8) plt.tight_layout() plt.savefig(OUTPUT_DIR / composition_stacked_bar.png, dpi150) 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_col0) group df[group] otu df.drop(columnsgroup) fiber_mask group fiber rows [] for taxon in otu.columns: fiber_vals