基于R語(yǔ)言的微生物群落組成多樣性分析——CCA分析

????之前文章中我們講到過(guò)冗余分析（redundancy analysis, RDA），今天我們來(lái)講另一種分析——典范對(duì)應(yīng)分析（canonical correspondence analysis, CCA），這是一種基于對(duì)應(yīng)分析（correspondence analysis, CA）發(fā)展而來(lái)的排序方法，又稱多元直接梯度分析，是研究?jī)山M變量之間相關(guān)關(guān)系的一種多元統(tǒng)計(jì)方法，它能夠揭示出兩組變量之間的內(nèi)在聯(lián)系。
????講到這兒也許很多同學(xué)會(huì)有疑問(wèn)：我怎么知道自己到底是選擇RDA分析還是CCA分析？其實(shí)RDA 和CCA 模型的選擇原則很簡(jiǎn)單，RDA分析一般是基于線性模型的，而CCA分析是基于單峰模型的。當(dāng)拿到數(shù)據(jù)之后呢，我們可以先對(duì)數(shù)據(jù)做DCA（detrendedcorrespondence analysis） 分析，然后我們根據(jù)DCA分析結(jié)果中DCA1的Axis Lengths值的大小進(jìn)行選擇，如果該值大于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、加載數(shù)據(jù)

#OTU表格
df <- read.table("otu.txt",sep="\t",header = T,row.names = 1,check.names = F)
df <-data.frame(t(df))
#環(huán)境因子數(shù)據(jù)
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()函數(shù)進(jìn)行CCA分析
df_otu_cca <- cca(df~., env)
#查看CCA結(jié)果信息，以 I 型標(biāo)尺為例，具體見參考文章
df_otu_cca.scaling1 <- summary(df_otu_cca, scaling = 1)

4、R2及P值校正、約束軸置換檢驗(yàn)

1）R2值校正

R2 <- RsquareAdj(df_otu_cca)
df_otu_cca_noadj <- R2$r.squared  #原R2
df_otu_cca_adj <- R2$adj.r.squared  #校正R2
#計(jì)算校正 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）約束軸的置換檢驗(yàn)及P值校正

## 置換檢驗(yàn)##
# 所有約束軸的置換檢驗(yàn)，即全局檢驗(yàn)，基于 999 次置換，詳情 ?anova.cca
df_otu_cca_test <- anova.cca(df_otu_cca, permutations = 999)
# 各約束軸逐一檢驗(yàn)，基于 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、提取作圖數(shù)據(jù)

###提取作圖數(shù)據(jù)
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")
#將繪圖數(shù)據(jù)和分組合并
df_otu_cca_sites <- merge(df_otu_cca_sites,group,by="samples")

6、繪圖

1）CCA散點(diǎn)圖繪制

color=c("#1597A5","#FFC24B","#FEB3AE") #顏色變量
p1<-ggplot(data=df_otu_cca_sites,aes(x=CCA1,y=CCA2,
                           color=group))+#指定數(shù)據(jù)、X軸、Y軸，顏色
  theme_bw()+#主題設(shè)置
  geom_point(size=3,shape=16)+#繪制點(diǎn)圖并設(shè)定大小
  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)+#添加數(shù)據(jù)點(diǎn)的標(biāo)簽
  # guides(color=guide_legend(title=NULL))+#去除圖例標(biāo)題
  labs(x=CCA1,y=CCA2)+#將x、y軸標(biāo)題改為貢獻(xiàn)度
  stat_ellipse(data=df_otu_cca_sites,
               level=0.95,
               linetype = 2,size=0.8,
               show.legend = T)+
  scale_color_manual(values = color) +#點(diǎn)的顏色設(shè)置
  scale_fill_manual(values = c("#1597A5","#FFC24B","#FEB3AE"))+
  theme(axis.title.x=element_text(size=12),#修改X軸標(biāo)題文本
        axis.title.y=element_text(size=12,angle=90),#修改y軸標(biāo)題文本
        axis.text.y=element_text(size=10),#修改x軸刻度標(biāo)簽文本
        axis.text.x=element_text(size=10),#修改y軸刻度標(biāo)簽文本
        panel.grid=element_blank())#隱藏網(wǎng)格線
p1
image.png

2）添加環(huán)境因子數(shù)據(jù)

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、環(huán)境因子與群落結(jié)構(gòu)差異性分析

1）顯著性計(jì)算

#描述統(tǒng)計(jì)
data<-summary(df_otu_cca)
#檢驗(yàn)環(huán)境因子相關(guān)顯著性（Monte Carlo permutation test）
df_permutest <- permutest(df_otu_cca,permu=999) # permu=999是表示置換循環(huán)的次數(shù)
#每個(gè)環(huán)境因子顯著性檢驗(yàn)
df_envfit <- envfit(df_otu_cca,env,permu=999)
#數(shù)據(jù)處理
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標(biāo)記為FALSE，p>0.05標(biāo)記為TRUE，使用此數(shù)據(jù)繪制柱形圖。

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, # 坐標(biāo)軸的粗細(xì)
              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

源碼及作圖數(shù)據(jù)可在后臺(tái)回復(fù)\color{red}{“CCA”}獲?。。?！

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容