GSA、GSEA、ssGSEA、GSVA的算法原理及它们的联系与区别

大家好，欢迎来到IT知识分享网。

一、 简述

    为了从基因的表达水平中得到更加具体直观的生物学功能变化的信息，多种基于已知的基因集的分析方法应运而生。其中，基因集分析（Gene Set Analysis）、基因集富集分析（Gene Set Enrichment Analysis）、单样本基因集富集分析（Single Sample Gene Set Enrichment Analysis, ssGSEA）、基因集变异分析（Gene Set Variation Analysis, GSVA）是我目前遇到比较多的基因集分析方法，在这里对四者做一个详细的阐述和比较，以加深对这些方法的理解。本文章不介绍代码，只阐述算法原理，因为现在这几种方法的代码的分享已经足够多了，而对于原理的详细阐述几乎没有。

二、GSA与GSEA

1. GSA

    GSA属于ORA法（Over-representation Analysis） ，过表达分析。统计原理其实我们在高中就已经学了。一个袋子里有100个球，其中红球20个，白球80个，现在随机不放回地从袋子中拿出10个球，其中有n个红球的概率？在这里，总共100个球对应总共2w多个基因，红球对应基因集中的基因，从袋子中拿球的数目对应DEG的数目。

   具体步骤如下：

   a) 通过差异分析（limma或DESeq2等R包）基于不同的阈值设定（p-value以及log2FC）得到不同组间的差异表达基因DEG。

   b) 将DEG与感兴趣的基因集做交集（KEGG、GO、MSigDB等数据库），得到一些共同的基因。

   c) 基于超几何分布的Fisher检验来评估，抽到这些共同的基因的的计数值是否显著高于随机，即待测功能集在基因列表中是否显著富集。

    本方法优缺点明显：摘自史上最通俗 Gene enrichment analysis 之 over representation analysis （ORA） 原理解释；

优点：基于完备的统计学理论, 结果稳健、可靠。
缺点：(1)仅使用了基因数目信息，而没有利用基因表达水平或表达差异值，为了获得感兴趣或者差异表达基因，需要人为的设置阈值；(2)ORA法通常仅使用最显著的基因，而忽略差异不显著的基因。在获得感兴趣的基因时, 往往需要选取合适的阈值, 有可能会丢失显著性较低但比较关键的基因, 导致检测灵敏性的降低；(3)将基因同等对待，ORA法假设每个基因都是独立的，忽视了基因在通路内部生物学意义的不同(如调控和被调控基因的不同)及基因间复杂的相互作用；(4)ORA假设通路与通路间是独立的，但这个前提假设是错误的。

    补充几点：(1)ORA方法只关心差异表达基因而不关心其上调、下调的方向，也许同一条通路里既有显著高表达的基因，也有显著低表达的基因，因此最后得到的通路结果很难结合表型进行分析。(2)根据上下调把差异表达基因分开，并分别进行GSA，似乎可以分别得到上下调的通路。但是人为对差异基因的范围进行调整后再进行GSA是否会引入误差呢？(3)针对上下调的差异表达基因分别进行GSA，也许会同时富集到一条通路上，这怎么解释该通路的变化呢？我在分析过程中的确遇到过这样的情况。

    与GSA不同，GSEA方法让通路的上下调分析成为可能：简单来说，先计算基因在不同组间表达量的log2FC，并以此从大到小进行排序，这样就得到了一个基因从上调到不变到下调的基因列表。然后，对于我们感兴趣的基因集，我们只需考察它们的成员是否富集在这个基因列表的上方（上调部分）或者下方（下调部分）即可判断这些基因集的上下调情况。

2. GSEA

    需要注意的是，在GSEA中，不仅仅可以用log2FC对基因进行排序，根据自己的研究内容，可以使用不同的参数比如相关性系数作为排序依据。实际上，GSEA原文是这么说的：“Genes are ranked based on the correlation between their expression and the class distinction by using any suitable metric”。在我来看目前应用得比较多的是不同组间比较得到的log2FC，因此在后面的阐述中我都将以log2FC做为排序的依据。

    在讲解GSEA、ssGSEA、GSVA前，需要了解几个统计学知识，因为它们与这三个方法都息息相关。

注释1. 分布函数（Cumulative Distribution Function，CDF）与概率密度函数（probability density function，PDF）

转存失败重新上传取消

概率密度函数与分布函数

注释2. 经验累计分布函数（Empirical Cumulative Distribution Function，eCDF）

    简而言之，这是一个根据抽样样本数据来近似总体分布函数的方法。我们从总体中抽出n个样本{x1,x2,x3…xn}，对于这个样本我们可以画一个频率密度直方图，并且我们设定每个样本的概率是1/n，于是据此可以画一个该抽样样本的分布函数。因为样本数目有限，样本变量为离散的，所以这个分布函数是阶梯函数（step function），每一步阶梯的高度都是1/n，代表每个对应的样本数据的概率为1/n，将所有的样本数据爬完后最终到达1。如果抽样样本足够多的话，eCDF也就越接近总体的CDF。
    在后面GSEA、ssGSEA、GSVA的详细介绍中，对于这种阶梯式的step function，被描述为random walk，也就是随机游走，随机游走也是一个统计学概念，在这里，我们考虑一个点从原点出发向右行走，当遇到抽样分布的样本点（数据点）时（对应的横坐标），就向上走1/n，如果没遇到就平行x轴行走。在后面的GSEA、ssGSEA、GSVA中，它们的随机游走可能有着不同的策略，比如说遇到在基因集中的基因就向上走一个特定的与表达量相关的数值，没有遇到就向下走一个特定的数值。后面我们会详细讲解。

