前面，我們仔細學習了Alpha多樣性的一些常見指數(shù)，今天我們嘗試用R來計算這些常用指數(shù)。
library(vegan)
library(picante)
library(reshape2)
library(ggplot2)
library(tidyverse)
library(ggpubr)
我們今天用一個簡單的OTU的豐度數(shù)據(jù)來做測試。
#讀入物種數(shù)據(jù)
otu <- read.delim("otu_table.tsv", row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
##物種豐富度 Richness 指數(shù)
richness <- rowSums(otu > 0)
#或
richness <- estimateR(otu)[1, ]
##Shannon（以下 Shannon 公式的對數(shù)底數(shù)均設(shè)使用 e，在 R 中即表示為 exp(1)）
#Shannon 指數(shù)
shannon_index <- diversity(otu, index = 'shannon', base = exp(1))
#Shannon 多樣性
shannon_diversity <- exp(1)^shannon_index
#Shannon 均勻度（Pielou 均勻度）
pielou <- shannon_index / log(richness, exp(1))
##Simpson
#Gini-Simpson 指數(shù)（一般用這個）
gini_simpson_index <- diversity(otu, index = 'simpson')
#經(jīng)典 Simpson 指數(shù)
simpson_index <- 1 - gini_simpson_index
#Invsimpson 指數(shù)（Gini-Simpson 的倒數(shù)）
invsimpson_index <- 1 / gini_simpson_index
#或
invsimpson_index <- diversity(otu, index = 'invsimpson')
#Simpson 多樣性
simpson_diversity <- 1 / (1 - gini_simpson_index)
#Simpson 均勻度（equitability 均勻度）
equitability <- 1 / (richness * (1 - gini_simpson_index))
##Chao1 & ACE
#Chao1 指數(shù)
chao1 <- estimateR(otu)[2, ]
#ACE 指數(shù)
ace <- estimateR(otu)[4, ]
##goods_coverage 指數(shù)
goods_coverage <- 1 - rowSums(otu == 1) / rowSums(otu)
##譜系多樣性（與上述相比，還需指定進化樹文件）
tree <- read.tree("rooted_tree.tre")
pd_whole_tree <- pd(otu, tree, include.root = FALSE)
為了方便，我們現(xiàn)在把這些指數(shù)封裝起來，封裝到一個文件中，一起計算：
alpha <- function(x, tree = NULL, base = exp(1)) {
? ? ? ? est <- estimateR(x)
? ? ? ? Richness <- est[1, ]
? ? ? ? Chao1 <- est[2, ]
? ? ? ? ACE <- est[4, ]
? ? ? ? Shannon <- diversity(x, index = 'shannon', base = base)
? ? ? ? Simpson <- diversity(x, index = 'simpson')? ? #Gini-Simpson 指數(shù)
? ? ? ? Pielou <- Shannon / log(Richness, base)
? ? ? ? goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
? ? ? ? result <- data.frame(Richness, Shannon, Simpson, Pielou, Chao1, ACE, 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
}
#所以，我們現(xiàn)在使用封裝好的函數(shù)，就可以一步計算所有樣品的alpha指數(shù)了
#加載 OTU 豐度表和進化樹文件
otu <- read.delim("otu_table.tsv", row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
tree <- read.tree("rooted_tree.tre")
samples_df <- read.delim("group.xls",header = T,sep="\t")? #這個不是必須的，我加載了樣品的分類文件，主要是為了后面畫圖用的。文件內(nèi)容很簡單
colnames(samples_df)<- c("ID","group")
#不包含譜系多樣性，無需指定進化樹；Shannon 公式的 log 底數(shù)我們使用 2
alpha_all <- alpha(otu, base = 2)
#包含譜系多樣性時，指定進化樹文件；Shannon 公式的 log 底數(shù)我們使用 2
alpha_all <- alpha(otu, tree, base = 2)
alpha_all$ID <- rownames(alpha_all)? ?#為了后面合并用
通過調(diào)用自己定義的alpha函數(shù)，我們便可以得到所有樣品的常用的alpha指數(shù)值。
#輸出保存在本地
write.csv(alpha_all, 'alpha.csv', quote = FALSE)
#下面，我們對每一個指數(shù)在樣品組內(nèi)進行對比，看一下有沒有顯著性差異。其實，也可以分面的，但是在分面的過程中發(fā)現(xiàn)，他們共用一個縱坐標，但是參數(shù)的scale差別太大，不適合共用一個縱坐標，所以寫了個for循環(huán)，單獨生成了一個小提琴圖，然后拼在了一起。
alpha_all2 <- melt(alpha_all,id="ID")
alpha_all3 <- merge(alpha_all2, samples_df, by = "ID")
class <- unique(alpha_all3$variable)
p_list <- list()
index <- 1
for(i in class)
{
? cat("i",i,"\n")
? cat("index",index,"\n")
? alpha_tmp <- alpha_all3 %>% filter(variable == i)
? fit <- aov(value~group, alpha_tmp)
? p_value <- round(summary(fit)[[1]][1,5],3)
? p <- ggplot(data = alpha_tmp, aes(x = group, y = value)) +
? ? ? ? ? ? ? geom_violin(aes(fill=group),trim=F,show.legend = FALSE)+
? ? ? ? ? ? ? geom_boxplot(width=0.03,fill="white")+
? ? ? ? ? ? ? labs(title = paste0(i,":"," p = ",p_value))+
? ? ? ? ? ? ? xlab(NULL)+ylab(NULL)+
? ? ? ? ? ? ? theme(text = element_text(size = 15))+
? ? ? ? ? ? ? theme(plot.title = element_text(hjust = 0.5))
? p_list[[index]] <- p
? index <- index+1
}
library(cowplot)
p <- plot_grid(p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]], p_list[[5]], p_list[[6]],p_list[[7]],p_list[[8]], ncol = 2)
ggsave("plot.pdf", plot = p, height = 10, width = 9)
從圖中，我們便可以明顯的對比不同指數(shù)在不同樣品組之間的差異情況了。
