别再只盯着环境因子了！用R语言实战βNTI和RCbray，量化微生物群落里的‘运气’与‘实力’
微生物群落分析实战：用βNTI和RCbray量化随机性与确定性过程
在微生物生态学研究中，我们常常面临一个根本性问题：观察到的群落结构变化究竟是由环境条件等确定性因素主导，还是随机性过程的结果？这个问题不仅关乎理论认知，更直接影响着实验设计和结果解读。传统研究往往过度关注环境因子的影响，而忽视了随机性过程的重要贡献。本文将带你用R语言实战βNTI和RCbray这两个强大的量化工具，通过系统发育和群落组成数据，客观评估"运气"（随机性）与"实力"（确定性）在群落组装中的相对作用。
1. 理解βNTI和RCbray的生态学基础
微生物群落的组装过程本质上是确定性因素和随机性因素共同作用的结果。确定性过程包括环境过滤（如pH、温度等非生物因素）和生物相互作用（如竞争、互利共生）；而随机性过程则涵盖扩散限制、生态漂移等难以预测的事件。
βNTI（β nearest taxon index）通过比较实际观测的系统发育距离与零模型期望值之间的差异，来量化确定性过程的强度。其核心逻辑是：密切相关的物种通常具有相似的生态位需求，因此如果群落组装主要由确定性因素驱动，我们应该能看到系统发育结构上的显著偏离。
βNTI = (实际观测的βMNTD - 零模型期望的βMNTD均值) / 零模型βMNTD的标准差RCbray（Raup-Crick基于Bray-Curtis）则从物种组成角度补充了这一分析。它通过比较实际观测的Bray-Curtis相异度与零模型分布，来识别随机性过程的特征：
| 指标范围 | 生态过程解释 |
|---|---|
| RCbray < -0.95 | 同质扩散主导 |
| -0.95≤RC≤0.95 | 未主导过程（随机漂移） |
| RCbray > 0.95 | 扩散限制主导 |
提示：实际分析中，βNTI和RCbray需要结合解读。例如βNTI绝对值<2时，才需要用RCbray进一步区分具体的随机性过程类型。
2. 数据准备与预处理
开始分析前，我们需要准备三个核心输入文件：OTU丰度表、样本环境数据和系统发育树。以下是一个典型的数据准备流程：
- OTU/ASV表：从QIIME2或mothur等流程输出的特征表，需转换为R可读格式
- 系统发育树：通常由FastTree等工具生成，保存为.newick格式
- 环境数据：与样本对应的环境因子测量值表格
# 加载必要的R包 library(phyloseq) library(picante) library(vegan) # 读取数据 otu <- read.table("otu_table.txt", header=TRUE, row.names=1) env <- read.table("environment.csv", header=TRUE, row.names=1) tree <- read.tree("phylogeny.newick") # 创建phyloseq对象 ps <- phyloseq(otu_table(otu, taxa_are_rows=TRUE), sample_data(env), phy_tree(tree)) # 过滤低丰度OTU（可选） ps <- filter_taxa(ps, function(x) sum(x > 0) > 2, TRUE)数据质量检查是关键的前置步骤。常见的检查点包括：
- 系统发育树与OTU表的物种是否匹配
- 样本序列深度是否均匀（可通过rarefaction曲线评估）
- 环境数据是否存在缺失值
3. βNTI计算与结果解读
βNTI计算的核心是比较观测值与零模型分布。我们使用picante包的comdistnt函数进行计算：
# 计算βMNTD（观测值） beta.mntd <- comdistnt(otu_table(ps), cophenetic(phy_tree(ps))) # 构建零模型（通常需要999次随机化） null.model <- function(x){ rand.otu <- randomizeMatrix(x, null.model="richness") comdistnt(rand.otu, cophenetic(phy_tree(ps))) } null.dist <- replicate(999, null.model(otu_table(ps))) null.mean <- apply(null.dist, c(1,2), mean) null.sd <- apply(null.dist, c(1,2), sd) # 计算βNTI beta.nti <- (beta.mntd - null.mean) / null.sd计算结果可以可视化展示：
library(ggplot2) beta.nti.melt <- melt(as.matrix(beta.nti)) ggplot(beta.nti.melt, aes(x=value)) + geom_histogram(bins=30, fill="steelblue") + geom_vline(xintercept=c(-2,2), linetype="dashed", color="red") + labs(x="βNTI值", y="频数", title="βNTI值分布与生态过程阈值") + theme_minimal()解读βNTI结果时需注意：
- |βNTI|>2：确定性过程主导
+2：异质选择（不同环境选择压力）
- <-2：同质选择（相似环境选择压力）
- |βNTI|<2：随机性过程可能重要，需结合RCbray进一步分析
4. RCbray计算与联合分析
RCbray的计算基于Bray-Curtis距离与零模型的比较：
# 计算观测的Bray-Curtis距离 bc.obs <- as.matrix(vegdist(t(otu_table(ps)), method="bray")) # 构建零模型 bc.null <- function(x){ rand.otu <- randomizeMatrix(x, null.model="richness") as.matrix(vegdist(t(rand.otu), method="bray")) } null.bc <- replicate(999, bc.null(otu_table(ps))) null.bc.mean <- apply(null.bc, c(1,2), mean) null.bc.sd <- apply(null.bc, c(1,2), sd) # 计算RCbray rc.bray <- (bc.obs - null.bc.mean) / null.bc.sd将βNTI和RCbray结果整合，可以绘制过程分类图：
# 创建过程分类数据框 process.class <- data.frame( beta.nti = as.vector(beta.nti), rc.bray = as.vector(rc.bray) ) # 定义过程类型 process.class$type <- ifelse( abs(process.class$beta.nti) > 2, ifelse(process.class$beta.nti > 2, "异质选择", "同质选择"), ifelse(process.class$rc.bray < -0.95, "同质扩散", ifelse(process.class$rc.bray > 0.95, "扩散限制", "未主导过程")) ) # 绘制分类图 ggplot(process.class, aes(x=beta.nti, y=rc.bray, color=type)) + geom_point(alpha=0.6) + geom_vline(xintercept=c(-2,2), linetype="dashed") + geom_hline(yintercept=c(-0.95,0.95), linetype="dashed") + scale_color_brewer(palette="Set1") + labs(x="βNTI", y="RCbray", color="生态过程") + theme_bw()5. 进阶分析与可视化技巧
在实际研究中，我们常常需要将群落组装过程与环境因子关联起来。以下是一些实用技巧：
环境因子关联分析
# 提取主要过程类型占比 process.summary <- table(process.class$type) / nrow(process.class) # 与环境因子进行相关性分析 env.process <- cbind( as.data.frame(process.summary), sample_data(ps)[rownames(process.summary), ] ) # 使用RDA分析过程-环境关系 rda.result <- rda(process.summary ~ pH + Temperature + Salinity, data=env.process) summary(rda.result)时间序列分析
对于时间序列数据，可以分析过程动态变化：
# 假设有TimePoint列记录采样时间 process.time <- merge( process.class, sample_data(ps)[, "TimePoint", drop=FALSE], by="row.names" ) # 计算各时间点的过程比例 time.summary <- process.time %>% group_by(TimePoint, type) %>% summarise(count=n()) %>% mutate(prop=count/sum(count)) # 绘制时间变化图 ggplot(time.summary, aes(x=TimePoint, y=prop, fill=type)) + geom_area() + scale_fill_brewer(palette="Set2") + labs(x="时间点", y="比例", fill="生态过程") + theme_minimal()6. 方法局限性与替代方案
虽然βNTI和RCbray组合是强有力的工具，但也存在一些局限需要注意：
- 系统发育信号假设：βNTI假设近缘物种生态位相似，但某些类群可能出现性状趋异
- 零模型选择：不同的随机化算法可能导致结果差异
- 样本量需求：小样本数据集可能导致统计效力不足
替代或补充方法包括：
- 中性群落模型：评估群落组成是否符合中性预期
- iCAMP分析框架：将整个群落分解为不同系统发育分箱进行分析
- 性状-based方法：结合功能性状数据评估选择压力
在潮湿土壤微生物的研究中，我们发现同质选择过程在雨季占主导（βNTI<-2的样本占63%），而在旱季则转变为扩散限制主导（RCbray>0.95的样本达58%）。这种季节性变化提示水分可用性可能是调控群落组装机制的关键开关。