738. 微生物组驱动种鉴定（Driver Species Identification）¶

# ===== 用SparCC构建微生物共现网络 =====
import pandas as pd  # 导入pandas
import numpy as np  # 导入numpy
import networkx as nx  # 导入网络分析包
import matplotlib.pyplot as plt  # 导入绑图

# 读取物种丰度表
abundance = pd.read_csv("species_abundance.tsv", sep="\t", index_col=0)
# 行=样本，列=物种

# 过滤低丰度物种（至少在20%的样本中出现）
prevalence = (abundance > 0).sum() / len(abundance)  # 计算每个物种的出现频率
abundant_species = prevalence[prevalence > 0.2].index  # 保留出现率>20%的物种
abundance_filtered = abundance[abundant_species]  # 过滤后的丰度表
print(f"保留物种数: {len(abundant_species)}")

# ===== 使用SparCC计算相关性 =====
# SparCC专门为组成型数据设计，比Pearson/Spearman更准确

# 安装
pip install sparcc  # 或从GitHub安装

# 运行SparCC
# 输入：物种丰度计数表（行=样本，列=物种）
python -m sparcc \
  abundance_counts.tsv \  # 输入：原始计数表（不是相对丰度！）
  --cor_file sparcc_corr.tsv \  # 输出：SparCC相关系数矩阵
  --cov_file sparcc_cov.tsv \  # 输出：协方差矩阵
  -i 20 \  # 迭代次数
  -x 10  # 排除迭代次数

# 计算p值（通过bootstrapping）
python -m sparcc.makebootstraps \
  abundance_counts.tsv \  # 输入
  --n_boot 100 \  # bootstrap次数
  --outdir bootstraps/  # 输出目录

# 对每个bootstrap计算SparCC
for i in $(seq 0 99); do
    python -m sparcc \
      bootstraps/boot_${i}.tsv \
      --cor_file bootstraps/boot_cor_${i}.tsv \
      --cov_file bootstraps/boot_cov_${i}.tsv \
      -i 20 -x 10
done

# 计算p值
python -m sparcc.pvalues \
  sparcc_corr.tsv \  # 真实相关系数
  bootstraps/ \  # bootstrap结果目录
  --outfile pvalues.tsv \  # 输出：p值矩阵
  --n_boot 100  # bootstrap次数

# ===== 构建网络并计算中心性指标 =====
import pandas as pd
import numpy as np
import networkx as nx

# 读取SparCC结果
corr_matrix = pd.read_csv("sparcc_corr.tsv", sep="\t", index_col=0)
pval_matrix = pd.read_csv("pvalues.tsv", sep="\t", index_col=0)

# 构建网络（只保留显著的强相关）
G = nx.Graph()  # 创建无向图
species_list = corr_matrix.index.tolist()

for i in range(len(species_list)):
    for j in range(i+1, len(species_list)):
        r = corr_matrix.iloc[i, j]  # 相关系数
        p = pval_matrix.iloc[i, j]  # p值

        if abs(r) > 0.3 and p < 0.05:  # 阈值：|r|>0.3且p<0.05
            G.add_edge(
                species_list[i],
                species_list[j],
                weight=abs(r),  # 边权重=相关系数绝对值
                correlation=r,  # 保存原始相关系数（正/负）
                sign="positive" if r > 0 else "negative"
            )

print(f"节点数: {G.number_of_nodes()}")
print(f"边数: {G.number_of_edges()}")

# 计算中心性指标
degree_cent = nx.degree_centrality(G)  # 度中心性：连接数/最大可能连接数
between_cent = nx.betweenness_centrality(G)  # 中介中心性：经过该节点的最短路径比例
close_cent = nx.closeness_centrality(G)  # 接近中心性：到其他节点的平均距离

# 汇总中心性指标
centrality_df = pd.DataFrame({
    "Species": list(degree_cent.keys()),
    "Degree": list(degree_cent.values()),  # 度中心性
    "Betweenness": [between_cent[s] for s in degree_cent.keys()],  # 中介中心性
    "Closeness": [close_cent[s] for s in degree_cent.keys()],  # 接近中心性
    "N_connections": [G.degree(s) for s in degree_cent.keys()],  # 连接数
})

# 综合排名：取三个指标的平均排名
for col in ["Degree", "Betweenness", "Closeness"]:
    centrality_df[f"{col}_rank"] = centrality_df[col].rank(ascending=False)

centrality_df["Avg_rank"] = centrality_df[
    ["Degree_rank", "Betweenness_rank", "Closeness_rank"]
].mean(axis=1)  # 平均排名

# 按综合排名排序，前10就是候选驱动种
drivers = centrality_df.sort_values("Avg_rank").head(10)
print("\n=== Top 10 候选驱动种 ===")
print(drivers[["Species", "Degree", "Betweenness", "N_connections", "Avg_rank"]])

# ===== 绘制驱动种网络图 =====
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

plt.figure(figsize=(14, 14))

# 布局
pos = nx.spring_layout(G, k=2, seed=42)  # 弹簧布局

# 节点大小按度中心性
node_sizes = [degree_cent[n] * 3000 + 100 for n in G.nodes()]

