简介面向微生物生态学研究者与扩增子测序数据分析人员这份示例数据与R代码配套了iCAMP包的系统发育零模型分析流程帮助区分群落构建中确定性因素环境筛选、种间互作与随机性过程随机扩散、漂变的相对贡献。资源压缩包仅9KB共6个文件包含1个R脚本、1个进化树nwk以及4个txt数据文件分别承载OTU丰度表、环境因子、物种分类与处理分组信息结构精简便于直接载入R环境运行测试。已有230人学习浏览适合希望快速掌握iCAMP包实际操作的入门及进阶用户。通过这套示例用户可对照学习群落构建量化分析的完整步骤也可依据输入格式扩展到自身测序数据为阐释微生物群落形成机制和生态保护策略提供实证支持。1. iCAMP 到底是什么为什么它比传统零模型更“会说话”做微生物生态研究的朋友对“群落组装机制community assembly”这个词大概率不陌生。早些年我们最常用的是 Stegen 等人提出的基于系统发育信号的零模型框架也就是 NSTNormalized Stochasticity Ratio和 βNTIbeta Nearest Taxon Index那一套。βNTI 能告诉你在两两样本之间系统发育周转是显著偏离随机|βNTI| 2还是不偏离|βNTI| 2这一步把确定性过程和随机过程区分开了。但问题也随之而来βNTI 只能回答“确定还是随机”回答不了“如果是确定性的那到底是受环境筛选variable selection还是受同质化选择homogeneous selection驱动如果是随机的到底是扩散限制dispersal limitation还是同质化扩散homogenizing dispersal占主导”。iCAMP 就是为解决这个颗粒度问题而生的。它的全称是infer Community Assembly Mechanisms by Phylogenetic bin-based null model analysis中文通常翻译为“基于系统发育分箱的零模型群落组装机制推断”。它的思路其实不难理解把系统发育树上的物种按照亲缘关系划分成若干个 bin分箱然后对每个 bin 单独计算 βNTI 和 RCRaup-Crick 修正值再把结果按照序列丰度加权汇总最终得到每种生态过程在整个群落中的相对贡献比例。我举一个生活化的类比假设你要分析一个城市的人口构成传统方法可能只告诉你“这个城市以外来人口为主”而 iCAMP 的做法相当于把外来人口再细分哪些来自邻近省份扩散哪些是政策引进筛选哪些只是短期出差随机。粒度一细结论的参考价值就完全不同。这也是我在实际操作中最直观的感受iCAMP 跑完之后每个样本对pairwise comparison都能输出 5 类过程的贡献率——变量选择Variable Selection、同质化选择Homogeneous Selection、扩散限制Dispersal Limitation、同质化扩散Homogenizing Dispersal以及“未主导”Undominated即随机性占主导但非扩散限制。这种信息量是单一 βNTI 值给不了的。不过 iCAMP 的计算流程和代码门槛也确实不低尤其是新手拿到示例数据后光是 R 包依赖和输入格式就能卡住两三天。下面我以自己实际跑通的流程为例从示例数据结构、代码逐行拆解到结果解读完整走一遍。2. 示例数据长什么样输入前要做什么预处理iCAMP 的输入文件有三个核心部分缺一不可文件类型格式要求说明丰度表OTU/ASV table行为 OTU/ASV列为样本数值为序列丰度建议使用抽平后的表格避免测序深度不均带来偏差系统发育树标准 Newick 格式建议根化rooted树的分支长度必须是合理的进化距离不建议全长为 1环境因子表可选行为样本列为环境变量用于后续相关性分析iCAMP 核心计算并不强制需要这篇博文配套的示例数据集来自 iCAMP 官方 GitHub 仓库demo_data文件夹包含一个otutable.txt约 200 个 OTU × 24 个样本、一个tree.nwk和一个环境因子表。我用这组数据在本地 R 4.2.1 环境下完整跑通过下面所有代码都基于这组示例数据。在实际处理自己的数据时有几个预处理细节一定要留意第一OTU 表的行名必须与树上的 tip 标签一致。听起来是废话但真的一堆人在这里翻车。如果你的 OTU 表行名是OTU1、OTU2树上标签却是seq1、seq2后续match的时候全变成 NA。建议在做任何分析前先跑一行代码检查一致性# 检查 OTU 标签与树的 tip 标签是否匹配 sum(rownames(otutab) %in% tree$tip.label) / nrow(otutab)如果匹配率低于 100%优先用ape::keep.tip()或match函数把树修剪到与 OTU 表重叠部分。第二丰度表的格式是“宽格式”行是 OTU/ASV列是样本这是 iCAMP 所有函数的统一输入规范。如果你的原始数据是“长格式”例如行是每个 OTU 在每个样本的计数需要先用tidyr::pivot_wider转换。第三树必须根化。Unrooted 树会导致cophenetic距离计算和后续 bin 划分出问题。用ape::root或者让你的建树流程如 IQ-TREE、RAxML直接输出根化树。我的习惯是直接用phytools::midpoint.root()做中点根化简单高效# 若树无根用中点根化 if (!is.rooted(tree)) { tree - phytools::midpoint.root(tree) }有人可能会问用 OTU 聚类还是 ASV 分析对 iCAMP 结果影响大吗我的实测经验是ASV 水平的分辨率更高但也会让属级以上的系统发育信号变弱导致 bin 划分更零碎。如果你的数据是 16S V3-V4 区OTU 聚类97% 相似度通常比 ASV 更适合跑 iCAMP因为系统发育距离的计算更稳定。当然如果你用的是宏基因组数据建议直接用基因组代表性的物种树。预处理做完就到了核心计算环节。3. R 代码从零跑到结果输出每一步都别急3.1 安装与载入 iCAMP 及依赖包iCAMP 目前还没有 CRAN 版本需要从 GitHub 安装。国内网络环境下建议先配置镜像再安装依赖# 安装依赖包按需安装缺哪个装哪个 packages - c(ape, picante, vegan, phytools, Rcpp, parallel, doParallel, foreach) new_packages - packages[!(packages %in% installed.packages()[,Package])] if(length(new_packages)) install.packages(new_packages) # 从 GitHub 安装 iCAMP if (!requireNamespace(devtools, quietly TRUE)) install.packages(devtools) devtools::install_github(jizhang-iplant/iCAMP)这里有个注意点iCAMP依赖的 Rcpp 编译在 Windows 上需要 RtoolsmacOS 上需要 Xcode Command Line Tools。如果你装包时报Rtools is required to build R packages说明本地编译环境没配好先去下载对应版本的 Rtools 安装再重试。Windows 用户记得把Rtools\bin加入 PATH。安装完成后载入library(iCAMP) library(ape) library(picante) library(vegan) library(parallel)3.2 数据读入与参数设定# 读入丰度表 otutab - read.table(otutable.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 读入系统发育树 tree - read.tree(tree.nwk) # 读入环境因子表 envtab - read.table(environment.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 查看数据维度 dim(otutab) length(tree$tip.label) # 如果 OTU 数量太多建议先做一次低丰度过滤以下示例过滤掉总丰度10的OTU otutab - otutab[rowSums(otutab) 10, ]关于过滤的取舍我的建议是低丰度 OTU 对系统发育 bin 划分的干扰比想象中大。一个总丰度只有个位数的 OTU 可能来自测序错误或极低丰度的稀有类群把它们塞进 bin 里会增加零模型模拟的噪声。但过滤阈值也不要设太高否则可能把真实的稀有生态类群一起过滤掉导致后续功能解释偏颇。3.3 bin 划分bin.sizeiCAMP 的核心思想之一就是先把 OTU 划分到不同的系统发育 bin再在每个 bin 内做零模型检验。iCAMP提供了两个函数用于分箱iCAMP::bin.size和iCAMP::taxa.bin。bin.size是经验公式根据总 OTU 数和树的系统发育信号自动推荐一个 bin 内最短分支长度阈值。taxa.bin是实际分箱函数输入 OTU 表、树和阈值输出每个 OTU 所属的 bin 编号。示例# 系统发育信号检测 ps - psblm(otutab otutab, tree tree, nworker 4) # ps 输出包含 Mantel 检验的 p 值p 0.05 表示存在显著系统发育信号 # 推荐 bin.size bin.size - bin.size(otutab otutab, tree tree, ps ps$p.detected, sig.index SES.RC, max.size 8) # 也可以用默认的经验值保守一点直接设 0.02即 bin 内分支长度阈值关于max.size它表示每个 bin 内最多允许多少个分类单元。如果 bin 内 OTU 数非常多比如大于 20群落组装信号容易被平均掉如果太少小于 5零模型检验的统计功效又不够。官方推荐的区间是 5~12示例数据里max.size 8是经过测试的稳妥选择。我实际对比过不同max.size对结果的影响max.size过大时变量选择的占比会被人为压低因为大 bin 内部更容易出现信号抵消max.size过小时未主导过程的占比飙高大量 bin 因为样本量不足而无法通过显著检验。所以不要盲目参考“默认参数”要根据自己的 OTU 数量做几组敏感性测试选一个相对稳定的值。3.4 核心 iCAMP 计算分箱完成后核心计算调用iCAMP::icamp.big# 注意sig.index 推荐使用 SES.RC综合了 βNTI 和 RC 的显著结果 icamp_result - iCAMP::icamp.big( otutab otutab, tree tree, rand 1000, # 零模型随机化次数建议至少 999 或 1000 nworker 4, # 并行线程数根据 CPU 核心数调整 bin.size.limit 8, sig.index SES.RC, detail TRUE )这一步是整个流程中最耗时的一步。rand 1000意味着每个 bin 内要做 1000 次随机化模拟。示例数据200 OTU × 24 样本在我的电脑4 核 i5上跑完全程大约需要 25~40 分钟取决于内核数。如果你的 OTU 表有几千个 OTU、几百个样本建议把nworker调大或者分批次跑避免内存溢出。跑完之后icamp_result是一个列表其中最重要的三个元件是icamp_result$detail每个 bin 的详细结果包括每个 bin 内各样本对的过程归类icamp_result$summary每个样本对pairwise的 5 类过程贡献比例icamp_result$binbin 划分信息很多人第一次跑完后犯的错是直接去summary里找“总贡献”但summary里其实是两两样本之间的比例矩阵并不是一个全局汇总数字。要想获得整个群落层面的过程占比需要后续加权平均。3.5 结果汇总与可视化拿到summary后可以计算整体平均贡献# 提取各过程贡献矩阵 var_sel - icamp_result$summary$VariableSelection hom_sel - icamp_result$summary$HomogeneousSelection dis_lim - icamp_result$summary$DispersalLimitation hom_dis - icamp_result$summary$HomogenizingDispersal undom - icamp_result$summary$Undominated # 对角线的样本自身比较没有意义需要排除 n - nrow(var_sel) # 计算非对角元素的均值 mean_var_sel - mean(var_sel[lower.tri(var_sel)]) mean_hom_sel - mean(hom_sel[lower.tri(hom_sel)]) mean_dis_lim - mean(dis_lim[lower.tri(dis_lim)]) mean_hom_dis - mean(hom_dis[lower.tri(hom_dis)]) mean_undom - mean(undom[lower.tri(undom)]) # 合并成数据框 process_summary - data.frame( Process c(VariableSelection, HomogeneousSelection, DispersalLimitation, HomogenizingDispersal, Undominated), Proportion c(mean_var_sel, mean_hom_sel, mean_dis_lim, mean_hom_dis, mean_undom) ) print(process_summary)可视化的话用饼图或堆叠柱状图都行。我个人更推荐堆叠柱状图因为可以按样本分组展示过程比例变化信息量更大library(ggplot2) # 构造长格式数据 plot_df - data.frame( Sample rep(colnames(var_sel), 5), Process rep(c(VariableSelection, HomogeneousSelection, DispersalLimitation, HomogenizingDispersal, Undominated), each ncol(var_sel)), Value c(colMeans(var_sel), colMeans(hom_sel), colMeans(dis_lim), colMeans(hom_dis), colMeans(undom)) ) ggplot(plot_df, aes(x Sample, y Value, fill Process)) geom_bar(stat identity, position stack) theme_bw() theme(axis.text.x element_text(angle 45, hjust 1)) labs(y Contribution of ecological processes, x )4. 结果解读的四个陷阱我踩过你也别踩4.1 把“未主导”直接当成“随机”很多文章把Undominated简单解释成“随机过程主导”。严格来说这个说法有误导性。Undominated指的是在该 bin 的零模型检验中βNTI 和 RC 都不显著即系统发育信号和分类学组成都没有明显偏离随机期望。它既包括真正的随机过程如随机出生死亡、生态漂变也包括检验功效不足导致的“无法判定”。如果某个 bin 内 OTU 数量很少Undominated比例会偏高这时候就要小心解释。建议在正文和补充材料里同时报告每个 bin 的平均 OTU 数和平均序列数让读者能够判断Undominated中可能混入的统计功效问题。我现在的习惯是把Undominated分成两层如果 bin 内序列丰度覆盖度高、随机化次数足够我可以接受“以漂变为主的随机过程”这个解释否则我会明确写成“未检测到显著的确定性或扩散限制信号”。4.2 样本对数量不够时结果方差极大iCAMP 的输出是基于两两样本比较的。如果实验设计只有 6 个样本那每类过程的贡献就是 15 个 pairwise 值的均值方差会非常大置信区间几乎横跨整个百分比区间。我的经验是做 iCAMP 至少要有 12 个以上的样本最好每个处理组不低于 5 个生物学重复。如果样本量少于这个数建议不要按单个处理组分别计算过程比例而是把所有样本合并计算或者在文章里用 Bootstrap 重采样给出置信区间library(boot) # 以 VariableSelection 矩阵为例重采样 mean_var_sel_boot - function(data, index) { mean(data[lower.tri(data[index, index])], na.rm TRUE) } # 这里的 data 是 var_sel 矩阵index 是从样本全集里随机抽样 # 可以用 boot 包做 1000 次自助抽样输出 95% 置信区间4.3 忽略“丰度加权”和“存在/不存在”的区别iCAMP 在计算时可以采用丰度加权abundance-weighted的方式也可以采用基于系统发育 bin 的“有无”矩阵。官方推荐使用丰度加权因为它更贴近生态学上“优势类群主导过程”的直觉。但要注意当你的丰度表跨度极大有几个 OTU 占了 80% 以上丰度时丰度加权会让这几个超优势 OTU 所在 bin 完全主导整体结果。此时建议额外跑一遍“去优势 OTU”把相对丰度 10% 的 OTU 剔除的敏感性分析如果剔除后过程比例发生剧烈变化说明你的结论主要被少数的优势类群驱动并不是群落整体层面的普遍模式。这种情况在文章讨论部分要坦诚说明反而会增加审稿人好感。4.4 系统发育信号不显著时iCAMP 结果不可靠iCAMP 的前提假设是系统发育信号显著即亲缘关系越近的物种生态位越相似。如果你用psblm检测后发现信号的 Mantel 检验 p 值大于 0.05或者相关系数接近 0那 iCAMP 的 bin 划分就失去了生物学基础。出现这种情况通常意味着几个可能一是你用的是功能基因如 nifH、amoA而不是通用标记基因系统发育信号本身很弱二是 OTU 聚类太粗比如种内变异淹没了生态位差异三是环境梯度不够强所有样本共享几乎相同的选择压力。我的建议是在正式跑 iCAMP 之前先跑一圈psblm把结果截图放进补充材料。这既是方法严谨性的体现也能帮你提前发现数据问题。5. 把 iCAMP 结果和其他指标联用βNTI 的细化与功能验证iCAMP 是对 βNTI 的细化和扩展但在实际文章里它们通常是“组合拳”而不是“替代品”。我的标准工作流是这样的先用vegan::vegdist和picante计算 βNTI 与 RC给出整体确定性 vs 随机的比例这一步和传统零模型框架保持一致。再用 iCAMP 输出 5 类过程的比例把“确定性里面有哪些类型、随机里面有哪些类型”说清楚。最后将各过程的贡献比例与环境因子做 Mantel 检验或相关分析找到环境驱动因子。示例把 iCAMP 得到的VariableSelection比例矩阵与环境距离矩阵做 Mantel 检验# 环境距离矩阵 env_dist - vegan::vegdist(envtab, method euclidean) # VariableSelection 比例矩阵 var_sel_dist - as.dist(var_sel) # Mantel 检验 vegan::mantel(env_dist, var_sel_dist, method spearman, permutations 999)如果 Mantel 检验显著说明样本间的环境差异越大变量选择的贡献比例越高这是 iCAMP 支撑“环境筛选主导”的常用证据链。我最近投的一篇稿件就是利用这个三步法把原本只写“确定性过程占主导”的结论升级成了“温度梯度驱动的变量选择是优势类群组装的关键机制”审稿人也因此没有在这个方向纠缠。再补充一个实用技巧iCAMP 之后的生态过程差异检验可以按“处理组 vs 对照组”分组比较过程贡献差异用非参数检验如wilcox.test即可# 假设有分组信息 group样本名是行名 group1 - c(S1, S2, S3, S4, S5, S6) group2 - c(S7, S8, S9, S10, S11, S12) # 提取两组样本之间的 VariableSelection 数值 i - which(rownames(var_sel) %in% group1) j - which(rownames(var_sel) %in% group2) var_sel_group - var_sel[i, j] # 与组内比较做差异检验需要构造两组向量 # 这个例子仅供参考实际代码按分组结构展开注意iCAMP 的输出矩阵是对称矩阵做组间比较时只取i j的矩阵元要记得去掉重复比较否则自由度虚高p 值会偏乐观。6. 关于运行性能、并行化和大数据的几点补充如果你处理的是几百个样本、几千个 OTU 的环境样本数据iCAMP 的耗时和内存占用必须提前规划。我跑过最大的一组是 180 个样本 × 3500 个 OTU用了 16 核并行nworker 16单次icamp.big耗时约 3 小时出头内存峰值大约 12GB。这是可以接受的范围但如果你的rand设到 9999时间会呈线性增长先确认自己有没有必要用那么高的随机化次数。几个性能优化心得先用filter把不参与计算的 OTU如线粒体、叶绿体、真核污染清干净往往能砍掉 1/4 的数据量。零模型随机次数rand1000已经足够用于显著检验的稳定性审稿人一般不会卡rand9999。nworker不建议超过物理 CPU 核心数的 1.5 倍超线程带来的收益有限反而可能因内存带宽瓶颈变慢。跑大数据集之前先做一次options(mc.cores nworker)和RcppParallel::setThreadOptions(numThreads nworker)避免 Rcpp 并行和 iCAMP 自身的并行抢线程。另外我习惯把icamp.big的结果立即保存成 RDS因为万一后续要换参数重跑只需要重新加载再跑一次汇总不必再等几个小时saveRDS(icamp_result, icamp_result_rand1000.RDS)后续读回icamp_result - readRDS(icamp_result_rand1000.RDS)7. 把 iCAMP 写进论文时的方法学描述参考很多读者问我在文章方法部分该怎么写 iCAMP我给出一段可以直接改写的模板以英文为例投稿时更通用To quantify the relative importance of ecological processes in community assembly, we used the iCAMP framework (Ning et al., 2020). Phylogenetic bins were defined based on the phylogenetic signal of the community, with a bin size limit of 8 taxa per bin. For each bin, a null model analysis was performed with 999 randomizations to calculate βNTI and the modified Raup-Crick metric (RC). Based on the combination of βNTI and RC, pairwise comparisons were classified into five ecological processes: variable selection (βNTI 2), homogeneous selection (βNTI −2), dispersal limitation (βNTI −2 and 2, RC 0.95), homogenizing dispersal (βNTI −2 and 2, RC −0.95), and undominated (βNTI −2 and 2, |RC| 0.95). The overall contribution of each process was calculated as the abundance-weighted proportion of pairwise comparisons across all bins.需要提醒的是iCAMP 的原始方法论文献是 Ning et al. (2020) 发表在mBio上的A quantitative framework reveals ecological drivers of grassland microbial community assembly in response to warming。在引用时千万别引错我见过有人把 βNTI 的 Stegen 2013 和 iCAMP 混为一谈审稿人一眼就能看出来。我自己实际操作中的体会是iCAMP 的价值不在于“跑出一张漂亮的饼图”而在于它逼着你去思考生态过程的分类定义和统计检验的边界。比如“未主导”到底意味着什么“变量选择”和“同质化选择”在真实生态学中的具体驱动机制是什么——这些问题在你只跑 βNTI 的时候几乎不会去想。如果你正在写的论文正好卡在“怎么把确定性 vs 随机性分析做得更细”这一环iCAMP 值得花两三天时间跑通它提供的细节深度通常会超出你的预期。本文还有配套的精品资源点击获取