基于R语言的微生物群落组成多样性分析——CCA分析
2022-07-20 本文已影响0人
科研那点事儿
之前文章中我们讲到过冗余分析（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
源码及作图数据可在后台回复获取！！！
