681 微生物组功能贡献度分析（SHAP/Shapley）¶

# SHAP分析实战
import shap           # SHAP库
import numpy as np    # 数值计算
import pandas as pd   # 数据处理
from sklearn.ensemble import RandomForestClassifier  # 随机森林
from sklearn.model_selection import train_test_split  # 数据分割
import matplotlib.pyplot as plt  # 绑图

# ============ 1. 准备数据和训练模型 ============
# 读取微生物组数据
abundance = pd.read_csv("species_abundance.tsv", sep='\t',
                         index_col=0)    # 物种丰度表
metadata = pd.read_csv("metadata.tsv", sep='\t',
                        index_col=0)     # 元数据

X = abundance.T                          # 特征矩阵（样本×物种）
y = (metadata["Group"] == "Disease").astype(int)  # 标签（0/1）

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# 训练随机森林模型
rf = RandomForestClassifier(
    n_estimators=500,                    # 500棵树
    max_depth=10,                        # 最大深度
    min_samples_leaf=5,                  # 叶节点最小样本数
    random_state=42,                     # 随机种子
    n_jobs=-1                            # 使用所有CPU
)
rf.fit(X_train, y_train)                 # 训练模型
print(f"训练集准确率: {rf.score(X_train, y_train):.3f}")
print(f"测试集准确率: {rf.score(X_test, y_test):.3f}")

# ============ 2. 计算SHAP值 ============
# 使用TreeSHAP（针对树模型的快速算法）
explainer = shap.TreeExplainer(rf)       # 创建解释器
shap_values = explainer.shap_values(X_test)  # 计算SHAP值
# shap_values是一个列表：[class_0的SHAP值, class_1的SHAP值]
# 我们关注class_1（疾病）的SHAP值

shap_disease = shap_values[1]            # 疾病类别的SHAP值
print(f"SHAP值矩阵维度: {shap_disease.shape}")
# shape = (n_samples, n_features)
# 每个值表示该特征对该样本预测的贡献

# ============ 3. 全局特征重要性（哪些物种最重要？） ============
# SHAP summary plot——最核心的图
plt.figure(figsize=(12, 8))
shap.summary_plot(
    shap_disease,                        # SHAP值
    X_test,                              # 特征值
    feature_names=X.columns.tolist(),    # 特征名
    max_display=20,                      # 显示前20个特征
    show=False
)
plt.tight_layout()
plt.savefig("shap_summary.png", dpi=150, bbox_inches='tight')
plt.close()
print("SHAP summary plot已保存")

# SHAP bar plot——平均绝对SHAP值排序
plt.figure(figsize=(10, 8))
shap.summary_plot(
    shap_disease,
    X_test,
    feature_names=X.columns.tolist(),
    plot_type="bar",                     # 柱状图模式
    max_display=20,
    show=False
)
plt.tight_layout()
plt.savefig("shap_importance.png", dpi=150, bbox_inches='tight')
plt.close()

# 提取特征重要性排名
mean_abs_shap = np.abs(shap_disease).mean(axis=0)  # 平均|SHAP|
importance_df = pd.DataFrame({
    "species": X.columns,
    "mean_abs_shap": mean_abs_shap
}).sort_values("mean_abs_shap", ascending=False)

print("\nTop 10 最重要的物种（SHAP）:")
for i, row in importance_df.head(10).iterrows():
    print(f"  {row['species']}: SHAP={row['mean_abs_shap']:.4f}")

# ============ 4. 单样本解释（为什么这个人被预测为疾病？） ============
# Waterfall plot——解释单个预测
sample_idx = 0                           # 第一个测试样本
plt.figure(figsize=(12, 6))
shap.waterfall_plot(
    shap.Explanation(
        values=shap_disease[sample_idx],  # 该样本的SHAP值
        base_values=explainer.expected_value[1],  # 基线值
        data=X_test.iloc[sample_idx],     # 该样本的特征值
        feature_names=X.columns.tolist()  # 特征名
    ),
    max_display=15,                      # 显示前15个特征
    show=False
)
plt.tight_layout()
plt.savefig("shap_waterfall_sample0.png", dpi=150, bbox_inches='tight')
plt.close()