转存失败重新上传取消

正态分布的经验分布函数

转存失败重新上传取消

样本数越多越近似总体分布函数

注释3 Kolmogorov–Smirnov test

    这是一个非参数检验，通过比较两个抽样样本的eCDF的形状，来检验它们是否来源于同一个总体分布。往往用于检验一个抽样分布是否属于正态分布。参考10: Kolmogorov-Smirnov test

    方法很简单，在同一个图上画出两个抽样样本的eCDF（样本数目分别为n、m），然后找到两条阶梯线最大的差距D，这个D就是我们需要的统计量，对于原假设H0：两抽样分布来源于同一总体，D就会很小，如果D过大（α=0.05）就可以拒绝原假设接受备择假设：两抽样分布来源于不同总体。

    D本身的分布是通过非常多次改变两样本在x轴上的排序从而计算得到的，每一次打乱样本顺序，都可以计算出相应的一个D，得到一个D的分布，这样就可以考察现在的D的水平是否满足p＜0.05。因为对于eCDF来说，确定样本量后，阶梯上升量就确定了，所以改变eCDF形状的因素就只剩下样本在x轴上的分布情况。在Myles Hollander等人的《Nonparametric Statistical Methods》的5.4的注释38里介绍了统计量D的显著性检验（书中用的统计量是J=（nm/(n和m的最大公约数)）*D）。其中介绍了两个抽样样本X、Y且n=1、m=3时的情况。当排序为XYYY时，J=3，排序为YXYY，J=2。

 

转存失败重新上传取消

https://www.youtube.com/watch?v=ZO2RmSkXK3c&ab_channel=MatthewE.Clapham

 

转存失败重新上传取消

《Nonparametric Statistical Methods》的5.4的注释38

    好了，现在我们可以详细具体介绍GSEA了。

GSEA的具体步骤

    为了方便讲解，我们假设现在有一个包含肿瘤样本和正常样本的基因表达矩阵，包含26个基因，我们感兴趣的基因集S包含其中9个基因（行为基因，列为样本）。

   第一步：计算得到与表型相关的参数。在这里，我们计算的是肿瘤样本和正常样本的每个基因表达量的log2FC值（肿瘤/正常）。

   第二步：将所有基因按照log2FC从大到小排列，这样就会形成一个基因list：L，顶部的基因有较大的log2FC，表示在肿瘤中相对高表达；底部的基因有较低的log2FC（负数），表示在肿瘤中相对低表达；中间的基因则比较接近于0，表示两样本中差距不大。

   第三步：计算富集分数（Enrichment Score，ES）：k-s like statistic：
先看看原文的解释：

转存失败重新上传取消

ES

转存失败重新上传取消

k-s like statistic：x轴上的红色线条代表在基因集S中的基因，蓝色表示不在基因集S中的基因。红色线条上升高度为1/9，蓝色线条上升高度为1/17。绿色线条表示两条eCDF在不同基因处的差距，等于红色减去蓝色。

    上图是我随手画的。类似于k-s检验，我们已经有了一个基因list L，我们直接把它当作x轴，于是x轴变量就是所有的基因，并且是按照log2FC从大到小排列的。对于一个我们感兴趣的基因集S（包含9个基因），y轴为累计概率（cumulative probability），我们可以画出它的eCDF：红色线条。同样我们可以画出不在基因集S中的17个基因的eCDF：蓝色线条。

转存失败重新上传取消

A GSEA overview illustrating the method

第四步：检验ES(S)的显著性（计算ES(S)的p-value）：
    “Step 2: Estimation of Significance Level of ES. We estimate the statistical significance (nominal P value) of the ES by using an empirical phenotype-based permutation test procedure that preserves the complex correlation structure of the gene expression data. Specifically, we permute the phenotype labels and recompute the ES of the gene set for the permuted data, which generates a null distribution for the ES. The empirical, nominal P value of the observed ES is then calculated relative to this null distribution. Importantly, the permutation of class labels preserves gene-gene correlations and, thus, provides a more biologically reasonable assessment of significance than would be obtained by permuting genes.“

    采用置换检验（empirical phenotype-based permutation test），参考Permutation test（置换检验）以及在R中的应用。在这里，我们可以重新随机分配基因的log2FC（原文中的phenotype label），然后再按照上述过程计算得到新的ES(S)，这样反复多次进行，次数足够大时，我们可以依据大量的新的ES(S)得到ES(S)的概率分布，然后我们在这个分布中考察之前计算的真正的ES(S)的p值（可以看出来与标准k-s test是一样的方法）。

第五步：多重假设检验：

