微生物组网络推断¶

一句话概述¶

利用SpiecEasi、FlashWeave和CoNet等工具从微生物组丰度数据推断物种间的共现/互斥关系网络，揭示微生物群落中的生态互作模式，是微生物组生态学研究的高级分析方法。

核心知识点表格¶

知识点	关键内容	常用工具
组成数据问题	相对丰度的组成性导致虚假相关	CLR转化, ALR
稀疏网络推断	假设生态网络稀疏，正则化方法	SpiecEasi (MB/glasso)
条件独立检测	区分直接关联和间接关联	FlashWeave, SPIEC-EASI
集成方法	多指标投票确定稳健边	CoNet, CCREPE
网络拓扑分析	度分布/模块化/hub鉴定	igraph, NetworkX
关键种鉴定	高连接度/介数中心性的hub物种	NetShift, NetCoMi
差异网络	条件间网络结构比较	NetCoMi, SPIEC-EASI

各步骤详解¶

第一步：理解组成数据问题¶

白话解释： 16S或宏基因组数据给出的是相对丰度（所有物种占比加起来等于100%）。这意味着如果A增多，其他所有物种的比例都会相对下降，即使它们的绝对量没变。直接对相对丰度计算相关性会产生大量虚假的负相关（称为"组成偏差"）。因此不能简单用Pearson/Spearman相关做网络推断。

技术细节： - 组成数据(compositional data)：和为常数的向量 - 简单相关会导致虚假负相关（Aitchison, 1986） - 解决方案：CLR转化（中心对数比）、ALR转化（加性对数比）、组成数据专用模型 - 零值处理：微生物数据含大量零值，需要适当处理 - 稀疏性假设：真实生态网络中大部分物种对不直接互作

代码示例：

# 演示组成偏差
set.seed(42)
# 模拟绝对丰度（物种间独立）
absolute <- matrix(rpois(100*5, lambda=50), ncol=5)
# 转为相对丰度
relative <- absolute / rowSums(absolute)

# 计算相关性
cor_abs <- cor(absolute, method="spearman")  # 接近0（无关联）
cor_rel <- cor(relative, method="spearman")  # 出现虚假负相关!

cat("绝对丰度相关:\n"); print(round(cor_abs, 2))
cat("相对丰度相关:\n"); print(round(cor_rel, 2))
# 相对丰度中出现强烈的负相关，但实际物种是独立的！

# CLR转化
library(compositions)
clr_data <- clr(relative + 0.5)  # 加伪计数处理零值
cor_clr <- cor(clr_data, method="spearman")

第二步：SpiecEasi网络推断¶

白话解释： SpiecEasi是最主流的微生物组网络推断工具。核心思路：先用CLR转化处理组成数据，再用稀疏正则化方法（LASSO/graphical lasso）推断条件独立性网络。"条件独立"意味着它找的是物种间的直接关联——如果A和C的相关性完全由B介导，那么SpiecEasi不会在A和C之间画边。

技术细节： - MB (Meinshausen-Bühlmann) 方法：对每个节点做LASSO回归 - glasso (Graphical LASSO) 方法：估计稀疏精度矩阵 - StARS (Stability Approach to Regularization Selection)：选择正则化参数 - 输出无权重、无方向的条件独立性网络 - 适合样本数>物种数的情况（但通过正则化也可处理p>n的高维情况）

代码示例：

# 安装SpiecEasi
# devtools::install_github("zdk123/SpiecEasi")
library(SpiecEasi)
library(phyloseq)

# 准备数据（phyloseq对象或OTU矩阵）
# 假设已有phyloseq对象 ps
otu_table <- as(otu_table(ps), "matrix")

# 过滤低丰度OTU（保留在至少20%样本中出现的OTU）
prevalence <- colSums(otu_table > 0) / nrow(otu_table)
otu_filtered <- otu_table[, prevalence >= 0.2]

# 方法1: MB (Meinshausen-Bühlmann)
se_mb <- spiec.easi(otu_filtered,
    method = "mb",
    lambda.min.ratio = 1e-2,
    nlambda = 20,
    sel.criterion = "stars",
    pulsar.params = list(
        rep.num = 50,      # StARS稳定性检验重复次数（默认20，建议提高到50）
        thresh = 0.05      # 不稳定性阈值（默认0.05）
    ),
    verbose = TRUE)

# 方法2: glasso (Graphical LASSO)
se_gl <- spiec.easi(otu_filtered,
    method = "glasso",
    lambda.min.ratio = 1e-2,
    nlambda = 20,
    sel.criterion = "stars",
    pulsar.params = list(rep.num=50, thresh=0.05))