# Force plot——另一种单样本可视化
shap.force_plot(
    explainer.expected_value[1],         # 基线
    shap_disease[sample_idx],            # SHAP值
    X_test.iloc[sample_idx],             # 特征值
    feature_names=X.columns.tolist(),    # 特征名
    matplotlib=True,                     # 用matplotlib
    show=False
)
plt.savefig("shap_force_sample0.png", dpi=150, bbox_inches='tight')
plt.close()

# ============ 5. SHAP依赖图（非线性关系） ============
# 查看单个物种丰度与其SHAP值的关系
top_species = importance_df.iloc[0]["species"]  # 最重要的物种
plt.figure(figsize=(8, 6))
shap.dependence_plot(
    top_species,                         # 目标特征
    shap_disease,                        # SHAP值
    X_test,                              # 特征值
    feature_names=X.columns.tolist(),    # 特征名
    interaction_index="auto",            # 自动选择交互特征
    show=False
)
plt.title(f"{top_species} 的SHAP依赖图")
plt.tight_layout()
plt.savefig(f"shap_dependence_{top_species}.png", dpi=150)
plt.close()

# ============ 6. SHAP交互效应 ============
# 计算特征间的交互SHAP值
shap_interaction = explainer.shap_interaction_values(X_test)
# shape = (n_samples, n_features, n_features)

# 提取最强的交互对
interaction_matrix = np.abs(shap_interaction[1]).mean(axis=0)
np.fill_diagonal(interaction_matrix, 0)  # 去除自身交互
top_pairs = []
n = interaction_matrix.shape[0]
for i in range(n):
    for j in range(i+1, n):
        top_pairs.append({
            "species_1": X.columns[i],
            "species_2": X.columns[j],
            "interaction_strength": interaction_matrix[i, j]
        })
interaction_df = pd.DataFrame(top_pairs).sort_values(
    "interaction_strength", ascending=False
)
print("\nTop 10 物种交互:")
print(interaction_df.head(10).to_string(index=False))

# 比较SHAP和传统特征重要性方法
import matplotlib.pyplot as plt  # 绑图

# 1. 随机森林内置重要性（Gini/MDI）
rf_importance = pd.DataFrame({
    "species": X.columns,
    "gini_importance": rf.feature_importances_
}).sort_values("gini_importance", ascending=False)

# 2. 置换重要性
from sklearn.inspection import permutation_importance
perm_imp = permutation_importance(
    rf, X_test, y_test,
    n_repeats=30,                        # 30次重复
    random_state=42
)
perm_importance = pd.DataFrame({
    "species": X.columns,
    "perm_importance": perm_imp.importances_mean
}).sort_values("perm_importance", ascending=False)

# 3. SHAP重要性（已计算）
shap_importance = importance_df.copy()

# 4. 三种方法对比
fig, axes = plt.subplots(1, 3, figsize=(18, 8))
top_n = 15

# Gini重要性
axes[0].barh(rf_importance.head(top_n)["species"][::-1],
             rf_importance.head(top_n)["gini_importance"][::-1],
             color="#FF6B6B")
axes[0].set_title("Gini重要性（MDI）")
axes[0].set_xlabel("重要性")

# 置换重要性
axes[1].barh(perm_importance.head(top_n)["species"][::-1],
             perm_importance.head(top_n)["perm_importance"][::-1],
             color="#4ECDC4")
axes[1].set_title("置换重要性")
axes[1].set_xlabel("重要性")

# SHAP重要性
axes[2].barh(shap_importance.head(top_n)["species"][::-1],
             shap_importance.head(top_n)["mean_abs_shap"][::-1],
             color="#45B7D1")
axes[2].set_title("SHAP重要性")
axes[2].set_xlabel("平均|SHAP值|")

plt.suptitle("三种特征重要性方法对比", fontsize=14)
plt.tight_layout()
plt.savefig("importance_comparison.png", dpi=150)
print("对比图已保存: importance_comparison.png")

# 排名一致性检验（Spearman相关）
from scipy.stats import spearmanr
# 合并排名
merged = rf_importance.merge(perm_importance, on="species")
merged = merged.merge(shap_importance, on="species")
r_gini_perm, _ = spearmanr(merged["gini_importance"],
                            merged["perm_importance"])
r_gini_shap, _ = spearmanr(merged["gini_importance"],
                            merged["mean_abs_shap"])
r_perm_shap, _ = spearmanr(merged["perm_importance"],
                            merged["mean_abs_shap"])
