gsea富集分析全面的实践指南与常见疑问解答

什么是GSEA富集分析？

GSEA，全称Gene Set Enrichment Analysis，即基因集富集分析。与传统的基于基因列表的富集分析（如ORA，Over-Representation Analysis）不同，GSEA不预先设定一个差异表达的阈值来筛选基因，而是基于整个排序的基因列表来进行分析。

它的核心思想是：
如果一个预先定义的基因集（比如某个通路或功能相关的基因集合）在基于某些指标（如差异表达倍数、t统计量等）排序的整个基因列表中，集中分布在列表的顶部或底部，那么就认为这个基因集在该条件下是显著富集的。

具体来说：

输入： 需要两个主要的输入文件：
- 一个基因排序列表：通常是基于不同实验条件（如疾病 vs 正常）下的基因表达数据计算得到的，衡量基因与表型的关联度。这个列表包含了实验中检测到的所有基因，并按照关联度从高到低（或从低到高）排序。
- 一个基因集数据库：这是一系列预先定义的基因集合，每个集合包含一组功能、通路或位置相关的基因（例如，KEGG通路、GO术语、微信号传导通路等）。
过程： GSEA算法沿着排序的基因列表“走”下去，计算一个富集分数（Enrichment Score, ES）。当遇到属于当前基因集的基因时，ES增加；遇到不属于该基因集的基因时，ES减少。ES的峰值反映了该基因集在列表顶部或底部的富集程度。
输出： 对于输入的每一个基因集，GSEA会计算其ES值、归一化ES（NES，用于跨分析比较）、名义p值以及经过多重检验校正的FDR q值。这些结果表明了哪些基因集在特定条件下显著富集。

为什么要使用GSEA？它有什么优势？

使用GSEA而非简单的差异表达基因列表富集，主要基于以下考虑和优势：

利用全量数据： GSEA考虑了所有检测到的基因及其表达变化趋势，而不仅仅是那些通过任意阈值筛选出的“显著”差异基因。这有助于捕捉那些变化幅度不大但协同变化的基因集。
检测细微但协调的变化： 即使某个通路中的基因个体差异表达不显著，但如果它们整体上都倾向于上调或下调，GSEA也能检测到这种通路层面的协同变化。
避免任意阈值问题： 传统的富集分析结果受差异表达基因筛选阈值影响很大，不同的阈值可能导致不同的富集结果。GSEA避免了这一主观选择。
提供更全面的生物学视角： 通过分析基因集，GSEA直接提供了与通路、功能或生物学过程相关的结果，帮助研究者从系统层面理解实验结果的生物学意义。
区分正向和负向富集： ES值的正负直观反映了基因集是在排序列表的顶部（通常对应于某个表型的正相关或上调）富集还是在底部（通常对应于负相关或下调）富集。

在哪里可以进行GSEA分析？

进行GSEA分析有几种常见的平台和工具：

Broad Institute GSEA Desktop Application： 这是GSEA算法最初开发者（Broad Institute）提供的官方桌面应用程序。
- 特点： 功能全面，提供详细的参数设置和可视化结果，内置或方便加载MSigDB等常用基因集数据库。它是基于Java开发的，需要在本地安装。
- 获取： 可在Broad Institute GSEA官方网站免费下载（需注册）。
R语言包： 在生物信息学分析中，使用R语言进行GSEA非常流行，因为它灵活且易于整合到其他分析流程中。
- 常用的包：
  - clusterProfiler：一个功能强大的生物信息学富集分析R包，其中包含GSEA分析的功能（GSEA函数）。它支持多种物种，可以方便地调用GO、KEGG等数据库，也可以使用自定义或MSigDB的基因集。
  - fgsea：一个快速的GSEA实现R包，尤其适合处理大型基因集数据库或需要进行大量样本分析的场景。其速度通常比原始算法实现更快。