# 提取网络
adj_matrix <- getRefit(se_mb)    # 邻接矩阵
opt_beta <- as.matrix(symBeta(getOptBeta(se_mb)))  # 边权重

# 转为igraph对象
library(igraph)
net <- adj2igraph(adj_matrix,
    vertex.attr = list(name = colnames(otu_filtered)))

# 网络基本统计
cat("节点数:", vcount(net), "\n")
cat("边数:", ecount(net), "\n")
cat("网络密度:", edge_density(net), "\n")
cat("平均度:", mean(degree(net)), "\n")

第三步：FlashWeave网络推断¶

白话解释： FlashWeave（当前版本v0.19.2）是一个基于条件互信息检验的工具，能处理异质性数据（如同时考虑环境因子）。它使用局部学习策略逐步构建网络，比全局方法更适合大规模数据。支持稀疏和高维数据。

技术细节： - 基于条件互信息(CMI)检验 - 支持混合数据类型（计数+连续+二分类） - FlashWeave-S (sensitive): 低假阴性，适合探索 - FlashWeave-HE (heterogeneous): 可整合元数据 - Julia语言实现，速度快

代码示例：

# 安装FlashWeave (Julia)
julia -e 'using Pkg; Pkg.add("FlashWeave")'

# 准备输入文件
# OTU表: samples × OTUs (TSV格式，第一行为OTU名)
# 元数据（可选）: samples × metadata variables

# Julia脚本
using FlashWeave

# 基本用法
netw_results = learn_network(
    "otu_table.tsv",
    sensitive = true,       # 高灵敏度模式(FlashWeave-S)
    heterogeneous = false,  # 纯微生物数据
    max_k = 3,             # 最大条件集大小
    alpha = 0.01           # 显著性阈值
)

# 整合元数据的异质模式
netw_results = learn_network(
    "otu_table.tsv",
    "metadata.tsv",
    sensitive = true,
    heterogeneous = true    # 支持异质数据
)

# 保存网络
save_network("flashweave_network.edgelist", netw_results)
save_network("flashweave_network.gml", netw_results)

# 命令行运行
julia -e '
using FlashWeave
netw = learn_network("otu_table.tsv", sensitive=true, heterogeneous=false)
save_network("network.gml", netw)
'

第四步：CoNet集成方法¶

白话解释： CoNet的策略是"集百家之长"：同时计算多种相关/距离指标（Pearson、Spearman、互信息、Bray-Curtis、KLD），只有被多数方法一致支持的边才保留。这样能大幅降低单一方法的假阳性。

技术细节： - 集成多种相关性/距离度量 - 使用permutation + bootstrap进行统计检验 - ReBoot (Renormalized Bootstrap) 处理组成偏差 - Brown's方法合并多指标的p-value - 可在Cytoscape中作为插件使用

代码示例：

# CoNet通过Cytoscape插件使用
# 或命令行版本:

# 准备输入（物种×样本矩阵）
# 运行CoNet
java -jar CoNet.jar \
    --input otu_table.txt \
    --output conet_output/ \
    --methods pearson,spearman,mutual_info,bray_curtis \
    --threshold top_bottom_100 \
    --permutations 1000 \
    --bootstraps 1000 \
    --reboot \            # ReBoot组成校正
    --merge brown \       # Brown's方法合并p-value
    --fdr bh \
    --significance 0.05

# R中使用ccrepe包（组成数据校正的相关性）
# 注意：ccrepe已从Bioconductor移除，可从GitHub安装：
# devtools::install_github("biobakery/ccrepe")
library(ccrepe)

# 计算组成校正的相关性
ccrepe_results <- ccrepe(
    x = otu_filtered,
    sim.score = cor,
    sim.score.args = list(method="spearman"),
    min.subj = 20,
    iterations = 1000,
    subset.cols.x = NULL)

# 提取显著关联
sig_pairs <- which(ccrepe_results$p.values < 0.05 & 
                   abs(ccrepe_results$sim.score) > 0.3,
                   arr.ind = TRUE)

第五步：网络拓扑分析与关键种鉴定¶

白话解释： 得到网络后，需要分析其结构特征：哪些物种是"中心节点"（连接很多其他物种的hub），网络中有没有模块（一群物种紧密互作的子群），网络是否具有"小世界"特征等。

代码示例：

library(igraph)
library(NetCoMi)

