美文网首页
基于OTU的alpha多样性指数计算和可视化

基于OTU的alpha多样性指数计算和可视化

作者: 吴十三和小可爱的札记 | 来源:发表于2020-04-01 14:17 被阅读0次

1. 简介

物种丰富度指数(Species richness)为群落中丰度大于0的物种数之和,一般用Observed OTU(observed species)表示,但只有物种种类信息,没有丰度信息,数值范围一般为几百至几千不等,范围很大,与研究对象有关;

香农指数(Shannon index)或称香农熵指数(Shannon entropy index)、香农-威纳指数(Shannon-Wiener index),大家最常用的Shannon index 数值为1-10左右的小数,是综合物种数量和丰度两个层面的结果。如果群落仅由单一物种组成(种群),那么我们确信随机选择的个体必定为那个唯一的物种,此时不确定性就为零;否则,我们将无法得知随机被选择的个体究竟属于什么物种,并且不确定性也会随着群落物种种类数的增多而增加。

辛普森指数(Simpson index)同样考虑了物种丰富度以及均匀度,但与Shannon指数相比,它更受均匀度的影响(Simpson,1949))。经典Simpson指数代表了在群落中两个随机选择的个体属于同一物种的概率,当群落物种丰富度增加时,这种概率降低,即Simpson指数随着物种丰富度的增加而降低。由于经典Simpson指数与物种丰富度相反的趋势不直观,如今常用演变而来的Gini-Simpson指数表示Simpson指数,即用1减去经典Simpson指数数值后得到,此时Simpson指数随着丰富度的增加而增加(二者保持一致的趋势)。

特别注意:通常情况下为保持物种丰度指数与Simpson指数的趋势一致,直接将“Gini-Simpson”定义为“Simpson”,即绝大多数情况下,我们在文献中看到的Simpson指数或者软件直接给出的Simpson指数结果,其实是Gini-Simpson指数,而并非经典Simpson指数。

均匀度(Evenness)用于度量群落中相对物种丰富度,Shannon均匀度(Shannon’s evenness),又称Pielou均匀度(Pielou’s evenness),为群落实际的Shannon指数与具有相同物种丰富度的群落中能够获得的最大Shannon指数的比值。如果所有物种具有相同的相对丰度,则该值为1。

Chao1是根据出现1/2次的OTU来估算总体;还有PD whole tree是考虑物种进化关系权重,认为分类学上非常上近的物种存在一定相关性;

物种丰富度指数(Species richness)为群落中丰度大于0的物种数之和,值越大表明群落中物种种类越丰富。

ACE指数在生态学中同样作为度量物种丰富度的指标(Chao and Yang,1993),其值越高代表群落物种越丰富。

goods_coverage常用在微生物16S/18S/ITS测序中,作为反映测序深度的指标。其值越接近于1,说明测序深度越合理,测序深度已经基本覆盖到样品中所有的物种;反之,测序深度不高,许多物种仅被测到了一次,暗示着很多低丰度物种可能尚未被测序测到。

数据样式:每一行为一个样本,每一列为一种OTU,交接位置为相应丰度计数。

2. alpha多样性指数计算

library(picante)
library(ggplot2)
otu <- matrix(sample(c(0:1000), 1200, replace = TRUE), ncol = 12, nrow = 100, 
              dimnames = list(row_names = paste0("OTU",seq(1:100)),
                              col_names = paste0("sample",seq(1:12))))
head(otu)
# 样本的OTU丰度/ 该样本OTU总丰度 = 相对丰度
otu <- apply(otu, 2, function(x){x/sum(x)})

# 如果计算无误,每个sample的相对丰度和为1
apply(otu, 2, sum)

# 转制 - 变成行代表样品,列代表OTU,交叉区域为丰值度
otu <- t(otu)

# 物种丰富度 Richness 指数(observed species)
richness <- rowSums(otu > 0)

# 
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))

# 转换组成数据集
shannon_diversity <- data.frame(shannon_diversity)
shannon_diversity$sample <- rownames(shannon_diversity)

