CCA分析图如何解读_微生物群落与环境因子关联：全自动的CCA/RDA分析流程！！…

环境因子关联

在微生物群落研究中，我们通常都希望知道微生物群落变化是会受到哪些因素的影响，在这些因素中，什么因素是主要的影响因素？

回答这一问题需要将微生物群落数据与其对应的环境因子进行关联分析，在这一项分析中，使用最早也是使用频率最高的分析就是CCA/RDA。

CCA/RDA

CCA的全称是典范对应分析 (Canonical correspondence analysis)，RDA的全称是冗余分析 (Redundancy analysis)。

CCA和RDA都属于排序分析，理论上其实和降维分析有点类似，都是尽量的将变量进行负责的整合，是的排在前几位的变量组合能够尽可能的解释更多的信息。

排序分析可以单独分析群落之间的关系，也可以分析群落与其环境之间的关系，只使用群落物种组成的排序分析成为间接排序，而使用物种和环境因子进行的排序分析成为直接排序。

CCA和RDA都属于直接排序分析，区别是使用的排序方法不一样，CCA使用的是单峰模型，而RDA使用的是线性模型。

CCA对应的间接排序方法是对应分析 (Correspondence analysis, CA)，RDA对应的间接排序方式是主成分分析 (Principal components analysis, PCA)。

在间接排序中的结果中，坐标轴是分析变量也就是物种的复杂函数组合，而在直接排序中结果中，坐标轴是环境因子的复杂函数组合。

方法选择

我们在阅读文献的时候，会发现有的文章使用的是CCA，而有的文章使用的是RDA，那么对于初学者来说就会产生一种困惑，到底我是使用CCA还是RDA呢？

在进行排序分析之前，我们可以先对物种群落数据进行去趋势的对应分析 (Detrended correspondence analysis, DCA)，根据结果中Lengths of gradient的数值来进行判断。

结果会给出4个Lengths of gradient的数值，如果其中最大的数值大于4，则应选择CCA，如果最大的数值小于3，则选择RDA，如果最大的数值在3-4之间，则两种分析方法都可以。

但是这种标准并不是100%合适，在实际的使用中，我们最好是同时进行CCA和RDA，根据结果进行选择。

在进行结果判断时，最主要是要看做出来的结果是否出现了“弓形效应”，比如下面这个图，样本点的分布想不想一个拉开的弓。(手头没有更好的弓型效应的结果了，大家凑合看一看吧)

一般来说如果环境因子的梯度范围较小，单峰模型和线性模型的结果差别不大，但如果环境因子的梯度范围较大，那线性模型就可能不太合适。

结果解释

CCA和RDA的结果图中使用点代表不同的样本，从原点发出的箭头代表不同的环境因子。

箭头的长度代表该环境因子对群落变化影响的强度，箭头的长度越长，表示环境因子的影响越大。

箭头与坐标轴的夹角代表该环境因子与坐标轴的相关性，夹角越小，代表相关性越高。

样本点到环境因子箭头极其延长线的垂直距离表示环境因子对样本的影响强度，样本点与箭头距离越近，该环境因子对样本的作用越强。

样本位于箭头同方向，表示环境因子与样本物种群落的变化正相关，样本位于箭头的反方向，表示环境因子与样本物种群落的变化负相关。

图像中坐标轴标签中的数值，代表了坐标轴所代表的环境因子组合对物种群落变化的解释比例。

使用蒙特卡洛置换检验，可以得到环境因子与物种群落是否显著相关，如果结果中p大于0.05，则表明该环境因子对物种群落变化的解释并不显著，也就是不是主要的影响因素。

分析过程

DCA

首先要进行的第一步就是DCA分析，以判断后续是使用CCA还是RDA，分析使用vegan包。

