GSEA简介

首先简单介绍一下GSEA，它是2005年在PNAS上发扬光大的方法，沿用至今，目的是看差异表达的基因在哪些基因集中富集。相比于Over-representation只关注显著差异表达的基因，GSEA分析纳入所有基因，将一些微弱但不显著的效应考虑在内。假设做了AHBA的10000个基因表达和大脑表征相关，其中300个基因经过检验是显著相关的，ORA富集分析这300个显著的基因，而GSEA分析10000个基因无论它们显著与否。

步骤1：计算富集得分(ES)

• 根据基因和表型的相关对其进行排序

• 计算排名基因的累积总和。

  • 当基因在集合S中时增加总和，否则减少。

  • 增量的大小取决于基因与表型的相关性。

• 记录与零的最大偏差作为富集得分

步骤2：评估显著性

• 将表型(phenotype)标签互换1000次(或更多)

• 按上述方法计算每个组合的ES得分，获得Null模型

• 将实际数据的ES得分与来自Null模型ES得分的分布进行比较

步骤3：多重比较矫正

• 根据每个基因集的大小对ES进行归一化处理获取归一化富集分数(NES)。

• 通过计算FDR来控制假阳性的比例。

注意事项：

在GSEA最初的版本中，每个基因使用相同的权重(此处对应的是clusterProfiler中的exponent=0的设置)。改善之后的方法是根据基因与表型的相关大小对其进行加权(对应exponent=1的设置)，这有可能导致Null模型中的ES分布不对称，需要分别对正负ES做标准化，评估显著性和做多重比较矫正。

好消息是，Y姥爷开发的clusterProfiler把计算和作图都搞定了。铺天盖地都是clusterProfiler的教程此处略过，下面是一些经验之谈。

基因集的问题

对于几个比较常用的库，比如GO，KEGG和DisGEnet，都有可观的基因集数量，例如GO中的Biological Process一类就有近1万多的类别。而对于GSEA分析来说，如果使用自己建立的Null模型来计算显著性，其富集分数的计算量就是: 类别数量*Null模型数量(一般是5000-10000)。因此根据基因集大小选择纳入分析的类别多少，就是一个需要考虑的问题。太多的话计算量大需要HPC才能完成，太少的话可能会错过一些重要的类别。

在Fulcher, B. D (2021)的文章中，他们用的基因集大小就是[10,200](这里的两个值指的都是与输入的基因集做交集之后的基因集大小的最小值和最大值，第一个值越小，对于的类别也就越细)，虽然Null模型数量高达40000，但是他们的方法可以向量化计算，并不需要考虑计算量。ClusterProfiler中默认的基因集大小是[10,500]。试过[30,200]的GO-BP+自定义的Null模型，还可以在个人电脑上完成。

GO和DisGeNET存在一个问题就是有很多非大脑相关的类别，如果做矫正的话，这些非大脑相关的类别可能大大降低类别过矫正的概率。另一个方面就是上面提到的，计算时间的问题。因此需要有人完成一些额外的工作，从中整理和大脑相关的GO和DisGeNET，比如做一个BrainGO。类似的工作已经开始出现👇，SynGO整理的就是Synapse相关的GO类别。

如何配置自己的基因集？

clusterProfiler的官方教程中讲的很清楚了:

https://yulab-smu.top/biomedical-knowledge-mining-book/universal-api.html

下面是几个有用的函数，不编程可以略过。

如果自己定义基因集，根据教程获取可以做出TERM2GENE和TERM2NAME，之后用build_Anno生成一个USER_DATA，可以直接传入内部函数GSEA_internal。

USER_DATA = build_Anno(TERM2GENE,TERM2NAME)

USER_DATA是一个environment类型的数，用getGeneSet获取到geneSets的数据，之后用geneSet_filter和自己的geneList做交集之后根据[minGSSize, maxGSSize]筛选纳入分析的基因集selected.gs。

USER_GS=getGeneSet(USER_DATA)
selected.gs=geneSet_filter(USER_GS, geneList, minGSSize,maxGSSize)

clusterProfiler已经整合了org.Hs.eg.db到了一站式分析中，其实它也是用的USER_DATA。

USER_DATA=get_GO_data('org.Hs.eg.db',ont2use,'ENTREZID')

NES计算的问题

为什么不直接用原始的ES? 

ES会受到基因集大小，以及基因集和表型相关的影响，因此在不同基因集之间ES没有可比性。同样的道理，基于ES获取的P值(称之为Nominal p value)只代表了单个基因组的富集分数的统计学意义。当评估多个基因集时，须对基因集的大小和多重假设检验进行修正。

NES计算

