RDA環(huán)境因子群落結構統(tǒng)計檢驗可視化環(huán)境因子的篩選及數(shù)據(jù)的轉化方面請參閱宏基因組公眾號之前的推文,，本文主要側重統(tǒng)計分析與可視化 看到師兄文章里的圖自己可能用到，想復現(xiàn)一下，于是就嘗試了一下，順便寫個推文記錄，在黃靜同學的幫助下完成 百度云鏈接：https://pan.baidu.com/s/1vNmxPV5kw51Ek_oXMmNycw
提取碼：no73 因公眾號文章不可修改，如以上鏈接失效，或想獲取代碼的更新版,，請在“宏基因組”公眾號后臺回復本文關鍵字“rdaenv”獲取最新下載地址。 
rm(list=ls())
library(pacman)
p_load(ggplot2,patchwork,vegan,geosphere,psych,corrplot,
permute,lattice,ggpubr,RColorBrewer,tidyverse,graphics)data=read.csv("sum_c.csv",row.names = 1)#讀入物種(以Phylum水平為例)矩陣表
head(data,n=3)
env=read.csv("env.csv",row.names = 1)#讀入環(huán)境因子數(shù)據(jù)(示例為隨機數(shù))
head(env,n=3)print(decorana(t(data)))
#DCA分析,，根據(jù)Axis lengths行的第一個值選擇排序分析模型
#Axis Lengths >4.0-CCA(基于單峰模型,，典范對應分析)；
#如果在3.0-4.0之間-RDA/CCA均可,；
#如果小于3.0-RDA(基于線性模型,，冗余分析)B.rda=rda(t(data),env[-1],scale = T)#RDA分析，如果沒有環(huán)境因子參數(shù),，就是PCA分析
#提取樣本得分
B.rda.data=data.frame(B.rda$CCA$u[,1:2],
env$Treat,rep(c("Mar","Apr","May","Jun","Jul","Aug"),each = 3))#為了仿師兄的圖添加的
colnames(B.rda.data)=c("RDA1","RDA2","group","Month")
head(B.rda.data,n=3)#提取物種得分
B.rda.spe=data.frame(B.rda$CCA$v[,1:2])
B.rda.spe=as.data.frame(B.rda.spe)
B.rda.spe$Species<-rownames(B.rda.spe)
head(B.rda.spe,n=3)#提取環(huán)境因子得分
B.rda.env <- B.rda$CCA$biplot[,1:2]
B.rda.env <- as.data.frame(B.rda.env)
head(B.rda.env,n=3)#帶有環(huán)境因子,，物種信息且進行不同月份不同處理標記的RDA圖(仿師兄)
yanse<-c("darkolivegreen3","gold","dodgerblue4","darkseagreen",
"chartreuse4","darkorange","burlywood2","brown3","#984EA3","cyan3")
p1=ggplot(data=B.rda.data,aes(RDA1,RDA2))+
geom_point(aes(color=Month,fill=Month,shape=group),size=5)+
scale_color_manual(values=yanse)+
scale_shape_manual(values = c(21,22,23))+
scale_fill_manual(values = yanse)+
geom_point(data=B.rda.spe,aes(RDA1,RDA2),pch=8,size=5)+
geom_text(data=B.rda.spe,aes(x=B.rda.spe[,1],y=B.rda.spe[,2],label=Species),
size=5.5,colour="black",hjust=0.5,vjust=1)+
labs(title = "B RDA plot",
x=paste("RDA1",round(B.rda$CCA$eig[1]/sum(B.rda$CCA$eig)*100,2)," %"),
y=paste("RDA2",round(B.rda$CCA$eig[2]/sum(B.rda$CCA$eig)*100,2)," %"))+
geom_hline(yintercept = 0,lty=3)+
geom_vline(xintercept = 0,lty=3)+
geom_segment(data=B.rda.env,aes(x=0,y=0,xend=B.rda.env[,1],yend=B.rda.env[,2]),
colour="blue",size=0.8,arrow=arrow(angle = 35,length=unit(0.3,"cm")))+
geom_text(data=B.rda.env,aes(x=B.rda.env[,1],y=B.rda.env[,2],
label=rownames(B.rda.env)),size=6.5,colour="blue",
hjust=(1-sign(B.rda.env[,1]))/2,angle=(180/pi)*atan(B.rda.env[,2]/B.rda.env[,1]))+
ggprism::theme_prism()
p1
#統(tǒng)計
B.sum=summary(B.rda)
B.sum$constr.chi/B.sum$tot.chi #constrained表示環(huán)境因子對群落結構差異的解釋度
B.sum$unconst.chi/B.sum$tot.chi#unconstrained表示環(huán)境因子對群落結構不能解釋的部分#環(huán)境因子對群落結構差異解釋量的餅圖繪制
cor_data=data.frame(row.names = c("Explained","Unexplained"),
B=c(B.sum$constr.chi/B.sum$tot.chi,B.sum$unconst.chi/B.sum$tot.chi))
cor_data$group=rownames(cor_data)
head(cor_data,n=3)
cor_data <- data.frame(cor_data)
cor_data=arrange(cor_data,B)
head(cor_data,n=3)
labs<-paste0(cor_data$group,"\n(",round(cor_data$B/sum(cor_data$B)*100,2),"%)")
pie(cor_data$B,labels=labs,init.angle = 90,col=brewer.pal(nrow(cor_data),"Reds"),boder="black")
#在R中手動導出，右側出圖區(qū)Export-PDF
#anova.cca檢驗
B.perm=permutest(B.rda,permu=999) # permu=999是表示置換999次
B.perm #是做環(huán)境因子整體與群落結構差異的相關性(解釋量),，anova.cca {vegan}#envfit檢驗 envfit函數(shù)跟mantel(見下文)的功能是一樣的
B.ef=envfit(B.rda,env[-1],permu=999)#是做每一個環(huán)境因子與群落結構差異的相關性(解釋量)
B.ef$vectors$r#R2值
B.ef$vectors$pvals#P值#每一個環(huán)境因子對群落結構差異解釋量的柱形圖繪制的數(shù)據(jù)整理
cor_com=data.frame(tax=rownames(B.rda.env),B.r=B.ef$vectors$r,B.p=B.ef$vectors$pvals)
cor_com=arrange(cor_com,B.r)
head(cor_com,n=3)
cor_com[c(3)]=cor_com[c(3)]>0.05
head(cor_com,n=3)
#將p<0.05標記為FALSE,，p>0.05標記為TRUE，使用此數(shù)據(jù)繪制柱形圖,，將其可視化#envfit檢驗可視化
cor_com$tax = factor(cor_com$tax,order = T,levels = row.names(cor_com))#按R2值排序
p2 <- ggplot(cor_com, aes(x =tax, y = B.r),size=2) +
geom_bar(stat = 'identity', width = 0.8,color="black",fill="red") +
scale_fill_manual(guide = FALSE)+
geom_text(aes(y = B.r+0.005, label = ifelse(B.p==TRUE,"","*")),
size = 5, fontface = "bold") +
xlab("Environmental factor")+
ylab(expression(r^"2"))+
scale_y_continuous(expand = c(0,0))+
ggprism::theme_prism()+
theme(axis.text.x = element_text(angle = 45))
p2
#mantel檢驗
data <- as.data.frame(t(data))
species.distance<-vegdist(data,method = 'bray')
soil <- NULL
for (i in 2:ncol(env)) {
dd <- mantel(species.distance, vegdist(env[,i], method = "euclidean"),
method = "pearson", permutations = 9999, na.rm = TRUE)
soil <- rbind(soil,c(colnames(env)[i],dd$statistic, dd$signif))
}
head(soil,n=3)
soil <- data.frame(aa=rownames(B.rda.env),M_r=soil[,2],M.p=soil[,3])
rownames(soil)=soil$aa
soil=arrange(soil,M_r)
soil[c(3)]=soil[c(3)]>0.05 # 將p<0.05標記為FALSE,，p>0.05標記為TRUE，同上
soil$aa = factor(soil$aa,order = T,levels = row.names(soil))
soil$M_r=round(as.numeric(soil$M_r),4)
head(soil,n=3)#mantel檢驗可視化
p3 <- ggplot(soil, aes(x =aa, y = M_r),size=2) +
geom_bar(stat = 'identity', width = 0.8,color="black",fill="red") +
scale_fill_manual(guide = FALSE)+
geom_text(aes(y = M_r+0.005, label = ifelse(M.p==TRUE,"","*")),
size = 5, fontface = "bold") +
xlab("Environmental factor")+
ylab(expression(r))+
ggprism::theme_prism()+
theme(axis.text.x = element_text(angle = 45))
p3
##################群落結構差異的統(tǒng)計檢驗(三種方法及可視化)
####Adonis
otu <- data.frame(data)#mantel檢驗時已經轉置,，勿要轉置
head(otu)
#樣本分組文件
group <- read.delim('metadata.txt', sep = '\t', stringsAsFactors = FALSE)
head(group)
#使用 Bray-Curtis 距離測度 unifrac_binary
adonis_result <- adonis(otu~Group, group, distance = 'Bray-Curtis', permutations = 999)
adonis_result$aov.tabgroup_name <- unique(group$Group)
adonis_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {
for (j in (i + 1):length(group_name)) {
group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))
otu_ij <- otu[group_ij$SampleID, ]
adonis_result_otu_ij <- adonis(otu_ij~Group, group_ij, permutations = 999, distance = 'bray')
adonis_result_two <- rbind(adonis_result_two, c(paste(group_name[i], group_name[j], sep = '_'), 'Bray-Curtis', unlist(data.frame(adonis_result_otu_ij$aov.tab, check.names = FALSE)[1, ])))
}
}
adonis_result_two <- data.frame(adonis_result_two, stringsAsFactors = FALSE)
names(adonis_result_two) <- c('group', 'distance', 'Df', 'Sums of squares', 'Mean squares', 'F.Model', 'R2', 'P')
adonis_result_two$R2<- as.numeric(adonis_result_two$R2)
adonis_result_two$P <- as.numeric(adonis_result_two$P)
#p值Benjamini校正 作圖時自己選擇用哪一個P值(下同)
adonis_result_two$P_adj_BH <- p.adjust(adonis_result_two$'P', method = 'BH')
head(adonis_result_two)
adonis_result_two=arrange(adonis_result_two,R2)
adonis_result_two
adonis_result_two$P=adonis_result_two$P>0.05#將p<0.05標記為FALSE,，p>0.05標記為TRUE，使用此數(shù)據(jù)繪制柱形圖(下同)
adonis_result_two$tax = factor(adonis_result_two$group,order = T,levels = adonis_result_two$group)
adonis <- ggplot(adonis_result_two, aes(x =tax, y = R2),size=2) +
geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+
scale_fill_manual(guide = FALSE)+
geom_text(aes(y = R2-0.02, label = ifelse(P==TRUE,"","*")),#選擇作圖的P值(P_adj_BH)(下同)
size = 5, fontface = "bold") +
xlab("")+
ylab(expression(r^"2"))+
scale_y_continuous(expand = c(0,0))+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
adonis ##圖片大小在保存時自己調合適
####Anosim
anosim_result <- anosim(otu, group$Group, distance = 'bray', permutations = 999)
anosim_result$signif #p 值
anosim_result$statistic #R 值group_name <- unique(group$Group)
anosim_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {
for (j in (i + 1):length(group_name)) {
group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))
otu_ij <- otu[group_ij$SampleID,]
anosim_result_otu_ij <- anosim(otu_ij, group_ij$Group, permutations = 999, distance = 'bray')
anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep = '_'), 'Bray-Curtis', anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif))
}
}
anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)
names(anosim_result_two) <- c('group', 'distance', 'R', 'P')
anosim_result_two$R<- as.numeric(anosim_result_two$R)
anosim_result_two$P <- as.numeric(anosim_result_two$P)
anosim_result_two$P_adj_BH <- p.adjust(anosim_result_two$P, method = 'BH')
head(anosim_result_two)
anosim_result_two=arrange(anosim_result_two,R)
anosim_result_two
anosim_result_two$P=anosim_result_two$P>0.05
anosim_result_two$tax = factor(anosim_result_two$group,order = T,levels = anosim_result_two$group)
anosim <- ggplot(anosim_result_two, aes(x =tax, y = R),size=2) +
geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+
scale_fill_manual(guide = FALSE)+
geom_text(aes(y =R-0.02, label = ifelse(P==TRUE,"","*")),
size = 5, fontface = "bold") +
xlab("")+
ylab(expression(r))+
scale_y_continuous(expand = c(0,0))+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
anosim
##MRPP
mrpp_result <- mrpp(otu, group$Group, distance = 'bray', permutations = 999)
mrpp_result$Pvalue # p 值
mrpp_result$A #相當于Anosim檢驗的R值roup_name <- unique(group$Group)
mrpp_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {
for (j in (i + 1):length(group_name)) {
group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))
otu_ij <- otu[group_ij$SampleID,]
mrpp_result_otu_ij <- mrpp(otu_ij, group_ij$Group, permutations = 999, distance = 'bray')
mrpp_result_two <- rbind(mrpp_result_two, c(paste(group_name[i], group_name[j], sep = '_'), 'Bray-Curtis', mrpp_result_otu_ij$A, mrpp_result_otu_ij$delta, mrpp_result_otu_ij$E.delta, mrpp_result_otu_ij$Pvalue))
}
}
mrpp_result_two <- data.frame(mrpp_result_two, stringsAsFactors = FALSE)
names(mrpp_result_two) <- c('group', 'distance', 'A', 'Observe_delta', 'Expect_delta', 'P')
mrpp_result_two$A<- as.numeric(mrpp_result_two$A)
mrpp_result_two$P <- as.numeric(mrpp_result_two$P)
mrpp_result_two$P_adj_BH <- p.adjust(mrpp_result_two$P, method = 'BH')
head(mrpp_result_two)
mrpp_result_two=arrange(mrpp_result_two,A)
mrpp_result_two
mrpp_result_two$P=mrpp_result_two$P>0.05
mrpp_result_two$tax = factor(mrpp_result_two$group,order = T,levels = mrpp_result_two$group)
MRPP <- ggplot(mrpp_result_two, aes(x =tax, y = A),size=2) +
geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+
scale_fill_manual(guide = FALSE)+
geom_text(aes(y =A-0.02, label = ifelse(P==TRUE,"","*")),
size = 5, fontface = "bold") +
xlab("")+
ylab(expression(r))+
scale_y_continuous(expand = c(0,0))+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
MRPP # 其余數(shù)據(jù)集重復繪制
#Output figure width and height
# Letter紙圖片尺寸為單欄89 mm,，雙欄183 mm,，頁面最寬為247 mm 推薦比例16：10，
# 即半版89 mm x 56 mm; 183 mm x 114 mm
#
##################保存
ggsave("./p1.pdf", p1, width = 183, height = 114, units = "mm")
ggsave("./p2.pdf", p2, width = 183, height = 114, units = "mm")
ggsave("./p3.pdf", p3, width = 183, height = 114, units = "mm")
ggsave("./p4.pdf", adonis, width = 183, height = 114, units = "mm")
ggsave("./p5.pdf", anosim, width = 183, height = 114, units = "mm")
ggsave("./p6.pdf", MRPP, width = 183, height = 114, units = "mm")參考資料R繪圖-RDA排序分析 R包vegan的置換多元方差分析（PERMANOVA）判斷群落結構差異 R包vegan的相似性分析（ANOSIM）判斷群落結構差異 R包vegan的MRPP分析判斷群落結構差異
|
