R语言计算α多样性指数与画图

Python022

R语言计算α多样性指数与画图,第1张

操作之前安装好ggplot2、vegan、ggpubr包。如下:

install.packages("ggplot2")

install.packages("ggpubr")

install.packages("vegan")

计算Shannon-香农指数和Simpson-辛普森指数的命令在vegan包中，计算各组显著性的命令在ggpubr包中；画图使用ggplot命令，在行使每个命令之前一定要加载相应的包，如下:

library(ggplot2)

library(ggpubr)

library(vegan)

拿到一个otu表格，要先计算香农指数和辛普森指数，操作如下:

otu=read.table('D:/r-working/feature-table.taxonomy.txt',row.names = 1,skip=1,header=T,comment.char ='',sep='\t')

#读取out表格

#'D:/feature table.taxonomy.txt'为文件路径，注意斜线方向

#row.names = 1指定第一列为行名

#skip=1跳过第一行不读

#header=T指定第一个有效行为列名

#sep='\t'表示指定制表符为分隔符

#comment.char=''表示设置注释符号为空字符‘’，这样#后面的内容就不会被省略

otu=otu[,-ncol(otu)]

#去除表格的最后一列，无用信息

otu=t(otu)

#表格转置，必须将样品名作为行名

shannon=diversity(otu,"shannon")

#计算香农指数,先加载vegan包

shannon

#查看香农指数

simpson=diversity(otu,"simpson")

#计算辛普森指数,先加载vegan包

simpson

#查看辛普森指数

alpha=data.frame(shannon,simpson,check.names=T)

#合并两个指数

write.table(alpha,"D:/r-working/alpha-summary.xls",sep='\t',quote=F)

#存储数据，注意路径使用反斜杠

将各样本进行分组，并进行画图，操作如下:

map<-read.table('D:/r-working/mapping_file.txt',row.names = 1,header=T,comment.char ='',sep='\t',check.names=F)

#读取分组表格

group<-map["Group1"]

#提取需要的分组,'Group1'是表中的分组列名，包括A,B,C三组

alpha<-alpha[match(rownames(group),rownames(alpha)),]

#重排alpha的行的顺序,使其与group的样本id(行名)一致

data<-data.frame(group,alpha,check.rows=T)

#合并两个表格.'<-'与'='同属赋值的含义.

p=ggplot(data=data,aes(x=Group1,y=shannon))+geom_boxplot(fill=rainbow(7)[2])

#data = data指定数据表格

#x=Group1指定作为x轴的数据列名

#y=shannon指定作为y轴的数据列名

#geom_boxplot()表示画箱线图

#fill=rainbow(7)[2]指定填充色

此处用到ggplot2包画箱线图，将画图函数赋值给p后，可以用‘+’不断进行图层叠加，给图片p增加新的特性

p

#查看p

mycompare=list(c('A','B'),c('A','C'),c('B','C'))

#指定多重比较的分组对

mycompare

p<-p+stat_compare_means(comparisons=mycompare,label = "p.signif",method = 'wilcox')

#添加显著性标记的第一种方法，在此之前先加载ggpubr包

p<-p+ylim(2,5.5)

#调整图像的外观

物种累积曲线( species accumulation curves)用于衡量和预测群落中物种丰富度随样本量扩大而增加的幅度，在生物多样性和群落调查中，被广泛用于判断样本量是否充分的并估计群落丰富度。

一般而言，在样本量较少的情况下，随着样本数量的增加，将有较大可能性发现大量新的物种，此时曲线呈急剧上升状态；当样本数量已经较大时，此时群落中的ASV/OTU总数将不再随着样本数量增加而显著增加，曲线将趋于平缓。

因此，通过物种累积曲线可以判断样本量是否充分：若曲线始终保持上升趋势，则表明样本量不足，反之，则表明样本量已足以反应群落的物种组成。在样本量充分的前提下，运用物种累积曲线还可以对物种丰富度进行预测。

下面我使用R脚本，vegan包对ASV/OTU丰度矩阵中每个样本所对应的SAV/OTU总数绘制specaccum物种累计曲线。

1.调用vegan包，读取物种数据；

2.使用 specaccum 函数用来计算物种的累计曲线；

3.作图展示。

如图，随着取样样本数量逐渐增加，所观测到的物种种类也不断增加。当曲线趋近平缓时，代表群落中的物种接近全部被观测到；反之不饱和，还需继续观测更多的样本数量。

alpha多样性指数的大小是与使用的ASV/OTU表的抽平深度有关，为探究样本alpha多样性随抽平深度的变化曲线，可绘制稀释曲线（rarefaction curve），这是生态领域的一种常用方法。 稀释曲线通过从每个样本中随机抽取一定数量的序列（即在不超过现有样本测序量的某个深度下进行重新抽样），可以预测在一系列给定的测序深度下，所可能包含的物种总数及其中每个物种的相对丰度。因此，通过绘制稀释曲线，还可以在相同测序深度下，通过比较样本中ASV/OTU数的多少，从而在一定程度上衡量每个样本的多样性高低。 Richness指数（物种丰富度指数）的Alpha多样性曲线在很多情况下等同于稀释曲线（rarefaction curve）。下面我开始绘制简单的稀释曲线和Richness指数曲线。 数据文件长这样： 1.调用相关R包，读取otu物种丰度表； 2.定义函数； 3.vegan包rarecurve()可以直接绘制稀释曲线； 4.下面绘制Richness指数曲线 稀释曲线可直接反映测序数据量的合理性，并间接反映样品中物种的丰富程度，当曲线趋向平坦时，说明测序数据量渐进合理，更多的数据量只会产生少量新OTUs（物种）；反之表明不饱和，增加数据量可以发现更多OTUs。