MSEA
什么是MSEA
MSEA，Metabolite Set Enrichment Analysis，代谢集富集分析，是代谢组数据分析方法的一种。MSEA是基因集富集分析的代谢物版本。
MESA分析流程（简略版）
- 设定一系列代谢集，每个代谢集代表某种生物学功能。
- 将自己的代谢组数据富集到这些代谢集。
- 统计找到有显著性差异的代谢集。
与KEGG Pathway富集分析的区别
- KEGG Pathway分析使用的是差异代谢物。
- MSEA分析使用的是所有代谢物。
- 组间差异不大的代谢物在KEGG Pathway分析中可能被忽略，而MSEA分析则不存在这种问题。
- MSEA分析中使用的代谢集可以根据试验需要进行设定，相对于KEGG Pathway分析更有针对性。
MSEA分析方法
MetaboAnalyst提供了在线MSEA分析工具，同时他们也提供了相应的R包可以进行本地分析
成品图预览
这是我们最终想要得到的图
图1
- 左边是代谢集名称，当然这里我瞎编了一些名字
- 条形的长度代表FE，富集程度
- 条形的颜色代表p值
- 所有代谢集按p值从小到大排序
OK，画图正式开始
数据准备
MetaboAnalyst在分析结束后会自动生成一个结果表，如下
MESA分析结果
对于我们来说，画图真正需要用到的共三列，分别是
- Set：代谢集的名称
- p：统计p值
- FE：Fold Enrichment，富集程度
注1：MetaboAnalyst分析结果中没有FE列，需要自行计算。
FE=StatisticQ / Expected Q
注2：MetaboAnalyst输出的分析结果表已经按p值从小到大进行了排序
基本思路
- 以Set为横坐标，以FE为纵坐标画柱状图
- 将p值映射到柱子的颜色
- 翻转柱状图，得到最终图像
R代码
############################################################
############# MSEA Plot ###############
############ Cdudu ##############
########### 2020.1.22 ##############
########################################################
library(readxl)
library(ggplot2)
library(export)
#按照前面的成品预览图，先进行一些背景设定
mytheme <- theme(axis.title=element_text(face="bold", size=10,colour = 'gray25'), #轴标题
axis.text=element_text(face="bold", size=10,colour = 'gray25'), #轴标签
axis.line = element_line(size=0.5, colour = 'black'), #轴线
axis.line.y = element_blank(), #关闭Y轴线
axis.ticks.y = element_blank(), #关闭Y轴刻度线
panel.background = element_rect(fill="white"), #背景色
panel.grid.major.y=element_blank(), #关闭Y轴主网格线
panel.grid.minor.y=element_blank(), #关闭Y轴次网格线
panel.grid.minor.x=element_blank()) #关闭X轴次网格线
MSEA<-read_excel('MSEA_result_sample.xlsx') #读入MSEA分析结果
p<-ggplot(MSEA,aes(x=Set,y=FE,fill=P)) + #以Set为横坐标，FE为纵坐标画柱状图，并用P填充颜色
geom_bar(stat='identity',color='black',width = 0.65) +
coord_flip() + #颠倒坐标轴
scale_fill_gradient(low='red', high='darkgoldenrod1') #设定柱子颜色变化范围，随着p值从低到高，柱子的颜色从红色向黄色渐变
p + labs(x='',y='Fold Enrichment',fill='p value') +
mytheme
输出结果
图2
很明显，图2和上面的图1不太一样，其中最大的问题是：
图2并没有按照p值从低到高的顺序对代谢集进行排序
这是因为Set列为字符型向量，ggplot2在画图时，会自动对Set列按字母顺序进行重排序。
为了防止ggplot2对Set列的重排序，我们需要将Set列转化为因子型，并设定因子水平的顺序，代码如下：
MSEA$Set<-factor(MSEA$Set,levels = rev(MSEA$Set)) #将Set转化为因子型
将这行代码放到刚才的代码中，再来看一下结果
############################################################
############# MSEA Plot NEW ###############
############ Cdudu ##############
########### 2020.1.21 ##############
########################################################
library(readxl)
library(ggplot2)
library(export)
mytheme <- theme(axis.title=element_text(face="bold", size=10,colour = 'gray25'), #轴标题
axis.text=element_text(face="bold", size=10,colour = 'gray25'), #轴标签
axis.line = element_line(size=0.5, colour = 'black'), #轴线
axis.line.y = element_blank(), #关闭Y轴线
axis.ticks.y = element_blank(), #关闭Y轴刻度线
panel.background = element_rect(fill="white"), #背景色
panel.grid.major.y=element_blank(), #关闭Y轴主网格线
panel.grid.minor.y=element_blank(), #关闭Y轴次网格线
panel.grid.minor.x=element_blank()) #关闭X轴次网格线
MSEA<-read_excel('MSEA_result_sample.xlsx')
MSEA$Set<-factor(MSEA$Set,levels = rev(MSEA$Set)) #将Set转化为因子型
p<-ggplot(MSEA,aes(x=Set,y=FE,fill=P)) +
geom_bar(stat='identity',color='black',width = 0.65) +
coord_flip() + #颠倒坐标轴
scale_fill_gradient(low='red', high='darkgoldenrod1')
p + labs(x='',y='Fold Enrichment',fill='p value') +
mytheme
输出结果
图3
现在看上去好多了，不过还需要一些微调
- 减少图例的刻度数量
scale_fill_gradient(breaks = c(0.01,0.05))
- 去掉图例的白色刻度线
guides(fill=guide_colorbar(ticks = F, frame.colour='black'))
- 调整FE轴刻度
scale_y_continuous(expand = c(0, 0),breaks = seq(0,6,1))
- 调整FE轴与柱子的距离
scale_x_discrete(expand = c(0,0.8))
最终代码
############################################################
############# MSEA Plot NEW ###############
############ Cdudu ##############
########### 2020.1.21 ##############
########################################################
library(readxl)
library(ggplot2)
library(export)
mytheme <- theme(axis.title=element_text(face="bold", size=10,colour = 'gray25'), #轴标题
axis.text=element_text(face="bold", size=10,colour = 'gray25'), #轴标签
axis.line = element_line(size=0.5, colour = 'black'), #轴线
axis.line.y = element_blank(), #关闭Y轴线
axis.ticks.y = element_blank(), #关闭Y轴刻度线
panel.background = element_rect(fill="white"), #背景色
panel.grid.major.y=element_blank(), #关闭Y轴主网格线
panel.grid.minor.y=element_blank(), #关闭Y轴次网格线
panel.grid.minor.x=element_blank()) #关闭X轴次网格线
MSEA<-read_excel('MSEA_result_sample.xlsx')
MSEA$Set<-factor(MSEA$Set,levels = rev(MSEA$Set)) #将Set转化为因子型
p<-ggplot(MSEA,aes(x=Set,y=FE,fill=P)) +
geom_bar(stat='identity',color='black',width = 0.65) +
coord_flip() + #颠倒坐标轴
scale_fill_gradient(low='red', high='darkgoldenrod1',breaks = c(0.01,0.05)) + #控制图例中的刻度
scale_x_discrete(expand = c(0,0.8)) # expand 调整Y轴与第一个柱子间的距离
p + labs(x='',y='Fold Enrichment',fill='p value') +
guides(fill=guide_colorbar(ticks = F, frame.colour='black')) + #关掉图例中的刻度线，图例边框色
scale_y_continuous(expand = c(0, 0),breaks = seq(0,6,1)) + #y轴刻度间隔
mytheme
最终结果
图4