print(f"\n排名一致性（Spearman）:")
print(f"  Gini vs 置换: r={r_gini_perm:.3f}")
print(f"  Gini vs SHAP: r={r_gini_shap:.3f}")
print(f"  置换 vs SHAP: r={r_perm_shap:.3f}")

# 按分类学级别汇总SHAP贡献度
import pandas as pd   # 数据处理
import numpy as np    # 数值计算

def summarize_shap_by_taxonomy(shap_values, feature_names, taxonomy_df,
                                level="Phylum"):
    """按分类学级别汇总SHAP贡献度"""
    # 创建SHAP数据框
    shap_df = pd.DataFrame(shap_values, columns=feature_names)

    # 按分类级别分组
    taxon_map = taxonomy_df.set_index("Species")[level]  # 物种→门映射
    grouped_shap = {}                    # 按门汇总
    for taxon in taxon_map.unique():
        species_in_taxon = taxon_map[taxon_map == taxon].index  # 该门的物种
        cols = [c for c in species_in_taxon if c in shap_df.columns]
        if cols:
            grouped_shap[taxon] = shap_df[cols].sum(axis=1)  # 求和

    grouped_df = pd.DataFrame(grouped_shap)

    # 平均绝对贡献
    mean_contrib = grouped_df.abs().mean().sort_values(ascending=False)
    print(f"\n{level}级别的平均|SHAP|贡献:")
    for taxon, val in mean_contrib.items():
        print(f"  {taxon}: {val:.4f}")

    return grouped_df, mean_contrib

# 功能贡献度分析
def functional_contribution(shap_values, feature_names, ko_annotation):
    """将SHAP值从物种级别映射到功能级别"""
    # ko_annotation: 物种→KO功能映射
    shap_df = pd.DataFrame(shap_values, columns=feature_names)

    # 按KO汇总
    ko_contrib = {}
    for ko, species_list in ko_annotation.items():
        cols = [s for s in species_list if s in shap_df.columns]
        if cols:
            ko_contrib[ko] = shap_df[cols].sum(axis=1).mean()

    contrib_df = pd.Series(ko_contrib).sort_values(ascending=False)
    print("\nTop 10 功能贡献（SHAP）:")
    print(contrib_df.head(10))
    return contrib_df

# Shapley群落贡献度（哪些物种"组合"最重要？）
def community_contribution_analysis(rf_model, X_test, top_n=5):
    """分析物种组合的协同贡献"""
    explainer = shap.TreeExplainer(rf_model)
    shap_values = explainer.shap_values(X_test)[1]  # 疾病类

    # 找到top_n个最重要的物种
    mean_shap = np.abs(shap_values).mean(axis=0)
    top_idx = np.argsort(mean_shap)[-top_n:]  # top物种索引
    top_species = X_test.columns[top_idx].tolist()

    # 计算这些物种的贡献占比
    top_shap_sum = np.abs(shap_values[:, top_idx]).sum(axis=1)
    total_shap_sum = np.abs(shap_values).sum(axis=1)
    contribution_pct = (top_shap_sum / total_shap_sum * 100)

    print(f"\nTop {top_n} 物种: {top_species}")
    print(f"贡献占比: {contribution_pct.mean():.1f}% ± {contribution_pct.std():.1f}%")

    return top_species, contribution_pct

# SHAP分析流程
1. 训练机器学习模型（随机森林/XGBoost）
2. 创建SHAP解释器（TreeExplainer最快）
3. 计算SHAP值
4. 全局重要性: summary_plot / bar_plot
5. 局部解释: waterfall_plot / force_plot
6. 非线性关系: dependence_plot
7. 交互效应: interaction_values（可选）

# SHAP vs 传统特征重要性
Gini重要性: 快但有偏（偏向高基数特征）
置换重要性: 无偏但慢
SHAP: 最完整（方向+大小+交互），有理论保证

# SHAP图表选择
summary_plot: 最核心，显示方向和大小
bar_plot: 简洁版，只显示重要性排序
waterfall_plot: 解释单个样本
dependence_plot: 特征的非线性效应
force_plot: 直观的单样本解释

# 在微生物组中的应用
疾病分类器解释: 哪些菌驱动了疾病预测
生物标志物发现: SHAP Top特征=候选标志物
功能解释: SHAP值按功能注释汇总
群落贡献度: 哪些物种组合最关键