GSEA——从原理到实战

大家好, 今天给大家介绍如何用clusterProfiler进行基因集富集分析。分为三个部分：原理，实战，总结。

原理部分主要是对GSEA作者们2005年文(https://www./content/102/43/15545)想法的解读，在实战部分，用GSEA软件进行基因富集分析，用clusterProfiler实现定制化的基因富集分析。

本文的新意是，我们会介绍如何从RNA-seq数据得到的表达谱进行GSEA分析，并且比较用GSEA和clusterProfiler对同样的数据进行基因集分析，看看其结果是否会有不同。

首先让我们从基因集富集分析原理开始讲起吧。

1. 基因集富集分析原理

1.1 背景

假定你是一个生物医学的研究人员，你对基因X的功能比较感兴趣，你读到一篇文献，里面的作者设计了一系列实验，其中一个实验是，knock down基因X，对WT和KD进行RNA-seq实验，比较WT和KD的表达谱的变化，从原始测序数据获得表达谱后，经过 差异表达分析之后，作者用了一种叫做基因集富集分析（Gene set enrichment analysis, GSEA）的方法，来推测在基因X knock down前后哪些通路里的相关基因出现了上调，哪些通路里相关基因出现了下调。

GSEA是如何做到这一点的呢？让我们从GSEA的原理讲起。

1.2 算法思想

从芯片时代就有的一个老问题是，如果我们对KD和WT做差异表达分析之后，软件会给出的差异表达基因list，按照某个统计量，比如fold change，也就是KD相较于WT的变化倍数，从小到大排序，得到一个rank list，记录为L。怎么从L中解读出有值得研究与挖掘的信息呢？

一种简单粗暴的方式是关注L的两端，一个个查这些上调或者下调的基因，如果出现预期的基因或者比较熟悉的基因就万事大吉了，不过这需要研究者有很强的背景或者先验知识。有没有什么通过计算和统计的方法，就能从L中做出一些有生物意义的推断呢。这正是GSEA算法的动机。

GSEA是Eric Lander教授（人类基因组计划的首席科学家，Broad Institute 的 director，GATK等软件也是Broad开发与维护的）带领的团队开发与维护的软件。GSEA的想法其实很朴素，就是看这些差异表达的基因在一些先验的通路中的富集情况。原假设是，某个通路的所有基因，在L中是随机的分布的，假如我们能观测到某个通路的所有基因突然富集与L中的一端，计算其富集程度，计算其统计显著性，如果小于某个cutoff，那么我们就可以拒绝原假设，认为该通路在L中富集，并且通过富集程度的打分，如果为正，则该通路倾向于在上调的基因中富集，如果为负，则该通路倾向于在下调的基因中富集。难点在于把这些想法转化为算法，代码与用户友好的软件。

GSEA流程示意图如下，来自文献(https://www./articles/ng1180)：

2003年最早版本的 GSEA workflow

精炼版本的流程图如下，来自文献(https://www./content/102/43/15545)

2005年的GSEA版本

在GSEA对应的软件中，用normalized enrichment score(NES)作为富集程度的度量，用p值和FDR作为统计显著性的度量。作者收集的基因集的数据库叫做MSigDB。

关于ES的计算，以及统计显著性的计算和矫正的细节请见原文献，虽然笔者很乐意将这些技术细节该大家讲清楚，但是笔者不愿意剥夺读者阅读原始文献的乐趣，所以此处从略啦。

1.3 差异基因排序度量的选取

做GSEA比较重要的一步是根据实验设计选择一个合适的统计量，作为差异基因排序的度量。GSEA中可以选如下几种统计量进行差异基因的排序，更细节的 说明见：GSEA User Guide(http://software./gsea/index.jsp)。

1.3.1 case-control 样本

对于case-control样本而言（每个样本至少三个重复），有如下几种度量方式，表示均值，表示标准差；

• Signal2Noise

• t-Test statistics:

• Ratio of Classes (也就是大家常说的fold change)

• log2_Ratio_of_Classes (  fold change)

1.3.2 连续表型的样本

对于连续表型的样本而言（例如时间序列表达谱），可以选如下统计量；

• Pearson 相关系数；
• Cosine
• Manhattan measure
• Euclidean measure

这些统计量是无监督学习常用的度量，不再赘述，可见李航的书《统计学习方法》，亦可见GSEA团队推荐的芯片时代的统计教材，Statistics reference: Statistics for Microarrays, Wit, E. and McClure J., John Wiley & Sons Ltd., 2004.对于bulk的RNA-seq的数据，在得到表达矩阵之后，在表达谱层面和芯片数据的数据处理方法并没有什么本质的差别，可以迁移。

2. GSEA分析实战

2.1 选择Cuffdiff2原始文献的结果作为我们的示例数据

GSAE原作者提供了基于Java写的进行GSEA分析的图形界面软件，下载以及使用说明请见GSEA User Guide(http://software./gsea/index.jsp)。

笔者并不想直接把GSEA团队的软件使用说明照着翻译一遍，因为那并没有什么意思，而且很枯燥乏味。在本文中，笔者尝试通过复现文献中的GSEA分析结果，来说明在高通量测序技术占主导的现在，如何从RNA-seq的结果进行GSEA分析。

我们规划如下的几步来进行我们的分析：

1. 找一个公共的RNA-seq数据集，并且该数据集已经有对应的GSEA分析结果，可以检验我们的分析结果；
2. 将该数据集进行预处理，转换为GSEA软件的输入；
3. 设置参数，运行；
4. 读取结果，与原结果对照。

数据集我们选择Cole Trapnell（Tophat2, cufflinks, cuffdiff, Monocle, Cicero等常用软件的作者，涵盖了bulkRNA-seq，scRNA-seq，scATAC-seq数据的分析）等人，发表于2012年的描述Cuffdiff2对应的经典文献(https://www./articles/nbt.2450)。文献中有一个结果(fig.5a)是对人肺纤维细胞的进行HOXA1转录因子的Knowkdown，对case-control的表达量计算的fold change，之后进行GSEA分析可以看到与细胞分裂相关的通路基因的下调：

Tranell_2012_fig5a

我们接下来将演示如何从作者提供的表达谱的信息中，用GSEA软件和clusterProfiler两种不同的方式复现这张图的结果。

2.2 用GSEA软件进行GSEA

2.2.1 所需输入

运行软件前，我们需要搞清楚需要输入的文件有哪些。GSEA软件需要3个输入文件，表达谱文件（.gct文件或者符合要求的.txt文件），样本信息文件（.cls文件），以及基因集文件（.gmt文件）。基因集文件可从官网下载。文件格式的细节见官网的格式说明(http://software./cancer/software/gsea/wiki/index.php/Data_formats)。

那么原作者是如何进行GSEA分析的呢，我们可以看文献中方法的描述：

Gene set enrichment analysis.
Enrichment for up- or downregulation sets of genes from the REACTOME pathway database was computed by running GSEA against the fold-change ranked list of genes in the experiment. Ranking was based on Cuffdiff 2–derived fold change. REACTOME gene sets with between 12 and 200 members in the MSigDB package “c2.all.v3.0.symbols.gmt” were downloaded from ftp://gseaftp./.

从这个描述中我们知道，我们需要得到原始RNA-seq处理之后的表达谱数据(文献作者是提出FPKM的人，所以这里用来衡量基因表达量的FPKM)，在GSEA的输入中用log2_Ratio_of_Classes，基因集选择从GSEA官网上下载v3.0版本的C2基因集，选择 REACTOME gene sets 作为输入基因集。

注意到作者提供的GEO号是GSE37704，在浏览器里输入：

https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE37704

挨个点开后，我们发现该数据集下还有子数据集，其中GSE37003是RNA-seq的数据集。进入GSE37003的页面后，我们发现作者提供了经过Tophat-Cuffdiff流程处理得到RNA-seq表达谱以及差异基因分析的文件（图中圈出来的两个文件，作者正文里放出来的图应该只用了Illumina Hiseq平台测出来的数据）, 这意味着我们可以不用再跑Tophat-Cuffdiff流程了，直接用其结果就可以了。

GSE37003的补充材料

2.2.2 安装GSEA软件与下载基因集

最简单的一步是安装软件与下载基因集。

进入GSEA官网的下载页面我们会发现需要登录：

GSEA login

也不用担心，我们找着要求，比如用学校邮箱注册之后登录就好了。进去之后呢，我们可以看到这个页面：

GSEA software download page

选择合适自己电脑版本的GSEA软件下载后安装。

这个页面的下半部分是基因集，但你会看到这个版本的基因集是最新的

MsigDB 7.0

我们需要的旧版的基因集，所以选择Archive，点击链接进入后，按照作者的说明，选择3.0版本的数据集下载。

MSigDB 3.0

下载压缩包解压，按照作者方法里的说明，选择其中的c2.cp.reactome.v3.0.symbols.gmt文件。该文件包涵了v3.0版本的MsigDB数据发布那年，所有REACTOME PATHWAY数据库中收集的通路相关的基因。为了后面用clusterProfiler，我们还需要选择其中的c2.cp.reactome.v3.0.entrez.gmt文件，该文件收录的通路和基因与上个文件没有不同，只是其中的基因用NCBI的ENTREZ ID代替。

2.2.3 获取基因表达数据并且进行预处理

接下来的工作可以用R做啦。打开R，设置选项默认不进行数据框里的字符串向量转因子的操作，并且选择bioconductor的镜像为清华的镜像。

options(BioC_mirror='https://mirrors.tuna./bioconductor')
options(stringsAsFactors = F)


安装之后，导入如下的包：

library(GEOquery) ###download GEO supplementary files
library(dplyr)
library(ggplot2)
library(ggthemes)
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)

我们从GEO上下载GSE37704的数据集，将其存于当前目录下建好的Trapnell2012文件夹中。

downloadGSE <- function(GSEid,baseDir){
attemptsLeft <- 20
geoFile <- NULL
while(attemptsLeft>0)
{
geoFile <- tryCatch(GEOquery::getGEOSuppFiles(GSEid, baseDir = baseDir,makeDirectory=FALSE), error=identity)
if('error' %in% class(geoFile))
{
attemptsLeft <- attemptsLeft-1
Sys.sleep(5)
}
else
attemptsLeft <- 0
}
}
downloadGSE(GSEid = 'GSE37703',baseDir = 'Trapnell2012/')



解压我们需要的两个文件备用

gunzip('Trapnell2012/GSE37703_hiseq_cuffdiff_genes.fpkm_tracking.gz',
       overwrite = T)
gunzip('Trapnell2012/GSE37703_hiseq_cuffdiff_gene_exp.diff.gz',
       overwrite = T)

打开这两个文件之后，进行对比后，我们发现作者有两个条件，每个条件有三个重复，作者选择了每个条件的三个重复的表达量的中位数来进行fold change的计算。所以我们接下来只需要分析和处理GSE37703_hiseq_cuffdiff_genes.fpkm_tracking文件就可以了。

将表达谱文件GSE37703_hiseq_cuffdiff_genes.fpkm_tracking进行基因名的去重，以及除去没有表达的基因，将data clean后的表达谱文件输出 ，并且导出响应的样本信息文件。

###export expression file
gene.exp <- gene.exp[!duplicated(gene.exp$gene_short_name),]
index <- apply(gene.exp[,-1],1,function(x) x[1]>0 & x[4] >0)
gene.exp <- gene.exp[index,]
write.table(gene.exp,file = 'HOXA1_KD.exp.txt',sep = '\t',quote = F,row.names = F)
###export phenotype file
cat(c(6,2,1,'\n'),file = 'HOXA1_KD.cls',sep = ' ',append = F)
cat(c('#','CTRL','KD','\n'),file = 'HOXA1_KD.cls',sep = ' ',append = T)
cat(c(rep(c(0,1),c(3,3)),'\n'),file = 'HOXA1_KD.cls',sep = ' ',append = T)


2.2.4 运行GSEA软件

打开GSEA软件，将我们准备好的文件导入后，调整为如下参数，以及结果输出路径。

run GSEArun GSEA advanced options

点击run，之后大概等3到5分钟就会出结果。输出结果的文件夹，打开其中的index.html页面：我们会看到一个结果汇总的网页：

GSEA result index page

第一个二级标题下的为Knock down敲除后上调基因倾向于富集的通路，第二个二级标题下的为Knock down敲除后下调基因倾向于富集的通路。分别点开，能看到原作者结果里出现的上调基因富集的通路（例如REACTOME_GPCR_LIGAND_BINDING）以及下调基因富集的通路（例如REACTOME_CELL_CYCLE_MITOTIC）：

up-regulated pathwaysdown-regulated pathways

点击其中的Details，例如我们点击REACTOME_CELL_CYCLE_MITOTIC的Details按钮，我们可以看见文献里经常出现的GSEA结果展示的时候常常展示的Enrichment plot:

2.2.5 后期处理

当然我们也可以根据这个输出的文件夹，复现原作者在正文里展示的气泡图：

###loading GSEA
GSEA_term <- readr::read_csv('GSEA_term.csv')
head(GSEA_term)
GSEA_result_up <- readr::read_delim('HOXA1_KD.Gsea.1574314099254/gsea_report_for_1_1574314099254.xls',delim = '\t')
GSEA_result_up <- GSEA_result_up %>%
  select(NAME,NES,SIZE,`FDR q-val`)
GSEA_result <- merge(GSEA_result_up,GSEA_term,by.x='NAME','ID')

GSEA_result_down <- readr::read_delim('HOXA1_KD.Gsea.1574314099254/gsea_report_for_0_1574314099254.xls',delim = '\t')
GSEA_result_down <- GSEA_result_down %>%
  select(NAME,NES,SIZE,`FDR q-val`)

GSEA_result <- rbind(GSEA_result,
                     merge(GSEA_result_down,GSEA_term,by.x='NAME','ID'))

mid <- mean(GSEA_result$`FDR q-val`)
ggplot(GSEA_result,aes(x=NES,y=name,color=`FDR q-val`,size=SIZE))+
  geom_point()+
  scale_size_continuous(breaks = c(40,80,120,180),
                        name = 'No. of\nsignificant genes')+
  scale_color_gradientn(colours = rainbow(4),limits=c(0,1),
                        breaks=c(0,0.25,0.5,0.75,1),
                        labels=c(0,0.25,0.5,0.75,1))+
  xlim(c(-3,3))+
  geom_vline(xintercept=0)+
  xlab('Normalized enrichment score')+
  ylab('')+
  theme_bw()
ggsave('Trapnell_2012_fig5a_re.pdf',
       width = 8,height = 6)

最后结果如下，和原文献的图比，内容差不多哈，只是排版后期还得用Adobe Illustrator处理一下 legend 和 color bar ：

Trapnell的fig5a复现

2.3 用clusterProfiler进行GSEA分析

南方医科大的余光创教授（人称Y叔）开发出了特别好用的进行各种富集分析的R包clusterProfiler，这个包里也集成了GSEA分析的功能，有两种算法可选，一种是Y叔自己写的DOSE包里的，根据GSEA原始文献的算法实现的函数，一种是俄国人开发的fgsea包实现的快速GSEA算法，默认为fgsea算法。为了节省时间，我们选择默认的fgsea算法进行。

注意clusterProfiler做富集的时候，需要ENTREZID，所以我们选择基因集的时候，用c2.cp.reactome.v3.0.entrez.gmt文件，并且需要把我们自己计算的log2 fold chagne进行转ID的操作。

代码如下：

h_c2_REACTOME_v3 <- read.gmt('msigdb_v3.0_files_to_download_locally/c2.cp.reactome.v3.0.entrez.gmt')

names(geneRanklist)[duplicated(names(geneRanklist))]

geneRanklist.df <- data.frame(geneName=gene.exp$gene_short_name,
log2FC=log2(gene.exp$Lung_A1KD_FPKM)-log2(gene.exp$Lung_Scramble_FPKM))
x <- geneRanklist.df$geneName
eg = bitr(x, fromType='SYMBOL',
toType='ENTREZID',
OrgDb=org.Hs.eg.db)
geneRanklist.df <- merge(eg,geneRanklist.df,
by.x='SYMBOL',
by.y='geneName')
geneRanklist <- geneRanklist.df$log2FC
names(geneRanklist) <- geneRanklist.df$ENTREZID
geneRanklist <- sort(geneRanklist,decreasing = T)
em2 <- GSEA(geneRanklist,pvalueCutoff = 1, TERM2GENE = h_c2_REACTOME_v3)
em2 <- setReadable(em2,OrgDb = org.Hs.eg.db,
keyType = 'ENTREZID')
em2.df <- as.data.frame(em2)
em2.df <- em2.df %>%
select(ID,NES,setSize,p.adjust)
GSEA_result2 <- merge(GSEA_term,em2.df,by='ID')


运行如下代码，我们也可以复现出Tranell文献里的结果：

ggplot(GSEA_result2,aes(x=NES,y=name,color=p.adjust,
                        size=setSize))+
  geom_point()+
  scale_size_continuous(breaks = c(40,80,120,180),
                        name = 'No. of\nsignificant genes')+
  scale_color_gradientn(colours = rainbow(4),limits=c(0,1),
                        breaks=c(0,0.25,0.5,0.75,1),
                        labels=c(0,0.25,0.5,0.75,1))+
  xlim(c(-3,3))+
  geom_vline(xintercept=0)+
  xlab('Normalized enrichment score')+
  ylab('')+
  theme_bw()
ggsave('Trapnell_2012_fig5a_re_clusterProfiler.pdf', width = 8,height = 6)

复现结果如下：

此外，配合Y叔的enrichplot包也可以做出和GSEA软件一样好看的Enrichment plot：

###choose Cell Cycle Mitotic Pathway
id <- 1
anno <- em2[id, c('NES', 'pvalue', 'p.adjust')]
lab <- paste0(names(anno), '=', round(anno, 3), collapse='\n')
gseaplot2(em2,geneSetID = id,
title = em2$Description[id])+
annotate('text', 0.8, 0.8412375,
label = lab, hjust=0, vjust=0)
ggsave('enplot_illustrate.pdf', width = 8,height = 6)



结果如下

3. 总结与讨论

• 富集分析是常用的一种解读与发现组学实验中蕴含的生物信息的统计方法。
• GSEA软件以及clusterProfiler都可以进行GSEA分析。结果相差不大。
• 如果进行人之外的模式动物研究，如小鼠或者斑马鱼，可用msigdbr提供的根据MSigDB中相应的人的基因进行同源转换的基因进行分析。
• 更多富集分析的知识请见Y叔的clusterProfiler book