R语言绘制物种稀释曲线及其它多种Alpha多样性曲线

R语言绘制Alpha多样性曲线

在微生物16S/ITS/18S测序分析及相关的文章中，我们经常可以看到这样的图，称为稀释曲线（Rarefaction curves）。

稀释曲线是从样品中随机抽取一定测序量的数据（序列条数），统计它们所对应的OTUs种类（代表物种），并以抽取的测序数据量与对应的代表OTUs来构建曲线。横坐标代表随机抽取的序列数量，纵坐标代表观测到的OTUs种类数量，样本曲线的延伸终点的横坐标位置为对应样本的测序数量。稀释曲线可直接反映测序数据量的合理性，并间接反映样品中物种的丰富程度，当曲线趋向平坦时，说明测序数据量渐进合理，更多的数据量只会产生少量新OTUs（物种）；反之表明不饱和，增加数据量可以发现更多OTUs。

这种稀释曲线实质上反映了Alpha多样性中的物种丰富度指数（Richness）信息。除此之外，在相关的文章中，我们有时还可以看到以Shannon指数、Chao1指数等绘制的稀释曲线，如下所示（注：在16S/ITS/18S测序分析中，常见observed species，其实就是Richness指数的另一种名称）。无论这些稀释曲线代表了物种丰富度（累计数量）、Shannon指数、Chao1指数还是其它，总之均代表了Alpha多样性指数的信息，因此它们均可以统称为Alpha多样性曲线。

在前篇，白鱼小编已经为大家介绍了群落多样性分析中常见Alpha多样性指数概念，及其在R语言中的计算方法。与之对应，本篇将简介如何使用R作统计并绘制Alpha多样性曲线。

示例数据、R脚本等，已上传至百度盘（提取码488s）。

https://pan.baidu.com/s/1HInewdVJ-T15hPHxI8gmYg

网盘文件“otu_table.txt”为某16S扩增子测序所得OTU丰度表格，可以理解为物种丰度表。表中每一列为一个样本，每一行为一种OTU，交叉区域为每种OTU在各样本中的丰度。

 另一文件“otu_tree.tre”为使用各OTU代表序列构建的进化树文件（即OTU水平的16S进化树，已存储为nwk格式）。纯文本模式打开查看该进化树文件，内容长这样。

接下来首先使用R计算微生物群落的多种Alpha多样性指数，每一种Alpha多样性指数都会用到OTU丰度表；对于进化树文件，将用于计算一种特殊的多样性指数，谱系多样性。之后，再通过计算结果绘制Alpha多样性曲线。

R语言绘制Alpha多样性曲线

Alpha多样性曲线作图思路：将原始的OTU丰度表读到R中；然后对该OTU表进行稀释抽样，并统计每种抽样深度下的Alpha多样性指数；最后将不同抽样深度下的数据合并，绘制曲线图就可以了。

定义统计函数

为了方便后续运行，首先对使用到的函数打个包。常见Alpha多样性指数统计函数主要包括（详情见前文）：

vegan包diversity()，计算Shannon指数、Simpson指数；estimateR()，计算Chao1指数、ACE指数；picante包pd()，计算谱系多样性（即PD_whole_tree）。 

library(vegan)
library(picante)

#读取 OTU 丰度表
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)

##定义函数
#计算多种 Alpha 多样性指数，结果返回至向量
alpha_index <- function(x, method = 'richness', tree = NULL, base = exp(1)) {
if (method == 'richness') result <- rowSums(x > 0) #丰富度指数
else if (method == 'chao1') result <- estimateR(x)[2, ] #Chao1 指数
else if (method == 'ace') result <- estimateR(x)[4, ] #ACE 指数
else if (method == 'shannon') result <- diversity(x, index = 'shannon', base = base) #Shannon 指数
else if (method == 'simpson') result <- diversity(x, index = 'simpson') #Gini-Simpson 指数
else if (method == 'pielou') result <- diversity(x, index = 'shannon', base = base) / log(estimateR(x)[1, ], base) #Pielou 均匀度
else if (method == 'gc') result <- 1 - rowSums(x == 1) / rowSums(x) #goods_coverage
else if (method == 'pd' & !is.null(tree)) { #PD_whole_tree
pd <- pd(x, tree, include.root = FALSE)
result <- pd[ ,1]
names(result) <- rownames(pd)
}
result
}

#根据抽样步长（step），统计每个稀释梯度下的 Alpha 多样性指数，结果返回至列表
alpha_curves <- function(x, step, method = 'richness', rare = NULL, tree = NULL, base = exp(1)) {
x_nrow <- nrow(x)
if (is.null(rare)) rare <- rowSums(x) else rare <- rep(rare, x_nrow)
alpha_rare <- list()

for (i in 1:x_nrow) {
step_num <- seq(0, rare[i], step)
if (max(step_num) < rare[i]) step_num <- c(step_num, rare[i])

alpha_rare_i <- NULL
for (step_num_n in step_num) alpha_rare_i <- c(alpha_rare_i, alpha_index(x = rrarefy(x[i, ], step_num_n), method = method, tree = tree, base = base))
names(alpha_rare_i) <- step_num
alpha_rare <- c(alpha_rare, list(alpha_rare_i))
}

names(alpha_rare) <- rownames(x)
alpha_rare
}

