680 微生物组共丰度网络分析¶

# SparCC分析
# SparCC专门为组成数据设计，解决假负相关问题

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

# 方法1: 命令行版本
# sparcc.py otu_table.tsv -i 20 --cor_file sparcc_correlation.tsv
# sparcc.py otu_table.tsv -i 20 --cov_file sparcc_covariance.tsv

# 方法2: Python实现SparCC核心思想
import numpy as np   # 数值计算
import pandas as pd  # 数据处理
from scipy.stats import pearsonr  # 相关性

def sparcc_correlation(counts, n_iter=20, threshold=0.1):
    """
    SparCC相关性计算（简化版）
    核心思想：迭代地从组成数据中估计真实的方差和协方差
    """
    n_otus = counts.shape[0]             # 物种数
    # 对数变换
    log_data = np.log(counts + 1)        # 加1避免log(0)
    # 初始化方差-协方差矩阵
    var_log = np.var(log_data, axis=1)   # 对数方差
    # 迭代估计
    for iteration in range(n_iter):
        # 计算对数比的方差
        T = np.zeros((n_otus, n_otus))   # 对数比方差矩阵
        for i in range(n_otus):
            for j in range(n_otus):
                if i != j:
                    T[i, j] = np.var(log_data[i] - log_data[j])
        # 从T矩阵估计方差和协方差
        # V_i = median_j((T_ij + V_i + V_j) / 2 - V_j)
        # ... 迭代收敛
    # 转换为相关系数
    # cor_ij = cov_ij / sqrt(var_i * var_j)
    print(f"SparCC迭代完成（{n_iter}轮）")
    return T  # 实际中用完整版SparCC

# 推荐：使用已有的SparCC实现
# 如fastspar (C++加速版，推荐)

# FastSpar (SparCC的C++实现，推荐)
conda install -c bioconda fastspar  # 安装

# 1. 计算SparCC相关性
fastspar \
  --otu_table otu_table.tsv \    # OTU表（行=物种，列=样本）
  --correlation sparcc_cor.tsv \ # 相关性矩阵输出
  --covariance sparcc_cov.tsv \  # 协方差矩阵输出
  --iterations 50 \              # SparCC迭代次数
  --threads 8                    # 线程数

# 2. 置换检验（评估显著性）
# 生成置换数据
mkdir bootstrap_counts
fastspar_bootstrap \
  --otu_table otu_table.tsv \    # 原始OTU表
  --number 1000 \                # 1000次置换
  --prefix bootstrap_counts/perm  # 输出前缀

# 对每个置换数据计算相关性
mkdir bootstrap_correlations
parallel fastspar \
  --otu_table {} \
  --correlation bootstrap_correlations/{/} \
  --covariance /dev/null \
  --iterations 10 \
  ::: bootstrap_counts/perm_*    # 并行处理

# 3. 计算p值
fastspar_pvalues \
  --otu_table otu_table.tsv \    # 原始OTU表
  --correlation sparcc_cor.tsv \ # 原始相关性
  --prefix bootstrap_correlations/perm \  # 置换相关性
  --permutations 1000 \          # 置换次数
  --outfile sparcc_pvalues.tsv   # p值输出

# SPIEC-EASI分析
# 安装
# devtools::install_github("zdk123/SpiecEasi")
library(SpiecEasi)    # SPIEC-EASI包
library(phyloseq)     # 数据管理
library(igraph)       # 网络分析

# 1. 运行SPIEC-EASI
se_result <- spiec.easi(
  ps,                                    # phyloseq对象
  method = "mb",                         # 方法：mb(Meinshausen-Bühlmann)
  # method = "glasso" 或 "mb"
  # mb: 更稀疏，推荐
  # glasso: 更密集
  lambda.min.ratio = 1e-2,               # 正则化参数范围
  nlambda = 20,                          # lambda值数量
  pulsar.params = list(
    rep.num = 50,                        # 子采样次数
    ncores = 8                           # 并行核数
  )
)

# 2. 提取网络
# 获取邻接矩阵
adj_matrix <- getRefit(se_result)        # 稀疏邻接矩阵
# 获取加权矩阵（包含关联强度）
opt_cov <- getOptCov(se_result)          # 最优协方差矩阵
weight_matrix <- as.matrix(
  symBeta(getOptBeta(se_result), mode = "maxabs")
)

# 3. 构建igraph网络
net <- adj2igraph(adj_matrix,
                   vertex.attr = list(
                     name = taxa_names(ps)  # 物种名
                   ))

