# 首先要安装devtools包，仅需安装一次
install.packages("devtools")
# 加载devtools包
library(devtools)
# 下载ggvegan包
devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
otu.tab <- read.csv("otutab.txt", row.names = 1, header=T, sep="\t")
env.data <- read.csv("new_meta.txt", row.names = 1, fill = T, header=T, sep="\t")
#transform data
otu <- t(otu.tab)
#data normolization (Legendre and Gallagher,2001)
##by log 环境因子进行log转化，以减少同一种环境因子之间本身数值大小造成的影响。
env.data.log <- log1p(env.data)##
##delete NA 删掉空值行
env <- na.omit(env.data.log)
###hellinger transform OTU进行hellinger转化，使数据均一性更好。
otu.hell <- decostand(otu, "hellinger")
#DCA analysis DCA分析看分析结果中Lengths of gradient的第一轴中的Axis lengths大小，
#如果大于4.0，就应该选CCA，如果在3.0-4.0之间，选RDA和CCA均可，如果小于3.0，RDA的结果要好于CCA。
#本例中，我们数据的Axis lengths 大小为0.63859，所有我们应该选择RDA进行分析
sel <- decorana(otu.hell)
sel
otu.tab.0 <- rda(otu.hell ~ 1, env) #no variables
#Axis 第一项大于四应该用CCA分析
otu.tab.1<- rda(otu.hell ~ ., env)
#我们在筛选完RDA和CCA分析后，我们需要对所有环境因子进行共线性分析，利用方差膨胀因子分析
vif.cca(otu.tab.1)
#删除掉共线性的环境因子，删掉最大的变量，直到所有的变量都小于10
otu.tab.1 <- rda(otu.hell ~ N+P+K+Ca+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)
otu.tab.1
vif.cca(otu.tab.1)
#进一步筛选
otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)
vif.cca(otu.tab.1)
#test again
otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env.data.log)
#方差膨胀因子分析,目前所有变量都已经小于10
vif.cca(otu.tab.1)
##用step模型检测最低AIC值
#mod.u <- step(otu.tab.0, scope = formula(otu.tab.1), test = "perm")
# "perm"增加P值等参数
mod.d <- step(otu.tab.0, scope = (list(lower = formula(otu.tab.0),
upper = formula(otu.tab.1))))
mod.d
##本处筛选的结果，找到一个Mg环境因子适合模型构建，为了下一步画图，我们
#保留上述筛选的结果所有非共线性的环境因子
#choose variables for best model and rda analysis again#
(otu.rda.f <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env))
anova(otu.rda.f)
anova(otu.rda.f, by = "term")
anova(otu.rda.f, by = "axis")
#计算db-rda
otu.tab.bray <- vegdist(otu.hell, "bray")
#dbrda() for dbRDA # cca() for CCA
otu.tab.b<- capscale(otu.tab.bray ~ ., env)
##绘制db-RDA图
plot(otu.tab.b ,display=c("si","bp","sp"))
## 用ggvegan绘制RDA图
p1<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = "text",
layers = c("sites"),
legend.position = "right", title = "db-RDA")
p1
p2<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = "text",
layers = c("biplot"),
legend.position = "right", title = "db-RDA")
p2
p3<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = "text",
layers = c("sites","biplot"),
legend.position = "right", title = "db-RDA")
p3
p4<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2),
geom = "text", layers = c( "species","sites",
"biplot", "centroids"), legend.position = "right",
title = "db-RDA")
p4
## 添加图层
p5<-p4 + theme_bw()+theme(panel.grid.major =element_blank(),
panel.grid.minor = element_blank())
p5
###################
###################
library(vegan)
library(tidyverse)
library(ggrepel)
genus <- read.delim("cpeBA.xls",check.names=F,row.names = 1) %>%
t() %>%
decostand(method = "hellinger") %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
arrange(sample) %>% column_to_rownames("sample")
envi <- read.delim("envir.xls",check.names=F) %>%
arrange(sample) %>%
column_to_rownames("sample") %>%
mutate(Depth=as.numeric(Depth)) %>% select(-Chla)
grp <- read.delim("group.txt",check.names = F)
p <- rda(genus~.,data=envi) %>% summary()
p
sp <- as.data.frame(p$species[,1:2])*2
st <- as.data.frame(p$sites[,1:2])
yz <- as.data.frame(p$biplot[,1:2])
colors <- c("#FED439FF","#709AE1FF","#8A9197FF","#D2AF81FF",
"#D5E4A2FF","#197EC0FF","#F05C3BFF","#46732EFF",
"#71D0F5FF","#370335FF","#075149FF","#C80813FF",
"#91331FFF","#1A9993FF","#FD8CC1FF")
ggplot()+
geom_point(data= sp %>% rownames_to_column("genus") %>% as.data.frame(),
aes(RDA1,RDA2,fill=genus,color=genus),
size=4,pch=21,color="white")+
scale_fill_manual(values = colors)+
geom_text_repel(data = sp,aes(RDA1,RDA2),label=row.names(sp))+
geom_segment(data=yz,aes(x=0,y =0,xend=RDA1,yend=RDA2),
arrow=arrow(angle=30,length=unit(0.25,"cm"),
type="closed"),linetype=1,
size=0.8,colour = "#00A08A")+
geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
labs(x=paste("RDA1(",format(100*p$cont[[1]][2,1],
digits=4),"%)",sep=""),
y=paste("RDA2(",format(100*p$cont[[1]][2,2],
digits=4),"%)",sep=""))+
geom_hline(yintercept=0,linetype=3,size=0.8)+
geom_vline(xintercept=0,linetype=3,size=0.8)+
theme_bw()+
theme(panel.grid=element_blank(),
legend.title =element_blank(),
plot.margin=unit(c(0.3,0.3,0.3,0.3),units=,"cm"),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_rect(colour="black",
size=1,fill="white"),
axis.title.y=element_text(size=11,color="black",
margin=margin(r=8),face="plain"),
axis.title.x=element_text(size=11,color="black",
margin=margin(t=8),face="plain"),
axis.text=element_text(size=11,color="black",face="plain"))+
annotate("text",x=0.2,y=0.6,label=expression(NO[3]^-""),
parse = TRUE,size=4,colour="black")+
annotate("text",x=0.2,y=0.4,label=expression(NH[4]^+""),
parse = TRUE,size=4,colour="black")+
annotate("text",x=0.6,y=0.4,label=expression(SiO[4]^4^- ""),
parse = TRUE,size=4,colour="black")+
annotate("text",x=0.7,y=0.5,label=expression(PO[4]^3^- ""),
parse = TRUE,size=4,colour="black")+
annotate("text",x=-0.6,y=0.5,label=expression(NO[2]^- ""),
parse = TRUE,size=4,colour="black")+
theme(legend.position=c(0,1),legend.justification = c(0,1),
legend.background=element_blank(),
legend.box.background=element_rect(colour = "black"),
legend.title = element_blank(), # 图例标题为空
legend.key=element_blank(), # 图例键为空
legend.text = element_text(color="black",size=10), # 定义图例文本
legend.spacing.x=unit(0,'cm'), # 定义文本书平距离
legend.key.width=unit(0.5,'cm'), # 定义图例水平大小
legend.key.height=unit(0.5,'cm'))