- 特点： 需要具备一定的R语言编程能力，但自动化和批量处理能力强，结果可视化选项丰富。
在线Web工具： 虽然有些在线工具提供富集分析功能，但需要注意它们是否执行的是真正的GSEA算法，还是基于基因列表的ORA。提供真正GSEA功能的在线工具相对较少，且可能对数据量有限制。不过，一些平台可能允许你上传GSEA桌面版或R包产生的富集结果进行进一步的可视化或解释。对于初学者或小规模数据，可以寻找提供GSEA功能的在线平台，但需仔细阅读其方法说明。

进行GSEA需要多少数据量和计算资源？

进行GSEA所需的数据量和计算资源取决于多个因素：

输入数据量：
- 样本量： 虽然理论上单个样本也能生成排序列表（例如基于与其他样本的相关性），但GSEA通常用于比较至少两个不同条件的样本组。样本量越大，差异表达分析结果越稳定，基因排序的可靠性也越高，从而提高GSEA结果的可信度。通常建议每组至少有3-5个生物学重复样本。
- 基因数量： 分析中包含的基因越多，计算量越大。
基因集数据库大小： 使用的基因集数据库越大（包含的基因集数量越多，每个基因集包含的基因数量越多），所需的计算资源和时间也越多。MSigDB等大型数据库包含数千甚至上万个基因集。
置换（Permutation）次数： GSEA通过对样本标签或基因进行置换来计算显著性p值。置换次数越多，p值计算越精确，但计算时间也越长。标准的GSEA分析通常建议至少进行1000次置换。
计算资源：
- 内存 (RAM)： GSEA，尤其是在加载大型基因表达矩阵和基因集数据库时，会消耗较多内存。使用Broad Institute桌面版或R包时，建议至少有8GB或更多的RAM，特别是处理大型数据集。
- 处理器 (CPU)： 置换过程是计算密集型的，多核处理器可以显著缩短运行时间。
- 存储空间： 输出文件可能包含详细的富集结果、基因列表和图表，会占用一定的硬盘空间。

总的来说，对于典型的转录组数据集（几万个基因，几十个样本）和常用的基因集数据库（如MSigDB Hallmark或C2），使用一台配置适中的台式电脑（8-16GB RAM，多核CPU）运行Broad Institute桌面版或R包，单次分析通常需要几分钟到几小时不等，具体时间取决于上述因素和设置的参数。

如何运行GSEA分析（以Broad Institute桌面版和R包为例）？

使用Broad Institute GSEA Desktop Application

这是最直观的运行方式，适合不熟悉编程的用户。

准备输入文件：
- 表达数据： 需要一个GCT或TXT格式的表达矩阵文件。行是基因，列是样本，包含表达值。第一列通常是基因ID。
- 表型标签： 需要一个CLS格式的文件，定义每个样本所属的实验组别（例如，Control, Treatment）。
- 基因集数据库： 可以使用内置的MSigDB数据库，或者提供一个GMT格式的基因集文件。
启动软件并加载数据： 打开GSEA软件，选择 “Load Data”，导入你的表达数据和表型文件。
运行GSEA分析：
- 选择 “Run GSEA”。
- 在 “Analysis & Reports” 区域：
  - 选择加载的表达数据文件。
  - 选择加载的表型文件，并选择你想要比较的两个表型组（例如，Treatment vs Control）。
  - 选择要使用的基因集数据库（Gene sets database）。
  - 设置输出目录。
- 在 “Parameters” 区域（重要设置）：
  - Number of permutations： 推荐设置为1000。
  - Phenotype labels： 选择 “phenotype”。
  - Collapse dataset to gene symbols： 如果你的表达数据使用探针ID或其他非基因Symbol的ID，勾选此项并提供平台注释文件进行ID转换（如果软件支持）。
  - Ranking metric： 选择用于基因排序的统计量。常用的有 Signal2Noise (适用于两个表型组比较)、t-Test (适用于两个表型组比较)、Log2 Ratio of Classes (简单倍数变化)、或Ranked list (如果你已经计算好并提供了排序列表)。通常推荐 Signal2Noise 或 t-Test。
  - Enrichment statistic： 默认 Weighted 优于 Classic。
  - Min/Max size exclusion list： 设置基因集包含的最少和最多基因数，通常默认即可（如 Min: 15, Max: 500）。过小或过大的基因集通常不太可靠。