# 添加分类学信息
V(net)$phylum <- as.character(
  tax_table(ps)[V(net)$name, "Phylum"])  # 门级分类

# 4. 网络分析
# 基本统计
cat("节点数:", vcount(net), "\n")         # 节点数
cat("边数:", ecount(net), "\n")           # 边数
cat("平均度:", mean(degree(net)), "\n")   # 平均度
cat("密度:", graph.density(net), "\n")    # 网络密度
cat("聚类系数:", transitivity(net, type = "global"), "\n")
cat("模块性:", modularity(cluster_louvain(net)), "\n")

# 度分布
deg <- degree(net)                       # 各节点的度
hist(deg, breaks = 30, main = "度分布",
     xlab = "度", ylab = "频率")

# 关键节点（hub）
hub_scores <- hub_score(net)$vector      # Hub分数
top_hubs <- sort(hub_scores, decreasing = TRUE)[1:10]
cat("\nTop 10 Hub物种:\n")
print(top_hubs)

# 5. 模块检测（社区发现）
modules <- cluster_louvain(net)          # Louvain算法
cat("\n检测到", length(modules), "个模块\n")
V(net)$module <- membership(modules)     # 模块标签

# 模块内物种组成
for (m in 1:length(modules)) {
  members <- V(net)$name[V(net)$module == m]
  phyla <- V(net)$phylum[V(net)$module == m]
  cat(sprintf("模块%d (%d个物种): %s\n",
              m, length(members),
              paste(table(phyla), collapse = ", ")))
}

# 6. 可视化
library(ggraph)       # ggplot风格的网络图
library(tidygraph)    # tidy网络数据

tg <- as_tbl_graph(net) %>%              # 转为tidy格式
  mutate(centrality = centrality_degree())  # 添加中心性

ggraph(tg, layout = "fr") +             # Fruchterman-Reingold布局
  geom_edge_link(alpha = 0.3) +          # 边（半透明）
  geom_node_point(aes(color = phylum,    # 按门着色
                      size = centrality)) +  # 大小=中心性
  scale_size(range = c(1, 8)) +          # 节点大小范围
  labs(title = "SPIEC-EASI微生物网络") +
  theme_graph()                          # 网络主题
ggsave("spieceasi_network.png", width = 12, height = 10, dpi = 150)

# FlashWeave安装（Julia语言）
# julia -e 'import Pkg; Pkg.add("FlashWeave")'

# Julia脚本运行FlashWeave
julia << 'EOF'
using FlashWeave  # 导入FlashWeave

# 读取数据
data_path = "otu_table.tsv"  # OTU表路径
meta_path = "metadata.tsv"  # 元数据路径

# 运行FlashWeave
netw_results = learn_network(
    data_path,
    meta_path,
    sensitive = true,              # 高灵敏度模式
    heterogeneous = false,         # 同质数据
    max_k = 3,                     # 最大条件集大小
    n_obs_min = 10,                # 最小观测数
    alpha = 0.01                   # 显著性水平
)

# 保存网络
save_network("flashweave_network.edgelist", netw_results)  # 边列表
save_network("flashweave_network.gml", netw_results)       # GML格式
EOF

# 网络拓扑分析和组间比较
import networkx as nx  # 网络分析
import numpy as np     # 数值计算
import pandas as pd    # 数据处理
import matplotlib.pyplot as plt  # 绑图

# 1. 从相关性矩阵构建网络
def build_network(cor_matrix, pval_matrix, cor_threshold=0.3,
                   pval_threshold=0.05):
    """从相关性矩阵构建网络"""
    G = nx.Graph()                       # 创建无向图
    n = len(cor_matrix)                  # 物种数
    species = cor_matrix.index.tolist()  # 物种名

    G.add_nodes_from(species)            # 添加节点

    for i in range(n):
        for j in range(i+1, n):
            cor = cor_matrix.iloc[i, j]  # 相关系数
            pval = pval_matrix.iloc[i, j]  # p值
            if abs(cor) >= cor_threshold and pval < pval_threshold:
                G.add_edge(species[i], species[j],
                           weight=cor,          # 权重=相关系数
                           sign="positive" if cor > 0 else "negative")
    return G