# 网络拓扑分析
# 基本指标
degree_dist <- degree(net)
betweenness_centrality <- betweenness(net, normalized=TRUE)
closeness_centrality <- closeness(net, normalized=TRUE)
eigenvector_centrality <- eigen_centrality(net)$vector

# Hub鉴定（多指标）
hub_score <- hub_score(net)$vector

# 综合hub定义：度、介数、特征向量中心性都高的节点
node_stats <- data.frame(
    taxon = V(net)$name,
    degree = degree_dist,
    betweenness = betweenness_centrality,
    eigenvector = eigenvector_centrality,
    hub = hub_score
)

# 鉴定hub：各指标z-score > 1.96
node_stats$is_hub <- with(node_stats,
    scale(degree) > 1.96 | scale(betweenness) > 1.96)

# 模块检测
modules <- cluster_louvain(net)
modularity(modules)
membership(modules)

# 小世界性检验
# 小世界指数 sigma > 1
transitivity_real <- transitivity(net, type="global")
path_length_real <- mean_distance(net)

# 与随机网络比较
random_nets <- replicate(100, {
    rg <- erdos.renyi.game(vcount(net), ecount(net), type="gnm")
    c(transitivity(rg, type="global"), mean_distance(rg))
})

sigma <- (transitivity_real / mean(random_nets[1,])) / 
         (path_length_real / mean(random_nets[2,]))
cat("小世界指数 sigma:", sigma, "\n")

# 网络可视化
# 先用getOptBeta的权重赋给边（SpiecEasi默认输出无权重）
beta_matrix <- as.matrix(symBeta(getOptBeta(se_mb)))
edge_weights <- beta_matrix[igraph::get.edgelist(net, names=FALSE)]
E(net)$weight <- edge_weights

V(net)$color <- membership(modules)
V(net)$size <- sqrt(degree_dist) * 3
E(net)$color <- ifelse(E(net)$weight > 0, "green", "red")

pdf("microbial_network.pdf", width=12, height=12)
plot(net, 
    layout = layout_with_fr(net),
    vertex.label.cex = 0.6,
    edge.width = abs(E(net)$weight) * 2,
    main = "Microbial Co-occurrence Network")
dev.off()

第六步：差异网络比较¶

白话解释： 比较不同条件（如健康vs疾病）下微生物网络结构的差异：哪些物种的连接关系发生了变化？哪些模块重组了？哪些hub物种消失或出现了？

代码示例：

# 使用NetCoMi进行差异网络分析
library(NetCoMi)
library(phyloseq)

# 分组构建网络
net_compare <- netConstruct(
    data = ps,
    group = "condition",  # phyloseq中的分组变量
    measure = "spieceasi",
    measurePar = list(method="mb", 
                      pulsar.params=list(rep.num=50)),
    filtTax = "highestFreq",
    filtTaxPar = list(highestFreq = 100),
    zeroMethod = "pseudo",
    normMethod = "clr",
    verbose = 3)

# 网络分析
net_anal <- netAnalyze(net_compare,
    centrLCC = TRUE,
    avDissIgnoreInf = TRUE,
    sPathNorm = TRUE,
    clustMethod = "cluster_fast_greedy",
    hubPar = "eigenvector",
    hubQuant = 0.95)

# 可视化差异
plot(net_anal,
    sameLayout = TRUE,
    layoutGroup = "union",
    rmSingles = "inboth",
    nodeSize = "eigenvector",
    cexNodes = 1.5,
    nodeSizeSpread = 3,
    edgeTranspLow = 80,
    edgeTranspHigh = 50,
    groupNames = c("Healthy", "Disease"),
    hubBorderCol = "red")

# 统计检验网络差异
net_comp <- netCompare(net_anal,
    permTest = TRUE,
    nPerm = 1000,
    adjust = "BH",
    seed = 42)

summary(net_comp)

实战命令¶

#!/bin/bash
# === 微生物组网络推断流程 ===

# R脚本完整流程
Rscript << 'EOF'
library(SpiecEasi)
library(phyloseq)
library(igraph)

# 1. 数据预处理
ps <- readRDS("phyloseq_object.rds")
ps_filtered <- filter_taxa(ps, function(x) sum(x > 0) > 0.2*nsamples(ps), TRUE)
ps_filtered <- prune_taxa(names(sort(taxa_sums(ps_filtered), TRUE))[1:200], ps_filtered)

