导读

宏基因组分析分为物种分析和功能分析两大块。物种组成分析是物种分析中最基本最常见的分析方法。利用R语言堆叠图，我们可以将一个项目中所有样品的物种组成展示出来。下面介绍如何利用R语言进行物种组成分析和可视化。过程分为以下几步：1）模拟丰度矩阵；2）模拟分组；3）标准化丰度；4）调整格式；5）ggplot2绘制堆叠图、冲积图、分面、分组、堆叠面积图。

1 模拟丰度矩阵

set.seed(1995)
# 随机种子
data=matrix(abs(round(rnorm(200, mean=1000, sd=500))), 20, 10)
# 随机正整数，20行，20列
colnames(data)=paste("Species", 1:10, sep=".")
# 列名
rownames(data)=paste("Sample", 1:20, sep=".")
# 行名
# 得到样品物种丰度矩阵，如下：

2 模拟分组

group=c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B")
sample_id=rownames(data)
data_group=data.frame(sample_id, group)
# 得到分组文件，如下：

3 标准化丰度

data_norm=data
for(i in 1:20){
sample_sum=apply(data, 1, sum)
    # 统计每个样品的总细菌数量
    for(j in 1:10){
        data_norm[i,j]=data[i,j]/sample_sum[i]
        # 将每个样品的总细菌数量控制为1
    }
}

4 调整格式

library(reshape2)
# 加载用于处理数据格式的reshape2包
Taxonomy=colnames(data)
# 从data矩阵中提取物种分类信息
data_frame=data.frame(t(data_norm), Taxonomy)
# 新建数据框

data_frame=melt(data_frame, id='Taxonomy')
# 根据Taxonomy和Sample将所有丰度竖着排列

names(data_frame)[2]='sample_id'
# 重命名variable为sample_id，保持与data_group的样品变量名一致
data_frame=merge(data_frame, data_group, by='sample_id')
# 根据样品变量名，给data_frame添加分组信息，如下：

5 ggplot2绘制堆叠图

1 普通堆叠图

geom_col(position = 'stack')”，y轴展示原始计数
geom_col(position = 'fill')，y轴展示菌丰度除以其在各样本中的菌总丰度

library(ggplot2)

stack_plot=ggplot(data_frame, aes(x=sample_id, fill=Taxonomy, y=value*100))+
# 数据输入：样本、物种、丰度
geom_col(position='stack') +
# stack：堆叠图
labs(x='Samples', y='Relative Abundance (%)')+
# 给xy轴取名
scale_y_continuous(expand=c(0, 0))+
# 调整y轴属性
theme(axis.text.x=element_text(angle=45, hjust=1))
# angle：调整横轴标签倾斜角度
# hjust：上下移动横轴标签

ggsave(stack_plot, filename="stack_plot.pdf")

2 拆成柱形图

geom_bar()和geom_col()都可以完成堆叠图和柱形图
position=position_dodge(0)默认值为0，即默认绘制堆叠图，如果position_dodge > width则能拆开堆叠图得到分组柱形图。

stack_plot = ggplot(data_frame, aes(x=sample_id, fill=Taxonomy, y=value))+
# 数据输入：样本、物种、丰度
geom_bar(stat="identity", position=position_dodge(0.75), width=0.5) +
# geom_col(position=position_dodge(0.75), width=0.5) +
# stack：堆叠图
labs(x='Samples', y='Relative Abundance (%)')+
# 给xy轴取名
scale_y_continuous(expand=c(0, 0))+
# 调整y轴属性
theme_classic() +
theme(axis.text.x=element_text(angle=45, hjust=1))

ggsave(stack_plot, filename="stack_plot.pdf", width=14)

3 添加冲积图

geom_bar(stat='identity') # 同样可以做堆叠图
geom_alluvium() # 添加冲积图
geom_stratum(width=0.45, size=0.1) # 添加阶层，下图中的黑线

安装依赖：

install.packages("ggalluvial")
library("ggalluvial")

install.packages("rlang", version="0.4.7")
packageVersion("rlang")

绘制冲积图：

stack_plot=ggplot(data_frame,
  aes(x=sample_id,
  y=value*100,
  fill=Taxonomy,
  stratum = Taxonomy,
  alluvium = Taxonomy)) +
  geom_bar(stat='identity', width=0.45) +
  geom_alluvium() +
  geom_stratum(width=0.45, size=0.1) +
  labs(x='Samples', y='Relative Abundance (%)')+
  scale_y_continuous(expand=c(0, 0))+
  theme(axis.text.x=element_text(angle=45, hjust=1))

ggsave(stack_plot, filename="stack_plot.pdf")

4 添加facet_wrap分面

facet_wrap(~group, scales = 'free_x', ncol = 2) # 按group组，X轴，分2面

