准备数据,即OTU table,原始的或者rare的均可。
下面是绘图代码:
```
library(tidyverse)
library(ggthemes)
library(ggsci)
#加载数据
otutab <- read_tsv("otutab.txt")
# 定义rank函数(需要排除为0的OTU)
myrank <- function(x){
res <- length(x)-rank(x, ties.method = "first")+1
res[x==0] <- NA
return(res)
}
# 生成rank矩阵
otutab_rank <- map_dfr(otutab, myrank)
rownames(otutab_rank) <- rownames(otutab)
# 重塑数据
otutab2 <- otutab %>% rownames_to_column("FeatureID") %>% gather("Group", "Count", -FeatureID)
otutab_rank2 <- otutab_rank %>% rownames_to_column("FeatureID") %>%
gather("Group", "Rank", -FeatureID)
# 合并otutab数据和ranks数据
data <- left_join(otutab2, otutab_rank2) %>% group_by(Group) %>%
mutate(Abundance = Count/sum(Count)) %>% na.omit()
# 绘制 Rank-Abundance,注意y轴数据进行对数变化
ggplot(data, aes(x = Rank, y = Abundance, color = Group)) + geom_line() +
scale_y_log10() +
scale_color_jco(name = NULL) +
labs(title="Rank-Abundance Curves") + ylab("% Abundance")+
theme_bw() +
theme(panel.grid=element_blank(), plot.title = element_text(lineheight=2.5, face="bold",hjust=0.5))
```
数据概览和结果
原始数据otutab:
otutabotutab_rank:
otutab_rank绘图数据data:
绘图数据dataRank-Abundance曲线
Rank-Abundance曲线
网友评论