ggplot2 | RDA 分析

library(vegan)

library(ggplot2)

library(rdacca.hp)

library(RColorBrewer)

#读取物种或者其他什么信息表，t()是由表格形式决定的，如果列标题为处理名，就需要这个，行名是处理名就不需要

species <- read.table("D:/aSbreviflora/data/community/Bacteria/species.txt",

                      header=T,row.names = 1,sep = '\t')%>% t()

#输入环境数据

env <- read.table("D:/aSbreviflora/data/env.txt",header=T,row.names = 1,sep = '\t')

#两个表的行号对应

env <-env[match(row.names(species),row.names(env)),]

#验证两个表的行号是否对应

print(rownames(env)==rownames(species))

#自相关性分析，自相关大于10的一般为环境因子自相关，建议去掉，如果很重要，可以选择性保留，具体可以参照大佬的解析

A <- cca(species~W+pH+SOC+TN+AP+NH4+NO3,data=env)

vif.cca(A)

#筛选环境因子

env=subset(env, select=c(W,pH,SOC,AP,NH4,NO3))

#两个表的行号对应

env<-env[match(row.names(species),row.names(env)),]

#判断适用于CCA还是RDA分析

decorana(species)#分析结束后Axis lengths 大于3适合用CCA分析，小于3适合用RDA分析

mycolor <- brewer.pal(12, "Set3")

spp <- decostand(species,method = "hellinger")

rdacca.hp(spp,env, method="RDA",add=F, n.perm=999)

uu=rda(spp~.,env)

ii=summary(uu)

head(ii)

r2<-RsquareAdj(uu)

r2

rda_noadj<-r2$r.squared

rda_adj <- r2$adj.r.squared #校正R2

rda_adj

# 所有约束轴的置换检验，即全局检验，基于 999 次置换,确定p值

df <- anova.cca(uu,by = "margin", permutations = 999)#by='margin'是决定这个检验值是分项的还是总项，具体自己看看大佬的解析

df

sp=as.data.frame(ii$species[,1:2])*5

st=as.data.frame(ii$sites[,1:2])

yz=as.data.frame(ii$biplot[,1:2])

#导入分组数据

design = read.table("D:/aSbreviflora/data/beta/group.txt", header=T, row.names=1, sep="\t")

design$group = factor(design$group, levels=c("North","West","Undisturbed"))

group = design$group

library(ggrepel)#运行geom_text_repel()

ggplot() +

  geom_point(data = st,aes(RDA1,RDA2,colour=group,shape = group),size=4,alpha=1)+

  scale_color_manual(values = c("#ff4400","#55aa00","#00AAFF"))+

  scale_shape_manual(values = c(16,16,16))+

  geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2),

              arrow = arrow(angle=22.5,length = unit(0.35,"cm"), type = "closed"),

              linetype=1, linewidth=0.8,colour = "#FFaa00")+

  geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+

  labs(x=paste("RDA 1 (", format(100 *ii$cont[[1]][2,1], digits=4), "%)", sep=""),

      y=paste("RDA 2 (", format(100 *ii$cont[[1]][2,2], digits=4), "%)", sep=""))+

  geom_hline(yintercept=0,linetype=3,linewidth=1) +

  geom_vline(xintercept=0,linetype=3,linewidth=1)+

  theme_bw()+theme(panel.grid=element_blank())···```