stack_plot=ggplot(data_frame, aes(x=sample_id, 
  fill=Taxonomy, 
  y=value*100,
  stratum = Taxonomy,
  alluvium = Taxonomy))+
  geom_col(position='stack') +
  geom_alluvium() +
  geom_stratum(width=0.45, size=0.1) +
  labs(x='Samples', y='Relative Abundance (%)')+
  scale_y_continuous(expand=c(0, 0))+
  theme(axis.text.x=element_text(angle=45, hjust=1))+
  facet_wrap(~group, scales = 'free_x', ncol = 2)

ggsave(stack_plot, filename="stack_plot.pdf")

5 添加geom_segment分组标记

数据准备：准备geom_segment需要的x、x_end值

x_start = c()
x_end = c()
for(i in 1:nrow(data_frame))
{
    tmp = unlist(strsplit(as.character(data_frame[,1])[i], split="\\."))
    x_start = c(x_start, as.numeric(tmp[2]) - 0.5)
    x_end = c(x_end, as.numeric(tmp[2]) + 0.5)
}
data_frame = data.frame(data_frame, x_start, x_end)

绘图：

stack_plot = ggplot(data=data_frame, mapping=aes(x=sample_id, 
  fill=Taxonomy, 
  y=value*100,
  stratum = Taxonomy,
  alluvium = Taxonomy)) +
  geom_col(position='stack') +
  geom_alluvium() +
  geom_stratum(width=0.45, size=0.1) +
  labs(x='Samples', y='Relative Abundance (%)') +
  theme_classic() +
  theme(axis.text.x=element_text(angle=45, hjust=1)) +
  scale_y_continuous(limits=c(0, 115), 
    # 定义y轴范围
    expand = c(0, 0), 
    # 定义y轴外展范围
    breaks = c(0, 20, 40, 60, 80, 100)) +
    # 定义y轴展示的每个刻度
  geom_segment(mapping=aes(
      x = x_start, 
      y = 105, 
      xend = x_end,
      yend = 105,
      color = group
    ), size = 5)

ggsave(stack_plot, filename="stack_plot.pdf")

6 翻转90度

facet_wrap(~group, scales = 'free_y', ncol = 2) # 按group组，Y轴，分2面
coord_flip() # 旋转90度

stack_plot=ggplot(data_frame, aes(x=sample_id, 
  fill=Taxonomy, 
  y=value*100,
  stratum = Taxonomy,
  alluvium = Taxonomy))+
  geom_col(position='stack') +
  geom_alluvium() +
  geom_stratum(width=0.45, size=0.1) +
  labs(x='Samples', y='Relative Abundance (%)')+
  scale_y_continuous(expand=c(0, 0))+
  theme(axis.text.x=element_text(angle=45, hjust=1))+
  facet_wrap(~group, scales = 'free_y', ncol = 2) +
  coord_flip()

ggsave(stack_plot, filename="stack_plot.pdf")

7 绘制堆叠面积图

数据准备：给每个样品按数字编号

id=rep(1:20, each=10)
data_frame=data.frame(data_frame, id)
# 给每个样品重新编号

绘图：

stack_plot=ggplot(data_frame, aes(id, fill=Taxonomy, value*100))+
geom_area() +
# 堆叠面积图
labs(x='Samples', y='Relative Abundance (%)')+
scale_x_continuous(breaks=1:20, labels=as.character(1:20), expand=c(0, 0))+
scale_y_continuous(expand=c(0, 0))+
# 调整x轴刻度和坐标轴属性
theme(panel.grid=element_blank(), panel.background=element_rect(color='black', fill='transparent'))
# 调整背景

ggsave(stack_plot, filename="stack_plot.pdf")

这配色似乎还可以

ggplot(input, aes(x=name, y=value, fill=variable)) +
  geom_col(position="stack") +
  theme_classic() +
  scale_fill_manual(values=brewer.pal(6, "Set2")) +
  theme(legend.text=element_text(size=15),
        legend.title=element_text(face='bold', size=20)) +
  theme(axis.title = element_text(size = 20),
        axis.text = element_text(size = 18),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1)) +
  theme(text=element_text(family="serif")) +
  labs(x="Phyla", y="CAZyme genes per genome", fill="CAZymes") +
  coord_flip() 

一组好看的堆叠图参数：

ggplot(input, aes(x=variable, y=value*100, fill=Genus)) +
  geom_col(position="stack") +
  theme_classic() +
  scale_fill_manual(values = colors) +
  theme(legend.text=element_text(size=15),
        legend.title=element_text(face='bold', size=20)) +
  labs(x="", 
       y="Relative abundance", 
       fill="Genus") +
  theme(title = element_text(size = 15, face="bold")) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.title = element_text(size = 25),
        axis.text.y = element_text(size = 18),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1)) +
  theme(axis.text.x = element_text(angle = 60, 
                                 hjust = 1, 
                                 size = 20,
                                 color = col_text))

ggsave(result, file="phylum_stack.png", width = 10)
ggsave(result, file="phylum_stack.pdf", width = 10)

参考：
R语言ggplot2绘制分组箱型图和分组柱状图
Make Grouped Boxplots with ggplot2