#载入样本物种群落数据，sampledata "otu_table.txt", head = TRUE, row.names=1,sep="\t")#分析要求行为样本，列为物种sampledata #载入分析包library(vegan)#进行DCA分析dca #调出Lengths of gradient值并保存dca1 $rproj[,1])dca2 $rproj[,2])dca3 $rproj[,3])dca4 $rproj[,4])GL rownames(GL) "Gradient length")write.csv(GL, file = "dca.csv")

生成的dca.csv文件中即为Lengths of gradient的结果数值。

CCA

CCA同样使用vegan包进行分析，此时需要载入环境因子数据文件。

#环境因子同样要求行为样本，列为因子env "environment.txt", header=TRUE, row.names=1,sep="\t")#进行CCA分析，并保存结果cca TRUE)ccascore write.csv(ccascore$sites, file = "cca.sample.csv")write.csv(cca$CCA$biplot, file = "cca.env.csv")write.csv(ccascore$species, file = "cca.species.csv")

这三个结果文件分别为样本、环境因子和物种在不同坐标轴上对应的具体数值。

得到CCA的结果之后，使用ggplot2包进行结果图的绘制，此时需要载入样本分组文件。

⚠️在进行cca分析时，设置了scale = TRUE，这是为了消除不同环境因子数值之间的固有差异，如果环境因子由于单位不同导致其数值之间差别很大，就一定要设置scale = TRUE，对其进行标准化。

#载入分组文件group "groups.txt", header = FALSE,  colClasses=c("character"))group as.list(group)#载入绘图数据包library(ggrepel)library(ggplot2)#制作绘图所需文件及相关参数文件pch 21:24)col "#E69F00", "#56B4E9", "#009E73", "#F0E442")CCAE as.data.frame(cca$CCA$biplot[,1:2])CCAS1 1]CCAS2 2]plotdata colnames(plotdata) "sample","CCAS1","CCAS2","group")cca1 1]/sum(cca$CCA$eig)*100,2)cca2 2]/sum(cca$CCA$eig)*100,2)#绘制图像plot_CCA  geom_label_repel(aes(CCAS1, CCAS2, label = sample), fill = "white", color = "black", box.padding = unit(0.6,"lines"), segment.colour = "grey50", label.padding = unit(0.35,"lines")) +  geom_point(aes(fill = group, shape = group),size = 8) +  scale_shape_manual(values = pch)+ scale_fill_manual(values = col)+ labs(title = "CCA Plot") +  xlab(paste("CCA1 ( ",cca1,"%"," )", sep = "")) +  ylab(paste("CCA2 ( ",cca2,"%"," )", sep = "")) + geom_segment(data = CCAE, aes(x = 0, y = 0, xend = CCAE[,1], yend = CCAE[,2]), colour = "black", size = 1.2, arrow = arrow(angle = 30, length = unit(0.4, "cm"))) + geom_label_repel(data = CCAE, fill = "grey90",segment.colour = "black", aes(x = CCAE[,1], y = CCAE[,2], label = rownames(CCAE))) + geom_vline(aes(xintercept = 0), linetype = "dotted") + geom_hline(aes(yintercept = 0), linetype = "dotted") + theme(panel.background = element_rect(fill = "white", colour = "black"),  panel.grid = element_blank(), axis.title = element_text(color = "black", size = 18), axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color = "black"), axis.line = element_line(colour = "black"), axis.title.x = element_text(colour = "black", size = 18), axis.title.y = element_text(colour="black", size = 18), axis.text = element_text(colour = "black", size = 18), legend.title = element_blank(), legend.text = element_text(size = 18), legend.key = element_blank(), plot.title = element_text(size = 22, colour = "black",  face = "bold", hjust = 0.5)) #保存图像png(filename = "CCA.png", res = 600, height = 5400, width = 7200)plot_CCAdev.off()

在得到结果图之后，我们来进行最后一步，使用蒙特卡洛置换检验评估环境因子与物种群落组成相关的显著性。

envfit 999) r as.matrix(envfit$vectors$r) p as.matrix(envfit$vectors$pvals)env.p colnames(env.p) "r2","p-value")write.csv(as.data.frame(env.p), file = "ccaenvfit.csv")

