跳转至

微生物组网络分析:SpiecEasi/SPRING/FlashWeave

一句话概述

微生物网络分析就是找出"谁和谁经常一起出现"或"谁和谁互相排斥"——把微生物之间的共现关系画成一张关系网,发现群落中的关键物种和互作模式。


核心知识点表格

知识点说明
共现网络基于微生物丰度相关性构建的网络,节点=物种,边=相关关系
直接关联 vs 间接关联传统相关分析无法区分直接互作和间接关联(被第三者调控)
SpiecEasi基于图模型的条件独立网络推断,推断直接关联
FlashWeave基于概率图模型,处理大规模异质数据集最快
SPRING半参数方法,处理组成性和零膨胀
组成性问题相对丰度数据不能直接算相关,需要专门处理

白话解释

想象一个班级的社交关系: - 相关分析(Spearman):小明和小红总是一起出现 → "他们相关" - 问题:可能只是因为他们都和小李是朋友,小李在哪他们就在哪 - 条件独立分析(SpiecEasi):在控制了小李的影响后,小明和小红还相关吗?→ 推断"直接关系"

这就是SpiecEasi等方法比简单相关分析更好的地方——它们试图找出直接互作而非间接关联


SpiecEasi

白话解释: SpiecEasi(SParse InversE Covariance Estimation for Ecological ASsociation Inference)先用CLR变换处理组成性数据,再用稀疏图模型推断条件独立网络。就像在一群人中找出"真正的朋友关系"而非"朋友的朋友"。

# 安装SpiecEasi
# install.packages("devtools")
devtools::install_github("zdk123/SpiecEasi")  # 从GitHub安装

library(SpiecEasi)  # 加载SpiecEasi
library(phyloseq)   # 加载phyloseq

# 准备数据(需要phyloseq对象)
# 过滤低丰度物种(建议保留在>20%样本中出现的物种)
ps_filtered <- filter_taxa(ps, function(x) sum(x > 0) > 0.2 * nsamples(ps), TRUE)

# ====== 方法1:Meinshausen-Bühlmann (MB) 方法 ======
# 特点:速度快,适合大数据集
se_mb <- spiec.easi(
    ps_filtered,                     # 输入phyloseq对象
    method = "mb",                   # 使用MB方法
    lambda.min.ratio = 1e-2,         # 正则化参数范围
    nlambda = 20,                    # 正则化参数个数
    pulsar.params = list(
        rep.num = 50,                # StARS重复次数
        ncores = 8                   # 并行核心数
    )
)

# ====== 方法2:Graphical Lasso (glasso) 方法 ======
# 特点:估计偏相关系数,可以得到关联强度
se_glasso <- spiec.easi(
    ps_filtered,
    method = "glasso",               # 使用glasso方法
    lambda.min.ratio = 1e-2,
    nlambda = 20,
    pulsar.params = list(rep.num = 50, ncores = 8)
)

# 提取邻接矩阵
adj_matrix <- getRefit(se_mb)        # 获取稀疏邻接矩阵(0/1)

# 转换为igraph对象进行可视化和分析
library(igraph)
ig <- adj2igraph(adj_matrix,
    vertex.attr = list(name = taxa_names(ps_filtered)))  # 创建igraph对象

# 计算网络拓扑属性
cat("节点数:", vcount(ig), "\n")           # 节点数
cat("边数:", ecount(ig), "\n")             # 边数
cat("网络密度:", edge_density(ig), "\n")    # 网络密度
cat("平均度:", mean(degree(ig)), "\n")      # 平均度
cat("聚类系数:", transitivity(ig), "\n")    # 聚类系数
cat("模块度:", modularity(cluster_louvain(ig)), "\n")  # 模块度

# 可视化网络
plot(ig,
    vertex.size = degree(ig) * 2,     # 节点大小按度值
    vertex.label = NA,                 # 不显示标签(太多会乱)
    edge.width = 0.5,                  # 边宽度
    layout = layout_with_fr(ig),       # 力导向布局
    main = "微生物共现网络"            # 标题
)

