基于OTU的稀释曲线（Rarefaction curves） +

1. 简介

稀释曲线（Rarefaction curves）是从样品中随机抽取一定测序量的数据（序列条数），统计它们所对应的OTUs种类（代表物种），并以抽取的测序数据量与对应的代表OTUs来构建曲线。

一般情况下，横坐标代表随机抽取的序列数量，纵坐标代表观测到的OTUs种类数量，样本曲线的延伸终点的横坐标位置为对应样本的测序数量，反映了Alpha多样性中的物种丰富度指数（Richness）信息。

稀释曲线可直接反映测序数据量的合理性，并间接反映样品中物种的丰富程度，当曲线趋向平坦时，说明测序数据量渐进合理，更多的数据量只会产生少量新OTUs（物种）；反之表明不饱和，增加数据量可以发现更多OTUs。

2. 实战

数据样式：常见OTU丰度表，列为样本，行为OTU，交叉区域为每种OTU在各样本中的丰度。

library(picante)
library(ggplot2)
rm(list = ls())
file <- file.path

genes_abundance <- read.delim(file = paste0(file, "otu_table.txt"),
                              row.names = 1, sep = '\t', 
                              stringsAsFactors = FALSE, 
                              check.names = FALSE)
# 去掉"#"注释开头,和物种注释结果
genes_abundance <- genes_abundance[-c(1:4), ]
genes_abundance <- genes_abundance[-ncol(genes_abundance)]

# 行列转换
otu <- t(genes_abundance)

# step 表示抽样步长
rare_otu <- rarecurve(otu, step = step, label = FALSE)

names(rare_otu) <- rownames(otu)

plot_data <- mapply(FUN = function(x, y) {
  mydf <- as.data.frame(x)
  colnames(mydf) <- "value"
  mydf$sample_name <- y
  mydf$subsample <- attr(x, "Subsample")
  mydf
}, x = rare_otu, y = as.list(rownames(otu)), SIMPLIFY = FALSE)

xy <- do.call(rbind, plot_data)
rownames(xy) <- NULL  # pretty
head(xy)

ggplot(xy, aes(x = subsample, y = value, color = sample_name)) +
  theme_bw() +
  scale_color_discrete(guide = FALSE) +  # turn legend on or off
  geom_line() + geom_point(