其他
基于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就行。
加载包
rm(list=ls())#clear Global Environmentsetwd("D:\\桌面\\CCA")#安装包install.packages("vegan")install.packages("ggplot2")install.packages("ggprism")install.packages("ggpubr")#加载包library(vegan)library(ggplot2)library(ggprism)library(ggpubr)加载数据
#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)CCA分析
#使用vegan包中的cca()函数进行CCA分析df_otu_cca <- cca(df~., env)
#查看CCA结果信息，以 I 型标尺为例，具体见参考文章df_otu_cca.scaling1 <- summary(df_otu_cca, scaling = 1)R2及P值校正、约束轴置换检验
1）R2值校正
R2 <- RsquareAdj(df_otu_cca)df_otu_cca_noadj <- R2$r.squared #原R2df_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.ccadf_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')提取作图数据
###提取作图数据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")绘图
1）散点图绘制
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())#隐藏网格线p12）添加环境因子数据
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环境因子与群落结构差异性分析
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")
p33）合并图形
ggarrange(p2,p3,ncol = 2,align="none",heights = c(1,1),widths = c(1,1))4）AI美化
我就知道你“在看”