741. 微生物组多样性指数深度解读¶

# ===== Alpha多样性指数的数学计算与解读 =====
import numpy as np  # 导入numpy
import pandas as pd  # 导入pandas
from scipy.stats import entropy  # 导入熵函数

# 模拟一个样本的物种丰度
sample_counts = np.array([100, 80, 50, 30, 20, 10, 5, 3, 1, 1])
# 10个物种，丰度从100到1不等

# ===== 1. 丰富度（Richness）=====
# 最简单：数有多少种
richness = np.sum(sample_counts > 0)  # 丰度>0的物种数
print(f"观察到的丰富度(Observed): {richness}")
# 问题：受测序深度影响大！深度越高发现的物种越多

# ===== 2. Chao1估计器 =====
# 估算"真实"丰富度（包括没测到的物种）
singletons = np.sum(sample_counts == 1)  # 只出现1次的物种数
doubletons = np.sum(sample_counts == 2)  # 只出现2次的物种数

if doubletons > 0:
    chao1 = richness + (singletons ** 2) / (2 * doubletons)
else:
    chao1 = richness + singletons * (singletons - 1) / 2
print(f"Chao1估计的总丰富度: {chao1:.1f}")
# Chao1远大于观察值 → 还有很多未发现的物种

# ===== 3. Shannon指数 (H') =====
# 同时考虑丰富度和均匀度
# H' = -sum(pi * ln(pi))
rel_abund = sample_counts / sample_counts.sum()  # 相对丰度
shannon = -np.sum(rel_abund * np.log(rel_abund))  # Shannon公式
print(f"Shannon H': {shannon:.4f}")
# H'高 → 物种多且丰度分布均匀
# H'低 → 物种少或分布极不均匀（一种菌占主导）

# Shannon的等效物种数（将H'转为"相当于多少个均匀分布的物种"）
effective_species = np.exp(shannon)  # e^H'
print(f"有效物种数: {effective_species:.1f}")
# 更直观！"这个样本相当于有X个均匀分布的物种"

# ===== 4. Simpson指数 =====
# 随机取两个个体属于不同物种的概率
# D = 1 - sum(pi^2)
simpson = 1 - np.sum(rel_abund ** 2)  # Simpson多样性
inv_simpson = 1 / np.sum(rel_abund ** 2)  # 逆Simpson
print(f"Simpson D: {simpson:.4f}")
print(f"Inverse Simpson: {inv_simpson:.4f}")
# Simpson对优势种更敏感，Shannon对稀有种更敏感

# ===== 5. Pielou均匀度 =====
# 丰度分布有多"均匀"
pielou_e = shannon / np.log(richness)  # J = H'/ln(S)
print(f"Pielou均匀度 J: {pielou_e:.4f}")
# J接近1 = 非常均匀（每种菌一样多）
# J接近0 = 极不均匀（一种菌独大）

# ===== 6. Faith系统发育多样性 =====
# 所有物种在进化树上覆盖的总枝长
# 需要系统发育树，这里给出计算逻辑
"""
Faith's PD = sum(所有存在物种的祖先节点到叶子节点的枝长)
用skbio计算：
from skbio.diversity import alpha_diversity
pd_value = alpha_diversity('faith_pd', counts, tree=tree, otu_ids=otu_ids)
"""

# ===== 对比不同场景下各指数的表现 =====
import matplotlib.pyplot as plt

# 场景1：均匀分布
uniform = np.array([20, 20, 20, 20, 20])  # 5种菌各20个

# 场景2：不均匀分布
uneven = np.array([96, 1, 1, 1, 1])  # 1种菌占主导

# 场景3：多物种均匀
diverse = np.array([10]*10)  # 10种菌各10个

# 场景4：多物种不均匀
diverse_uneven = np.array([50, 20, 10, 5, 5, 3, 3, 2, 1, 1])

scenarios = {
    "5sp均匀": uniform,
    "5sp不均匀": uneven,
    "10sp均匀": diverse,
    "10sp不均匀": diverse_uneven
}

results = []
for name, counts in scenarios.items():
    p = counts / counts.sum()
    h = -np.sum(p * np.log(p))
    d = 1 - np.sum(p**2)
    j = h / np.log(len(counts))

    results.append({
        "场景": name,
        "丰富度": len(counts),
        "Shannon": round(h, 3),
        "Simpson": round(d, 3),
        "Pielou": round(j, 3),
        "有效种数": round(np.exp(h), 1)
    })

result_df = pd.DataFrame(results)
print(result_df.to_string(index=False))

# 关键洞察：
# 1. "5sp均匀" vs "5sp不均匀"：丰富度相同但Shannon/Simpson差异大 → 均匀度的影响
# 2. "5sp均匀" vs "10sp均匀"：均匀度相同但Shannon不同 → 丰富度的影响
# 3. Shannon对稀有种更敏感，Simpson对优势种更敏感

# ===== Beta多样性距离度量对比 =====
from scipy.spatial.distance import braycurtis, jaccard
import numpy as np

# 两个样本的丰度
sample_a = np.array([50, 30, 20, 0, 0])
sample_b = np.array([0, 0, 20, 30, 50])
sample_c = np.array([40, 25, 20, 10, 5])

# 1. Jaccard距离：只看"有没有"，不看"有多少"
# J = 1 - (共有种数 / 总种数)
a_binary = (sample_a > 0).astype(int)
b_binary = (sample_b > 0).astype(int)
c_binary = (sample_c > 0).astype(int)

j_ab = jaccard(a_binary, b_binary)
j_ac = jaccard(a_binary, c_binary)
print(f"Jaccard(A,B) = {j_ab:.3f}")  # A和B只共享1种 → 距离大
print(f"Jaccard(A,C) = {j_ac:.3f}")  # A和C共享3种 → 距离小