# 识别关键节点(Hub物种)
hub_scores <- hub_score(ig)$vector     # 计算hub得分
top_hubs <- sort(hub_scores, decreasing = TRUE)[1:10]  # Top10 hub物种
cat("关键物种(Hub):\n")
print(top_hubs)

FlashWeave

白话解释: FlashWeave用Julia语言编写,特点是特别快、能处理大规模异质数据集,还能同时把元数据(如pH、温度、分组)纳入网络中。

# FlashWeave用Julia运行
# 安装Julia和FlashWeave
# 下载Julia: https://julialang.org/downloads/
# 在Julia REPL中安装FlashWeave
julia -e 'using Pkg; Pkg.add("FlashWeave")'  # 安装FlashWeave包

# 准备输入文件
# 1. OTU表(TSV格式,行=物种,列=样本)
# 2. 元数据文件(可选)

# 运行FlashWeave(Julia脚本)
julia << 'JLEOF'
using FlashWeave  # 加载FlashWeave包

# 从文件推断网络
netw_results = learn_network(
    "otu_table.tsv",              # 输入OTU表
    "metadata.tsv",               # 输入元数据(可选)
    sensitive = true,             # 使用FlashWeave-S(更敏感)
    heterogeneous = true,         # 异质数据模式
    n_obs_min = 10,               # 最少观测数
    normalize = true              # 自动标准化
)

# 保存网络结果
save_network(
    "flashweave_network.edgelist",  # 输出边列表
    netw_results                    # 网络结果对象
)

# 保存为GraphML格式(可用Cytoscape打开)
save_network(
    "flashweave_network.gml",
    netw_results
)
JLEOF

在R中使用FlashWeave结果

# 读取FlashWeave的边列表
edges <- read.delim("flashweave_network.edgelist", header = FALSE)
colnames(edges) <- c("source", "target", "weight")  # 命名列

# 创建igraph对象
library(igraph)
ig_fw <- graph_from_data_frame(edges, directed = FALSE)  # 无向图

# 正边和负边分开
positive_edges <- E(ig_fw)[E(ig_fw)$weight > 0]  # 正相关边
negative_edges <- E(ig_fw)[E(ig_fw)$weight < 0]  # 负相关边
cat("正相关边数:", length(positive_edges), "\n")
cat("负相关边数:", length(negative_edges), "\n")

SPRING

白话解释: SPRING用半参数秩相关来处理组成性和零膨胀问题,介于简单相关和复杂图模型之间。

# 安装SPRING(在SpiecEasi包内)
# SPRING是SpiecEasi包的一部分,已包含在内

library(SpiecEasi)

# 运行SPRING
spring_result <- spiec.easi(
    ps_filtered,
    method = "mb",                   # 基础方法
    sel.criterion = "stars",         # 模型选择标准
    lambda.min.ratio = 1e-2,
    nlambda = 20,
    pulsar.params = list(rep.num = 50, ncores = 8)
)

网络分析与可视化

使用microeco包(一站式)

# microeco包整合了多种网络推断方法
# install.packages("microeco")
library(microeco)

# 创建microtable对象
dataset <- microtable$new(
    otu_table = otu_table,       # OTU丰度表
    sample_table = metadata,      # 样本元数据
    tax_table = taxonomy           # 分类信息
)

# 创建网络对象
network <- trans_network$new(
    dataset = dataset,
    cor_method = "sparcc",         # 相关方法:sparcc/spearman/spieceasi
    filter_thres = 0.001,          # 过滤阈值
    use_NetCoMi_env = FALSE        # 不使用NetCoMi环境
)

# 计算相关网络
network$cal_network(
    COR_p_thres = 0.01,           # p值阈值
    COR_cut = 0.6                  # 相关系数阈值
)

# 计算网络属性
network$cal_network_attr()         # 计算网络拓扑属性
network$cal_node_type()            # 计算节点类型(模块内/间连接者)