##测试
#统计 OTU 丰度表中各样本的 Shannon 指数，对数底数使用 e
shannon_index <- alpha_index(otu, method = 'shannon', base = exp(1))
#以 1000 条序列为抽样步长，依次对 OTU 表稀释抽样，直到最大序列深度；并统计各抽样梯度下的 OTU 丰度表中各样本的 Shannon 指数，对数底数使用 e
shannon_curves <- alpha_curves(otu, step = 1000, method = 'shannon', base = exp(1))

第一个函数alpha_index()用于计算Alpha多样性指数。其中，x为输入的物种丰度表；method选项对应具体的Alpha多样性指数（见注释）；若计算PD_whole_tree，则还需指定进化树文件tree，计算其它指数时可缺省进化树；base用于设置Shannon指数的对数底数，默认使用自然对数底数e（在R中表示为exp(1)）。

第二个函数套用了第一个函数，首先对丰度表进行稀释抽样，并使用alpha_index()计算各抽样深度下的Alpha多样性指数。其中，x为输入的物种丰度表；step为抽样步长，例如当设置step=1000时，则以1000为一个梯度，抽样1000、2000、3000……条序列，直到最大深度或者达到设定值；rare，设定一个值，抽样到此深度下停止，否则将默认抽样至各样本的最大深度；其余三个参数method、tree、base和alpha_index()一致。

接下来，我们使用这两个函数去完成以上目的。

稀释曲线

这个是最常见的，也是最简单的。如本篇一开始时提到的，我们平时常说的稀释曲线（rarefaction curve），绝大多数情况下等同于Richness指数（物种丰富度指数）的Alpha多样性曲线。当我们统计不同抽样深度下的物种丰富度时，所绘制的Richness稀释曲线即反映了不同个体数量下的物种类型数量信息。

 其实，通过vegan包rarecurve()可以直接绘制稀释曲线。

#以 2000 条序列为一抽样深度（步长）
rarecurve(otu, step = 2000, col = c('red', 'green', 'blue', 'orange', 'purple', 'black'))

考虑到rarecurve()只能绘制稀释曲线，本文为了方便和其它Alpha多样性指数结合在一起计算，以及更好地实现可视化，未再考虑rarecurve()。下文开始统一使用自定义的函数执行Alpha多样性统计，并使用ggplot2作图展示。

Richness指数曲线

以下示例，设定以2000条序列为一抽样深度，依次对OTU表中的OTU序列条数进行抽样并统计各深度下的Richness指数，一直统计到各样本的最大深度。“richness_curves”以列表类型存储了统计结果，然后我们继续读取列表“richness_curves”中的内容并转化为数据框类型，构建ggolot2作图文件，Alpha多样性曲线使用ggplot2绘制。

##以下以物种丰富度指数为例绘制 Alpha 多样性曲线
#以 2000 步长（step=2000）为例统计
richness_curves <- alpha_curves(otu, step = 2000, method = 'richness')

#获得 ggplot2 作图文件
plot_richness <- data.frame()
for (i in names(richness_curves)) {
richness_curves_i <- (richness_curves[[i]])
richness_curves_i <- data.frame(rare = names(richness_curves_i), alpha = richness_curves_i, sample = i, stringsAsFactors = FALSE)
plot_richness <- rbind(plot_richness, richness_curves_i)
}

rownames(plot_richness) <- NULL
plot_richness$rare <- as.numeric(plot_richness$rare)
plot_richness$alpha <- as.numeric(plot_richness$alpha)

#ggplot2 作图，更好的可视化效果，还请自行修改 ggplot2 作图细节
library(ggplot2) #若未加载时先加载

ggplot(plot_richness, aes(rare, alpha, color = sample)) +
geom_line() +
labs(x = 'Number of sequences', y = 'Richness', color = NULL) +
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'), legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = min(rowSums(otu)), linetype = 2) +
scale_x_continuous(breaks = seq(0, 30000, 5000), labels = as.character(seq(0, 30000, 5000)))

示例结果如下所示，图中的虚线为具有最低OTU序列条数的样本中的OTU序列条数。

如果觉得基于一次抽样的结果可能存在误差，希望在同深度下多次抽样并统计均值和标准差，并最终绘制成带有误差棒的曲线图，做法如下（以重复抽样5次为例计算）。

##多计算几次以获取均值 ± 标准差，然后再展示出也是一个不错的选择
#重复抽样 5 次
plot_richness <- data.frame()

for (n in 1:5) {
richness_curves <- alpha_curves(otu, step = 2000, method = 'richness')

for (i in names(richness_curves)) {
richness_curves_i <- (richness_curves[[i]])
richness_curves_i <- data.frame(rare = names(richness_curves_i), alpha = richness_curves_i, sample = i, stringsAsFactors = FALSE)
plot_richness <- rbind(plot_richness, richness_curves_i)
}
}