- 点击 “Run” 开始分析。
查看结果： 分析完成后，软件会生成一个HTML格式的报告文件。打开报告，可以查看显著富集的基因集列表、富集分数、p值、FDR q值以及每个富集基因集的详细结果页面和富集图。

使用R包 (以`clusterProfiler`为例)

R语言提供了更大的灵活性，尤其适合自动化和批量分析。

安装并加载包：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db") # 示例：人类基因注释包
BiocManager::install("enrichplot") # 可选，用于可视化

library(clusterProfiler)
library(org.Hs.eg.db) # 示例
library(enrichplot) # 可选

准备输入数据：

基因排序列表： 通常从差异表达分析结果中获得。你需要一个命名向量，向量的值是用于排序的统计量（如logFC乘以-log10(p值)），向量的名称是基因ID（建议使用Symbol或Entrez ID）。

# 假设你有一个差异表达结果数据框 diff_results 包含 gene_id, logFC, pvalue
# 生成排序列表 (示例使用 logFC * -log10(pvalue))
geneList <- diff_results$logFC * -log10(diff_results$pvalue)
names(geneList) <- diff_results$gene_id
geneList <- sort(geneList, decreasing = TRUE) # 从高到低排序

基因集数据库： 可以使用clusterProfiler内置支持的数据库（如GO、KEGG），或者从MSigDB下载GMT文件后读取。

# 使用内置GO BP数据库 (需要ENTREZID)
# 对 geneList 进行ID转换，如果原始ID不是ENTREZID
# geneList_entrez <- bitr(names(geneList), fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
# names(geneList_entrez) <- bitr(names(geneList), fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$SYMBOL # 或者保留SYMBOL作为name，ENTREZID在向量值中，取决于GSEA函数的要求

# 或者读取GMT文件
# msd_geneset <- read.gmt("path/to/your/msigdb.gmt")

运行GSEA分析： 使用GSEA函数。

# 使用GO数据库示例
gsea_result_go <- GSEA(geneList, # 基因排序列表
                     OrgDb = org.Hs.eg.db, # 物种注释数据库
                     ont = "BP", # 选择本体 (BP, CC, MF)
                     nPerm = 1000, # 置换次数
                     minGSSize = 15, # 最小基因集大小
                     maxGSSize = 500, # 最大基因集大小
                     pvalueCutoff = 0.05, # p值阈值
                     verbose = FALSE)

# 使用GMT文件示例
# gsea_result_msigdb <- GSEA(geneList,
#                          TERM2GENE = msd_geneset, # 提供GMT数据
#                          nPerm = 1000,
#                          minGSSize = 15,
#                          maxGSSize = 500,
#                          pvalueCutoff = 0.05,
#                          verbose = FALSE)

查看和可视化结果： GSEA函数返回一个GSEA结果对象，可以像数据框一样查看，也可以使用enrichplot包进行可视化。

# 查看结果摘要
head(gsea_result_go)

# 查看所有显著结果
# View(gsea_result_go[gsea_result_go$qvalues < 0.05,])

# 绘制富集图 (比如显示第一个显著富集通路)
# gseaplot2(gsea_result_go, geneSetId = 1, title = gsea_result_go$Description[1])

# 绘制多个通路的富集图
# gseaplot2(gsea_result_go, geneSetId = 1:5)

# 绘制核心富集基因图 (可选)
# heatplot(gsea_result_go, foldChange=geneList) # 需要安装并加载DOSE包

如何解读GSEA分析结果？

解读GSEA结果报告或R包输出的关键在于理解以下几个核心指标：