得到的结果文件中包含的每一个环境因子的相关系数R2和显著性p-value。

RDA

RDA的整个流程与CCA基本上一致，这里就直接给出完整的代码了。

rda TRUE)rdascore write.csv(rdascore$sites,file="rda.sample.csv")write.csv(rda$CCA$biplot,file="rda.env.csv")write.csv(rdascore$species,file="rda.species.csv")RDAE as.data.frame(rda$CCA$biplot[,1:2])RDAS1 1]RDAS2 2]plotdata colnames(plotdata) "sample","RDAS1","RDAS2","group")rda1 1]/sum(rda$CCA$eig)*100,2)rda2 2]/sum(rda$CCA$eig)*100,2)#RDA plot plot_RDA geom_label_repel(aes(RDAS1, RDAS2, label = sample), fill = "white", color = "black", box.padding = unit(0.6,"lines"), segment.colour = "grey50", label.padding = unit(0.35,"lines")) +  geom_point(aes(fill = group, shape = group),size = 8) +  scale_shape_manual(values = pch)+ scale_fill_manual(values = col)+ labs(title = "RDA Plot") +  xlab(paste("CCA1 ( ",rda1,"%"," )", sep = "")) +  ylab(paste("CCA2 ( ",rda2,"%"," )", sep = "")) + geom_segment(data = RDAE, aes(x = 0, y = 0, xend = RDAE[,1], yend = RDAE[,2]), colour = "black", size = 1.2, arrow = arrow(angle = 30, length = unit(0.4, "cm"))) + geom_label_repel(data = RDAE, fill = "grey90",segment.colour = "black", aes(x = RDAE[,1], y = RDAE[,2], label = rownames(RDAE))) + geom_vline(aes(xintercept = 0), linetype = "dotted") + geom_hline(aes(yintercept = 0), linetype = "dotted") + theme(panel.background = element_rect(fill = "white", colour = "black"),  panel.grid = element_blank(), axis.title = element_text(color = "black", size = 18), axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color = "black"), axis.line = element_line(colour = "black"), axis.title.x = element_text(colour = "black", size = 18), axis.title.y = element_text(colour="black", size = 18), axis.text = element_text(colour = "black", size = 18), legend.title = element_blank(), legend.text = element_text(size = 18), legend.key = element_blank(), plot.title = element_text(size = 22, colour = "black",  face = "bold", hjust = 0.5))png(filename="RDA.png",res=600,height=5400,width=7200)plot_RDAdev.off()#Permutation test envfit 999)r as.matrix(envfit$vectors$r)p as.matrix(envfit$vectors$pvals)env.p colnames(env.p) "r2","p-value")write.csv(as.data.frame(env.p),file="rdaenvfit.csv")

从结果可以看出，这套数据使用RDA并不合适，一方面是弓型效应，另一方面是代表环境因子的箭头太短了。

写在后面

CCA和RDA无法给出单个环境因子对物种群落变化的解释比例，要想知道单个因子对物种群落变化的解释，需要使用pCCA/pRDA分析，或者叫VPA，这个方法的实现我在之后的推文中会介绍。

不过我目前只是能够实现这个分析，但是对具体的原理还是不是很清楚，经常会得到一些解释比例为负值的情况，所以什么时候更新这一项分析还不清楚😄

有对这个分析原理清楚的大神，希望能加我微信教教我😊

关注公众号“红皇后学术”，后台回复“CCA”或“RDA”获取示例文件和完整代码。

扩展阅读

• 高通量测序基础知识
• 转录组测序技术和结果解读
• 红皇后学术文献解读列表
• 基本分子生物学实验
• PAST：最简便易用的统计学分析软件教程目录
• 每天学习一点R系列
• 微生物群落数据分析系列教程
• 微生物研究相关工具
• 微生物研究投稿期刊简介