主编寄语
微生物组发展到现在，相关分析技术和体系已经十分成熟，而微生物组数据的分析方法更是各式各样，我们借助“R语言在微生物组的最佳实践”这篇之前我们发表的论文带领大家从简单到困难，从基础到流程，从单一功能到组合功能，从上游分析到下游分析的全部代码。全套流程全部借助于Rstudio作为平台，在win和linux和mac下都可以完善运行，大家可以自己准备相关硬件资源，跟着我们这批教程进行全面学习。
这篇综述是我们小组在2023年发表在protein cell 期刊，整体内容几乎囊括了微生物组数据分析的全部内容，代码脚本有2w+行，全部属于公开分享阶段，但由于内容很多，为了方便大家理解和使用学习，这里选择分期带大家运行。
微生信生物分享的这批内容经过如下步骤：
同学A进行代码整理运行，基础注释书写，形成分享稿件的基本样貌；
同学B对代码复现一遍，解决其中无法顺利运行的地方，以及再次注释内容无法让自己明白的地方，完善分享稿件；
同学C再次运行一遍，进行全部流程代码学习，继续注释自己不明白的地方，最终形成完善分享稿件。
我们希望通过这个教程带领需要学习的小伙伴开启全民免费学习微生物组联合挖掘时代。给大家尽量做高质量的教程；
写在前面
“R语言在微生物组的最佳实践”的综述流程，我们已经学习了七节，有很多朋友问怎么做群落构建过程，那我们今天就先跳过前面的一些基础分析，先介绍一下群落构建过程部分。
今天继续带大家学习一下：群落构建过程（Neutral model（中性模型）、RCbray calculation（计算RCbray）、βNTI-nearest taxon index（最近种间亲缘关系指数）
需要注意的是：每次分析之前，需要加载”数据分析准备“ 这一节的代码，才可以开始下游分析
这次的分享内容希望能够帮助大家更加深入的挖掘自己的数据，给大家带来分析数据的便捷，同时也欢迎大家可以引用我们的文章。
下面是文章的具体信息：
Wen, T., et al., The best practice for microbiome analysis using R. Protein & Cell, 2023. 14(10): p. 713-725.
文章链接：https://doi.org/10.1093/procel/pwad024
流程实践数据：https://github.com/taowenmicro/EasyMicrobiomeR
数据代码获取地址
数据地址：https://github.com/taowenmicro/EasyMicrobiomeR/tree/master/data
全套数据代码整套流程下载地址：https://github.com/taowenmicro/EasyMicrobiomeR/tree/master
实例介绍
以下图片来自“Specific metabolites drive the deterministic assembly of diseased rhizosphere microbiome through weakening microbial degradation of autotoxin”文章，于2022年发表于Microbiome，查看链接https://link.springer.com/article/10.1186/s40168-022-01375-z#/Abs1
图例书写参考
c Contributions of deterministic and stochastic processes on community assembly in the sampling data. βNTI calculations of two soils conditioned with metabolites indicate that variable selection has greater effects on disease than health. d The relative influence of each community assembly process of two soils conditioned with metabolites as defined by the percentage of site pairs governed by each process. S1C1 means soil1 conditioned by metabolites at a concentration of 1 μM, S1C2 means soil1 conditioned by metabolites at a concentration of 100 μM, S2C1 means soil2 conditioned by metabolites at a concentration of 1 μM, and S2C2 means soil2 conditioned by metabolites at a concentration of 100 μM. Different letters mean significant difference among groups (p < 0.05, Wilcoxon t-test)以下图片来自“Specific metabolites drive the deterministic assembly of diseased rhizosphere microbiome through weakening microbial degradation of autotoxin”文章，于2022年发表于Microbiome，查看链接https://link.springer.com/article/10.1186/s40168-022-01375-z#/Abs1
图例书写参考
b Contributions of deterministic and stochastic processes in community assembly within collected rhizosphere soil samples. βNTI calculation of phylogenetic turnover among diseased and healthy samples indicates that variable selection was more consistent in diseased soils. c The relative influence of each community assembly process among diseased and healthy samples was defined by the percentage of site pairs governed by each process. BD, diseased banana; BH, healthy banana (from Hainan); CD, diseased cucumber; CH, healthy cucumber (from Guangdong); WD, diseased watermelon; WH, healthy watermelon (from Beijing); LD, diseased lily; LH, healthy lily (from Hunan). e Contributions of deterministic and stochastic processes on community assembly within diseased and healthy soils of collected metadata. βNTI calculations of phylogenetic turnover between diseased and healthy soils indicate that variable selection has greater effects on disease than health. f The relative influence of each community assembly process between diseased and healthy soils as defined by the percentage of site pairs governed by each process.以下图片来自“Specific metabolites drive the deterministic assembly of diseased rhizosphere microbiome through weakening microbial degradation of autotoxin”文章，于2022年发表于Microbiome，查看链接https://link.springer.com/article/10.1186/s40168-022-01375-z#/Abs1
图例书写参考
c Scatter plots illustrating the correlation between βNTI and the relative abundance of the five featured metabolites. p-values were evaluated and significant correlations were determined at p < 0.05
群落构建过程
Neutral model（中性模型）
中性模型（Neutral model）是生态学中用来描述和理解物种多样性和群落动态的一个理论框架。它的核心概念是基于中立生态理论（Neutral Theory of Biodiversity），该理论认为在某些生态系统中，物种之间的竞争作用和生态位分化可能不如随机过程对物种多样性和群落结构的影响显著。
基本概念：
随机性：中性模型强调物种的出生、死亡、迁移和扩散等过程是随机的，而不是由物种之间的竞争或适应性差异所决定。这意味着物种在群落中的丰度和分布模式主要是由随机事件和概率决定的。
等效性：在中性模型中，所有物种被假设为在生态位上是等效的，即它们对资源的利用效率和生态策略没有差异。这与传统的生态学理论不同，后者认为物种之间的竞争和适应性差异是群落结构和动态的主要驱动力。
物种多样性：中性模型提供了一种解释物种多样性的机制，即物种多样性可以通过随机过程在没有竞争压力的情况下自然产生和维持。
RCbray calculation（计算RCbray）
RCbray是一种衡量微生物群落在空间分布上的均匀性的指标。它通过计算微生物在不同空间位置的分布变异系数来评估群落的空间异质性。RCbray的值越低，表示微生物在空间上的分布越均匀；值越高，表示空间分布的异质性越大。这对于理解微生物如何在不同生境中分布和迁移非常重要。
βNTI-nearest taxon index（最近种间亲缘关系指数）
βNTI是一种衡量微生物群落组成差异的指标。它通过比较不同样本之间的微生物种类和相对丰度来评估群落的异质性。βNTI的值越高，表示不同样本间的微生物组成差异越大，即群落结构的异质性越高。这可以反映环境因素、宿主差异或其他生态压力对微生物群落结构的影响。
βNTI and RCbray combination analysis（βNTI和RCbray联合出图）
结合βNTI和RCbray分析，研究者可以更全面地了解微生物群落的多样性和空间分布模式。例如，通过这种分析，科学家可以识别出影响微生物群落结构和功能的关键环境因子，评估不同生态系统中微生物群落的稳定性和恢复力，以及预测和解释微生物群落对环境变化的响应。
下面开始实践操作
Neutral model（中性模型）
简化假设：中性模型通过简化的假设减少了生态系统复杂性，使得模型更易于处理和分析。
数学表达：中性模型通常通过一系列数学方程式来描述物种的随机出生、死亡和迁移过程，这些方程式可以用来模拟群落的动态变化。
预测能力：中性模型可以用来预测群落结构的变化，例如物种丰度分布的模式和物种多样性的变化趋势。
应用场景：
微生物群落分析：中性模型特别适用于微生物群落的研究，因为微生物生态系统通常具有高度的物种多样性和复杂的相互作用，这些特性使得随机过程在群落动态中起到重要作用。
生态系统管理：中性模型可以帮助生态学家和环境管理者理解生态系统在受到干扰后如何恢复，以及物种多样性如何随时间变化。
生物多样性保护：通过中性模型，可以预测不同保护措施对物种多样性的影响，从而为生物多样性保护提供理论支持。
#数据准备library(picante)library(ape)library(vegan)library(FSA)library(eulerr)library(grid)library(gridExtra)require(minpack.lm)require(Hmisc)require(stats4)library(parallel)env = read.csv("./data/dataNEW/env.csv")envRDA = envhead(env)row.names(envRDA) = env$IDenvRDA$ID = NULLhead(envRDA)psphy = filter_taxa(ps, function(x) sum(x ) > 1000 , TRUE);psphyres1path = "./result_and_plot/"phypath = paste(res1path,"/Phylogenetic_analyse_spacies/",sep = "")dir.create(phypath)map = sample_data(ps)n = map$Group %>% unique() %>%length()n#Neutral model（中性模型）neutralModel = function(otu = NULL,tax = NULL,map = NULL,tree = NULL,ps = NULL,group = "Group",ncol = 3,nrow = 1){# 抽平，默认使用最小序列抽平ps = inputMicro(otu,tax,map,tree,ps,group = group)psset.seed(72) #设置随机种子，保证结果可重复psrare = rarefy_even_depth(ps)# 标准化ps.norm = transform_sample_counts(psrare, function(x) x/sum(x))#------------------------------------------开始计算中性模型----------------------------------------------------------map = as.data.frame(sample_data(psrare))aa = levels(map$Group)aamap$ID = row.names(map)plots = list()dat1 = list()dat2 = list()i =1for (i in 1:length(aa)) {maps<- dplyr::filter(as.tibble(map),Group %in%aa[i])maps = as.data.frame(maps)row.names(maps) = maps$IDps_sub = psraresample_data(ps_sub) =maps ;ps_sub# 提取OTU表格OTU.table = t(otu_table(ps_sub))head(OTU.table )# 将整个群落看做一个整体，计算每个样本的序列数，并求取均值Calculate the number of individuals in the meta community (Average read depth)N <- mean(apply(OTU.table, 1, sum))#计算每个OTU的的平均序列数 Calculate the average relative abundance of each taxa across communitiesp.m <- apply(OTU.table, 2, mean)#去除OTU序列数为0的OTUp.m <- p.m[p.m != 0]p <- p.m/Np.df = data.frame(p) %>%rownames_to_column(var="OTU")# Calculate the occurrence frequency of each taxaOTU.table.bi <- 1*(OTU.table>0)freq.table <- apply(OTU.table.bi, 2, mean)freq.table <- freq.table[freq.table != 0]freq.df = data.frame(OTU=names(freq.table), freq=freq.table)#CombineC <- inner_join(p.df,freq.df, by="OTU") %>%arrange(p)# Remove rows with any zero (absent in either source pool or local communities). You already did this, but just to make sure we will do it again.C.no0 <- C %>%filter(freq != 0, p != 0)#Calculate the limit of detectiond <- 1/N##Fit model parameter m (or Nm) using Non-linear least squares (NLS)p.list <- C.no0$pfreq.list <- C.no0$freqm.fit <- nlsLM(freq.list ~ pbeta(d, N*m*p.list, N*m*(1-p.list), lower.tail=FALSE), start=list(m=0.1))m.ci <- confint(m.fit, 'm', level=0.95)m.sum <- summary(m.fit)m.coef = coef(m.fit)freq.pred <- pbeta(d, N*coef(m.fit)*p.list, N*coef(m.fit)*(1-p.list), lower.tail=FALSE)Rsqr <- 1 - (sum((freq.list - freq.pred)^2))/(sum((freq.list - mean(freq.list))^2))# Get table of model fit statsfitstats <- data.frame(m=m.coef, m.low.ci=m.ci[1], m.up.ci=m.ci[2],Rsqr=Rsqr, p.value=m.sum$parameters[4], N=N,Samples=nrow(OTU.table), Richness=length(p.list),Detect=d)# Get confidence interval for predictionsfreq.pred.ci <- binconf(freq.pred*nrow(OTU.table), nrow(OTU.table), alpha=0.05, method="wilson", return.df=TRUE)# Get table of predictionspred.df <- data.frame(metacomm_RA=p.list, frequency=freq.pred,frequency_lowerCI=freq.pred.ci[,2],frequency_upperCI=freq.pred.ci[,3]) %>%unique()# Get table of observed occupancy and abundanceobs.df = C.no0 %>%dplyr::rename(metacomm_RA = p, frequency=freq)head(obs.df)p = ggplot(data=obs.df) +geom_point(data=obs.df, aes(x=log10(metacomm_RA), y=frequency),alpha=.3, size=2, color="#8DD3C7") +geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency), color="#FFFFB3") +geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_lowerCI), linetype=2, color="#FFFFB3") +geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_upperCI), linetype=2, color="#FFFFB3") +# geom_text(data=fitstats, aes(label = paste("R^2 == ", round(Rsqr, 3))),# x=1, y=0.75, size=4, parse=TRUE) +# geom_text(data=fitstats, aes(label = paste("italic(m) ==", round(m, 3))),# x=-1, y=0.85, size=4, parse=TRUE) +labs(x="Log10 abundance in\nmetacommunity", y="Frequency detected",title = paste(aa[i],paste("R^2 == ", round(fitstats$Rsqr, 3)),paste("italic(m) ==", round(fitstats$m, 3)))) +theme_bw() +theme(axis.line = element_line(color="black"),legend.position = "none",axis.title = element_text(size=14),axis.text = element_text(size=12))pplots[[aa[i]]] = pdat1[[aa[i]]] = obs.dfdat2[[aa[i]]] = pred.df}# plots$ABCD# library(ggpubr)# nrow=2,,ncol=4p = ggpubr::ggarrange(plotlist = plots,common.legend = TRUE, legend="right",ncol = ncol,nrow = nrow)preturn(list(p,plots,dat1,dat2))}result = neutralModel(ps = ps,group = "Group",ncol = 3)#--合并图表p1 = result[[1]]p1FileName <- paste(otupath,"1_neutral_modelCul", ".pdf", sep = "")ggsave(FileName, p1,width = 12,height = 4)FileName <- paste(otupath,"1_neutral_modelCul", ".png", sep = "")ggsave(FileName, p1,width = 12,height = 4)#--系统发育信号# 将环境因子数据和样本信息数据进行合并，形成一个新的数据框mapE。# 从合并后的数据框中提取出对应的样本信息，用于后续的处理和分析。phyloSignal = function(otu = NULL,tax = NULL,map = NULL,tree = NULL ,ps = NULL,env = env,group = "Group",path = "./"){# 抽平，默认使用最小序列抽平ps = inputMicro(otu,tax,map,tree,ps,group = group)ps# 抽平set.seed(72)psrare = rarefy_even_depth(ps)# 标准化ps.norm = transform_sample_counts(psrare, function(x) x/sum(x))map = as.data.frame(sample_data(psrare))mapE =merge(map,env,by = "row.names",all= TRUE)row.names(mapE) = mapE$Row.namesmapE$Row.names = NULLmapE$ID = row.names(mapE)sample_data(ps.norm) = mapEaa = levels(mapE$Group)dir.create(path)#----------分组计算门特尔相关,将结果保存，因为计算时间很长，只需计算一个就好了#-------eco = "Endosp."for (eco in as.character(unique(mapE$Group))){# Subset dataprint(paste("Now running", eco))# sub.physeq = phyloseq::subset_samples(ps.norm , Group == eco)sub.physeq = ps.normotu = as.data.frame(vegan_otu(ps.norm))head(otu)map = as.data.frame(sample_data(ps.norm))mapsub <- map[map$Group == eco,]sample_data(sub.physeq) = mapsub# Remove OTUs not found in at least 3 samplesOTU.table = otu_table(sub.physeq)OTU.table[OTU.table > 0] = 1OTU.freq = rowSums(OTU.table)OTU.freq = OTU.freq[OTU.freq > 2]sub.physeq = prune_taxa(names(OTU.freq), sub.physeq)sub.physeq# get phylogenetic distancestree = phy_tree(sub.physeq)phylo.dist = cophenetic(tree)sample_OTUs = tree$tip.labelsam.phylo.dist = phylo.dist[sample_OTUs, sample_OTUs]sam.phylo.dist[upper.tri(sam.phylo.dist, diag=TRUE)] = NA# Generate dataframe of niche preference for pH, SOC and CN# site.chem.mat = data.frame(sample_data(sub.physeq)) %>%# # mutate(CN = percent_C / percent_N) %>%# dplyr::select(ID, colnames(env))site.chem.mat = env[row.names(env) %in% row.names(mapsub),]# rownames(site.chem.mat) = site.chem.mat$ID# site.chem.mat$ID = NULLsite.chem.mat = as.matrix(site.chem.mat)otu.table = t(otu_table(sub.physeq))# head(otu.table)match(row.names(otu.table),row.names(site.chem.mat))OTU.niche = wascores(site.chem.mat, otu.table)OTU.niche.df = data.frame(OTU.niche)head( OTU.niche.df)# i =1for (i in 1:dim(OTU.niche.df)[2]) {pH.pref = OTU.niche.df[[i]]names(pH.pref) = rownames(OTU.niche.df)pH.dist = as.matrix(dist(pH.pref), labels=TRUE)sam.pH.dist = pH.dist[sample_OTUs, sample_OTUs]sam.pH.dist[upper.tri(sam.pH.dist, diag=TRUE)] = NAsam.pH.crlg = mantel.correlog(sam.pH.dist, sam.phylo.dist)# ?mantel.correlogfilename = paste(path,eco,colnames(OTU.niche.df[i]), "_crlg.rds", sep="_")saveRDS(sam.pH.crlg, file=filename)}}}# env参数作为环境因子的数据框传递给函数。函数的主要目的是根据不同的生态群落（通过Group参数确定）和不同的环境因子（通过colnames(env)确定）绘制Mantel相关性的散点图，并标记显着性。phySigPlot = function(otu = NULL,tax = NULL,map = NULL,tree = NULL,ps = NULL,group = "Group",env = env,path = "./"){# 抽平，默认使用最小序列抽平ps = inputMicro(otu,tax,map,tree,ps,group = group)psmapE = as.data.frame(sample_data(ps))for (eco in levels(mapE$Group)) {# eco = "KO"# i = 1for (i in 1:length(colnames(env))) {ag.pH.crlg = data.frame(readRDS(file=paste(path,eco,colnames(env[i]), "_crlg.rds", sep="_"))$mantel.res) %>%mutate(Group = eco, property = colnames(env)[i])if (i == 1) {data = ag.pH.crlg}if (i != 1) {data = rbind(data,ag.pH.crlg )}}if (eco == levels(mapE$Group)[1]) {data2 = data}if (eco != levels(mapE$Group)[1]) {data2 = rbind(data2, data)}}dim(data2)eco.crlg = data2 %>%mutate(sig = ifelse(Pr.corrected. <= 0.05, "significant", "non-significant")) %>%filter(!(is.na(Pr.corrected.)))eco.crlg$Group= factor(eco.crlg$Group)p = ggplot(data=eco.crlg, aes(x=class.index, y=Mantel.cor)) +geom_point(data=eco.crlg[eco.crlg$sig=="significant",], color = "black", size=2, shape=16) +geom_point(data=eco.crlg[eco.crlg$sig=="non-significant",], color = "black",size=2, shape=1) +geom_line(data=eco.crlg, aes(color=property)) +geom_hline(yintercept = 0, linetype=2) +labs(x = "Phylogenetic distance class", y="Mantel correlation", color="property") +# facet_grid(~Group)facet_wrap(~Group,scales="free_y",ncol = 4)return(list(p,eco.crlg,data2))}# 从文件中读取环境因子数据，并将其存储在名为env的数据框中。# 复制env数据框到名为envRDA的新数据框中，并对envRDA进行一些处理：将ID列的值作为行名，并将ID列从envRDA中移除。# 创建一个文件夹phypath2，用于存储生成的生物学信号分析结果。# 调用phyloSignal函数对生态群落数据进行生物学信号分析。该函数需要传递生态群落数据（通过ps参数传递），分组信息（通过group参数传递），以及处理后的环境因子数据（通过env参数传递）。# 调用phySigPlot函数绘制Mantel相关性的散点图，其中包括生态群落数据和环境因子数据。该函数需要传递生态群落数据（通过ps参数传递），分组信息（通过group参数传递），以及处理后的环境因子数据（通过env参数传递）。env = read.csv("./data/dataNEW/env.csv")head(env)envRDA = envhead(env)row.names(envRDA) = env$IDenvRDA$ID = NULLhead(envRDA)#phypath2 = paste(otupath,"/phyloSignal/",sep = "")dir.create(phypath2)phyloSignal(ps = ps %>% filter_OTU_ps(400),group = "Group",env = envRDA[,2:3],path = phypath2)result = phySigPlot(ps = ps,group = "Group",env = envRDA[,2:3],path = phypath2)##提取图片p2 = result[[1]] + mytheme1p2#-提取作图数据data = result[[2]]head(data)FileName <- paste(phypath2,"2_phySigPlot", ".pdf", sep = "")ggsave(FileName, p2,width = 15,height = 6)FileName <- paste(phypath2,"2_phySigPlot", ".csv", sep = "")write.csv(data,FileName)
系统发育（Phylogenetic）距离分类的图表。图表展示了不同分类水平上的系统发育距离，并对比了不同环境因子（SOC和TN）对这些分类的影响。以下是对图表核心内容的整理：
系统发育距离分类:
图表的横轴表示系统发育距离的分类，通常这些分类代表物种之间的进化关系差异。
纵轴可能表示环境因子的浓度或某种生物学上的重要性，但具体含义未在图表中给出。
环境因子:
SOC和TN代表土壤有机碳（Soil Organic Carbon）和总氮（Total Nitrogen）含量，这两个因子是土壤生态系统中的重要参数。
图表中展示了SOC和TN对不同系统发育距离分类的影响，但具体的数值和单位未在图表中明确。
数据分布:
图表中的线型可能表示随着系统发育距离的增加，SOC和TN的浓度或影响如何变化。
可以看到，对于某些系统发育距离分类，SOC和TN的影响呈现出不同的模式，有的分类显示出正相关（例如，随着SOC的增加，某个指标也增加），而有的分类显示出负相关或无明显关系。
中性模型（neutral model）的图表，展示了三个不同群体（Group1, Group2, Group3）在元群落（metacommunity）中的对数丰度（Log10 abundance）分布情况。图表中的每个群体都有五个数据点，分别对应不同的对数丰度值。同时，图表提供了每个群体的R^2值和m值，这些值是用来评估数据拟合优度的统计参数。
以下是对图表核心内容的整理：
Group1:
R^2值为0.658，表示模型对数据的解释程度较高。
m值为0.077，这是与R^2值一起使用的统计参数，但在此上下文中没有具体解释。
对数丰度分布在0到1的范围内，其中最高丰度值为0.75。
Group2:
R^2值为0.589，略低于Group1，但仍表示模型有较高的解释能力。
m值为0.072，与Group1的m值相近。
对数丰度分布同样在0到1的范围内，最高丰度值也为0.75。
Group3:
R^2值为0.774，是三个群体中最高的，表明模型对数据的拟合最佳。
m值为0.299，显著高于前两个群体的m值。
对数丰度分布的波动范围更大，从0.25到1.00，最高丰度值为1。
从这些数据可以看出，Group3的模型拟合效果最好，而Group1和Group2的拟合效果相近但略低于Group3。对数丰度的分布范围在三个群体中有所不同，Group3显示出更广泛的丰度变化。
βNTI-nearest taxon index（最近种间亲缘关系指数）
特点：
量化群落组装的确定性与随机性：βNTI通过比较观察到的群落与通过零模型预期的群落之间的系统发育距离差异，来量化确定性过程（如生态选择）和随机过程（如漂移）在群落组装中的相对贡献。
系统发育信号的考量：βNTI考虑了系统发育信号，即亲缘关系密切的物种倾向于具有相似的生态位和环境偏好，这一假设在较短的系统发育距离内更为显著。
数值范围：βNTI的值范围通常在-2到+2之间。值大于+2或小于-2表示群落组装主要受确定性过程影响，正值与变量选择相关，负值与同质选择相关。值接近0则表明群落组装主要受随机过程影响。
对环境变化的敏感性：βNTI能够反映环境变化对群落组装的影响，如研究发现，随着环境条件（如土壤pH和稀释梯度）的变化，βNTI值也会发生变化，从而揭示确定性与随机性过程的相对重要性。
使用场景：
微生物群落研究：βNTI广泛应用于微生物群落的研究中，尤其是在分析微生物群落对环境变化的响应时。例如，在研究土壤微生物群落对不同土壤类型和pH值的响应时，βNTI可以帮助理解确定性过程和随机性过程的作用。
植物群落动态：在植物群落动态研究中，βNTI可以用来分析物种多样性、物种丰富度和群落结构随时间的变化趋势，以及这些变化背后的生态过程。
生态系统管理与保护：βNTI可以用于评估生态系统管理措施的效果，如在自然保护区中，通过监测βNTI的变化，可以了解保护措施是否有效地维护了群落的多样性和稳定性。
生态模型的验证：βNTI还可以用于验证和改进生态模型，通过比较模型预测的群落结构与实际观测到的群落结构，可以评估模型对生态过程的模拟精度。
bNTICul = function(otu = NULL,tax = NULL,map = NULL,tree = NULL ,ps = NULL,group = "Group",num = 99,thread = 6){ps = inputMicro(otu,tax,map,tree,ps,group = group)psps_sub <- ps# tree = phy_tree(ps)# tree#-------------调整map文件-----------------------------------------------------------------#添加一个ID列map = as.data.frame(sample_data(ps_sub))map$ID = row.names(map)sample_data(ps) = map#-----------准备OTU表格---------------------抽平-不设置抽平条数，默认按照最小序列数数目抽平set.seed(72) # setting seed for reproducibilitypsrare = rarefy_even_depth(ps_sub)#检查序列数量sample_sums(psrare)# 标准化数据ps.norm = transform_sample_counts(psrare, function(x) x/sum(x))# 计算βMNTD对每个随机零模型群落bMNTD_null_func <- function(i, OTU.table, tree){tree$tip.label = sample(tree$tip.label)bMNTD_s = comdistnt(OTU.table, cophenetic(tree), abundance.weighted = TRUE)A <- attr(bMNTD_s, "Size")B <- if (is.null(attr(bMNTD_s, "Labels"))) sequence(A) else attr(bMNTD_s, "Labels")if (isTRUE(attr(bMNTD_s, "Diag"))) attr(bMNTD_s, "Diag") <- FALSEif (isTRUE(attr(bMNTD_s, "Upper"))) attr(bMNTD_s, "Upper") <- FALSEbMNTD_s.df = data.frame(Sample_1 = B[unlist(lapply(sequence(A)[-1], function(x) x:A))],Sample_2 = rep(B[-length(B)], (length(B)-1):1),bMNTD = as.vector(bMNTD_s),rep=i)return(bMNTD_s.df)}# 计算βNTIPhylo_turnover <- function(physeq, reps, nproc){# Extract OTU tableOTU.table = t(otu_table(physeq))# Extract phylogenetic treetree = phy_tree(physeq)# Get βMNTD between all communitiesbMNTD_o = comdistnt(OTU.table, cophenetic(tree), abundance.weighted = TRUE)A <- attr(bMNTD_o, "Size")B <- if (is.null(attr(bMNTD_o, "Labels"))) sequence(A) else attr(bMNTD_o, "Labels")if (isTRUE(attr(bMNTD_o, "Diag"))) attr(bMNTD_o, "Diag") <- FALSEif (isTRUE(attr(bMNTD_o, "Upper"))) attr(bMNTD_o, "Upper") <- FALSEbMNTD_o.df = data.frame(Sample_1 = B[unlist(lapply(sequence(A)[-1], function(x) x:A))],Sample_2 = rep(B[-length(B)], (length(B)-1):1),bMNTD = as.vector(bMNTD_o))# Get βMNTD for randomized null communitiesrep.list = seq(1, reps)bMNTD_s.df.list = mclapply(rep.list, bMNTD_null_func, OTU.table=OTU.table, tree=tree, mc.cores=nproc)# Combine all data together and calculate βNTI for each sample pairbMNTD_s.df <- do.call("rbind", bMNTD_s.df.list)bMNTD_s.means.df = bMNTD_s.df %>%group_by(Sample_1, Sample_2) %>%dplyr::summarize(mean_bMNTD = mean(bMNTD),sd_bMNTD = sd(bMNTD))bMNTD_o.df = inner_join(bMNTD_o.df, bMNTD_s.means.df, by=c("Sample_1", "Sample_2")) %>%mutate(bNTI = (bMNTD - mean_bMNTD)/sd_bMNTD)return(bMNTD_o.df)}#========这里一把单核就真实数据而言需要超过10个小时，跑999次，所以需要多核# 计算bnti，这里可以设置线程数量，是第三个参数，我们在linux下面可以设置，30个线程# 第二个参数设置迭代数量，这里文献一般999嘛。bNTI = Phylo_turnover(psrare, num, thread)return(list(bNTI))}result = bNTICul(ps = psphy,group = "Group",num = 100,thread = 1)bNTI = result[[1]]head(bNTI)filename = paste(phypath,"/4_bNTI.csv",sep = "")write.csv(bNTI, filename)RCbray calculation（计算RCbray）区分生态过程：RCbray可以帮助研究者区分影响群落结构的生态过程。通过比较观察到的物种组成更替与通过零模型预期的更替，RCbray可以揭示是确定性过程还是随机性过程在群落构建中起主导作用。评估群落相似性：RCbray通过比较实际观察到的Bray-Curtis相似性指数与零模型预期的相似性指数之间的偏差，来评估不同群落之间的相似性。这个值的范围在-1到+1之间，正值表示观察到的群落比随机预期的更相似，而负值则表示观察到的群落比随机预期的更不同。理解群落动态：RCbray可以用于理解群落随时间或空间变化的动态。例如，在研究土壤微生物群落对环境变化的响应时，RCbray能够揭示环境因素如何影响群落的组成和结构。生态系统管理和保护：RCbray的计算结果可以为生态系统管理和保护提供重要信息。通过识别影响群落构建的主要生态过程，管理者可以采取针对性的措施来维护或恢复生态系统的健康和多样性。微生物群落研究：在微生物生态学中，RCbray尤其有用，因为它可以帮助研究者理解微生物群落如何组装，以及环境因素（如土壤类型、pH值、有机质浓度等）如何影响这些过程。RCbary = function(otu = NULL,tax = NULL,map = NULL,tree = NULL ,ps = NULL,group = "Group",num = 99,thread = 6){ps_sub <- ps#----------------整理map文件map = as.data.frame(sample_data(ps_sub))map$ID = row.names(map)sample_data(ps) = map#-------------------准备OTU表格#-----------------抽平-不设置抽平条数，默认按照最小序列数数目抽平set.seed(72) # setting seed for reproducibilitypsrare = rarefy_even_depth(ps_sub )#检查序列数量sample_sums(psrare)# 标准化数据ps.norm = transform_sample_counts(psrare, function(x) x/sum(x))#--------------两个函数# 对模拟群落计算距离RCbray_null_func <- function(i, freq.abd.df, alpha1, alpha2, N){# Get simulated communities and distance## initally select OTUs weighted by their frequency. The number of OTUs selected should equal the richness of the samples.simcom1 = data.frame(table(sample(freq.abd.df$OTU, size=alpha1, replace=FALSE, prob=freq.abd.df$freq)), stringsAsFactors = F)colnames(simcom1) = c("OTU","simcom1")simcom1$OTU = as.character(simcom1$OTU)simcom1 = inner_join(simcom1, freq.abd.df, by="OTU")simcom2 = data.frame(table(sample(freq.abd.df$OTU, size=alpha2, replace=FALSE, prob=freq.abd.df$freq)), stringsAsFactors = F)colnames(simcom2) = c("OTU","simcom2")simcom2$OTU = as.character(simcom2$OTU)simcom2 = inner_join(simcom2, freq.abd.df, by="OTU")## Now recruit OTUs based on their abundance in the metacommunitysimcom1.abd = data.frame(table(sample(simcom1$OTU, size=N-alpha1, replace=TRUE, prob=simcom1$p)), stringsAsFactors = F)colnames(simcom1.abd) = c("OTU","simcom1.abd")simcom1.abd$OTU = as.character(simcom1.abd$OTU)simcom1 = full_join(simcom1, simcom1.abd, by="OTU") %>%mutate(simcom1.abd = ifelse(is.na(simcom1.abd), 1, simcom1.abd)) %>%select(OTU, simcom1.abd)simcom2.abd = data.frame(table(sample(simcom2$OTU, size=N-alpha2, replace=TRUE, prob=simcom2$p)), stringsAsFactors = F)colnames(simcom2.abd) = c("OTU","simcom2.abd")simcom2.abd$OTU = as.character(simcom2.abd$OTU)simcom2 = full_join(simcom2, simcom2.abd, by="OTU") %>%mutate(simcom2.abd = ifelse(is.na(simcom2.abd), 1, simcom2.abd)) %>%select(OTU, simcom2.abd)simcom = full_join(simcom1, simcom2, by="OTU")simcom[is.na(simcom)] = 0rownames(simcom) = simcom$OTUsimcom$OTU = NULLnull.dist = vegdist(t(simcom), method="bray")[1]return(null.dist)}# 计算RCbray的主功能Calc_RCbray <- function(physeq, reps, nproc){# Get OTU table from phyloseq objectotu.table = otu_table(physeq)# Get alpha diversity for each sampleotu.PA.table = otu.tableotu.PA.table[otu.PA.table > 0] = 1alpha.df = data.frame(Sample_ID = colnames(otu.PA.table), OTU.n = colSums(otu.PA.table), stringsAsFactors = F)# Get beta diversity matrixbeta.table = as.matrix(vegdist(t(otu.PA.table), method="bray", diag=TRUE, upper=TRUE))## Get metacommunity# Calculate the number of individuals in the meta community (Average read depth)N <- mean(apply(t(otu.table), 1, sum))# Calculate the average relative abundance of each taxa across communitiesp.m <- apply(t(otu.table), 2, mean)p.m <- p.m[p.m != 0]p <- p.m/N# Calculate the occurrence frequency of each taxa across communitiesotu.table.bi <- 1*(t(otu.table)>0)freq <- apply(otu.table.bi, 2, mean)freq <- freq[freq != 0]# Combinefreq.abd.df = data.frame(p=p, freq=freq) %>%tibble::rownames_to_column(var="OTU") %>%filter(p != 0, freq != 0) %>%arrange(p)# For each pair of samples run the RCbray analysiscomps = combn(alpha.df$Sample_ID, m=2, simplify = F)RCb.df = data.frame(Site1 = character(), Site2 = character(), RCb = numeric(), stringsAsFactors = F)for (j in seq(1, length(comps))){sam = comps[[j]]alpha1 = alpha.df[alpha.df$Sample_ID == sam[1],]$OTU.nalpha2 = alpha.df[alpha.df$Sample_ID == sam[2],]$OTU.n# Permute "reps" many timesrep.list = seq(1, reps)null.list = mclapply(rep.list, RCbray_null_func, freq.abd.df=freq.abd.df, alpha1=alpha1, alpha2=alpha2, N=N, mc.cores=nproc)RCb = (length(null.list[null.list > beta.table[sam[1], sam[2]]]) + (0.5*length(null.list[null.list == beta.table[sam[1], sam[2]]])))/repsRCb = (RCb - 0.5)*2RCb.df = rbind(RCb.df, data.frame(Site1=sam[1], Site2=sam[2], RCb=RCb, stringsAsFactors = F))}RCb.dfreturn(RCb.df)}# 运行RCbray的计算，这个运算再5个小时左右999重复RCb = Calc_RCbray(psrare, num, thread)head(RCb)return(list(RCb))}result = RCbary(ps = psphy ,group = "Group",num = 10,thread = 1)RCbary = result[[1]]head(RCbary)filename = paste(phypath,"/5_RCb.csv",sep = "")write.csv(RCbary,filename)
βNTI and RCbray combination analysis（βNTI和RCbray联合出图）
特点：
多维度分析：
βNTI和RCbray联合出图允许研究者同时考虑微生物群落的组成差异（βNTI）和空间分布均匀性（RCbray），提供了一个多维度的视角来理解微生物群落的复杂性。
直观性：
通过图表的形式，研究者可以直观地看到不同样本或条件下微生物群落的变化，包括物种丰富度、组成差异以及空间分布模式。
比较性：
联合出图使得不同样本、时间点或实验条件下的微生物群落特征可以直接比较，有助于识别关键的环境或生物学因素。
解释性：
通过分析βNTI和RCbray的变化，研究者可以解释微生物群落对环境变化的响应，以及这些响应如何影响生态系统的功能和服务。
使用场景：
环境监测：
在评估环境变化对微生物群落的影响时，βNTI和RCbray联合出图可以帮助研究者理解污染、气候变化或其他环境压力如何改变微生物群落的结构和功能。
生态系统管理：
在生态系统恢复和管理项目中，这种分析可以帮助制定策略，通过优化微生物群落的结构和功能来提高生态系统的恢复力和稳定性。
农业应用：
在农业土壤健康和作物生产力的研究中，βNTI和RCbray联合出图可以用来评估不同农业实践对土壤微生物群落的影响，从而指导更可持续的农业管理方法。
生物技术：
在开发新的生物制品或生物技术应用时，了解微生物群落的多样性和分布特征对于设计有效的生物反应器和生物处理系统至关重要。
bNTIRCPlot = function(otu = NULL,tax = NULL,map = NULL,tree = NULL ,ps = NULL,RCb = RCb,bNTI = bNTI,group = "Group"){ps = inputMicro(otu,tax,map,tree,ps,group = group)pspsrare <- psmap = as.data.frame(sample_data(psrare))map$ID = row.names(map)sample_data(psrare) = map# Get habitat metadata and add it to the βNTI then merge with the RCbray dataseteco.meta1 = data.frame(sample_data(psrare)) %>%select(ID, Group) %>%dplyr::rename(Sample_1 = ID, Group_1 = Group)eco.meta2=data.frame(sample_data(psrare)) %>%select(ID, Group) %>%dplyr::rename(Sample_2 = ID, Group_2 = Group)# bNTI 匹配第一列和第二列的分组信息bNTI.df = inner_join(bNTI, eco.meta1) %>%inner_join(eco.meta2)# 合并两个数据turnover.df = inner_join(bNTI.df, RCb)head(turnover.df)dim(turnover.df)#--------------合并文件保存# write.csv(turnover.df,"./Result/bNTI//bNTI_RCbray.csv")#-----按照分组统计作图#------------bNIT作图dim(bNTI.df)within.bNTI.df = bNTI.df %>%filter(Group_1 == Group_2) %>%mutate(Group = Group_1)head(within.bNTI.df )eco.bNTI.plot <- ggplot(within.bNTI.df, aes(x=Group, y=bNTI)) +geom_jitter(alpha = 0.1,color ="#984EA3") +geom_boxplot(outlier.shape=1,outlier.alpha = 0,fill = "#984EA3") +geom_hline(yintercept = 2, linetype=2, size=0.5) +geom_hline(yintercept = -2, linetype=2, size=0.5) +labs(x="", y="bNTI") +theme_classic() +theme(legend.position = "none",axis.text = element_text(size=12),axis.text.x = element_text(angle=45, hjust=1),axis.title = element_text(size=14))# 现在按照RCbray进行分开标记系统发育过程eco.turnover.df = turnover.df %>%filter(Group_1 == Group_2) %>%mutate(Group = Group_1)head(eco.turnover.df )## Calculate the relative influence of each processeco.turnover.df = eco.turnover.df %>%mutate(process = ifelse(abs(bNTI) < 2,ifelse(abs(RCb) < 0.95, "Drift",ifelse(RCb >= 0.95, "Dispersal Limited",ifelse(RCb <= -0.95, "Homogenizing Dispersal", "ERROR"))),ifelse(bNTI >= 2, "Variable Selection",ifelse(bNTI <= -2, "Homogeneous Selection", "ERROR"))))eco.turnover.df$process = factor(eco.turnover.df$process, levels = c("Drift","Dispersal Limited", "Homogenizing Dispersal","Variable Selection", "Homogeneous Selection"))head(eco.turnover.df)#------计算每个组的系统发育过程中五个部分分别占有的比例pre = eco.turnover.df %>%dplyr::group_by(Group, process) %>%dplyr::summarize(n_sites = n(),perc=(n()/45)*100) %>%as.data.frame# head(numeco )numeco <- pre %>% dplyr::group_by(Group) %>%dplyr::summarise(num = sum(n_sites))alleco <- pre %>% dplyr::left_join(numeco,by = "Group")alleco$perc = alleco$n_sites/ alleco$num * 100sum.eco.turnover.df = allecoeco.turnover.plot = ggplot(sum.eco.turnover.df, aes(x=Group, y=perc, fill=process)) +geom_bar(stat="identity", color="black") +# scale_fill_manual(values = c("white", "grey75", "grey50", "black")) +labs(x="", y="Percent of pairs (%)", fill="Process") +theme_bw() +theme(panel.grid = element_blank(),axis.text = element_text(size=12),axis.text.x = element_text(angle=45, hjust=1),axis.title = element_text(size=14),legend.key.size = unit(10, "mm"),legend.text = element_text(size=12),legend.title = element_text(size=14))eco.turnover.plot# Merge the plotseco.plot = cowplot::plot_grid(eco.bNTI.plot, eco.turnover.plot,rel_widths=c(0.6, 1), labels=c("A", "B"))eco.plotreturn(list( eco.bNTI.plot, eco.turnover.plot,eco.plot,turnover.df,sum.eco.turnover.df))}bNTI = read.csv(paste(phypath,"/4_bNTI.csv",sep = ""),row.names = 1)head(bNTI)# RCbray 数据读入，修改列名RCb = read.csv(paste(phypath,"/5_RCb.csv",sep = ""),row.names = 1) %>%dplyr::mutate(Sample_1 = Site2, Sample_2 = Site1)head(RCb)result = bNTIRCPlot(ps = psphy ,RCb = RCb,bNTI = bNTI,group = "Group")#--bNTI出图片p3 <- result[[1]]p3#RCbary可视化p4 <- result[[2]]p4#组合图片BNTI，RCbrayp5 <- result[[3]]p5plotdata = result[[4]]head(plotdata)dat = result[[5]]head(dat)filename = paste(phypath,"/6_bNTI_RCbray.csv",sep = "")write.csv(plotdata,filename)FileName <- paste(phypath,"6_bNTI", ".pdf", sep = "")ggsave(FileName, p3,width =8,height = 6)FileName <- paste(phypath,"6_RCbary", ".pdf", sep = "")ggsave(FileName, p4,width = 6,height = 6)FileName <- paste(phypath,"6_BNTI_RCbray", ".pdf", sep = "")ggsave(FileName, p5,width = 12,height = 8)FileName <- paste(phypath,"6_bNTI", ".png", sep = "")ggsave(FileName, p3,width =8,height = 6)FileName <- paste(phypath,"6_RCbary", ".png", sep = "")ggsave(FileName, p4,width = 6,height = 6)FileName <- paste(phypath,"6_BNTI_RCbray", ".png", sep = "")ggsave(FileName, p5,width = 12,height = 8)FileName <- paste(phypath,"6_RCbray.percent.csv", sep = "")write.csv(dat,FileName, quote = F)
βNTI 图形解读：
βNTI值的范围通常在-2到+2之间。
当βNTI值接近0时，表明群落构建主要受随机过程（如生态漂变）的影响，没有显著的选择或扩散限制。
当βNTI值小于-2时，表明群落构建可能主要受到同质选择的影响，即群落中的物种更倾向于在相似的环境条件下被选择。
当βNTI值大于+2时，表明群落构建可能主要受到变量选择的影响，即不同环境下物种的适应性差异导致了群落组成的显著变化。
图中通常会展示βNTI的分布或其与某些环境因子的关系，例如，βNTI与土壤有机质(SOM)浓度的关系。
通过观察βNTI值的变化，可以推断群落构建过程中随机性和确定性（选择）的相对贡献。
RCbray 图形解读：
RCbray是基于Bray-Curtis dissimilarity计算的，用于量化两个群落之间的物种组成差异。
RCbray值越接近0，表示两个群落之间的物种组成越相似；值越大，表示物种组成差异越大。
在图形中，RCbray值可以用来展示不同群落之间的相似性或差异性，或者与βNTI值一起展示，以揭示群落构建过程中的随机性和选择性。
通过比较不同环境条件下的RCbray值，可以分析环境因素如何影响群落的构建和物种的分布。
运行时的错误
Error in nlsLM(freq.list ~ pbeta(d, N * m * p.list, N * m * (1 - p.list), :
could not find function "nlsLM"
Error in wascores(site.chem.mat, otu.table) :
could not find function "wascores"
R包没有加载library()加载对应R包即可
library(minpack.lm)
library(Hmisc)
library(stats)
library(vegan)根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣院士土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 Nature Communications，ISME J，Microbiome，SCLS，New Phytologist，iMeta，Fundamental Research， PCE，SBB，JAFC（封面），Horticulture Research，SEL（封面），BMC plant biology等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。
撰写：牛国庆 文浩仰
修改：文涛
排版：杨雯儀
审核：袁军 团队工作及其成果 （点击查看）
了解 交流 合作
小组负责人邮箱 袁军：junyuan@njau.edu.cn；
小组成员文涛：taowen@njau.edu.cn等
团队公众号：微生信生物 添加主编微信，或者后台留言。
1.仅限相关专业或研究方向人员添加，必须实名，不实名则默认忽略。
2.非相关专业的其他人员及推广宣传人员禁止添加。
3.添加主编微信需和简单聊一聊专业相关问题，等待主编判断后，可拉群。
4 微生信生物VIP微信群不受限制，给微生信生物供稿一次即可加入（群里发送推文代码+高效协助解决推文运行等问题+日常问题咨询回复）。
加主编微信 加入群聊
目前营销人员过多，为了维护微生信生物3年来维护的超5500人群聊，目前更新进群要求：
