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

做完抽平以后就是我们最关心的多样性分析，其实所有过程都能在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等方法。如有错误欢迎大家指出。谢谢。

用R做群落α、β多样性分析相关推荐

在线作图|在线做Unifrac PCoA分析
Unifrac PCoA分析 UniFrac分析利用各样品序列间的进化信息来比较环境样品在特定的进化谱系中是否有显著的微生物群落差异.UniFrac 可用于beta 多样性的评估分析,即对样品两两之间 ...
使用R语言包clusterProfiler做KEGG富集分析时出现的错误及解决方法
使用enrichKEGG做通路富集分析时,一直报错:显示No gene can be mapped.... k <- enrichKEGG(gene = gene, organism = &qu ...
用R做中文LDA主题模型可视化分析
LDA主题模型在2002年被David M. Blei.Andrew Y. Ng(是的,就是吴恩达老师)和Michael I. Jordan三位第一次提出,近几年随着社会化媒体的兴起,文本数据成为越来 ...
R数据分析：生存分析的做法和结果解释
今天给大家写写生存分析: Survival analysis corresponds to a set of statistical approaches used to investigate th ...
Ecol. Lett.：写给实践生态学家的β多样性分析指南 | 朝花夕拾
本文转载自"生态学文献分享",已获授权原文信息: Navigating the multiple meanings of β diversity: a roadmap for t ...
图表复现|PRD地下水微生物群落的多样性分析文献
前言之前有小伙伴问道,能否复现一篇关于微生物多样性文献里的图片,今天小编给大家分享一下微生物分析相关的画图方法.这篇文献<Diversity and predictive metabolic ...
R软件中文本分析安装包 Rjava 和 Rwordseg 傻瓜式安装方法四部曲
这两天,由于要做一个文本分析的内容,所以搜索了一天R语言中的可以做文本分析的加载包,但是在安装包的过程,真是被虐千百遍,总是安装不成功.特此专门写一篇博文,把整个心塞史畅快的释放一下. ------- ...
中科院单细胞分析算法开发博士带你做单细胞转录组分析
" 福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组和Python课程的线上直播课.报名参加线上直播课的老师可在1年内选择参加同课程的一次 ...
《数据分析实战》--用R做交叉列表
<数据分析实战>–用R做交叉列表本文参考的是<数据分析实战>第四章. 背景:针对某公司的产品,发现当月的用户使用量减少了很多,但是和上月相比,本月的商业宣传和月度活动并无大的 ...

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

用R做群落α、β多样性分析相关推荐

最新文章

热门文章