# 节点颜色：驱动种标红
driver_species = set(drivers["Species"].head(5).tolist())  # Top5驱动种
node_colors = ["#e74c3c" if n in driver_species else "#3498db" for n in G.nodes()]

# 边颜色：正相关=绿色，负相关=红色
edge_colors = ["#2ecc71" if G[u][v]["sign"] == "positive" else "#e74c3c"
               for u, v in G.edges()]

# 绘制
nx.draw_networkx_nodes(G, pos, node_size=node_sizes, 
                       node_color=node_colors, alpha=0.8)
nx.draw_networkx_edges(G, pos, edge_color=edge_colors, 
                       alpha=0.3, width=0.5)

# 只给驱动种加标签
labels = {n: n.split("__")[-1][:15] for n in G.nodes() if n in driver_species}
nx.draw_networkx_labels(G, pos, labels, font_size=8, font_weight="bold")

# 图例
legend_elements = [
    mpatches.Patch(color="#e74c3c", label="Driver Species"),
    mpatches.Patch(color="#3498db", label="Other Species"),
]
plt.legend(handles=legend_elements, loc="upper left", fontsize=12)
plt.title("Microbial Co-occurrence Network with Driver Species", fontsize=14)
plt.axis("off")
plt.tight_layout()
plt.savefig("driver_species_network.png", dpi=300, bbox_inches="tight")
plt.show()

# ===== 用随机森林找驱动种 =====
from sklearn.ensemble import RandomForestClassifier  # 导入随机森林
from sklearn.model_selection import cross_val_score  # 导入交叉验证
import pandas as pd
import numpy as np

# 准备数据
X = abundance_filtered.values  # 特征：物种丰度
y = metadata["Group"].values  # 标签：分组（如Healthy vs T2D）

# 训练随机森林
rf = RandomForestClassifier(
    n_estimators=500,  # 500棵树
    max_depth=None,    # 不限深度
    random_state=42,   # 随机种子
    n_jobs=-1           # 使用所有CPU
)

# 10折交叉验证评估模型
cv_scores = cross_val_score(rf, X, y, cv=10, scoring="accuracy")
print(f"交叉验证准确率: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# 训练最终模型并提取特征重要性
rf.fit(X, y)
importances = pd.Series(
    rf.feature_importances_,  # 特征重要性（基于基尼不纯度减少）
    index=abundance_filtered.columns  # 物种名
).sort_values(ascending=False)

# Top 20重要物种
print("\n=== Top 20 Feature Importance ===")
print(importances.head(20))

# 绘制特征重要性条形图
plt.figure(figsize=(10, 8))
top20 = importances.head(20)
plt.barh(range(len(top20)), top20.values, color="#3498db")
plt.yticks(range(len(top20)), top20.index)
plt.xlabel("Feature Importance (Gini)")
plt.title("Top 20 Driver Species by Random Forest")
plt.gca().invert_yaxis()  # 最重要的在上面
plt.tight_layout()
plt.savefig("rf_driver_species.png", dpi=300)
plt.show()

# ===== 使用LEfSe鉴定差异物种 =====
# LEfSe是最经典的微生物组差异分析工具

# 安装
conda install -c bioconda lefse  # 安装LEfSe

# 第一步：格式化输入文件
# LEfSe需要特定格式：第一行是分组，后面每行是一个物种的丰度
lefse_format_input.py \
  lefse_input.tsv \  # 输入：格式化的丰度表
  lefse_formatted.in \  # 输出：LEfSe格式
  -c 1 \  # 第1行是类别（分组）
  -u 2 \  # 第2行是子分组（可选）
  -o 1000000  # 归一化到百万

# 第二步：运行LEfSe分析
lefse_run.py \
  lefse_formatted.in \  # 输入
  lefse_result.res \  # 输出：结果文件
  -l 2.0 \  # LDA阈值（默认2.0）
  -a 0.05 \  # KW检验的p值阈值
  -w 0.05  # Wilcoxon检验的p值阈值

# 第三步：绘制结果
lefse_plot_res.py \
  lefse_result.res \  # 输入：结果文件
  lefse_barplot.png \  # 输出：条形图
  --dpi 300  # 分辨率

# 绘制进化分支图
lefse_plot_cladogram.py \
  lefse_result.res \  # 输入
  lefse_cladogram.png \  # 输出：进化分支图
  --dpi 300 \
  --format png

# ===== 驱动种鉴定速查 =====

# 方法一：网络分析
# SparCC → 网络构建 → 度/中介/接近中心性 → Top10

# 方法二：随机森林
# RF分类 → feature_importances_ → Top20

# 方法三：LEfSe
lefse_format_input.py input.tsv formatted.in -c 1
lefse_run.py formatted.in result.res -l 2.0
lefse_plot_res.py result.res barplot.png

# 综合策略：取三种方法的交集

# 网络中心性指标含义：
# 度中心性 = 连了多少条边（人缘好）
# 中介中心性 = 多少最短路径经过它（信息枢纽）
# 接近中心性 = 到其他节点平均距离（位置核心）