GSEA

GSEA背景

比如你读到一篇文献，里面的作者设计了一系列实验，其中一个实验是，knock down基因X，对WT和KD进行RNA-seq实验，比较WT和KD的表达谱的变化，从原始测序数据获得表达谱后，经过差异表达分析之后，作者用了GSEA的方法，来推测在基因X knock down前后哪些通路里的相关基因出现了上调，哪些通路里相关基因出现了下调。GSEA是如何做到这一点的呢？让我们从GSEA的原理讲起。

GSEA的优势和必要性

一般的差异分析（GO和KEGG pathway）往往侧重于比较两组间(比如处理组和对照组)的基因表达差异，根据差异表达的基因得到一个差异基因列表L (gene list)。针对这个基因列表L，一般会关注少数几个top基因(上调/下调基因)，利用差异倍数(FC，fold change)变化大的基因进行GO 富集分析，Pathway富集分析以及各种图片的绘制。但是差异基因是通过人为定义的阈值得到的，而这种一刀切的阈值，会遗漏掉那些表达差异不显著但有重要生物学意义的基因，因为实际通过芯片观测到的RNA 表达变化，往往是层层的负反馈调控后的结果(个人理解：即由小差异累积起来的大差异)，并且不同组织对于表达差异的敏感度是不同的，在神经递质系统内，一个1.2 倍的表达差异就可产生极其显著的效应。而GSEA不需要指定明确的差异基因阈值就能把表达谱芯片数据与生物学意义很好地衔接起来。

另外，假如一般的差异分析（GO和KEGG pathway）富集到的某个通路，既有上调差异基因，也有下调差异基因，那么该通路是被抑制还是被激活？或者更直观点说，这条通路下的基因表达水平在实验处理后是上升还是下降? 也有人灵光一闪，想出一个解决方案，在进行传统的富集分析时，假如每次只提取上调或者下调的差异基因来进行分析，那么如果上调基因和下调基因分开富集，然后富集到了同一条通路，又该怎么解释？传统的富集分析只能定位到功能，即这些差异基因与哪些功能相关，而不能回答这条通路是被抑制还是被激活的问题。而GSEA可以通过预定义的基因集在排序好的基因列表中的分布回答某通路被抑制还是被激活！(预定义的基因集通常来自功能注释或先前实验的结果，可以是GO注释、MsigDB的注释或其它符合格式的基因集定义；排序好的基因列表是将基因按照在两类样本中的差异表达程度排序；查看预定义基因集在排序好的基因列表中的分布是指查看预定义基因集是否在这个排序表的顶端或者底端富集。)所以GSEA对于GO和KEGG是一个补充，因为GSEA检测基因集而不是单个基因的表达变化，因此可以包含细微的表达变化，预期得到更为理想的结果。

GSEA有三个特点：分析的基因集合而不是单个基因；将基因与预定义的基因集进行比较；富集分析；

GSEA其实就是得到一堆基因，想看这些基因在哪些功能上有富集，David和Kobas也可以做这个工作，但是它们需要先找差异基因，然后将差异基因输入，最后得到GO term 或者KEGG；而GSEA不需要先找差异基因。GSEA能反映一批基因的微量变化的累积所造成的显著功能差异。

对于背景中的例子，就是对KD和WT做差异表达分析之后，软件会给出的差异表达基因list，按照某个统计量，比如fold change，也就是KD相较于WT的变化倍数，从小到大排序，得到一个rank list，记录为L。怎么从L中解读出信息呢？一种简单粗暴的方式是关注L的两端，一个个查这些上调或者下调的基因，如果出现预期的基因或者比较熟悉的基因就万事大吉了，不过这需要研究者有很强的背景或者先验知识。而通过GSEA就能够从L中做出一些有生物意义的推断。

GSEA是Broad 研究所的Eric Lander教授在2005年的时候发表的。GSEA就是看这些差异表达的基因在一些先验的通路中的富集情况。原假设是，某个通路的所有基因，在L中是随机的分布的，假如我们能观测到某个通路的所有基因突然富集与L中的一端，计算其富集程度，计算其统计显著性，如果小于某个cutoff，那么我们就可以拒绝原假设，认为该通路在L中富集，并且通过富集程度的打分，如果为正，则该通路倾向于在上调的基因中富集，如果为负，则该通路倾向于在下调的基因中富集。

GSEA软件以及clusterProfiler都可以进行GSEA分析。结果相差不大。

GSEA分析结果解读

tvwj8x.jpg

因为GSEA需要一个排序的基因表L和一个预定义的基因集S (比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因，或同一GO注释下的基因)。

GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部或底部。

这些基因排序的依据是其在不同表型状态下的表达差异，若研究的基因集S的成员显著聚集在L的顶部或底部，则说明此基因集成员对表型的差异有贡献，也是我们关注的基因集。

这张原理图，GSEA的输入是一个基因表达量矩阵，其中的样本分成了A和B两组(可以认为处理组和对照组)，首先对所有基因进行排序，排序的标准即foldchange, (基因在两组间表达量的变化趋势)。排序之后的基因列表其顶部可以看做是上调的差异基因，其底部是下调的差异基因。GSEA分析的是一个预定义的基因集下的所有基因是否在这个排序列表的顶部或者底部富集，如果在顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势。

tv0PVH.jpg

分成3个部分，第一部分为基因Enrichment Score的折线图，横轴为该预定义基因集下的每个基因，纵轴为对应的Running ES, 在折线图中有个峰值，该峰值就是该基因集的Enrichemnt score，峰值之前的基因就是该基因集下的核心基因。这个峰出现在左侧，Enrichment Score都是正数，可看做该基因集是上调趋势；若峰出现在右侧，则Enrichment Score都是负数，可看做该基因集是下调趋势。

第二部分为hit，每条竖线都是该基因集下的一个基因，

第三部分为所有基因的rank值分布图, 默认采用Signal2Noise算法，对应了纵轴的标题。

tv0nsS.jpg

一篇Cell文章(Disease Model of GATA4 Mutation Reveals Transcription Factor Cooperativity in Human Cardiogenesis)

与心脏发育有关的基因集在iwt组中普遍表达更高，而在G296S组中表达更低；而对于参与内皮或内膜发育的基因集，在iwt组中表达更低，在G296S组中表达更高。根据这个图和其它证据推测iwt组的心脏发育更加完善，而G296S组更倾向于心脏内皮或内膜的发育，即GATA基因的这种突变可能导致心脏内皮或内膜的过度发育而导致心脏相关疾病的产生。

tv0MZQ.jpg

图中间，就是我们每个gene set里面的基因在所有的2万个排序好基因的位置，如果gene set里面的基因集中在2万个基因的前面部分，就是在case里面富集，如果集中在后面部分，就是在control里面富集着。