通过元数据（vegan 包）（phyloseq 包）为 Rarefaction 曲线线着色

这是如何使用 ggplot 绘制稀疏曲线的示例。我使用了bioconductor提供的phyloseq包中的可用数据。

安装 phyloseq：

source('http://bioconductor.org/biocLite.R')
biocLite('phyloseq')
library(phyloseq)

需要其他库

library(tidyverse)
library(vegan)

数据：

mothlist <- system.file("extdata", "esophagus.fn.list.gz", package = "phyloseq")
mothgroup <- system.file("extdata", "esophagus.good.groups.gz", package = "phyloseq")
mothtree <- system.file("extdata", "esophagus.tree.gz", package = "phyloseq")
cutoff <- "0.10"
esophman <- import_mothur(mothlist, mothgroup, mothtree, cutoff)

提取OTU表，转置并转换为数据帧

otu <- otu_table(esophman)
otu <- as.data.frame(t(otu))
sample_names <- rownames(otu)

out <- rarecurve(otu, step = 5, sample = 6000, label = T)

现在你有一个列表，每个元素对应一个样本：

稍微清理一下列表：

rare <- lapply(out, function(x){
  b <- as.data.frame(x)
  b <- data.frame(OTU = b[,1], raw.read = rownames(b))
  b$raw.read <- as.numeric(gsub("N", "",  b$raw.read))
  return(b)
})

标签列表

names(rare) <- sample_names

转换为数据框：

rare <- map_dfr(rare, function(x){
  z <- data.frame(x)
  return(z)
}, .id = "sample")

让我们看看它的外观：

head(rare)
  sample       OTU raw.read
1      B  1.000000        1
2      B  5.977595        6
3      B 10.919090       11
4      B 15.826125       16
5      B 20.700279       21
6      B 25.543070       26

用 ggplot2 绘图

ggplot(data = rare)+
  geom_line(aes(x = raw.read, y = OTU, color = sample))+
  scale_x_continuous(labels =  scales::scientific_format())

素食情节：

rarecurve(otu, step = 5, sample = 6000, label = T) #low step size because of low abundance

可以根据这一点制作额外的分组和颜色列。

这是如何添加另一个分组的示例。让我们假设你有一个表格：

groupings <- data.frame(sample = c("B", "C", "D"),
                       location = c("one", "one", "two"), stringsAsFactors = F)
groupings
  sample location
1      B      one
2      C      one
3      D      two

其中样本根据另一个特征进行分组。您可以使用lapply或map_dfr来groupings$sample标记rare$location。

rare <- map_dfr(groupings$sample, function(x){ #loop over samples
  z <- rare[rare$sample == x,] #subset rare according to sample 
  loc <- groupings$location[groupings$sample == x] #subset groupings according to sample, if more than one grouping repeat for all
  z <- data.frame(z, loc) #make a new data frame with the subsets
  return(z)
})

head(rare)
  sample       OTU raw.read loc
1      B  1.000000        1 one
2      B  5.977595        6 one
3      B 10.919090       11 one
4      B 15.826125       16 one
5      B 20.700279       21 one
6      B 25.543070       26 one

让我们用这个做一个像样的情节

ggplot(data = rare)+
  geom_line(aes(x = raw.read, y = OTU, group = sample, color = loc))+
  geom_text(data = rare %>% #here we need coordinates of the labels
              group_by(sample) %>% #first group by samples
              summarise(max_OTU = max(OTU), #find max OTU
                        max_raw = max(raw.read)), #find max raw read
              aes(x = max_raw, y = max_OTU, label = sample), check_overlap = T, hjust = 0)+
  scale_x_continuous(labels =  scales::scientific_format())+
  theme_bw()