归一化富集分数（NES）是考察基因集(geneset)富集结果的主要统计指标。通过归一化富集分数，GSEA考虑了基因集大小的差异以及基因集和表型数据集之间的相关性；因此，归一化富集分数（NES）可用于比较不同基因组的分析结果。计算方法是:

https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideTEXT.htm#_Nominal_P_Value

btw

第一，实际计算时，正负ES是分开计算的，也就是说

Pos_ES/mean(Pos_null_ES)

Neg_ES/mean(Neg_null_ES)

第二，这样用NES得到的p值和原始ES算出来的nominal P value一样，因为所有ES值都除以了同一个值。

Null模型的问题

在步骤2评估显著性时，是对表型而不是基因标签进行排列组合，这样可以保持基因表达数据的复杂相关结构。clusterProfiler作为最受欢迎的富集分析工具包之一，似乎也只提供了置换基因标签的方法评估显著性。

差不多在20年前，生信人就已意识到了要置换表型。最近的Neuroimaging领域的文章也开始注意到基因表达数据的复杂结构，比如co-expression。相比于其他领域，Neuroimaging领域还需要注意空间自相关的问题。

👇下面这篇文章讲了co-expression和spatial-autocorrelation对富集结果的影响，之前介绍过这一篇文章(传送门)，现在看懂了他们用的并不是富集分数，而是类别分数----类别内的平均相关值。

Fulcher, B. D., Arnatkeviciute, A., & Fornito, A. (2021). Overcoming false-positive gene-category enrichment in the analysis of spatially resolved transcriptomic brain atlas data. Nature communications, 12(1), 1-13.

👇另外一篇2022年发在HBM的文章，虽然讲的是感兴趣基因集和脑表型做相关，细节上和富集分析有诸多不同，但思想是一样的。他们提出了在创建Null模型时匹配co-expression的方法，相比置换大脑表型的方法效率不高。

Wei, Y., de Lange, S. C., Pijnenburg, R., Scholtens, L. H., Ardesch, D. J., Watanabe, K., ... & van den Heuvel, M. P. (2022). Statistical testing in transcriptomic‐neuroimaging studies: A how‐to and evaluation of methods assessing spatial and gene specificity (Vol. 43, No. 3, pp. 885-901). Human Brain Mapping.

此外，该用什么统计量在实际模型和Null模型之间做比较，也是需要考虑的问题。在对比真实模型和Null模型时，我们可以对比以下统计量:

• 平均值

• 绝对值的平均

• 中位数
  

• 平方的平均
  

• 显著相关的数量
  

• 富集分数

• 等等等
  

平均值作为在置换检验中的首选，在某些情景下并非最优。有研究表明，用显著相关的数量作为统计量时，似乎受Null模型的类型影响较小。例如，即使在用较宽松的Null模型(只置换基因的标签)时，用显著相关的数量作为统计量(黄色)，比使用平均数(紫色)和绝对值的平均数(绿色)作为统计量更好地控制假阳性的概率。

其实，这是一个trade-off，它在大大降低假阳性率的同时又可能过于严格而有较高的假阴性率。由于目前缺乏grand truth，无法对这一点进行定量分析。此外，根据个人经验，使用富集分数也可能会有类似的表现。

其他关于GSEA

[官方文档]

https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideTEXT.htm

[相关Slides]

https://bioinformatics.mdanderson.org/MicroarrayCourse/Lectures09/gsea1_bw.pdf

计算ES得分的函数

👉DOSE中的gseaScores

https://github.com/YuLab-SMU/DOSE/blob/52613be3ae/R/gsea.R#L441:1

👉fgsea中的 calcGseaStat

https://rdrr.io/github/ctlab/fgsea/src/R/fgsea.R

👉GSEA.1.0.R

https://rdrr.io/bioc/EnrichmentBrowser/src/R/GSEA.1.0.R#sym-GSEA.EnrichmentScore [有FWE矫正的算法]

btw，代码中的if_else用的贼溜，相比而言，matlab的if_else就稍微麻烦一点，需要自己定义一个匿名函数，然后就可以在cellfun中愉快地使用if条件了。

if_else= @(varargin) varargin{3-(varargin{1}>0)}
% 用法:
if_else(condition, true_value, false_value)
% 好处在于可以在cellfun中增加if条件
numeric_array = cellfun( @(x) if_else(isnumeric(x) & ~isempty(x),x,NaN), Ra);

👉clusterProfiler

https://yulab-smu.top/biomedical-knowledge-mining-book/index.html

集成度和效率都很高，中文教程也很多，但似乎需要编程置换phenotype。对于GO的富集分析，可以一站式对BP/MF/CC做富集。查代码得知，矫正是根据这三个类别分开做的。