# 2. Bray-Curtis距离：看"有多少"的差异
# BC = sum|ai-bi| / sum(ai+bi)
bc_ab = braycurtis(sample_a, sample_b)
bc_ac = braycurtis(sample_a, sample_c)
print(f"Bray-Curtis(A,B) = {bc_ab:.3f}")  # 组成完全不同
print(f"Bray-Curtis(A,C) = {bc_ac:.3f}")  # 组成相似

# 3. UniFrac（需要进化树）
# Weighted UniFrac: 考虑丰度和进化关系
# Unweighted UniFrac: 只考虑"有无"和进化关系
# 区别：Weighted对优势种敏感，Unweighted对稀有种敏感

# ===== R语言PERMANOVA检验 =====
library(vegan)

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

# 计算Bray-Curtis距离
bc_dist <- vegdist(species, method="bray")

# PERMANOVA检验
result <- adonis2(bc_dist ~ Group + Age + Sex,  # 分组+协变量
                  data=metadata,
                  permutations=9999,  # 置换次数
                  method="bray")
print(result)
# R² = 解释的变异比例
# p值 < 0.05 → 组间差异显著

# 注意！PERMANOVA假设组内方差齐性
# 用betadisper检验
bd <- betadisper(bc_dist, metadata$Group)
permutest(bd)  # p>0.05 → 方差齐性 → PERMANOVA可靠
               # p<0.05 → 方差不齐 → PERMANOVA结果需谨慎解读

# ===== PCoA可视化 =====
from sklearn.manifold import MDS
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

# 读取数据
abundance = pd.read_csv("species_abundance.tsv", sep="\t", index_col=0)
metadata = pd.read_csv("metadata.tsv", sep="\t", index_col=0)

# 计算Bray-Curtis距离矩阵
dist_matrix = squareform(pdist(abundance.values, metric="braycurtis"))

# PCoA（用MDS with precomputed实现）
pcoa = MDS(n_components=2, dissimilarity="precomputed", random_state=42)
coords = pcoa.fit_transform(dist_matrix)

# 绘图
fig, ax = plt.subplots(figsize=(10, 8))
groups = metadata["Group"].unique()
colors = {"Healthy": "#2ecc71", "T2D": "#e74c3c"}

for group in groups:
    mask = metadata["Group"] == group
    ax.scatter(coords[mask, 0], coords[mask, 1],
               c=colors[group], label=group, s=60, alpha=0.7)

    # 添加95%置信椭圆
    from matplotlib.patches import Ellipse
    group_coords = coords[mask]
    mean_x, mean_y = group_coords.mean(axis=0)
    cov = np.cov(group_coords.T)
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    angle = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]))
    width, height = 2 * np.sqrt(eigenvalues) * 2.447  # 95%置信椭圆
    ellipse = Ellipse(xy=(mean_x, mean_y), width=width, height=height,
                      angle=angle, alpha=0.15, color=colors[group])
    ax.add_patch(ellipse)

ax.set_xlabel("PCoA1")
ax.set_ylabel("PCoA2")
ax.set_title("PCoA of Microbial Community (Bray-Curtis)")
ax.legend()
plt.tight_layout()
plt.savefig("pcoa_plot.png", dpi=300)
plt.show()

# ===== 稀释曲线（判断测序是否足够）=====
def rarefaction_curve(counts, n_steps=50):
    """绘制稀释曲线"""
    total = counts.sum()
    steps = np.linspace(100, total, n_steps).astype(int)
    richness_at_depth = []

    for depth in steps:
        # 随机抽样多次取平均
        n_repeats = 10
        rich_values = []
        for _ in range(n_repeats):
            # 按概率抽样depth条reads
            probs = counts / counts.sum()
            subsample = np.random.multinomial(depth, probs)
            rich_values.append(np.sum(subsample > 0))
        richness_at_depth.append(np.mean(rich_values))

    return steps, richness_at_depth

# 绘制多个样本的稀释曲线
plt.figure(figsize=(10, 6))
for sample_name in abundance.index[:10]:  # 前10个样本
    counts = abundance.loc[sample_name].values
    counts = counts[counts > 0]
    depths, richness = rarefaction_curve(counts)
    plt.plot(depths, richness, alpha=0.7, label=sample_name[:10])

plt.xlabel("Sequencing Depth")
plt.ylabel("Observed Species")
plt.title("Rarefaction Curves")
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
plt.tight_layout()
plt.savefig("rarefaction_curves.png", dpi=300)
plt.show()
# 曲线趋平 → 测序够了
# 曲线还在上升 → 测序不够

# ===== 多样性指数速查 =====

# Alpha多样性：
# 丰富度 = 物种数（受测序深度影响大）
# Chao1 = 估计真实丰富度（考虑未发现物种）
# Shannon H' = -Σ(pi·ln(pi))，同时反映丰富度和均匀度
# Simpson D = 1-Σ(pi²)，对优势种敏感
# Pielou J = H'/ln(S)，纯均匀度
# Faith PD = 进化树枝长之和

# Beta多样性：
# Jaccard = 共有种比例（定性：有/无）
# Bray-Curtis = 丰度差异（定量）
# UniFrac = 考虑进化关系的差异
#   Weighted = 考虑丰度+进化
#   Unweighted = 只考虑有无+进化

# 统计检验：
# Alpha差异: Wilcoxon/Kruskal-Wallis
# Beta差异: PERMANOVA (adonis2)
# 方差齐性: betadisper

# R关键代码：
# diversity(x, "shannon")     # Shannon
# vegdist(x, "bray")          # Bray-Curtis
# adonis2(dist~Group, data=meta)  # PERMANOVA