转存失败重新上传取消

多重假设检验

转存失败重新上传取消

Leading-Edge Subset

    综上我们就介绍完了GSEA。总结一下，相比起GSA，GSEA不再关注于差异基因，因此不受p-value以及log2FC阈值的影响，可以获得更多生物学功能变化的信息。富集分数ES，实际上是k-s like test的统计量，所以ES主要表示基因集S的基因的log2FC的分布与不在基因集S的其他基因的log2FC的分布是否一致，当ES大于0并且具有统计学意义时，那我们可以说基因集S内基因相比其他基因表达上调。

三、ssGSEA

转存失败重新上传取消

文章method部分：Signature projection method（part 1）

转存失败重新上传取消

文章method部分：Signature projection method（part 2）

    第一， 基因的排序方法不同。在GSEA中，基因的是按照log2FC从大到小排列，而在ssGSEA中，对于单个样本，将基因按照其绝对表达量从大到小排序。

   第二， K-s like test中的eCDF阶梯上升高度所依赖的值不同。
   在GSEA中，最后基因集S中基因的阶梯上升高度依赖于log2FC的加权值。而在ssGSEA中，为了消除异常值对结果的影响，在第一步排序基因后，会将表达量用秩次（rank）替换。比如说一个样本有100个基因的信息，先按照基因表达量排好序后，这个genelist对应的值是基因的表达量，然后经过秩次标准化转换后，这个genelist已经由表达量的降序排列变成了100，99，98，…，3，2，1。也就是原文L中的

转存失败重新上传取消, , ,…,。所以最后基因集S中的基因的阶梯上升高度依赖于秩次的加权值。
    并且，ssGSEA默认加权值为0.25（），

   但是，不在基因集S中的其他基因的eCDF上升高度与GSEA一致的（就是GSEA中的），为：

    另外，后面大家会知道，GSEA、ssGSEA、GSVA三者的不在基因集S中的其他基因的eCDF上升高度是相同的。

    第三， ES(S)的计算方式不同。计算两条eCDF在每个基因处的D后，在GSEA中，ES(S)等于绝对值最大的D；而在ssGSEA中，ES(S)等于所有D的加和（积分）。

   最后，对ssGSEA的计算细节感兴趣的话，可以逐行运行此代码：A simple implementation of ssGSEA (single sample gene set enrichment analysis)。

四、GSVA

    第一：基因集的排序方法不同。竟然又有不同，这次GSVA能玩出什么花？我们看看原文：

转存失败重新上传取消

基因排序part 1

转存失败重新上传取消

基因排序part 2

转存失败重新上传取消在n个样本中的表达量，通过核密度估计（kernel estimation）得到基因的概率密度函数，然后积分得到不同样本的CDF值。当所有基因的所有表达值都转换成CDF值后，再在每一个样本中对基因按照CDF值从大到小排序并最终如ssGSEA一样用秩次代替CDF值。这么听着似乎有点疑惑对吧，下面再仔细点解释。

注释4：核密度估计（kernel density estimation）

    可仔细参考什么是核密度估计？如何感性认识？以及Intro to Kernel Density Estimation (KDE)

    回到我们的基因数据来，我们得到的是一系列样本的

genej <- c(6, 3, 2, 7, 4, 2, 7, 3, 1, 5, 3, 2, 5, 7) hist(genej,xlim=c(-2,10),xlab = "",ylab = "",xaxt="n",yaxt="n",main = "") gene_kernel <- density(genej)#kernel estimation par(new=T) plot(density(genej),xlim=c(-2,10),ylim=c(0,0.2)) p <- data.frame(approx(x=gene_kernel$x,y=gene_kernel$y,xout = genej)) plot(x=p$x,y=p$y,xlim=c(-2,10),ylim=c(0,0.2),xlab = "",ylab = "") #get the CDF KPDF <- approxfun(x=gene_kernel$x,y=gene_kernel$y,yleft=0, yright=0) cdf<-sapply(Fh$x,function(probability){integrate(KPDF, -Inf,probability )}) 

转存失败重新上传取消

基因表达数据

转存失败重新上传取消

表达数据的频率直方图

转存失败重新上传取消

表达数据的核密度估计

转存失败重新上传取消

计算基因表达值对应点的CDF值（积分）

转存失败重新上传取消

用CDF值替换基因表达值

    第二：K-s like test中的eCDF阶梯上升高度所依赖的值不同：

转存失败重新上传取消

eCDF阶梯上升高度

转存失败重新上传取消的基因的rank值为，基因总数为，则用于计算阶梯上升高度的中心化后的rank值为：

    如此操作后可以提升基因列表头部及尾部基因的阶梯上升高度。因此，对于样本，基因集S内的基因其eCDF上升高度依赖于标准化后的rank值的加权值，默认权重与GSEA相同（）：

    至于不在基因集S的其他基因eCDF上升高度就不赘述了，如前所述，与GSEA、ssGSEA一致。

    第三：ES(S)的计算方式不同。

转存失败重新上传取消

GSVA score计算方法 part 1

转存失败重新上传取消

GSVA score计算方法 part 2

总结