前面，我们仔细学习了Alpha多样性的一些常见指数，今天我们尝试用R来计算这些常用指数。
library(vegan)
library(picante)
library(reshape2)
library(ggplot2)
library(tidyverse)
library(ggpubr)
我们今天用一个简单的OTU的丰度数据来做测试。
#读入物种数据
otu <- read.delim("otu_table.tsv", row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
##物种丰富度 Richness 指数
richness <- rowSums(otu > 0)
#或
richness <- estimateR(otu)[1, ]
##Shannon（以下 Shannon 公式的对数底数均设使用 e，在 R 中即表示为 exp(1)）
#Shannon 指数
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 指数（一般用这个）
gini_simpson_index <- diversity(otu, index = 'simpson')
#经典 Simpson 指数
simpson_index <- 1 - gini_simpson_index
#Invsimpson 指数（Gini-Simpson 的倒数）
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 指数
chao1 <- estimateR(otu)[2, ]
#ACE 指数
ace <- estimateR(otu)[4, ]
##goods_coverage 指数
goods_coverage <- 1 - rowSums(otu == 1) / rowSums(otu)
##谱系多样性（与上述相比，还需指定进化树文件）
tree <- read.tree("rooted_tree.tre")
pd_whole_tree <- pd(otu, tree, include.root = FALSE)
为了方便，我们现在把这些指数封装起来，封装到一个文件中，一起计算：
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 指数
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
}
#所以，我们现在使用封装好的函数，就可以一步计算所有样品的alpha指数了
#加载 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") #这个不是必须的，我加载了样品的分类文件，主要是为了后面画图用的。文件内容很简单
colnames(samples_df)<- c("ID","group")
#不包含谱系多样性，无需指定进化树；Shannon 公式的 log 底数我们使用 2
alpha_all <- alpha(otu, base = 2)
#包含谱系多样性时，指定进化树文件；Shannon 公式的 log 底数我们使用 2
alpha_all <- alpha(otu, tree, base = 2)
alpha_all$ID <- rownames(alpha_all) #为了后面合并用
通过调用自己定义的alpha函数，我们便可以得到所有样品的常用的alpha指数值。
#输出保存在本地
write.csv(alpha_all, 'alpha.csv', quote = FALSE)
#下面，我们对每一个指数在样品组内进行对比，看一下有没有显著性差异。其实，也可以分面的，但是在分面的过程中发现，他们共用一个纵坐标，但是参数的scale差别太大，不适合共用一个纵坐标，所以写了个for循环，单独生成了一个小提琴图，然后拼在了一起。
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)
从图中，我们便可以明显的对比不同指数在不同样品组之间的差异情况了。
网友评论