基于R语言的微生物群落组成多样性分析——CCA分析

    之前文章中我们讲到过冗余分析（redundancy analysis, RDA），今天我们来讲另一种分析——典范对应分析（canonical correspondence analysis, CCA），这是一种基于对应分析（correspondence analysis, CA）发展而来的排序方法，又称多元直接梯度分析，是研究两组变量之间相关关系的一种多元统计方法，它能够揭示出两组变量之间的内在联系。
    讲到这儿也许很多同学会有疑问：我怎么知道自己到底是选择RDA分析还是CCA分析？其实RDA 和CCA 模型的选择原则很简单，RDA分析一般是基于线性模型的，而CCA分析是基于单峰模型的。当拿到数据之后呢，我们可以先对数据做DCA（detrendedcorrespondence analysis） 分析，然后我们根据DCA分析结果中DCA1的Axis Lengths值的大小进行选择，如果该值大于4.0就选CCA；如果该值在3.0-4.0 之间，选RDA 和CCA都可以；如果该值小于3.0，选择RDA就行。

1、加载包

rm(list=ls())#clear Global Environment
setwd("D:\\桌面\\CCA")
#安装包
install.packages("vegan")
install.packages("ggplot2")
install.packages("ggprism")
install.packages("ggpubr")
#加载包
library(vegan)
library(ggplot2)
library(ggprism)
library(ggpubr)

2、加载数据

#OTU表格
df <- read.table("otu.txt",sep="\t",header = T,row.names = 1,check.names = F)
df <-data.frame(t(df))
#环境因子数据
env <- read.table("env.txt",sep="\t",header = T,row.names = 1,check.names = F)
head(df)
head(env)

image.png
image.png

3、CCA分析

#使用vegan包中的cca()函数进行CCA分析
df_otu_cca <- cca(df~., env)
#查看CCA结果信息，以 I 型标尺为例，具体见参考文章
df_otu_cca.scaling1 <- summary(df_otu_cca, scaling = 1)

4、R2及P值校正、约束轴置换检验

1）R2值校正

R2 <- RsquareAdj(df_otu_cca)
df_otu_cca_noadj <- R2$r.squared  #原R2
df_otu_cca_adj <- R2$adj.r.squared  #校正R2
#计算校正 R2 后的约束轴解释率
df_otu_cca_exp_adj <- df_otu_cca_adj * df_otu_cca$CCA$eig/sum(df_otu_cca$CCA$eig)
CCA1 <- paste("CCA1 (",round(df_otu_cca_exp_adj[1]*100, 1),"%)")
CCA2 <- paste("CCA2 (",round(df_otu_cca_exp_adj[2]*100, 1),"%)")

2）约束轴的置换检验及P值校正

## 置换检验##
# 所有约束轴的置换检验，即全局检验，基于 999 次置换，详情 ?anova.cca
df_otu_cca_test <- anova.cca(df_otu_cca, permutations = 999)
# 各约束轴逐一检验，基于 999 次置换
df_otu_cca_test_axis <- anova.cca(df_otu_cca, by = 'axis', permutations = 999)
# p值校正（Bonferroni为例）
df_otu_cca_test_axis$`Pr(>F)` <- p.adjust(df_otu_cca_test_axis$`Pr(>F)`, method = 'bonferroni')

5、提取作图数据

###提取作图数据
df_otu_cca_sites <- data.frame(df_otu_cca.scaling1$sites)[1:2]
df_otu_cca_env <- data.frame(df_otu_cca.scaling1$biplot)[1:2]
#######添加分组信息
df_otu_cca_sites$samples <- rownames(df_otu_cca_sites)
#读入分组信息
group <- read.table("group.txt", sep='\t', header=T)
#修改列名
colnames(group) <- c("samples","group")
#将绘图数据和分组合并
df_otu_cca_sites <- merge(df_otu_cca_sites,group,by="samples")

6、绘图

1）CCA散点图绘制

