740. 环境微生物组生态位建模¶

# ===== 计算Levins生态位宽度 =====
import pandas as pd  # 导入pandas
import numpy as np  # 导入numpy

# 读取物种在不同环境梯度上的丰度
# 假设我们沿pH梯度采样了多个样本
abundance = pd.read_csv("species_by_pH.tsv", sep="\t", index_col=0)
# 行=pH区间(如 4-5, 5-6, 6-7, 7-8, 8-9)
# 列=物种
# 值=相对丰度

def levins_niche_width(abundances):
    """
    Levins生态位宽度指数
    B = 1 / sum(pi^2)
    标准化后: Ba = (B - 1) / (n - 1)
    """
    p = abundances / abundances.sum()  # 转为比例
    p = p[p > 0]  # 去掉0值
    B = 1.0 / np.sum(p ** 2)  # Levins B
    n = len(abundances)  # 环境类别数
    Ba = (B - 1) / (n - 1)  # 标准化到0-1
    return Ba

# 计算每个物种的生态位宽度
niche_widths = {}
for species in abundance.columns:  # 遍历每个物种
    niche_widths[species] = levins_niche_width(abundance[species].values)

niche_df = pd.DataFrame.from_dict(
    niche_widths, orient="index", columns=["Niche_Width"]
).sort_values("Niche_Width", ascending=False)

print("=== 生态位宽度排名 ===")
print(niche_df.head(20))
# Niche_Width接近1 = 广适种（generalist，什么环境都能活）
# Niche_Width接近0 = 特化种（specialist，只能在特定环境活）

# ===== 计算Pianka生态位重叠 =====
def pianka_overlap(sp1_abund, sp2_abund):
    """
    Pianka生态位重叠指数
    O = sum(pi * qi) / sqrt(sum(pi^2) * sum(qi^2))
    值域: 0-1, 越大重叠越多
    """
    p = sp1_abund / sp1_abund.sum()  # 物种1的比例分布
    q = sp2_abund / sp2_abund.sum()  # 物种2的比例分布

    numerator = np.sum(p * q)  # 分子
    denominator = np.sqrt(np.sum(p**2) * np.sum(q**2))  # 分母

    if denominator == 0:
        return 0
    return numerator / denominator

# 计算所有物种对之间的生态位重叠
species_list = abundance.columns.tolist()
n = len(species_list)
overlap_matrix = np.zeros((n, n))  # 初始化重叠矩阵

for i in range(n):
    for j in range(n):
        overlap_matrix[i, j] = pianka_overlap(
            abundance.iloc[:, i].values,
            abundance.iloc[:, j].values
        )

overlap_df = pd.DataFrame(overlap_matrix, 
                           index=species_list, columns=species_list)

# 绘制热图
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(12, 10))
sns.heatmap(overlap_df, cmap="YlOrRd", vmin=0, vmax=1,
            annot=False, square=True)
plt.title("Niche Overlap Between Species (Pianka Index)")
plt.tight_layout()
plt.savefig("niche_overlap_heatmap.png", dpi=300)
plt.show()

# ===== 用随机森林预测微生物分布 =====
from sklearn.ensemble import RandomForestRegressor  # 导入随机森林回归
from sklearn.model_selection import cross_val_score  # 导入交叉验证
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 读取数据
# 环境变量表（行=样本，列=环境因子）
env_data = pd.read_csv("environmental_variables.tsv", sep="\t", index_col=0)
# 列名示例: pH, Temperature, DO, Salinity, TotalN, TotalP

# 物种丰度表
species_data = pd.read_csv("species_abundance.tsv", sep="\t", index_col=0)

# 选择目标物种
target_species = "Lactobacillus"
y = species_data[target_species].values  # 目标变量：物种丰度
X = env_data.values  # 特征：环境变量
feature_names = env_data.columns.tolist()

# 训练随机森林
rf = RandomForestRegressor(
    n_estimators=500,  # 500棵树
    max_depth=10,      # 最大深度10
    random_state=42,
    n_jobs=-1
)

# 交叉验证评估
cv_scores = cross_val_score(rf, X, y, cv=5, scoring="r2")
print(f"R² = {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# 拟合最终模型
rf.fit(X, y)

# 特征重要性 = 哪些环境因子最影响这种菌
importances = pd.Series(rf.feature_importances_, index=feature_names)
importances = importances.sort_values(ascending=True)

plt.figure(figsize=(8, 5))
importances.plot(kind="barh", color="#3498db")
plt.xlabel("Feature Importance")
plt.title(f"Environmental Drivers of {target_species}")
plt.tight_layout()
plt.savefig(f"{target_species}_env_drivers.png", dpi=300)
plt.show()

# ===== 偏依赖图：展示环境因子的非线性效应 =====
from sklearn.inspection import PartialDependenceDisplay

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
top3_features = importances.index[-3:].tolist()  # 最重要的3个因子

for i, feature in enumerate(top3_features):
    feat_idx = feature_names.index(feature)
    PartialDependenceDisplay.from_estimator(
        rf, X, [feat_idx], feature_names=feature_names, ax=axes[i]
    )
    axes[i].set_title(f"Effect of {feature}")

plt.suptitle(f"Partial Dependence: {target_species}", fontsize=14)
plt.tight_layout()
plt.savefig(f"{target_species}_partial_dependence.png", dpi=300)
plt.show()

# ===== R语言：典范对应分析（CCA）=====
library(vegan)    # 群落生态分析
library(ggplot2)  # 绑图

# 读取数据
species <- read.table("species_abundance.tsv", header=TRUE,
                      row.names=1, sep="\t")
env <- read.table("environmental_variables.tsv", header=TRUE,
                  row.names=1, sep="\t")

# 选择CCA还是RDA？先做DCA看梯度长度
dca <- decorana(species)  # DCA分析
print(dca)
# 如果第一轴梯度长度 > 4 → 用CCA（单峰响应）
# 如果第一轴梯度长度 < 3 → 用RDA（线性响应）

# 运行CCA
cca_result <- cca(species ~ ., data=env)  # 用所有环境变量做CCA
summary(cca_result)

# 检验整体显著性
anova(cca_result, permutations=999)  # 全局检验

# 检验每个环境因子的显著性
anova(cca_result, by="terms", permutations=999)  # 逐项检验

# VIF检查多重共线性
vif.cca(cca_result)  # VIF>10的变量考虑去掉

# 方差分解
# 看各环境因子解释了多少变异
varpart_result <- varpart(species, 
                          ~pH, ~Temperature, ~DO,  # 分三组环境变量
                          data=env)
plot(varpart_result)  # 绘制Venn图

# 绘图
plot(cca_result, display=c("sites", "bp"), type="n")  # 创建空图
points(cca_result, display="sites", pch=19, 
       col=ifelse(env$Group=="Healthy", "blue", "red"))  # 添加样本点
text(cca_result, display="bp", col="darkgreen", cex=0.8)  # 添加环境变量箭头

# ===== 生态位分析速查 =====

# 生态位宽度（Levins）
# Ba = (1/sum(pi²) - 1) / (n - 1)
# 接近1=广适种, 接近0=特化种

# 生态位重叠（Pianka）
# O = sum(pi*qi) / sqrt(sum(pi²)*sum(qi²))
# 接近1=完全重叠, 接近0=不重叠

# R语言CCA/RDA选择：
# DCA梯度长度>4 → CCA
# DCA梯度长度<3 → RDA

# Python SDM：
# RandomForest/MaxEnt → 特征重要性 → 偏依赖图

# 关键R包：vegan（CCA, RDA, DCA, PERMANOVA）
# 关键Python包：sklearn（RF, PDPlot）