微生物组随机森林分类¶
一句话概述¶
利用随机森林算法对微生物组数据进行疾病分类预测，结合SHAP值进行模型解释找出关键微生物标志物，是微生物组生物标志物发现和疾病诊断模型构建的标准机器学习流程。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| 数据预处理 | CLR转化、过滤、标准化 | phyloseq, scikit-learn |
| 特征选择 | 方差过滤/Boruta/递归特征消除 | Boruta, sklearn.RFE |
| 随机森林 | 集成学习、超参数调优 | sklearn, ranger(R) |
| 交叉验证 | Nested CV、LOOCV、分层K-fold | sklearn.model_selection |
| 模型评估 | AUC-ROC、准确率、混淆矩阵 | sklearn.metrics |
| 模型解释 | SHAP值、特征重要性 | shap, TreeSHAP |
| 标志物验证 | 独立验证集、队列验证 | - |
| 多组学整合 | 微生物+代谢物联合预测 | mixOmics |
各步骤详解¶
第一步：理解随机森林在微生物组中的优势¶
白话解释： 随机森林是多棵决策树的"投票系统"。每棵树用部分样本和部分特征来训练，最终由所有树投票决定分类结果。它特别适合微生物组数据，因为：1）能处理高维数据（物种数>>样本数）；2）不怕特征间相关性；3）自带特征重要性评分；4）对异常值不敏感。
技术细节： - Bagging: 有放回抽样（~63%样本被选中） - Random subspace: 每个节点只考虑√p个候选特征 - OOB (Out-Of-Bag): 未被选中的样本用于内部验证 - 特征重要性: Mean Decrease Accuracy/Gini - 超参数: n_estimators, max_depth, min_samples_leaf, max_features
代码示例：
# 为什么随机森林适合微生物组？
"""
微生物组数据特点：
1. 高维低样本 (p >> n): 数百个OTU/ASV vs 几十-几百个样本
2. 零膨胀: 大量零值（物种缺失）
3. 组成性: 相对丰度和为1
4. 高共线性: 生态相关物种间强相关
5. 非线性关系: 物种与疾病间非线性关联
随机森林优势:
- 不需要假设数据分布
- 自动处理高维和共线性
- 捕获非线性和交互效应
- OOB估计避免过拟合
- 特征重要性直接指示标志物
"""
第二步：数据预处理¶
白话解释： 原始微生物组数据需要先做预处理：去除低丰度特征、处理零值、进行适当转化。CLR（中心对数比）转化是处理组成数据的标准方法。
代码示例：
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from scipy.stats import gmean
import warnings
warnings.filterwarnings('ignore')
# 读取OTU/ASV表和元数据
otu_table = pd.read_csv("otu_table.csv", index_col=0) # samples x OTUs
metadata = pd.read_csv("metadata.csv", index_col=0)
# 1. 过滤低丰度特征
# 保留在至少10%样本中出现且总丰度>0.01%的OTU
prevalence = (otu_table > 0).mean(axis=0)
abundance = otu_table.sum(axis=0) / otu_table.sum().sum()
keep = (prevalence >= 0.1) & (abundance >= 0.0001)
otu_filtered = otu_table.loc[:, keep]
print(f"过滤后保留 {otu_filtered.shape[1]} 个特征 (原始 {otu_table.shape[1]})")
# 2. CLR转化（处理组成数据）
def clr_transform(df):
"""Center Log-Ratio transformation"""
# 添加伪计数（处理零值）
df_pseudo = df + 0.5
# 计算几何均值
gm = gmean(df_pseudo, axis=1)
# CLR = log(x / geometric_mean)
clr = np.log(df_pseudo.divide(gm, axis=0))
return clr
# 演示用法；实际项目中 CLR 的几何均值应在每折训练集内部计算，避免信息泄露
X_clr = clr_transform(otu_filtered)
# 3. 准备标签
y = metadata.loc[X_clr.index, "disease_status"] # 'healthy' vs 'disease'
y_binary = (y == "disease").astype(int)
print(f"样本数: {len(y_binary)}")
print(f"类别分布:\n{y.value_counts()}")
第三步：特征选择¶
白话解释： 虽然随机森林能处理高维数据，但去除噪声特征能提高模型性能和可解释性。常用方法：Boruta（基于随机森林的全相关特征选择）、递归特征消除（RFE）、或基于方差的简单过滤。
代码示例：
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFECV
from boruta import BorutaPy
# 警告：BorutaPy 已停止维护，与 NumPy ≥1.24 不兼容（AttributeError: module 'numpy' has no attribute 'int'）
# 替代方案：pip install eboruta（eBoruta）或使用 sklearn 的 SequentialFeatureSelector
# 方法1: Boruta特征选择
rf_boruta = RandomForestClassifier(
n_estimators=500,
max_depth=7,
n_jobs=-1,
random_state=42)
boruta_selector = BorutaPy(
rf_boruta,
n_estimators='auto',
max_iter=100,
random_state=42,
verbose=2)
boruta_selector.fit(X_clr.values, y_binary.values)
# 查看选择的特征
selected_features = X_clr.columns[boruta_selector.support_].tolist()
tentative_features = X_clr.columns[boruta_selector.support_weak_].tolist()
print(f"确认重要的特征: {len(selected_features)}")
print(f"待定特征: {len(tentative_features)}")
# 方法2: RFECV (带交叉验证的递归特征消除)
from sklearn.model_selection import StratifiedKFold
rf_rfe = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
rfecv = RFECV(
estimator=rf_rfe,
step=5, # 每步删除5个特征
cv=cv,
scoring='roc_auc',
n_jobs=-1,
verbose=1)
rfecv.fit(X_clr.values, y_binary.values)
print(f"最优特征数: {rfecv.n_features_}")
# 使用选择的特征
X_selected = X_clr[selected_features]
第四步：模型训练与交叉验证¶
白话解释： 使用嵌套交叉验证(Nested CV)同时进行超参数调优和性能评估。外层循环评估模型泛化性能，内层循环选择最优超参数。这避免了"信息泄露"导致的过于乐观估计。
代码示例：
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import (
StratifiedKFold, GridSearchCV, cross_val_predict,
RepeatedStratifiedKFold
)
from sklearn.metrics import (
roc_auc_score, accuracy_score, classification_report,
roc_curve, confusion_matrix
)
import matplotlib.pyplot as plt
# 超参数搜索空间
param_grid = {
'n_estimators': [100, 300, 500, 1000],
'max_depth': [3, 5, 7, 10, None],
'min_samples_leaf': [1, 3, 5, 10],
'max_features': ['sqrt', 'log2', 0.3]
}
# 嵌套交叉验证
outer_cv = RepeatedStratifiedKFold(n_splits=5, n_repeats=10, random_state=42)
inner_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# 存储每折结果
auc_scores = []
all_y_true = []
all_y_prob = []
for fold, (train_idx, test_idx) in enumerate(outer_cv.split(X_selected, y_binary)):
X_train, X_test = X_selected.iloc[train_idx], X_selected.iloc[test_idx]
y_train, y_test = y_binary.iloc[train_idx], y_binary.iloc[test_idx]
# 内层CV调参
# 提示：全网格搜索在嵌套 CV 中极耗时（可达数万次拟合），实际项目建议改用 RandomizedSearchCV(n_iter=20)
grid_search = GridSearchCV(
RandomForestClassifier(random_state=42, n_jobs=-1),
param_grid,
cv=inner_cv,
scoring='roc_auc',
n_jobs=-1)
grid_search.fit(X_train, y_train)
# 在测试集评估
y_prob = grid_search.predict_proba(X_test)[:, 1]
auc = roc_auc_score(y_test, y_prob)
auc_scores.append(auc)
all_y_true.extend(y_test)
all_y_prob.extend(y_prob)
print(f"Nested CV AUC: {np.mean(auc_scores):.3f} ± {np.std(auc_scores):.3f}")
# 绘制ROC曲线
fpr, tpr, _ = roc_curve(all_y_true, all_y_prob)
plt.figure(figsize=(8, 6))
plt.plot(fpr, tpr, lw=2, label=f'AUC = {roc_auc_score(all_y_true, all_y_prob):.3f}')
plt.plot([0,1], [0,1], 'k--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve - Microbiome Random Forest Classifier')
plt.legend()
plt.savefig("roc_curve.pdf")
plt.close()
# 最终模型（使用全部数据训练）
final_rf = RandomForestClassifier(
**grid_search.best_params_,
random_state=42, n_jobs=-1)
final_rf.fit(X_selected, y_binary)
第五步：SHAP模型解释¶
白话解释： SHAP (SHapley Additive exPlanations) 基于博弈论的Shapley值，能精确量化每个特征对每个预测的贡献。比传统的feature importance更准确，还能看到每个特征的方向性效应（增加还是降低风险）。
代码示例：
import shap
# 计算SHAP值
explainer = shap.TreeExplainer(final_rf)
shap_values = explainer.shap_values(X_selected)
# 对于二分类，shap_values[1]是阳性类别的SHAP值
shap_values_disease = shap_values[1] # shap ≥ 0.40：shap_values 返回列表 [class_0_array, class_1_array]，[1] 取正类
# 1. Summary plot (全局特征重要性 + 方向)
plt.figure(figsize=(10, 8))
shap.summary_plot(shap_values_disease, X_selected,
max_display=20, show=False)
plt.tight_layout()
plt.savefig("shap_summary.pdf")
plt.close()
# 2. Bar plot (平均绝对SHAP值排序)
plt.figure(figsize=(10, 8))
shap.summary_plot(shap_values_disease, X_selected,
plot_type="bar", max_display=20, show=False)
plt.tight_layout()
plt.savefig("shap_importance.pdf")
plt.close()
# 3. 单样本解释（Force plot）
# 解释某一个预测为什么得到那个结果
idx = 0 # 第一个样本
# matplotlib=True 输出静态图（脚本环境必须加）；Jupyter 中可不加（默认交互式 JS 图）
shap.force_plot(explainer.expected_value[1],
shap_values_disease[idx],
X_selected.iloc[idx],
matplotlib=True)
plt.savefig("shap_force_single.pdf")
# 4. Dependence plot (单特征效应)
top_feature = X_selected.columns[
np.abs(shap_values_disease).mean(0).argmax()]
shap.dependence_plot(top_feature, shap_values_disease, X_selected)
plt.savefig("shap_dependence.pdf")
# 5. 提取Top标志物
mean_shap = np.abs(shap_values_disease).mean(0)
feature_importance = pd.DataFrame({
'feature': X_selected.columns,
'mean_abs_shap': mean_shap
}).sort_values('mean_abs_shap', ascending=False)
print("Top 10 微生物标志物:")
print(feature_importance.head(10))
feature_importance.to_csv("biomarker_ranking.csv", index=False)
第六步：独立验证与报告¶
白话解释： 模型在训练队列表现好不代表真的有效，必须在独立队列中验证。理想情况是使用不同地区、不同时间收集的样本。还需要与临床指标比较，评估其实用价值。
代码示例：
# 独立验证集评估
X_val = pd.read_csv("validation_cohort_otu.csv", index_col=0)
y_val = pd.read_csv("validation_cohort_labels.csv", index_col=0)["status"]
# 确保使用相同的特征和预处理
X_val_selected = clr_transform(X_val[selected_features])
# 预测
y_val_prob = final_rf.predict_proba(X_val_selected)[:, 1]
y_val_pred = final_rf.predict(X_val_selected)
# 评估
val_auc = roc_auc_score(y_val, y_val_prob)
print(f"验证集 AUC: {val_auc:.3f}")
print(classification_report(y_val, y_val_pred))
# 与临床指标比较
from sklearn.metrics import average_precision_score
# 假设有临床标志物（如CRP）
clin_auc = roc_auc_score(y_val, clinical_marker_values)
print(f"临床标志物 AUC: {clin_auc:.3f}")
print(f"微生物模型 AUC: {val_auc:.3f}")
# DeLong test比较两个AUC
# 使用scipy或专门的包
实战命令¶
#!/usr/bin/env python3
"""微生物组随机森林分类完整流程"""
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import RepeatedStratifiedKFold, GridSearchCV
from sklearn.metrics import roc_auc_score, classification_report
from scipy.stats import gmean
import shap
import matplotlib.pyplot as plt
# 1. 数据加载
otu = pd.read_csv("otu_table.csv", index_col=0)
meta = pd.read_csv("metadata.csv", index_col=0)
# 2. 预处理
prev = (otu > 0).mean(); otu = otu.loc[:, prev >= 0.1]
X = np.log((otu + 0.5).divide(gmean(otu + 0.5, axis=1), axis=0))
y = (meta.loc[X.index, "group"] == "disease").astype(int)
# 3. 嵌套交叉验证
outer_cv = RepeatedStratifiedKFold(5, 10, random_state=42)
params = {'n_estimators':[300,500], 'max_depth':[5,7,None], 'min_samples_leaf':[3,5]}
aucs = []
for tr, te in outer_cv.split(X, y):
gs = GridSearchCV(RandomForestClassifier(random_state=42, n_jobs=-1),
params, cv=5, scoring='roc_auc', n_jobs=-1)
gs.fit(X.iloc[tr], y.iloc[tr])
aucs.append(roc_auc_score(y.iloc[te], gs.predict_proba(X.iloc[te])[:,1]))
print(f"AUC: {np.mean(aucs):.3f} ± {np.std(aucs):.3f}")
# 4. SHAP解释
rf_final = RandomForestClassifier(n_estimators=500, max_depth=7, random_state=42, n_jobs=-1)
rf_final.fit(X, y)
shap_vals = shap.TreeExplainer(rf_final).shap_values(X)[1]
shap.summary_plot(shap_vals, X, max_display=15, show=False)
plt.savefig("shap_biomarkers.pdf", bbox_inches='tight')
面试常问点¶
Q1: 为什么选择随机森林而不是其他算法？ A: 随机森林对微生物组数据的优势：1）处理高维低样本（p>>n）时不容易过拟合；2）对数据分布无假设；3）能捕获非线性关系和交互效应；4）自带特征重要性，适合标志物发现；5）对离群值和噪声鲁棒。劣势是不如XGBoost精确，不如DNN处理复杂非线性。
Q2: SHAP值与传统特征重要性有什么区别？ A: 传统Gini importance有偏（倾向高基数分类特征）、只有全局排名。SHAP基于博弈论有理论保证（满足效率、对称、虚拟性公理），能给出：1）每个样本每个特征的贡献；2）方向性（正/负贡献）；3）交互效应；4）一致性（全局+局部解释）。
Q3: 如何避免过拟合？ A: 1）使用嵌套交叉验证（分离调参和评估）；2）控制树深度(max_depth=5-10)；3）增加min_samples_leaf；4）特征选择减少噪声；5）独立验证集验证；6）检查训练集和测试集AUC差距（<0.1为宜）。
Q4: 类别不平衡如何处理？ A: 1）使用class_weight='balanced'自动调权；2）SMOTE过采样少数类；3）使用AUC而非accuracy评估（accuracy在不平衡时misleading）；4）分层抽样保持训练/测试集比例一致；5）调整分类阈值。
Q5: 如何报告模型性能使文章可信？ A: 1）报告nested CV的mean±std AUC；2）展示完整ROC曲线；3）报告准确率、灵敏度、特异度、F1；4）独立验证队列验证；5）SHAP解释标志物生物学合理性；6）公开代码和数据；7）与已有文献方法比较。
易错点¶
信息泄露：在交叉验证之前对全部数据做特征选择/标准化，导致测试集信息泄露到训练过程。正确做法：特征选择必须在每折的训练集内部进行。
过拟合未检测：只报告训练集性能而不报告交叉验证性能。OOB估计虽可参考但不能完全替代。
特征数过多：数百个特征建模可能导致过拟合。建议先筛选到20-50个核心特征。
CLR零值处理不当：log(0)无意义，伪计数选择影响结果。常用方法：加0.5或用multiplicative replacement。
忽略混淆因素：年龄、性别、BMI、药物使用等可能是真正的分类驱动因素而非微生物组。需要在模型中控制或子集分析。
单一队列验证：只在同一队列中做CV不能证明模型泛化性。必须有独立验证队列。
补充知识¶
微生物组ML常用算法比较¶
| 算法 | 优势 | 劣势 | 适用场景 |
|---|---|---|---|
| Random Forest | 稳定，可解释 | 非最优性能 | 标志物发现 |
| XGBoost/LightGBM | 性能好 | 需调参 | 追求最优分类 |
| LASSO-LR | 简单，高可解释 | 只能线性 | 小样本/临床模型 |
| SVM | 高维表现好 | 不可解释 | 高维分类 |
| DNN | 复杂模式 | 需大样本 | 大队列(>1000) |