# 加入分组信息
group_data <- data.frame(
  sample = paste0("sample",seq(1:12)),
  group = rep(LETTERS[1:4], each = 3),
  sites = rep(c("soil", "root"), each = 6),
  stringsAsFactors = FALSE)

# 根据sample 合并
shannon_diversity <- merge(shannon_diversity, 
                           group_data, by = "sample")
str(shannon_diversity)
# 箱线图展示
ggplot(data = shannon_diversity) + geom_boxplot(aes(x = group, 
                                                    y = shannon_diversity,
                                                    fill = sites))

其他包括:Simpson指数(Gini-Simpson指数)、Pielou均匀度、Chao1指数、ACE指数、goods_coverage等可以换index进行计算,或者用estimateR()计算常见几种;谱系多样性指数需要引入进化树,用pd(){picante}计算。
如果用qiime2 分析平台,则可以用一句命令计算所有多样性指数。16S rRNA 流程目录 - qiime2

打包进function中,避免太多重复代码

library(vegan)
library(picante)


otu <- read_xlsx(path = file, col_names = TRUE)
otu <- t(otu)

#  大多数alpha多样性计算
# 公式的对数底数均设使用 e(exp(1))
# 给出的simpson 指数是Gini-Simpson 指数值与物种丰富度趋势一致
Alpha_diversity_index <- function(x, tree = NULL, base = exp(1)) {
  est <- estimateR(x)
  Obs <-  est[1, ]
  Shannon <- diversity(x, index = 'shannon', base = base)
  Simpson <- diversity(x, index = 'simpson')   
  Pielou <- Shannon / log(Obs, base)
  goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
  result <- rbind(est, Shannon, Simpson,
                  Pielou, goods_coverage)
  # 如果检测到树文件,将添加Pd指数计算
  if (!is.null(tree)) {
    # Pd 同时计算谱系多样性(PD)和物种丰富度(SR)
    # 如果此处删除[1], 最终结果有两个我中丰富度值
    Pd <- pd(x, tree, include.root = FALSE)[1]
    Pd <- t(Pd)
    result <- rbind(result, Pd)
  }
  result
}
alpha.png

3. 批量箱线图

# 读入实验设计和Alpha多样性值
alpha_file <-  "C:\\Users\\...\\total_data\\"
simple_site <- read_xls(path = paste0(alpha_file, "the_information_of_groups.xls"))
alpha_diversity <- read.table(file = paste0(alpha_file, "alpha_diversity_index.txt"), header = TRUE, sep = "\t")

theme <- theme(panel.background = element_blank(), # 去掉背景格子
               # 显示x平行网格线
               panel.grid.major.y = element_line(colour = "black"), 
               # 显示x轴坐标
               axis.line.x = element_line(colour = "black"),
               axis.line.y = element_line(colour = "black"),
               axis.title.x = element_blank()
               )

# 按Sample_name
data <- merge(simple_site, alpha_diversity, by = "Sample_name")

# Grouped boxplot  --------------------------------------------
p <- list()
for (i in colnames(alpha_diversity)[2:8]){
  p[[i]] <- ggplot(data = data, aes_string(x = "Group1",
                         y = i)) +  
    geom_jitter(aes(color = Group1)) +
    geom_boxplot(aes(fill = Group1), width = .2, position = position_nudge(x = - 0.01),
                 fill = "white", size = 0.8) + 
    theme
  }


plot = wrap_plots(p, nrow = 2) + plot_annotation(tag_levels = 'a')

file_save <- "C:\\Users\\Administrator\\Desktop\\total_data\\alpha_diversity\\"

ggsave(filename = paste0(file_save, "ABCD_alpha_diversity.tiff"),
       plot = plot )

棒棒糖图

当发现箱形图非常丑,如:箱体很长,离群值的点非常多,可以采用棒棒糖图直接可视化每个样点的多样性指数。

4. 统计组间是否显著差异

  1. anova对指数与分组统计
  2. 使用TukeyHSD对组间进行检验,效正pvalue
  3. R - t检验
  4. R - 非参数检验

相关文章

网友评论

      本文标题:基于OTU的alpha多样性指数计算和可视化

      本文链接:https://www.haomeiwen.com/subject/ungouhtx.html