技术贴 | R语言物种组成分析和绘图
本文由阿童木根据实践经验而整理，希望对大家有帮助。
原创微文，欢迎转发转载。
宏基因组分析分为物种分析和功能分析两大块。物种组成分析是物种分析中最基本最常见的分析方法。利用R语言堆叠图，我们可以将一个项目中所有样品的物种组成展示出来。下面介绍如何利用R语言进行物种组成分析和可视化。过程分为以下几步：1）模拟丰度矩阵；2）模拟分组；3）标准化丰度；4）调整格式；5）ggplot2绘制堆叠图；6）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绘制堆叠图
library(ggplot2)
stack_plot=ggplot(data_frame, aes(sample_id, fill=Taxonomy, 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：上下移动横轴标签
pdf('stack_plot.pdf')
stack_plot
dev.off()
# 保存结果，打开stack_plot.pdf，结果如下：
# 6 ggplot2绘制堆叠面积图
id=rep(1:20, each=10)
data_frame=data.frame(data_frame, id)
# 给每个样品重新编号
stack_area_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'))
# 调整背景
pdf('stack_area_plot.pdf')
stack_area_plot
dev.off()
# 保存结果，打开stack_area_plot.pdf，结果如下：