#计算均值 ± 标准差（doBy 包中的 summaryBy() 函数）
library(doBy)
plot_richness_stat <- summaryBy(alpha~sample+rare, plot_richness, FUN = c(mean, sd))
plot_richness_stat$rare <- as.numeric(plot_richness_stat$rare)
plot_richness_stat[which(plot_richness_stat$rare == 0),'alpha.sd'] <- NA

#ggplot2 作图，更好的可视化效果，还请自行修改 ggplot2 作图细节
ggplot(plot_richness_stat, aes(rare, alpha.mean, color = sample)) +
geom_line() +
geom_point() +
geom_errorbar(aes(ymin = alpha.mean - alpha.sd, ymax = alpha.mean + alpha.sd), width = 500) +
labs(x = 'Number of sequences', y = 'Richness', color = NULL) +
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'), legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = min(rowSums(otu)), linetype = 2) +
scale_x_continuous(breaks = seq(0, 30000, 5000), labels = as.character(seq(0, 30000, 5000)))

与上图相比，在每次抽样步长位置处，展示了多次抽样统计结果的均值±标准差。

其它Alpha多样性指数曲线

同样地，对于其它类型的Alpha多样性指数曲线，如Shannon指数、Simpson指数、Chao1指数、ACE指数等，方法一致。

不过由于不同指数的特征不同，统计作图时可能会存在一些需要调整的小细节。以下继续展示两种其它类型的Alpha多样性指数曲线示例。

Shannon指数曲线

如Shannon指数稀释曲线，方法和上述Richness稀释曲线（物种累计曲线）一致。

##对于 Shannon 指数等，方法类似
#以 2000 步长（step=2000）为例统计每个稀释梯度下的 Shannon 指数，Shannon 公式的对数底数默认为 e，若有需要可更改（例如 2）
shannon_curves <- alpha_curves(otu, step = 2000, method = 'shannon', base = 2)

#获得 ggplot2 作图文件（略，参见上述）
#ggplot2 作图（略，参见上述）

如果觉得折线图不美观（有时折线图的线条可能会出现很大的波动），期望使曲线更加平滑，我们可以考虑替换为平滑拟合曲线的样式，而不再直接使用直线段连接各点。做拟合曲线的方法有很多，以下方法只是其一，大家有兴趣还可自行了解更多方法。

##若简单的“geom_line()”样式波动幅度过大，不平滑等，可以尝试拟合曲线的样式
#获得作图数据。前面多生成一个点，目的使 Shannon 拟合曲线更加平滑
shannon_curves1 <- alpha_curves(otu, step = 200, rare = 200, method = 'shannon')
shannon_curves2 <- alpha_curves(otu, step = 2000, method = 'shannon')
shannon_curves <- c(shannon_curves1, shannon_curves2)

plot_shannon <- data.frame()
for (i in 1:length(shannon_curves)) {
shannon_curves_i <- shannon_curves[[i]]
shannon_curves_i <- data.frame(rare = names(shannon_curves_i), alpha = shannon_curves_i, sample = names(shannon_curves)[i], stringsAsFactors = FALSE)
plot_shannon <- rbind(plot_shannon, shannon_curves_i)
}

rownames(plot_shannon) <- NULL
plot_shannon$rare <- as.numeric(plot_shannon$rare)
plot_shannon$alpha <- as.numeric(plot_shannon$alpha)
plot_shannon <- plot_shannon[order(plot_shannon$sample, plot_shannon$rare), ]

#ggplot2 作图，更好的可视化效果，还请自行修改 ggplot2 作图细节
library(ggalt) #若未加载时先加载，使用 ggalt 包的 geom_xspline() 绘制平滑拟合线

ggplot(plot_shannon, aes(rare, alpha, color = sample)) +
geom_xspline() +
labs(x = 'Number of sequences', y = 'Shannon', color = NULL) +
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'), legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = min(rowSums(otu)), linetype = 2) +
scale_x_continuous(breaks = seq(0, 30000, 5000), labels = as.character(seq(0, 30000, 5000)))

同样地，对于拟合线，我们也可以基于多次抽样统计的结果绘制，以减小仅基于一次抽样而带来的误差。通过在不同深度下重复抽样后，不再根据多次的计算结果统计均值和标准差，而是直接使用这些点去做拟合曲线。这里就不再展示了，方法很简单。

PD_whole_tree曲线

PD_whole_tree指数（谱系多样性指数）则比较特殊，它除了需要使用到OTU丰度表外，还需指定进化树文件。对于PD_whole_tree稀释曲线的做法，只需在指数计算时添加进化树信息即可，其它过程与上述一致。

##对于 PD_whole_tree，除了 OTU 丰度表，还使用到进化树文件
#加载 OTU 丰度表和进化树文件
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
tree <- read.tree('otu_tree.tre')

#以 2000 步长（step=2000）为例统计
pd_curves <- alpha_curves(otu, tree = tree, step = 2000, method = 'pd')

#统计及做图方法同上述，不再展示

友情链接