# 2. 计算网络拓扑指标
def network_topology(G):
    """计算全局和局部网络拓扑指标"""
    if len(G.nodes()) == 0:
        return {}
    # 全局指标
    topology = {
        "nodes": G.number_of_nodes(),    # 节点数
        "edges": G.number_of_edges(),    # 边数
        "density": nx.density(G),        # 密度
        "avg_degree": np.mean([d for _, d in G.degree()]),  # 平均度
        "avg_clustering": nx.average_clustering(G),  # 平均聚类系数
        "transitivity": nx.transitivity(G),  # 传递性
    }
    # 连通分量
    components = list(nx.connected_components(G))
    topology["n_components"] = len(components)  # 连通分量数
    largest_cc = max(components, key=len)
    topology["largest_cc_size"] = len(largest_cc)  # 最大连通分量

    # 正/负边比例
    pos_edges = sum(1 for _, _, d in G.edges(data=True) if d.get("sign") == "positive")
    neg_edges = sum(1 for _, _, d in G.edges(data=True) if d.get("sign") == "negative")
    topology["positive_edges"] = pos_edges
    topology["negative_edges"] = neg_edges
    topology["pos_neg_ratio"] = pos_edges / max(neg_edges, 1)

    # 模块性
    from community import best_partition  # python-louvain
    if len(G.edges()) > 0:
        partition = best_partition(G)    # Louvain社区发现
        topology["modularity"] = nx.algorithms.community.quality.modularity(
            G, [{n for n, m in partition.items() if m == c}
                for c in set(partition.values())]
        )
        topology["n_modules"] = len(set(partition.values()))

    return topology

# 3. 比较两组网络
def compare_networks(G1, G2, name1="Group1", name2="Group2"):
    """比较两组的网络拓扑差异"""
    topo1 = network_topology(G1)         # 组1网络指标
    topo2 = network_topology(G2)         # 组2网络指标

    comparison = pd.DataFrame({
        name1: topo1,
        name2: topo2
    })
    print("网络拓扑比较:")
    print(comparison.to_string())

    # 可视化比较
    metrics = ["nodes", "edges", "avg_degree", "density",
               "avg_clustering", "modularity"]
    fig, axes = plt.subplots(2, 3, figsize=(15, 8))
    for idx, metric in enumerate(metrics):
        ax = axes[idx // 3, idx % 3]
        vals = [topo1.get(metric, 0), topo2.get(metric, 0)]
        ax.bar([name1, name2], vals,
               color=["#FF6B6B", "#4ECDC4"])  # 柱状图
        ax.set_title(metric)
        ax.set_ylabel("值")
    plt.suptitle("网络拓扑比较", fontsize=14)
    plt.tight_layout()
    plt.savefig("network_comparison.png", dpi=150)

    return comparison

# 4. 关键物种鉴定
def identify_keystone_species(G, top_n=10):
    """鉴定关键物种（基于多种中心性指标）"""
    centralities = pd.DataFrame({
        "degree": dict(G.degree()),                     # 度中心性
        "betweenness": nx.betweenness_centrality(G),    # 介数中心性
        "closeness": nx.closeness_centrality(G),        # 接近中心性
        "eigenvector": nx.eigenvector_centrality_numpy(G),  # 特征向量中心性
    })
    # 综合排名
    centralities["rank_sum"] = centralities.rank(ascending=False).sum(axis=1)
    centralities = centralities.sort_values("rank_sum")  # 按排名排序

    print(f"\nTop {top_n} 关键物种（Keystone species）:")
    print(centralities.head(top_n))
    return centralities

# 共丰度网络分析流程
1. 数据预处理：过滤低丰度物种（出现率>20%）
2. 选择相关性方法：
   - 组成偏差处理: SparCC/FastSpar
   - 直接关联推断: SPIEC-EASI (MB/glasso)
   - 大规模数据: FlashWeave
   - 纵向数据: SpeSpeNet (2025)
3. 置换检验评估显著性（≥1000次）
4. 构建网络（阈值过滤）
5. 拓扑分析（度、模块性、中心性）
6. 关键物种鉴定
7. 组间网络比较

# 方法选择
SparCC: 处理组成偏差，简单相关（可能含间接关联）
SPIEC-EASI: 推断直接关联（去间接效应），推荐
FlashWeave: 条件独立性检验，适合>1000个物种
CoNet: 多方法集成（取交集）

# 网络过滤标准
|r| > 0.3-0.6（根据数据调整）
p < 0.05（经FDR校正）
出现率 > 20-30%（去除稀有物种）

# 网络拓扑解读
高模块性(>0.4): 群落有明显的功能分组
正/负边比>1: 合作关系多于竞争
高Hub度: 该物种连接多个生态位

# 重要警告
共丰度 ≠ 互作（不能从网络直接推断生物学互作）
相关 ≠ 因果（需要实验验证）