NAME/Description： 基因集的名称或描述（如“KEGG_GLYCOLYSIS_GLUCONEOGENESIS”）。
SIZE： 该基因集中包含的基因数量。
ES (Enrichment Score)： 富集分数。正值表示基因集倾向于在排序列表的顶部富集（通常对应于比较组中的上调）；负值表示基因集倾向于在列表的底部富集（通常对应于比较组中的下调）。ES值本身受基因集大小和排名影响，主要看其正负和归一化后的NES。
NES (Normalized Enrichment Score)： 归一化富集分数。ES值经过多次置换计算得到的随机分布进行归一化，使得不同分析结果或不同基因集之间的富集程度可以进行比较。NES值越大（正向或负向），表示富集程度越高。
NOM p-val (Nominal p-value)： 名义p值。基于置换测试计算出的显著性p值。表示观察到的ES值是否显著高于或低于随机期望。
FDR q-val (False Discovery Rate q-value)： 经过多重检验校正后的p值。这是评估富集结果可靠性的最重要指标。通常会设定一个阈值（如FDR q值 < 0.25 或 < 0.05）来确定哪些基因集是统计学显著富集的。FDR q值越小，结果越可信。
FWER p-val (Family-wise Error Rate p-value)： 更严格的多重检验校正方法。如果FDR q值显著，通常FWER p值也会较小，但FDR更常用。
RANK AT MAX ES： 达到最大ES值时，当前基因集中的基因在排序列表中的位置。
LEADING EDGE： 核心富集基因集。这是基因集中对富集分数贡献最大的基因列表。这些基因位于排序列表的顶部（正ES）或底部（负ES），是驱动该基因集富集的主要因素。分析Leading Edge基因可以帮助深入理解富集的生物学机制。
富集图 (Enrichment Plot)：
- 图的上半部分显示了基因集中的基因在排序列表中的位置（垂直黑线）。
- 图的中间部分是ES曲线，曲线的峰值即为ES值。
- 图的下半部分显示了用于排序的基因列表的排名值（如logFC或t-statistic）的分布。
- 通过富集图可以直观地看到基因集中的基因是否集中分布在列表的顶端或底端。

解读步骤建议：

首先关注FDR q-val，筛选出统计学显著富集的基因集（例如，FDR q值 < 0.05 或 < 0.25）。
对于显著富集的基因集，查看其NES值和方向（正或负），判断它们在哪个条件下上调或下调。
查看富集图，确认基因集成员确实集中在列表的相应位置。
分析Leading Edge基因列表，识别出在该基因集富集中起关键作用的基因，结合这些基因的功能进一步解释生物学意义。
将显著富集的基因集进行分类或聚类（例如，使用Enrichment Map可视化），以发现更高级别的生物学主题或关联。

还有哪些与GSEA相关的实践问题？

基因集的选择： 选择合适的基因集数据库对结果至关重要。MSigDB是最常用的数据库，包含多种集合（如Hallmark专注于代表明确定义的生物状态或过程的基因集，C2包含来自各种来源的通路和基因集，C5是GO基因集等）。应根据研究问题和数据特性选择合适的基因集集合。
基因ID的匹配： 输入的基因排序列表中的基因ID必须与基因集数据库中的基因ID格式一致（如都使用Symbol，或都使用Entrez ID）。如果不一致，需要进行ID转换。
排序指标的选择： 不同的排序指标可能影响结果。对于差异表达分析，通常使用t统计量、Signal2Noise或logFC结合p值（如logFC * -log10(p值)）作为排序指标。选择能最好区分表型组的指标。
置换方法的选择： GSEA可以基于样本表型进行置换（Phenotype permutation，适用于比较不同样本组）或基于基因集进行置换（Gene Set permutation，适用于评估特定基因集的内部结构）。对于比较组别差异，通常使用Phenotype permutation。
小样本量问题： GSEA对样本量不是非常敏感，因为它利用了所有基因的信息。但样本量过小（如每组少于3个）会使得差异表达分析结果不稳定，从而影响基因排序的可靠性，间接影响GSEA结果。
结果的可视化进阶： 除了标准的富集图，还可以使用其他工具或R包（如Cytoscape的Enrichment Map插件，或clusterProfiler/enrichplot提供的多种绘图函数）来可视化多个富集基因集之间的关系，例如通过基因重叠来构建网络图。
如何处理多种比较： 如果有多个实验组别或多种比较需求，可以对每对比较独立运行GSEA，然后汇总和比较结果。

GSEA富集分析是一个强大的工具，能够从基因表达数据中挖掘通路层面的生物学信息。掌握其原理、操作流程和结果解读是进行高通量测序数据分析的重要技能。

gsea富集分析