|||
library(ggrepel) library(vegan) library(ggrepel) phylum=read.csv("phylum.csv",header = T,row.names = 1) env=read.csv("env.csv",header = T,row.names = 1) phylum1=decostand(phylum,method = "hellinger")###对响应变量做hellinger转化 env=log10(env) rda_analysis<-rda(phylum1~.,env) result<-summary(rda_analysis) result sp=as.data.frame(result$species[,1:2])*5###提取相应变量坐标,乘以5是使图美观,不影响分析 st=as.data.frame(result$sites[,1:2])####提取样方坐标,如果不想让样点名称显示可以st=as.data.frame(result$sites[,1:2],row.names = F) yz=as.data.frame(result$biplot[,1:2])###提取解释变量坐标 group=as.data.frame(c(rep("I",10),rep("II",10),rep("III",10)))####创建分组信息 colnames(group)="groups"####将分组列命名为groups p<-ggplot() +geom_text_repel(data = st,aes(RDA1,RDA2,label=row.names(st)),size=4)+geom_point(data = st,aes(RDA1,RDA2,shape=group$groups,fill=group$groups), size=3.5)+scale_shape_manual(values = c(21:23))+ geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), arrow = arrow(angle=22.5,length = unit(0.35,"cm"), type = "closed"),linetype=1, size=0.6,colour = "red")+ 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=22.5,length = unit(0.35,"cm"), type = "closed"),linetype=1, size=0.6,colour = "blue")+ geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+ labs(x="RDA1 58.76%",y="RDA2 3.22%")+ ##RDA1,RDA2的值需要从结果中获得 geom_hline(yintercept=0,linetype=3,size=1) + geom_vline(xintercept=0,linetype=3,size=1)+ theme_bw()+theme(panel.grid=element_blank())
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2023-3-23 20:02
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社