扩增子测序|中性模型（NCM）在R中快速复现之群落装配机制

1、引言

虽然群落装配过程，确定性和随机性研究已经进入了新的阶段，开发了许多新的分析方法，比如零模型。但是sloan在2006提出中性模型理论因其具有计算方便、简化了复杂生态过程等特点被研究者广泛认同和接受，目前不少nature系列期刊采用该工具分析群落构建。

1.1、NCM中关键参数及其意义说明

R² (拟合优度): 定义: 衡量中性模型与实际数据的匹配程度。

意义: R²值越高，说明实际群落的结构更符合中性模型，表示随机性对群落影响更大。

N (元群落规模): 定义: 群落中所有物种的总丰度。

意义: N值越大，群落规模越大，物种池越丰富。

m (迁移率): 定义: 物种在群落间迁移的能力。

意义: m值越高，物种迁移越容易，扩散受限越少；m值越低，则迁移受限越多。

Nm (迁移扩散量): 定义: 元群落规模与迁移率的乘积（Nm = N * m）。

意义: 衡量群落间物种扩散的总能力，高Nm值表示群落间物种分布更均匀。

随机性 vs 确定性: 通过模型拟合优度（R²）和迁移扩散量（Nm）的结果，可以确定随机性（即中性过程）与确定性（即竞争、适应等过程）在群落结构中的作用。

如果R²高且Nm高: 说明群落结构主要受随机过程影响。

如果R²低且Nm低: 说明群落结构受确定性过程影响较大。

1.2 文献美图展示

（Cheng et al., 2024, Nature Communications）

在水系生态领域，该研究使用了中性群落模型（NCM）来评估青藏高原水生态系统随机过程的作用。结果显示，大多数微生物群落的NCM拟合度较低（R² < 0.3），河流微生物群落的拟合度最低（R² = 0.021）。这表明，群落组装中的随机过程贡献有限，环境选择压力在塑造微生物群落多样性方面起主导作用。

（Liu et al., 2023, The ISME Journal）

在土壤系统，使用NCM确定决定性和随机过程在功能基因群落组装中的贡献，NS的R²（0.79）高于AS的R²（0.67），表明NS中的功能基因群落更符合中性模型。

（Zorea et al., 2024, Nature Communications）

人类肠道质粒探索，NCM指出质粒的扩散与中性模型拟合度较高（R² = 0.5，图3A），这表明质粒扩散主要是中性的。

2、示例数据与R代码

2.1、示例数据

准备数据，只需要准备高通量测序后获得的OTU（ASV表格）

2.2、NCM代码

# 清理工作环境中的所有对象
rm(list = ls())

# 安装并加载所需的包
if (!requireNamespace("devtools", quietly = TRUE)) {
  install.packages("devtools")
}
if (!requireNamespace("ggimage", quietly = TRUE)) {
  install.packages("ggimage")
}
if (!requireNamespace("plyr", quietly = TRUE)) {
  install.packages("plyr")
}

library(devtools)
library(ggplot2)
library(ggimage)
library(plyr)


# 加载 remotes 包（用于 install_github 函数）
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
# 检查并安装 MicEco 包
if (!requireNamespace("MicEco", quietly = TRUE)) {
  remotes::install_github("Russel88/MicEco")
}
library(MicEco)

# 读取数据
read_table <- read.csv(file = "read_table.csv", row.names = 1)

# 运行中性模型分析
res <- neutral.fit(t(read_table))

m <- res[[1]][1]
N <- res[[1]][4]
Nm<- N*m
r2 <- res[[1]][3]
out <- res[[2]]

# 处理数据
out$group <- with(out, ifelse(freq < Lower, "#509579",
                              ifelse(freq > Upper, "#cf9198", "#485970")))

# 绘制模型结果图
p1 <- ggplot(data = out) +
  geom_line(aes(x = log(p), y = freq.pred), size = 1.2, linetype = 1) +
  geom_line(aes(x = log(p), y = Lower), size = 1.2, linetype = 2) +
  geom_line(aes(x = log(p), y = Upper), size = 1.2, linetype = 2) +
  geom_point(aes(x = log(p), y = freq, color = group), size = 2) +
  xlab("log10(mean relative abundance)") +
  ylab("Occurrence frequency") +
  scale_colour_manual(values = c("#485970", "#cf9198", "#509579")) +  # 应用新的颜色方案
  annotate("text", x = -4, y = 0.10, label = paste("Nm = ", round(Nm, 3)), size = 7) +
  annotate("text", x = -4, y = 0.20, label = paste("R² = ", round(r2, 3)), size = 7) +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.line.x = element_line(size = 0.5, colour = "black"),
    axis.line.y = element_line(size = 0.5, colour = "black"),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black", size = 18),
    legend.position = "none",
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.text = element_text(size = 18),
    text = element_text(family = "sans", size = 18),
    panel.border = element_rect(color = "black", fill = NA, size = 1)  # 添加图表边框
     )


p1
###########要是需要计算各分类占比可以使用在图中加入饼图，代码如下#########

# 计算饼图数据(下面不清晰的话可以参考)
#low = nrow(out[out[,6]== "#509579",])
#med = nrow(out[out[,6]== "#485970",])
#high = nrow(out[out[,6]== "#cf9198",])
#type <- c('med','high','low')
#nums <- c(med,high,low)
#df <- data.frame(type = type, nums = nums)

data_summary <- as.data.frame(table(out$group))
colnames(data_summary) <- c("group", "nums")
data_summary$type<-(c("med","low","high"))
data_summary$percentage <- round(data_summary$nums / sum(data_summary$nums) * 100, 1)
data_summary$label <- paste(data_summary$type, paste(data_summary$percentage, "%", sep = ''))

# 绘制饼图
p2 <- ggplot(data_summary, aes(x = "", y = nums, fill = group)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar(theta = "y") +
  scale_fill_manual(
    values = c("#485970","#509579","#cf9198"),  # 应用新的颜色方案
    labels = data_summary$label
  ) +
  theme_void() +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.text = element_text(size = 18)
  )
p2

# 嵌入饼图到主图中
p_final <- p1 + geom_subview(subview = p2, x = -11.5, y = 0.95, w = 4.5, h = 4.5)

# 显示最终图形
print(p_final)

2.3、输出结果

3、参考文献

[1] Sloan, W.T., Lunn, M., Woodcock, S., Head, I.M., Nee, S. and Curtis, T.P. (2006), Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environmental Microbiology, 8: 732-740.

[2] Cheng, M., Luo, S., Zhang, P. et al. A genome and gene catalog of the aquatic microbiomes of the Tibetan Plateau. Nat Commun 15, 1438 (2024).

[3] Zorea, A., Pellow, D., Levin, L. et al. Plasmids in the human gut reveal neutral dispersal and recombination that is overpowered by inflammatory diseases. Nat Commun 15, 3147 (2024).

[4] 中性群落模型（neutral community model，NCM）分析及R语言计算，小白鱼的生统笔记。

4、相关信息

！！！有需要的小伙伴可以去评论区自取今天的测试代码和实例数据。

📌示例代码中提供了数据和代码，已经小编测试，可直接运行。

⚠️本账号会不时更新代码，代码失效可联系小编~

以上就是利用NCM中性模型全部内容。