737. 微生物群落稳定性分析¶
一句话概述：量化微生物群落抵抗干扰和恢复原状的能力——就像测试一棵树在大风中"抗风能力"有多强，风停后"恢复能力"有多快。
核心知识点速查表¶
| 概念 | 白话解释 | 计算方式 |
|---|---|---|
| 抵抗力(Resistance) | 干扰时群落变化有多小 | 干扰前后的距离 |
| 恢复力(Resilience) | 干扰后恢复到原状的速度 | 恢复速率曲线斜率 |
| 时间稳定性 | 群落组成随时间波动的程度 | 均值/标准差 |
| 物种周转率 | 物种来来去去的速度 | Jaccard时间序列 |
| 群落抗性指数 | 综合衡量群落稳定性 | 多维指标 |
| 吸引子 | 群落倾向回归的"平衡状态" | 动力学分析 |
一、原理（白话版）¶
1.1 稳定性的两个维度¶
- 抵抗力(Resistance)：面对干扰，群落组成改变得有多少
- 高抵抗力 = 抗生素吃下去，菌群几乎没变
低抵抗力 = 吃了抗生素，菌群大洗牌
恢复力(Resilience)：干扰去除后，多快回到原来的状态
- 高恢复力 = 停药1周就恢复了
- 低恢复力 = 停药3个月还没恢复
1.2 为什么要研究群落稳定性？¶
| 应用场景 | 具体问题 |
|---|---|
| 疾病诊断 | 不稳定的菌群是否是疾病前兆？ |
| 抗生素使用 | 用药后菌群能恢复吗？需要多久？ |
| 益生菌评估 | 益生菌能持久改变菌群吗？ |
| FMT效果 | 移植的菌群能稳定定植吗？ |
二、时间序列稳定性分析¶
2.1 数据准备¶
# ===== 时间序列微生物组数据的稳定性分析 =====
import pandas as pd # 导入数据处理包
import numpy as np # 导入数值计算包
from scipy.spatial.distance import braycurtis # 导入Bray-Curtis距离
import matplotlib.pyplot as plt # 导入绑图包
import seaborn as sns # 导入增强绑图
# 读取时间序列丰度表
# 行=样本(含时间信息)，列=物种
abundance = pd.read_csv("time_series_abundance.tsv", sep="\t", index_col=0)
# 读取样本元数据（含时间点和受试者信息）
metadata = pd.read_csv("metadata.tsv", sep="\t")
# Subject Timepoint Group Treatment
# S001 0 Control None
# S001 1 Control None
# S001 7 Control None
# S002 0 Treatment Antibiotics
# S002 1 Treatment Antibiotics
# S002 7 Treatment Antibiotics
2.2 计算时间稳定性¶
# ===== 方法一：变异系数法（Coefficient of Variation）=====
def temporal_stability(abundance_df, metadata_df, subject_col="Subject"):
"""
计算每个受试者的时间稳定性
稳定性 = 均值 / 标准差（越大越稳定）
"""
results = []
for subject in metadata_df[subject_col].unique(): # 遍历每个受试者
# 获取该受试者的所有时间点数据
samples = metadata_df[metadata_df[subject_col] == subject].index
sub_data = abundance_df.loc[samples] # 该受试者的丰度数据
# 对每个物种计算CV
means = sub_data.mean() # 各物种的时间平均丰度
stds = sub_data.std() # 各物种的时间标准差
# 群落整体稳定性：所有物种丰度的均值/标准差的中位数
cv_values = stds / (means + 1e-10) # CV = std/mean（加小数避免除0）
stability = 1 / cv_values.median() # 稳定性 = 1/CV的中位数
results.append({
"Subject": subject,
"Stability": stability,
"N_timepoints": len(samples),
"Mean_richness": (sub_data > 0).sum(axis=1).mean() # 平均物种丰富度
})
return pd.DataFrame(results)
stability_df = temporal_stability(abundance, metadata)
print(stability_df.sort_values("Stability", ascending=False))
# ===== 方法二：Bray-Curtis时间序列距离 =====
def bc_temporal_variability(abundance_df, metadata_df, subject_col="Subject"):
"""
计算连续时间点之间的Bray-Curtis距离
距离越小 = 越稳定
"""
results = []
for subject in metadata_df[subject_col].unique():
sub_meta = metadata_df[metadata_df[subject_col] == subject].sort_values("Timepoint")
samples = sub_meta.index.tolist()
bc_dists = [] # 存储相邻时间点的BC距离
for i in range(len(samples) - 1):
d = braycurtis(
abundance_df.loc[samples[i]].values,
abundance_df.loc[samples[i+1]].values
)
bc_dists.append(d)
results.append({
"Subject": subject,
"Mean_BC": np.mean(bc_dists), # 平均BC距离（越小越稳定）
"Max_BC": np.max(bc_dists), # 最大BC距离（最大波动）
"Variability": np.std(bc_dists) # BC距离的波动性
})
return pd.DataFrame(results)
variability_df = bc_temporal_variability(abundance, metadata)
2.3 抵抗力与恢复力¶
# ===== 计算抵抗力（Resistance） =====
def calculate_resistance(abundance_df, metadata_df,
baseline_tp=0, perturb_tp=1):
"""
抵抗力 = 1 - BC(baseline, perturbation)
值越接近1 = 抵抗力越强
"""
results = []
for subject in metadata_df["Subject"].unique():
sub_meta = metadata_df[metadata_df["Subject"] == subject]
# 获取基线和干扰时的样本
baseline_sample = sub_meta[sub_meta["Timepoint"] == baseline_tp].index[0]
perturb_sample = sub_meta[sub_meta["Timepoint"] == perturb_tp].index[0]
# 计算BC距离
bc = braycurtis(
abundance_df.loc[baseline_sample].values,
abundance_df.loc[perturb_sample].values
)
resistance = 1 - bc # 抵抗力
results.append({
"Subject": subject,
"Group": sub_meta["Group"].iloc[0],
"Resistance": resistance,
"BC_distance": bc
})
return pd.DataFrame(results)
# ===== 计算恢复力（Resilience） =====
def calculate_resilience(abundance_df, metadata_df,
baseline_tp=0, recovery_tps=[7, 14, 28]):
"""
恢复力 = 恢复曲线的斜率（BC距离随时间减小的速度）
"""
results = []
for subject in metadata_df["Subject"].unique():
sub_meta = metadata_df[metadata_df["Subject"] == subject]
baseline_sample = sub_meta[sub_meta["Timepoint"] == baseline_tp].index[0]
bc_over_time = []
for tp in recovery_tps:
tp_samples = sub_meta[sub_meta["Timepoint"] == tp]
if len(tp_samples) > 0:
recovery_sample = tp_samples.index[0]
bc = braycurtis(
abundance_df.loc[baseline_sample].values,
abundance_df.loc[recovery_sample].values
)
bc_over_time.append({"Timepoint": tp, "BC": bc})
if len(bc_over_time) >= 2:
bc_df = pd.DataFrame(bc_over_time)
# 用线性回归拟合恢复曲线
slope = np.polyfit(bc_df["Timepoint"], bc_df["BC"], 1)[0]
# 负斜率 = 恢复中（BC在减小）
resilience = -slope # 恢复力（正值=在恢复，负值=还在恶化）
results.append({
"Subject": subject,
"Group": sub_meta["Group"].iloc[0],
"Resilience": resilience,
"Final_BC": bc_df["BC"].iloc[-1]
})
return pd.DataFrame(results)
# 计算
resistance_df = calculate_resistance(abundance, metadata)
resilience_df = calculate_resilience(abundance, metadata)
# ===== 可视化 =====
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左图：抵抗力比较
sns.boxplot(data=resistance_df, x="Group", y="Resistance", ax=axes[0],
palette=["#2ecc71", "#e74c3c"])
axes[0].set_title("Resistance to Perturbation")
axes[0].set_ylabel("Resistance (1 - BC distance)")
# 右图：恢复力比较
sns.boxplot(data=resilience_df, x="Group", y="Resilience", ax=axes[1],
palette=["#2ecc71", "#e74c3c"])
axes[1].set_title("Resilience After Perturbation")
axes[1].set_ylabel("Resilience (recovery rate)")
plt.tight_layout()
plt.savefig("stability_analysis.png", dpi=300)
plt.show()
三、R语言实现¶
# ===== R语言：群落稳定性分析 =====
library(vegan) # 群落生态分析
library(ggplot2) # 绑图
# 读取数据
abundance <- read.table("abundance.tsv", header=TRUE,
row.names=1, sep="\t") # 丰度表
metadata <- read.table("metadata.tsv", header=TRUE, sep="\t") # 元数据
# 计算连续时间点间的BC距离
bc_dist <- vegdist(abundance, method="bray") # 计算BC距离矩阵
# 提取同一个体相邻时间点的距离
# （需要根据metadata中的Subject和Timepoint信息）
# NMDS可视化时间轨迹
nmds <- metaMDS(abundance, distance="bray", k=2) # NMDS降维
# 提取坐标
coords <- as.data.frame(scores(nmds, display="sites"))
coords$Subject <- metadata$Subject # 添加受试者信息
coords$Timepoint <- metadata$Timepoint # 添加时间点信息
coords$Group <- metadata$Group # 添加分组信息
# 绘制时间轨迹
ggplot(coords, aes(x=NMDS1, y=NMDS2, color=Group, group=Subject)) +
geom_path(alpha=0.5, arrow=arrow(length=unit(0.15, "cm"))) + # 轨迹线+箭头
geom_point(aes(size=Timepoint)) + # 点大小表示时间
theme_minimal() +
labs(title="Community Trajectory Over Time",
subtitle="Arrows show temporal direction") +
scale_color_manual(values=c("Control"="#2ecc71", "Treatment"="#e74c3c"))
ggsave("temporal_trajectory.png", width=10, height=8, dpi=300)
四、常见报错与解决¶
| 报错信息 | 原因 | 解决方案 |
|---|---|---|
Too few timepoints | 时间点<3个无法计算恢复力 | 设计实验时至少采集4个时间点 |
BC distance = NaN | 某个样本全为0 | 过滤掉reads数太少的样本 |
Negative resilience | 群落还在恶化/没开始恢复 | 延长观察时间或增加时间点 |
vegdist: negative values | 丰度表有负值 | 检查归一化步骤 |
NMDS: stress > 0.2 | 二维投影效果差 | 增加维度(k=3)或增加迭代 |
五、面试高频问题¶
Q1: 稳定的菌群就是健康的菌群吗？¶
A: 不一定。过于稳定（无法响应饮食变化）可能反映功能受限。健康菌群的特点是"适度稳定"——能抵抗病理性干扰，但能响应生理性变化。
Q2: 如何设计实验来研究菌群稳定性？¶
A: ①纵向采样（多个时间点）是必须的；②至少包括：基线期（干扰前多个点）、干扰期、恢复期；③建议每组至少10个受试者，每人≥5个时间点。
Q3: 物种多样性和群落稳定性的关系？¶
A: 经典的"多样性-稳定性"假说认为物种越多群落越稳定（保险效应）。但在微生物组中这个关系更复杂——关键是功能冗余而非单纯的物种数量。
六、速查表¶
# ===== 稳定性分析速查 =====
# 核心指标：
# 抵抗力 = 1 - BC(baseline, perturbation)
# 恢复力 = -slope(BC vs time)
# 时间稳定性 = mean(abundance) / std(abundance)
# 物种周转率 = 1 - Jaccard(t1, t2)
# Python关键包：
pip install pandas numpy scipy scikit-bio matplotlib seaborn
# R关键包：
install.packages(c("vegan", "ggplot2"))
# 实验设计最低要求：
# 时间点数 ≥ 5
# 受试者数 ≥ 10/组
# 包含：基线期 + 干扰期 + 恢复期