10X单细胞（10X空间转录组）富集分析GSEA、GSVA算法回顾

又是周五，又一周即将过去，人生总会失去很多东西，但是，想到自己的初心，会明白很多的事情。好了，这一篇我们来简单回顾一下常用的富集分析方法GSEA、GSVA的分析算法原理，如有不准确之处，希望大家指出，共同进步。

期待有缘人的相逢

GSEA部分

Gene Set Enrichment Analysis (GSEA，基因集富集分析)用来评估一个预先定义的基因集的基因在与给定按照一定分类标准的基因表（可以是某个功能相关的基因列表，也可以是某个信号通路相关的基因列表）的分布趋势，从而判断其对这个功能或者信号通路的贡献。

其输入数据包含两部分，一是已知功能的基因集(可以是GO注释、MsigDB的注释或其它符合格式的基因集定义)，一是表达矩阵 (也可以是排序好的列表)，软件会对基因根据其与功能基因集的关联度(可以理解为表达值的变化)从大到小排序，然后判断基因集内每条注释下的基因是否富集于该功能基因集相关度排序后基因表的上部或下部，从而判断此基因集内基因的协同变化对该基因集描述的功能变化的影响。

GSEA原理

给定一个排序的基因表L和一个预先定义的基因集S(比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因，或同一GO注释下的基因)，GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异，若研究的基因集S的成员显著聚集在L的顶部或底部，则说明此基因集成员对表型的差异有贡献，也是我们关注的基因集。

GSEA目的就是看样本差异表达基因在一些先验的基因通路或者给定的基因集合中的富集情况。原定假设是某个通路所有基因，在L中是随机分布的，假设我们能观测到某个通路的所有基因突然富集于L中的一端，计算其富集程度，计算其统计显著性，设定截断值，小于这个截断值，则拒绝原假设，认为该通路在L中被富集到，并进行富集程度打分，如果为正，则该通路倾向于在上调的基因中富集，如果为负，则该通路倾向于在下调的基因中富集。

1.对差异基因排序度量的选取

GSEA分析中，首先对样本检测的基因进行排序，用什么样的指标进行大小排序，这个非常关键，往往根据实验方案设计来进行选择。

对于实验vs对照的实验方案设计往往度量都是均值、标准差、log2FoldChange来进行排序。
对于连续性样本呢，往往可以使用Pearson相关系数、Cosine、Manhattan measure、Euclidean measure这些参数进行排序。

2. 计算富集得分（ ES, enrichment score）

然后根据基因集S与样本中基因排序L的顺序依次开始打分，计算ES分数。具体如下：

从基因集 L 的第一个基因开始，计算一个累计统计值（sum for ES）。
累加规则：当遇到一个落在 s 里面的基因，则增加统计值。遇到一个不在 s 里面的基因，则降低统计值。
每一步统计值增加或减少的幅度与基因的表达变化程度相关（统计使用的指标是第一步给定的指标进行）。
富集得分ES最后定义为最大的峰值。正值 ES表示基因集在列表的顶部富集，负值 ES表示基因集在列表的底部富集。

18814178-7794798d11b20ffc.png

3.计算富集得分的显著性

通过基于表型而不改变基因之间关系的置换检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少，也可基于基因集做置换检验 (permutation test)，计算 p-value 。

4. 多重假设检验校正

首先对每个基因子集 s 计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score ( NES )。随后针对NES计算假阳性率。（计算NES也有另外一种方法，是计算出的ES除以置换检验得到的所有ES的平均值）

5.Leading edge analysis and core enriched genes

Leading-edge subset，对富集得分贡献最大的基因成员。Tags说明对富集分数有贡献的基因的百分比，List指出在列表中获得富集分数的位置，Signal是富集信号的强度。获得有助于富集的核心富集基因也将是非常有趣的。

GSEA应用范围

GSEA能在两种不同的生物学状态中，判断某一组有特定意义的基因集合的表达模式更接近于其中哪一种。因此GSEA是一种非常常见且实用的分析方法，可以将数个基因组成的基因集与整个转录组、修饰组等做出简单而清晰的关联分析。

除了对特定gene set的分析，反过来GSEA也可以用于发现两组样本从表达或其它度量水平分别与哪些特定生物学意义的基因集有显著关联，或者发现哪些基因集的表达模式或其他模式更接近于表型A、哪些更接近于表型B。这些特定的基因集合可以从GO、KEGG、、hallmark或MSigDB等基因集中获取，其中MSigDB数据库整合了上述所有基因集。也可自定义gene set (即新发现的基因集或其它感兴趣的基因的集合，所以有时候也用GSEA做细胞定义)。