# 可视化
network$plot_network(
    node_size_range = c(2, 8)      # 节点大小范围
)

# Zi-Pi图(识别关键石物种)
network$plot_taxa_roles(
    roles_color_map = c("Peripherals" = "grey",
                         "Connectors" = "blue",
                         "Module hubs" = "red",
                         "Network hubs" = "darkred")
)

关键石物种(Keystone Species)识别

# 关键石物种特征:
# 1. 高度中心性(连接很多其他物种)
# 2. 高介数中心性(是网络中的"桥梁")
# 3. 在Zi-Pi图中属于module hub或network hub

# 计算多种中心性指标
degree_centrality <- degree(ig, normalized = TRUE)     # 度中心性
betweenness_cent <- betweenness(ig, normalized = TRUE)  # 介数中心性
closeness_cent <- closeness(ig, normalized = TRUE)      # 接近中心性
eigenvector_cent <- eigen_centrality(ig)$vector         # 特征向量中心性

# 综合排名
centrality_df <- data.frame(
    taxon = V(ig)$name,
    degree = degree_centrality,
    betweenness = betweenness_cent,
    closeness = closeness_cent,
    eigenvector = eigenvector_cent
)

# 按度中心性排序
top_keystone <- centrality_df[order(-centrality_df$degree), ][1:10, ]
print(top_keystone)  # 打印Top10关键石物种

常见报错与解决

报错原因解决方案
SpiecEasi pulsar收敛失败参数范围不对或数据太稀疏增大lambda.min.ratio,过滤更多低频物种
FlashWeave Julia not foundJulia未安装或不在PATH安装Julia并添加到PATH
网络边数为0阈值设太严降低p值阈值或相关系数阈值
内存不足物种太多过滤保留前200-500个物种
网络太密看不清阈值太松提高阈值或用模块化布局
igraph版本冲突R包版本不兼容更新所有相关包

速查表

# 方法选择
简单快速 → SparCC相关(处理组成性)
推断直接关联 → SpiecEasi (MB或glasso)
大数据+元数据 → FlashWeave
一站式分析 → microeco包

# 网络拓扑指标
节点数(Nodes) — 物种数量
边数(Edges) — 关联数量
密度(Density) — 实际边数/可能边数
平均度(Mean degree) — 每个节点平均连接数
聚类系数(Clustering coefficient) — 局部聚集程度
模块度(Modularity) — 模块化程度
平均路径长度(Average path length) — 节点间平均距离

# 关键石物种识别
度中心性(Degree) — 连接数量多
介数中心性(Betweenness) — 网络桥梁位置
Zi-Pi图 — Module hub/Network hub

# 物种过滤建议
至少在20%样本中出现
保留Top 200-500个物种
去除单例和极低丰度物种

面试高频问题

Q1:为什么不能直接用Spearman/Pearson相关做微生物网络?

A:两个原因:①微生物数据是组成性数据(相对丰度),直接算相关会产生假相关;②Spearman/Pearson只能检测成对相关,无法区分直接关联和间接关联(A和B相关可能只是因为都与C相关)。

Q2:SpiecEasi如何解决组成性问题?

A:SpiecEasi先对数据做CLR(中心对数比)变换消除组成性的影响,然后用稀疏逆协方差估计(Graphical Lasso)或邻域选择(MB方法)来推断条件独立关系,得到的是直接关联网络。

Q3:什么是关键石物种(Keystone Species)?

A:关键石物种是在微生物网络中处于关键位置的物种,移除它会导致网络结构崩塌。通常具有高度中心性和高介数中心性。在Zi-Pi图中表现为Module Hub或Network Hub。

Q4:不同网络推断方法为什么结果差异很大?

A:2025年研究表明,不同方法基于不同的数学假设和模型,产生的网络结构确实差异很大。这是该领域的核心挑战。建议:使用多种方法比较,关注多种方法一致的连接;或使用共识网络方法取多个方法的交集。