# 2. SpiecEasi网络推断
se_net <- spiec.easi(ps_filtered, method="mb",
    lambda.min.ratio=1e-2, nlambda=20,
    pulsar.params=list(rep.num=50, thresh=0.05))

# 3. 提取网络
adj <- getRefit(se_net)
ig <- adj2igraph(adj, vertex.attr=list(name=taxa_names(ps_filtered)))

# 4. 拓扑分析
cat("Nodes:", vcount(ig), "\n")
cat("Edges:", ecount(ig), "\n")
cat("Density:", edge_density(ig), "\n")
cat("Modularity:", modularity(cluster_louvain(ig)), "\n")

# 5. Hub鉴定
deg <- degree(ig)
hub_taxa <- names(deg[deg > quantile(deg, 0.95)])
cat("Hub taxa:", hub_taxa, "\n")

# 6. 导出
write_graph(ig, "network.graphml", format="graphml")
saveRDS(se_net, "spieceasi_result.rds")
EOF

面试常问点¶

Q1: 为什么不能直接对微生物丰度数据计算Pearson/Spearman相关来构建网络？ A: 两个原因：1）组成偏差——相对丰度数据的和为1，导致虚假负相关；2）间接关联——简单相关无法区分直接关联和由第三方介导的间接关联。需要使用条件独立性方法（如SpiecEasi）来推断直接互作。

Q2: SpiecEasi的MB和glasso方法有什么区别？ A: MB方法对每个节点单独做LASSO回归推断其邻居，最后取对称；glasso直接估计整体精度矩阵。MB更快、可并行化，但可能产生不对称的结果需要额外处理。glasso理论更优雅但对大网络较慢。一般推荐MB。

Q3: 共现网络中的edge代表什么生物学含义？ A: 正edge可能代表：互利共生、共同偏好相同环境、代谢互补。负edge可能代表：竞争、拮抗、生态位排斥。但注意：统计关联不等于因果关系，需要实验验证。

Q4: 网络推断需要多少样本量？ A: 经验法则：样本数至少是物种数的3-5倍。SpiecEasi在样本数少时（如n<p）通过正则化处理，但至少需要30-50个样本。FlashWeave论文建议至少50个样本。

Q5: 如何验证推断的网络是否可靠？ A: 1）StARS稳定性选择（SpiecEasi内置）；2）交叉验证：随机留出部分样本重新构建网络检查一致性；3）与已知互作文献比对；4）体外共培养实验验证；5）比较不同方法的一致边。

易错点¶

直接用相关系数：这是最常见错误。简单相关产生大量虚假关联，必须用专门处理组成数据的方法。
样本量不足：样本数远小于物种数时，网络推断结果不可靠。需要进行物种过滤或增加样本。
零值处理不当：微生物数据零值很多，不同方法对零值的处理策略不同。过多伪计数会扭曲结果。
过度解读edge含义：统计共现/互斥不能直接等同于生态互作。可能是共同环境偏好而非直接互作。
忽略批次效应：不同批次的样本如果混在一起分析，技术变异可能被误认为生物学关联。
网络可视化误导：force-directed layout可能将不连接的节点放在一起，造成视觉上的"模块"假象。

补充知识¶

网络推断方法选择指南¶

场景	推荐方法	原因
标准16S分析	SpiecEasi (MB)	处理组成偏差，成熟稳定
大规模数据(>1000样本)	FlashWeave	快速，可扩展
整合环境因子	FlashWeave-HE	支持异质数据
稳健性优先	CoNet集成	多方法投票
差异网络	NetCoMi	统计比较框架
纵向数据	LUPINE/GGM时间序列	考虑时间依赖

关键网络拓扑指标解读¶

度(Degree): 连接数 → 高度节点是hub/关键种
介数中心性(Betweenness): 桥梁作用 → 控制信息/物质流动
模块化(Modularity): 群落结构 → 功能群组/生态位
连通性(Connectivity): 网络脆弱性 → 去除hub后碎片化
正/负edge比: 群落稳定性 → 高比例负edge表示竞争主导

微生物组网络推断¶

一句话概述¶

核心知识点表格¶

各步骤详解¶

第一步：理解组成数据问题¶

第二步：SpiecEasi网络推断¶

第三步：FlashWeave网络推断¶

第四步：CoNet集成方法¶

第五步：网络拓扑分析与关键种鉴定¶

第六步：差异网络比较¶

实战命令¶

面试常问点¶

易错点¶

补充知识¶

网络推断方法选择指南¶

关键网络拓扑指标解读¶

📚 相关文章推荐