生物信息分析必备 | GSEA富集分析一文通透

最新推荐文章于 2025-02-23 04:22:12 发布

mutouyuaichi

最新推荐文章于 2025-02-23 04:22:12 发布

阅读量3.5k

点赞数 25

分类专栏： R进阶学习 R生信分析文章标签：数据分析

本文链接：https://blog.csdn.net/mutouyuaichi/article/details/144982294

版权

R进阶学习同时被 2 个专栏收录

3 篇文章

订阅专栏

R生信分析

3 篇文章

订阅专栏

生物信息分析必备 | GSEA富集分析一文通透

1. 富集分析概述

基因富集分析（Gene Enrichment Analysis）是一种常用的生物信息学方法，用于解释在基因组或基因集合中出现的显著富集的功能或特定特征。这种分析用于高通量基因表达数据的解释，比如基因芯片数据或RNA测序数据。其通常包含以下四种功能注释分析：

图1 富集分析

2. 富集分析数据库简介

不管想做哪种富集通路分析，均可去在线网站进行操作，下面列出了常用的富集通路数据库，本文主要介绍R代码实现通路富集分析，详细示例见目录章节 “5. 示例操作”

👉 温馨提示：R代码实操示例请直接下滑到文末 “5. 示例操作” 章节查看。

1. KEGG（Kyoto Encyclopedia of Genes and Genomes）

图2 KEGG

简介：是一个包含生物系统功能和信息的数据库资源，提供了关于基因组、化学物质、代谢途径、疾病和药物等方面的信息。KEGG数据库涵盖了多个生物学领域，包括基因组学、代谢组学、药物发现等，为研究人员提供了丰富的生物信息资源，用于理解生物系统的结构、功能和动态变化。
特点：
- 整合基因组、代谢通路、药物和疾病信息
- 覆盖几乎所有物种
- 提供详细的通路图
访问地址：https://www.kegg.jp/

2. GO（Gene Ontology）

图3 GO

简介：是一种用于描述基因功能和基因组学数据的标准化语言和分类体系。它是一个用于注释基因和蛋白质的功能、过程和细胞组件的分类系统。GO 分为三个主要的分类：分子功能（Molecular Function）、生物过程（Biological Process）和细胞组件（Cellular Component）。研究人员可以使用GO标准对基因和蛋白质进行分类和注释，从而更好地理解它们在生物学过程中的作用和相互关系。
三大分类：
- 生物学过程（Biological Process, BP）
- 分子功能（Molecular Function, MF）
- 细胞组件（Cellular Component, CC）
访问地址：http://geneontology.org/

3. Reactome（Reactome Pathway Database）

图4 Reactome

简介：Reactome 是一个开源、手工注释、经过同行评议的生物通路数据库。它提供了从小分子代谢到转录调控、信号转导和复杂的细胞过程等多个层面的生物学通路信息。
特点：
- 提供直观的通路可视化工具
- 包含30多个物种信息
- 定期更新维护
访问地址：https://reactome.org/

4. DO（Disease Ontology）数据库

图5 DO

简介：Disease Ontology 是一个标准化的疾病本体知识库，旨在提供人类疾病的统一命名和分类体系。它将人类疾病与基因组、变异信息、药物反应等多个生物医学领域关联起来。
特点：
- 标准化疾病术语
- 疾病表型特征描述
- 相关医学词汇整合
访问地址：http://www.disease-ontology.org/

3. GSEA的原理

GSEA的核心思想是将所有基因按照表达差异的显著性进行排序，然后检测特定基因集在这个排序中的富集程度。主要步骤包括：

基因排序：根据差异表达分析的统计量（如logFC或p值）对基因进行排序。
富集评分（Enrichment Score, ES）：计算每个基因集在基因排序中的富集分数。
统计显著性评估：使用置换检验（Permutation Test）来评估富集分数的显著性。
多重检验校正：对显著性结果进行FDR校正，减少假阳性。

4. GSEA的分析步骤

数据准备：整理表达矩阵和分组信息。
基因集选择：选择合适的基因集（如KEGG、GO、Hallmark等）。
运行GSEA：使用GSEA工具对数据进行富集分析。
结果可视化：生成富集曲线、热图等图表。
生物学解释：对显著富集的通路或功能进行解读。

5. 示例操作

# 加载必要的包
library(ReactomePA)
library(tidyverse)
library(data.table)
library(org.Hs.eg.db)    # 人类基因注释数据库
library(clusterProfiler) # GSEA核心包
library(enrichplot)      # 富集分析可视化包
library(dplyr)

# 读取数据
gene_data <- read.csv("APT_vs_AHa.csv", header = TRUE, row.names = 1)
# 注意：输入文件必须包含log2FoldChange列

# 准备基因列表
gene_list <- data.frame(
  Gene = rownames(gene_data),
  Score = gene_data$log2FoldChange  # 重要：这里使用log2FoldChange作为排序依据
) %>%
  filter(is.finite(Score))  # 过滤掉无穷大值

# 转换基因ID
gene_ids <- mapIds(org.Hs.eg.db,
                   keys = gene_list$Gene,
                   column = "ENTREZID",    # 转换为ENTREZID格式
                   keytype = "SYMBOL",     # 输入为基因符号格式
                   multiVals = "first")    # 对于多个匹配只取第一个

# 合并数据并创建排序的基因列表
gene_list$ENTREZID <- gene_ids
final_list <- gene_list %>%
  filter(!is.na(ENTREZID)) %>%
  arrange(desc(Score))

# 创建命名向量
gene_scores <- final_list$Score
names(gene_scores) <- final_list$ENTREZID

# 执行GSEA分析
# GO分析
go_result <- gseGO(gene_scores, 
                   OrgDb = org.Hs.eg.db,
                   ont = "ALL",           # 可选 "BP", "CC", "MF"
                   minGSSize = 10,        # 最小基因集大小
                   maxGSSize = 500,       # 最大基因集大小
                   pvalueCutoff = 0.05)   # p值阈值

# KEGG分析
kegg_result <- gseKEGG(gene_scores,
                       organism = "hsa",   # 人类的organism代码
                       minGSSize = 10,
                       maxGSSize = 500,
                       pvalueCutoff = 0.05)

# Reactome分析
reactome_result <- gsePathway(gene_scores,
                             organism = "human",
                             minGSSize = 10,
                             maxGSSize = 500,
                             pvalueCutoff = 0.05)

# 保存结果
write.csv(as.data.frame(go_result), "test_analysis_GO.csv", row.names = FALSE)
write.csv(as.data.frame(kegg_result), "test_analysis_KEGG.csv", row.names = FALSE)
write.csv(as.data.frame(reactome_result), "test_analysis_Reactome.csv", row.names = FALSE)

# 绘图
pdf("GSEA_plots.pdf", width = 10, height = 8)  # 可根据需要调整图片大小

# GO富集分析图
dotplot(go_result, showCategory = 20) + ggtitle("GO Enrichment Analysis")

# KEGG通路富集分析图
ridgeplot(kegg_result, showCategory = 15) + ggtitle("KEGG Pathway Enrichment")

# Reactome通路富集分析图
dotplot(reactome_result, showCategory = 20) + ggtitle("Reactome Pathway Enrichment")

dev.off()