Loading...
墨滴

fafu生信小蘑菇

2021/07/05  阅读:695  主题:自定义主题1

R语言实现Alpha多样性指数的计算

R语言实现Alpha多样性指数的计算

上次我们已经使用R语言来对OTU表的抽平分析,那么我们如何使用抽平后的OTU表来重新计算Alpha多样性呢?接下来我们就来学习一下。

不过你可能会说,这个不是测序公司都计算好了吗,而且还可以用qiime软件,为啥还要使用R,因为有些测序公司并没有帮你抽平,再一个我使用R语言比较多,如果你会使用其他软件当然更好。欢迎你与我分享一下。

1 所需的数据类型

这里我们需要使用到两个数据集:一个是经过抽平分析处理后的otu表(如果不会抽平分析的可以查看该文章

另外一个数据集是使用各OTU代表序列构建的进化树文件“otu_tree.tre”。

抽平分析后得到的otu表:(当然你也可以选择不抽平)

图1-1 otu表
图1-1 otu表
图1-1 otu_tree.tre文件
图1-1 otu_tree.tre文件

计算每一种Alpha多样性指数都会用到otu表,但是对于otu_tree.tre文件,只用于计算谱系多样性。

2 使用R语言计算常用的Alpha多样性指数

我们接下来会使用到两个包,一个是vegan包,另一个是picante包。如果没安装这两个包,需要提前安装好。

vegan包可以用来计算多种Alpha多样性指数,例如这次我们要学习计算的物种丰富度(Richness)、Chao 1指数、ACE指数、Shannon指数、Simpson指数等。

谱系多样性(即PD_whole_tree)需要使用picante包,该多样性除了物种丰富度数据外还需要进化树文件。

2.1 加载包以及数据集

#设置工作目录
setwd("D:/R_wenji/06-微信公众号/21_07_05")
#需要加载vegan包和picante包,没有安装需要先安装
library(vegan)
library(picante)

#读入抽平后的otu表
otu <- read.delim('otu.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#将otu数据转置
otu <- t(otu)
#加载进化树文件
tree <- read.tree('otu_tree.tre')
图2.1-1 读入抽平后的otu表
图2.1-1 读入抽平后的otu表

2.2 计算-物种丰富度 Richness 指数,又称observed species 指数

#计算方法一:
observed_species <- rowSums(otu > 0)
#计算方法二:
observed_species <- estimateR(otu)[1, ]

#输出查看结果
observed_species
图2.2 物种丰富度 Richness 指数
图2.2 物种丰富度 Richness 指数

2.3 计算 Chao 1指数

Chao1  <- estimateR(otu)[2, ]

Chao1
图2.3 Chao 1指数
图2.3 Chao 1指数

2.4 计算ACE 指数

ACE  <- estimateR(otu)[4, ]

ACE
图2.4 ACE指数
图2.4 ACE指数

2.5 计算Shannon指数

#Shannon 指数,通常使用2、e作为底数
#以e作为底数表示方法
Shannon <- diversity(otu, index = 'shannon', base = exp(1))
#以2作为底数表示方法
Shannon <- diversity(otu, index = 'shannon', base = 2)
#输出Shannon_index结果
Shannon
图2.5 Shannon指数
图2.5 Shannon指数

2.6 计算Simpson指数

#Simpson指数分为经典 Simpson 指数和Gini-Simpson 指数,不过平时常用的 Simpson 指数即为 Gini-Simpson 指数
#Gini-Simpson 指数代码
Gini_simpson  <- diversity(otu, index = 'simpson')
#经常使用
Gini_simpson
#经典 Simpson 指数
simpson_index <- 1 - Gini_simpson
图2.6 Simpson指数
图2.6 Simpson指数

2.7 计算goods_coverage 指数

goods_coverage <- 1 - rowSums(otu == 1) / rowSums(otu)

goods_coverage
图2.7 goods_coverage 指数
图2.7 goods_coverage 指数

2.8 计算谱系多样性(PD)

#除了otu文件,需要指定进化树文件

PD_whole_tree <- pd(otu, tree, include.root = FALSE)[1]
PD_whole_tree
图2.8 谱系多样性(PD)
图2.8 谱系多样性(PD)

单个计算Alpha多样性指数的方法已经讲完了,那么我该如何使用会比较方便呢?

那么请使用下面这个自定义函数,函数不够完美,你有需求可以自己修改。

3 使用自定义alpha_diversity函数来快速计算多种Alpha多样性指数

首先我们需要定义alpha_diversity函数:

library(vegan)
library(picante)      

alpha_diversity <- function(x, tree = NULL) {
  observed_species <- estimateR(x)[1, ]
  Chao1 <- estimateR(x)[2, ]
  ACE <- estimateR(x)[4, ]
  Shannon <- diversity(x, index = 'shannon',base = 2)
  Simpson <- diversity(x, index = 'simpson')    #注意,这里是Gini-Simpson 指数
  goods_Coverage <- 1 - rowSums(x == 1) / rowSums(x)
  
  #保留四位小数
  Shannon <- sprintf("%0.4f", Shannon)
  Simpson <- sprintf("%0.4f", Simpson)
  goods_Coverage <- sprintf("%0.4f", goods_Coverage)
  
  
  result <- data.frame(observed_species, ACE,Chao1, Shannon, Simpson, goods_Coverage)
  
  if (!is.null(tree)) {
    PD_whole_tree <- pd(x, tree, include.root = FALSE)[1]
    names(PD_whole_tree) <- 'PD_whole_tree'
    result <- cbind(result, PD_whole_tree)
 
    result <- data.frame(observed_species, ACE,Chao1, Shannon, Simpson,
                         PD_whole_tree ,goods_Coverage)
  }
  
  
  result
}

alpha_diversity函数定义好了,我们就可以导入数据进行计算了

#加载OTU 表
otu <- read.delim('otu.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
#加载进化树文件
tree <- read.tree('otu_tree.tre')

#如果不需要计算谱系多样性
alpha <- alpha_diversity (otu)

#需要计算谱系多样性时,需要指定进化树文件
alpha1 <- alpha_diversity (otu, tree)
图3-1 不包括谱系多样性的Alpha多样性指数
图3-1 不包括谱系多样性的Alpha多样性指数
图3-2 包括谱系多样性的Alpha多样性指数
图3-2 包括谱系多样性的Alpha多样性指数
#将结果输出,保存在本地
write.csv(alpha, 'alpha_diversity.csv', quote = FALSE)
write.csv(alpha1, 'alpha_diversity1.csv', quote = FALSE)
图3-3 你可以在文件夹中看到它们
图3-3 你可以在文件夹中看到它们

是不是非常简单,如果你需要获取示例数据及代码,可以给我留言,如果你觉得对你有帮助,记得点个赞。

让我们一起加油吧。

fafu生信小蘑菇

2021/07/05  阅读:695  主题:自定义主题1

作者介绍

fafu生信小蘑菇