color=c("#1597A5","#FFC24B","#FEB3AE") #颜色变量
p1<-ggplot(data=df_otu_cca_sites,aes(x=CCA1,y=CCA2,
                           color=group))+#指定数据、X轴、Y轴，颜色
  theme_bw()+#主题设置
  geom_point(size=3,shape=16)+#绘制点图并设定大小
  theme(panel.grid = element_blank())+
  geom_vline(xintercept = 0,lty="dashed",color = 'black', size = 0.8)+
  geom_hline(yintercept = 0,lty="dashed",color = 'black', size = 0.8)+#图中虚线
  geom_text(aes(label=samples, y=CCA2+0.1,x=CCA1+0.1,  vjust=0),size=3)+#添加数据点的标签
  # guides(color=guide_legend(title=NULL))+#去除图例标题
  labs(x=CCA1,y=CCA2)+#将x、y轴标题改为贡献度
  stat_ellipse(data=df_otu_cca_sites,
               level=0.95,
               linetype = 2,size=0.8,
               show.legend = T)+
  scale_color_manual(values = color) +#点的颜色设置
  scale_fill_manual(values = c("#1597A5","#FFC24B","#FEB3AE"))+
  theme(axis.title.x=element_text(size=12),#修改X轴标题文本
        axis.title.y=element_text(size=12,angle=90),#修改y轴标题文本
        axis.text.y=element_text(size=10),#修改x轴刻度标签文本
        axis.text.x=element_text(size=10),#修改y轴刻度标签文本
        panel.grid=element_blank())#隐藏网格线
p1

image.png

2）添加环境因子数据

p2<-p1+geom_segment(data=df_otu_cca_env,aes(x=0,y=0,xend=CCA1*3,yend=CCA2*3),
                color="red",size=0.8,
                arrow=arrow(angle = 35,length=unit(0.3,"cm")))+
  geom_text(data=df_otu_cca_env,aes(x=CCA1,y=CCA2,
                                label=rownames(df_otu_cca_env)),size=3.5,
            color="blue", 
            hjust=(1-sign(df_otu_cca_env$CCA1))/2,angle=(180/pi)*atan(df_otu_cca_env$CCA2/df_otu_cca_env$CCA1))+
  theme(legend.position = "top")
p2

image.png

7、环境因子与群落结构差异性分析

1）显著性计算

#描述统计
data<-summary(df_otu_cca)
#检验环境因子相关显著性（Monte Carlo permutation test）
df_permutest <- permutest(df_otu_cca,permu=999) # permu=999是表示置换循环的次数
#每个环境因子显著性检验
df_envfit <- envfit(df_otu_cca,env,permu=999)
#数据处理
cor_data<-data.frame(data$constr.chi/data$tot.chi, data$unconst.chi/data$tot.chi)
cor_com <- data.frame(tax=colnames(env),r=df_envfit$vectors$r,p=df_envfit$vectors$pvals)
cor_com[1:5,3]=cor_com[,3]>0.05 # 将p<0.05标记为FALSE，p>0.05标记为TRUE，使用此数据绘制柱形图。

2）绘图

p3 <- ggplot(cor_com,aes(x =tax, y = r),size=2) +
  geom_bar(aes(fill=tax),stat = 'identity', width = 0.8)+
  geom_text(aes(y = r+0.05, label = ifelse(p==T,"","*")),size = 5, fontface = "bold") +
  labs(x = '', y = '')+
  xlab("Environmental factor")+
  ylab(expression(r^"2"))+
  theme_prism(palette = "candy_bright",
              base_fontface = "plain", # 字体样式，可选 bold, plain, italic
              base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
              base_size = 16,  # 图形的字体大小
              base_line_size = 0.8, # 坐标轴的粗细
              axis_text_angle = 45)+ # 可选值有 0，45，90，270
  scale_fill_prism(palette = "candy_bright") 

p3

image.png

3）合并图形

ggarrange(p2,p3,ncol = 2,align="none",heights = c(1,1),widths = c(1,1))

image.png

4）AI美化

image.png

参考：

1）https://copyfuture.com/blogs-details/20200723174028438au5ftbuawf9sdyo
2）https://zhuanlan.zhihu.com/p/399810564?ivk_sa=1024320u

源码及作图数据可在后台回复获取！！！