用R做群落α、β多样性分析

最新推荐文章于 2024-12-26 16:42:01 发布

聂小七f

最新推荐文章于 2024-12-26 16:42:01 发布

阅读量2.2w

点赞数 48

文章标签： r语言开发语言

本文链接：https://blog.csdn.net/qq_35823062/article/details/123315832

版权

做完抽平以后就是我们最关心的多样性分析，其实所有过程都能在qiime2里面做。但是，当我们只看群落中某一科的多样性变化时，我发现qiime2就不再那么灵活。

α多样性是指一个特定区域或生态系统内的多样性，是反映丰富度和均匀度的综合指标。 α多样性主要与两个因素有关：一是种类数目，即丰富度；二是多样性，群落中个体分配上的均匀性。α多样性主要关注局域均匀生境下的物种数目，因此对它的分析也是最直观。平时在分析多样性时，通常是先分析α多样性，如果没有显著差异，我们在分析β多样性，或者将二者结合分析。

###α多样性计算
#清除所有变量
rm(list=ls())
#加载vegan包
library(vegan)
#读入物种数据
otu<-read.table('totalcp.txt',header = T,row.names = 1,check.names=F)
#Shannon 指数
Shannon<-diversity(otu, index = "shannon", MARGIN = 2, base = exp(1))
#Simpson 指数
Simpson<-diversity(otu, index = "simpson", MARGIN = 2, base = exp(1))
#Richness 指数
Richness <- specnumber(otu,MARGIN = 2)
#合并
index<-as.data.frame(cbind(Shannon,Simpson,Richness))
#接下来分析的多样性指数一般不作为重点分析对象，但既然要写，就整理的完整一些
#转置物种数据
totu<-t(otu)
totu<-ceiling(as.data.frame(t(otu)))
#多样性指数
obs_chao1_ace<-t(estimateR(totu))
obs_chao1_ace<-obs_chao1_ace[rownames(index),]
index$Chao1<-obs_chao1_ace[,2]
index$Ace<-obs_chao1_ace[,4]
index$Sobs<-obs_chao1_ace[,1]
index$Pielou <- Shannon / log(Richness, 2)
index$Goods_coverage <- 1 - colSums(otu == 1) / colSums(otu)
#合并、导出数据
write.table(cbind(sample=c(rownames(index)),index),'totalcp.alpha.txt',row.names = F,sep = '\t',quote = F)

β多样性又称生境间的多样性（between-habitat diversity），是指沿环境梯度不同生境群落之间物种组成的相异性或物种沿环境梯度的更替速率。控制β多样性的主要生态因子有土壤、地貌及干扰等。β多样性计算不再是直观地依据物种数量来，而是均以群落相似（或相异）程度为基础，即不同群落之间的距离。常用的距离指数包括Bray-Crutis距离、Unifrac距离等。而具体的分析方法，通常包括PCoA、NMDS等分析。

今天仅以基于Bray-Crutis距离的NMDS分析为例，计算β多样性。

###NMDS分析
#加载包，ape是计算bray距离的包
library(ggplot2)
library(vegan)
library(ape)
#导入群落数据
sample <- read.table("juke.txt",sep = "\t",header = T)
rownames(sample)<-sample[,1]
sample<-sample[,-1]
sample<-data.frame(t(sample))
#计算距离指数
bray<-vegdist(sample,method="bray")
bray<-as.matrix(bray)
#导入群落分组信息
group <- read.table("env.txt",sep = "\t",row.names=1,header = T)
nmds<-metaMDS(bray,k = 2)
summary(nmds)
###提取数据
#NMDS的压迫系数，大于0.2表明不适用NMDS分析
nmds.stress <- nmds$stress
nmds.point <- data.frame(nmds$point)
#导出NMDS的1、2轴，留作后续分析使用，这一步可不做
write.table (nmds.point, file ="totalcp_NMDS12.csv",sep =",", quote =FALSE) 
nmds.species <- data.frame(nmds$species)
sample_site <- nmds.point[1:2]
sample_site$names <- rownames(sample_site)
colnames(sample_site)[1:2] <- c('NMDS1', 'NMDS2')
#合并分组数据
sample_site <- cbind(sample_site,group)
#NMDS图绘制
P1<-ggplot(data=sample_site,aes(NMDS1,NMDS2))+theme_bw()+theme(panel.grid= element_line(color =NA),
                                                               panel.grid.minor = element_line(color = NA),
                                                               panel.border = element_rect(fill = NA, colour ="black"),
                                                               axis.text.x  = element_text(size=15,family="A", colour="black",hjust = 0.7),
                                                               axis.title.x = element_text(vjust=0.2, size = 15,family="A"),
                                                               axis.text.y  = element_text(size=15,family="A", colour="black",hjust = 0.7),
                                                               axis.title.y = element_text(vjust=1, size = 15,family="A", colour="black"))+
  geom_point(aes(color=Biodiversity),size=2,alpha=0.9)+theme(legend.position= "top")+theme(panel.grid=element_blank())+
#下面的#00FF00等分别是颜色的16进制代码，可自己百度，需要注意的是，点的颜色用scale_color_manual（）函数，置信椭圆用scale_fill_manual（）
  scale_color_manual(values=c(R1 = "#00FF00", R2 = "#F74ED6", R4 = "#AD07E3"))+
  stat_ellipse(data=sample_site,geom = "polygon",aes(fill=Biodiversity),alpha=0.35,level = 0.95)+
  scale_fill_manual(values=c(R1 = "#00FF00", R2 = "#F74ED6", R4 = "#AD07E3"))
P1

如果以上代码不能满足大家的需求，可以评论或给我私信。我再更新PCoA等方法。如有错误欢迎大家指出。谢谢。