R語言實現(xiàn)Alpha多樣性指數(shù)的計算

上次我們已經(jīng)使用R語言來對OTU表的抽平分析，那么我們?nèi)绾问褂贸槠胶蟮腛TU表來重新計算Alpha多樣性呢？接下來我們就來學習一下。

不過你可能會說，這個不是測序公司都計算好了嗎，而且還可以用qiime軟件，為啥還要使用R，因為有些測序公司并沒有幫你抽平，再一個我使用R語言比較多，如果你會使用其他軟件當然更好。歡迎你與我分享一下。

1 所需的數(shù)據(jù)類型

這里我們需要使用到兩個數(shù)據(jù)集：一個是經(jīng)過抽平分析處理后的otu表（如果不會抽平分析的可以查看該文章
；

另外一個數(shù)據(jù)集是使用各OTU代表序列構(gòu)建的進化樹文件“otu_tree.tre”。

抽平分析后得到的otu表：（當然你也可以選擇不抽平）

圖1-1 otu表

圖1-1 otu_tree.tre文件

計算每一種Alpha多樣性指數(shù)都會用到otu表，但是對于otu_tree.tre文件，只用于計算譜系多樣性。

2 使用R語言計算常用的Alpha多樣性指數(shù)

我們接下來會使用到兩個包，一個是vegan包，另一個是picante包。如果沒安裝這兩個包，需要提前安裝好。

vegan包可以用來計算多種Alpha多樣性指數(shù)，例如這次我們要學習計算的物種豐富度（Richness）、Chao 1指數(shù)、ACE指數(shù)、Shannon指數(shù)、Simpson指數(shù)等。

譜系多樣性（即PD_whole_tree）需要使用picante包，該多樣性除了物種豐富度數(shù)據(jù)外還需要進化樹文件。

2.1 加載包以及數(shù)據(jù)集

#設置工作目錄
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數(shù)據(jù)轉(zhuǎn)置
otu <- t(otu)
#加載進化樹文件
tree <- read.tree('otu_tree.tre')

圖2.1-1 讀入抽平后的otu表

2.2 計算-物種豐富度 Richness 指數(shù)，又稱observed species 指數(shù)

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

#輸出查看結(jié)果
observed_species

圖2.2 物種豐富度 Richness 指數(shù)

2.3 計算 Chao 1指數(shù)

Chao1  <- estimateR(otu)[2, ]

Chao1

圖2.3 Chao 1指數(shù)

2.4 計算ACE 指數(shù)

ACE  <- estimateR(otu)[4, ]

ACE

圖2.4 ACE指數(shù)

2.5 計算Shannon指數(shù)

#Shannon 指數(shù),通常使用2、e作為底數(shù)
#以e作為底數(shù)表示方法
Shannon <- diversity(otu, index = 'shannon', base = exp(1))
#以2作為底數(shù)表示方法
Shannon <- diversity(otu, index = 'shannon', base = 2)
#輸出Shannon_index結(jié)果
Shannon

圖2.5 Shannon指數(shù)

2.6 計算Simpson指數(shù)

#Simpson指數(shù)分為經(jīng)典 Simpson 指數(shù)和Gini-Simpson 指數(shù)，不過平時常用的 Simpson 指數(shù)即為 Gini-Simpson 指數(shù)
#Gini-Simpson 指數(shù)代碼
Gini_simpson  <- diversity(otu, index = 'simpson')
#經(jīng)常使用
Gini_simpson
#經(jīng)典 Simpson 指數(shù)
simpson_index <- 1 - Gini_simpson

圖2.6 Simpson指數(shù)

2.7 計算goods_coverage 指數(shù)

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

goods_coverage

圖2.7 goods_coverage 指數(shù)

2.8 計算譜系多樣性（PD）

#除了otu文件，需要指定進化樹文件

PD_whole_tree <- pd(otu, tree, include.root = FALSE)[1]
PD_whole_tree

圖2.8 譜系多樣性（PD）

單個計算Alpha多樣性指數(shù)的方法已經(jīng)講完了，那么我該如何使用會比較方便呢？

那么請使用下面這個自定義函數(shù)，函數(shù)不夠完美，你有需求可以自己修改。

3 使用自定義alpha_diversity函數(shù)來快速計算多種Alpha多樣性指數(shù)

首先我們需要定義alpha_diversity函數(shù)：

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 指數(shù)
  goods_Coverage <- 1 - rowSums(x == 1) / rowSums(x)
  
  #保留四位小數(shù)
  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函數(shù)定義好了，我們就可以導入數(shù)據(jù)進行計算了

#加載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多樣性指數(shù)

圖3-2 包括譜系多樣性的Alpha多樣性指數(shù)

#將結(jié)果輸出，保存在本地
write.csv(alpha, 'alpha_diversity.csv', quote = FALSE)
write.csv(alpha1, 'alpha_diversity1.csv', quote = FALSE)

圖3-3 你可以在文件夾中看到它們

是不是非常簡單，如果你需要獲取示例數(shù)據(jù)及代碼，可以給我留言，如果你覺得對你有幫助，記得點個贊。

讓我們一起加油吧。