GSVA部分

Gene Set Variation Analysis（GSVA）与 GSEA 的原理类似，也是计算每个基因集合在每个样本中的 enrichment statistic（ES，或 GSVA score），其算法流程如下：

企业微信截图_16451485352420.png

GSVA原理

不同于 GSEA 之处在于，对于不同的数据类型（只支持 log 表达值或原始的 read counts 值），假设了不同的累积密度函数（cumulative density function，CDF）。

芯片数据：正态分布密度函数 RNA-seq 数据：泊松分布密度函数（这个我们应该用不到）

而且，GSVA 是为每个样本的每个基因计算对应的 CDF 值，然后根据该值对基因进行排序，这样，每个样本都有一个从大到小排序的基因列表。

为什么这样做呢？因为要使得数据符合一定的数据分布，才能使用对应的统计检验进行p值计算，检验其显著性或准确性。

对于某一基因集合，计算其在每个样本中的 ES 值，也就是评估基因集合在基因列表中的富集情况。例如，我们有一个排序后的样本：

图片.png

基因集合包含：B、E、H，我们可以绘制这样一张 K-S 分布图

企业微信截图_16451487568443.png

x 轴为排序后的基因顺序，依照这一顺序，如果基因在集合内，则累积和会加上该基因的值（与基因的顺序有关，排名越靠前值越大），否则累积和不变。

将基因列表分为基因集内核基因集外两个集合，就可以绘制两个分布（红色和绿色曲线），分别计算两个分布之间的最大间距，以基因集内的分布值更大视为正间距（即红色曲线更高），两个最大间距之和即为该基因集的 ES 值。

这部分其实和GSEA计算ES得分有点相似。

这样，就把基因水平的表达矩阵转换成了基因集水平的评分矩阵，可以使用差异表达基因识别算法，寻找显著差异的基因集，从而达到类似功能富集的作用。

GSVA计算步骤

第一步：对于RNA-seq的原始count数据我们需要利用泊松核函数计算Fr：

企业微信截图_1645148882870.png

其中：

x_ij 是第 j 个sample中第 i 个基因的原始count数
n为总的sample数
k与sample数有关系，k=1...n
y=1...x_ij
x_ik 为第 k 个sample中第 i 个基因的count值
r = 0.5

对于芯片数据，首先得对芯片的表达矩阵取一个log2，然后利用高斯核函数计算F_hi(这部分了解即可)：

企业微信截图_16451490996731.png

其中：

xij 是第 j 个sample中第 i 个基因的log2表达值
n为总的sample数
k与sample数有关系，k=1...n
t 为形式变量
xik 为第 k 个sample中第 i 个基因的log2表达值
hi = si / 4，si 为第 i 个基因在所有sample中表达量的标准差；hi可以控制分辨率

第二步：

我们把计算出来的 Fr 和 Fhi 称为zij，那么每一个sample每一个基因都会有对应的 Fr 和 Fhi 的值，接下来对于每一个sample由zij的大小进行排序。

企业微信截图_16451492419776.png

即对于每一个sample，zij大的排前面，小的排后面，排完序后将zij改写为z(i)j，z(i)j定义为排过序的zij；并定义：rij = | p/2 - z(i)j |，其中 p 为表达矩阵中所有的基因数。

第三步：

这一步是利用 Kolmogorov-Smirnov 计算富集分数，首先计算vjk。

企业微信截图_1645149337877.png

其中：

τ 取 1 为权重参数
rij 为第二步所定义,rij = | p/2 - z(i)j |
p 为表达矩阵中所有的基因数
ℓ 代表的是排序后从i = 1...ℓ 的基因
γk 代表第 k 个gene set ，类似于GSEA，每个gene set 对应一个功能（pathway）
I 其实是判别函数，当g(i)∈γk，值为1；当g(i)不属于γk，值为0

那么富集分数怎么计算呢？类似于GSEA：对于某一个sample j，定义富集分数为vjk(ℓ)所构成的集合中（ℓ = 1...p）的最大值减去vjk(ℓ)所构成的集合中（ℓ = 1...p）的最小值，用它们的差值作为富集分数。其数学表达式为：