ixxmu
commented
1 month ago bulk RNA-seq | 下游分析 | 功能富集分析 ORA - clusterProfiler by 被炸熟的虾
在完成了bulk RNA-seq数据质控与差异分析后，通过设置不同的阈值，我们获得的差异基因数量可能从几个到上千个不等。
bulk RNA-seq | 下游分析 | 分析前评估 bulk RNA-seq | 下游分析 | 差异分析 DESeq2 bulk RNA-seq | 下游分析 | 差异分析 DESeq2-补1-统计结果缺失值
富集指的是使用某种方法使目的物质占比增大。如使用磁珠富集T细胞，最终的结果是使T细胞的占比增大。功能富集分析就是在目标基因集中寻找占比高的功能，从而推断目标基因集的关键功能。例如对差异基因集进行功能富集分析，寻找差异基因集“富集”的关键功能，进而推出样本间重要的功能差异。
严谨来说，基于数据来源和算法分类，富集分析方法主要有过度代表性分析（Over Representation Analysis, ORA）、功能分类打分（Functional Class Scoring，FCS）、通路拓扑结构（Pathway Topology, PT）和网络拓扑结构（Network Topology, NT）四种。
我们常听说的GO富集分析，KEGG富集分析，WikiPathways富集分析等直接以基因功能注释数据库取名的富集分析采用的都是ORA算法；而基因集富集分析（GSEA）、基因集变异分析（GSVA）、单样本基因集富集分析（ssGSEA）等，采用的是FCS算法。
很久前写过一篇笔记：记录 8 种常见的基因功能富集在线工具（Metascape、Enrichtool、DAVID、GeneTrail、G:profiler、WEBgestalt、KOBAS与modEnrichr）。曾有大佬比较不同工具的分析效能，强调数据集的更新时间、内置数据集种类是选择富集分析工具时的主要参考指标，尤其是前者，如果不及时更新，很容易遗漏重要信息（PMID: 27575621）。
clusterProfiler包可以抓取GO、KEGG等数据库的最新信息，因此可以说是R语言富集分析（ORA）的不二选择。今天以我们以DESeq2差异分析结果为例，记录一下基础流程。
library(clusterProfiler)
packageVersion("clusterProfiler")#[1]‘4.11.1’
提取差异表达基因：
setwd("D:/bioinformatics/2.NGS/1.bulk RNA-seq")
load("./data/Merge2.Rdata")
load("./data/DEG_DESeq2.Rdata")
# 设定阈值，筛选显著上下调差异基因
logFC <- 2
Pvalue <- 0.01
DEG_DESeq2$Group <- ifelse(DEG_DESeq2$log2FoldChange > logFC & DEG_DESeq2$padj < Pvalue,"Up",
ifelse(DEG_DESeq2$log2FoldChange < -logFC & DEG_DESeq2$padj < Pvalue,
"Down","Stable"))
table(DEG_DESeq2$Group)
#Down Stable Up
#1824 18306 506
一、ORA 原理
（想了一下还是记录了一点）
ORA算法的主要思想是统计评估目标基因集（例如差异基因集）中属于特定功能条目（term）的基因的比例。如果这种比例出现的可能性（期望）很低，远低于随机事件发生的概率，那么我们就认为我们观察到的富集情况是“显著”的。
有点拗口，我们套用一下经典的“摸球”问题：
蓝色方框中的球是本次测序获得的所有基因【共N个】，某个特定通路P包括基因共M个（红色小球）；绿色圆圈是一次摸球事件，用来表示做了一次差异分析得到的差异基因【共n个】，这些基因中，有属于通路P的（红色球）【共k个】，有不属于的（黑色球）。问这次分析结果中，通路P的基因占比是否显著过高？
无放回的抽样问题，至少有k个小球属于通路P的概率可以用超几何检验来计算。若p ≤0.05，说明在随机状态下，当前目标基因集中含有k个以上的P通路基因是极不可能的，但是当下却发生了，该事件就是一个“有意义”的事件，因此P通路显著富集。反之，富集不显著。
推荐看一下Y叔的讲解：富集分析的p值是怎么算出来的？
二、GO & KEGG
Gene Ontology (GO)是一个用于描述基因和蛋白质功能的标准化系统，GO数据库将基因/蛋白的属性分为3大类：
生物过程（Biological Process,BP）：描述基因产物参与的生物学过程，例如细胞周期、信号转导等。 分子功能（Molecular Function,MF）：描述基因产物的具体分子活动，例如酶活性、结合活动等。 细胞组分（Cellular Component,CC）：描述基因产物的细胞位置或组成，例如细胞膜、细胞核等。
KEGG (Kyoto Encyclopedia of Genes and Genomes)数据库的侧重点是描述生物体内的代谢通路，以及这些通路在不同生物体系中的变化和调控。
两者通常会结合使用，以全面理解基因和蛋白质的功能及其在生物体内的作用。对于这些最常用的注释数据库，clusterprofiler包分别提供了专用函数，函数语法很简单：
enrichKEGG(gene, #待富集的基因列表，只接受entrez gene id
organism = "hsa", #指定物种
keyType = "kegg",
pvalueCutoff = 0.05, #指定 p 值阈值（可指定 1 以输出全部）
qvalueCutoff = 0.2, #指定 q 值阈值（可指定 1 以输出全部）
pAdjustMethod = "BH",#指定 p 值校正方法，比如holm,hochberg,hommel,bonferroni,BH,BY,fdr,none
universe,
minGSSize = 10, #通路最少基因数量
maxGSSize = 500, #通路最大基因数量
use_internal_data = FALSE) #使用本地KEGG.db还是在线抓取
enrichGO(gene, #待富集的基因列表
OrgDb, #指定物种的基因数据库
keyType = "ENTREZID", #指定给定的基因名称类型
ont = "MF", #GO Ontology，可选 BP、MF、CC，也可以指定ALL同时计算3者
pvalueCutoff = 0.05, #指定 p 值阈值（可指定 1 以输出全部）
qvalueCutoff = 0.2, #指定 q 值阈值（可指定 1 以输出全部）
pAdjustMethod = "BH", #指定 p 值校正方法，比如holm,hochberg,hommel,bonferroni,BH,BY,fdr,none
universe,
minGSSize = 10, #Ontology最少基因数量
maxGSSize = 500, #Ontology最大基因数量
readable = FALSE, #显示结果是否将gene ID转为gene Name
pool = FALSE)KEGG富集分析
#BiocManager::install("org.Hs.eg.db")
#library(org.Hs.eg.db)
#packageVersion("org.Hs.eg.db")#[1]‘3.15.0’
#BiocManager::install("org.Mm.eg.db")
library(org.Mm.eg.db)
packageVersion("org.Mm.eg.db")#[1]‘3.15.0’
ID_entrez <- bitr(rownames(DEG_DESeq2[DEG_DESeq2$Group == "Up",]),#以上调基因为例
fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的
toType = "ENTREZID", #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = "org.Mm.eg.db")#Orgdb是指对应的注释包是哪个
#'select()' returned 1:1 mapping between keys and columns
#Warning message:
#In bitr(rownames(DEG_DESeq2[DEG_DESeq2$Group == "Up", ]), fromType = "SYMBOL", :
# 4.35% of input gene IDs are fail to map...
head(ID_entrez)
# SYMBOL ENTREZID
#1 Ltf 17002
#2 Lcn2 16819
#3 Farsb 23874
#4 Oas2 246728
#5 Clca3a2 80797
#6 Aqp5 11830
个人习惯，如果有ENSEMBL，可以也尝试一下，看看转换率会不会高一丢丢，本例就没有。
head(ids)
ID_ENSEMBL <- ids[ids$SYMBOL %in% rownames(DEG_DESeq2[DEG_DESeq2$Group == "Up",]),"ENSEMBL"]
ID_ENSEMBL <- sapply(strsplit(ID_ENSEMBL,"[.]"),"[",1)
library(org.Mm.eg.db)
ID_entrez2 <- bitr(ID_ENSEMBL,
fromType = "ENSEMBL", #fromType是指你的数据ID类型是属于哪一类的
toType = "ENTREZID", #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = "org.Mm.eg.db")#Orgdb是指对应的注释包是哪个
#'select()' returned 1:1 mapping between keys and columns
#Warning message:
#In bitr(rownames(DEG_DESeq2[DEG_DESeq2$Group == "Up", ]), fromType = "SYMBOL", :
# 6.32% of input gene IDs are fail to map...
head(ID_entrez2)
# ENSEMBL ENTREZID
#1 ENSMUSG00000023902 332221
#2 ENSMUSG00000027514 58203
#3 ENSMUSG00000018849 211652
#4 ENSMUSG00000035112 69847
#5 ENSMUSG00000040154 209232
#6 ENSMUSG00000017723 67701
packageVersion("clusterProfiler")#[1]‘4.11.1’
Enrichment_KEGG <- enrichKEGG(gene = unique(ID_entrez$ENTREZID),
organism = "mmu", #Mus musculus (house mouse)
pvalueCutoff = 0.5,
qvalueCutoff = 0.5,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 10000)
head(Enrichment_KEGG)
library(openxlsx)
write.xlsx(Enrichment_KEGG@result,file = "./Enrichment_KEGG.xlsx",colNames = TRUE,rowNames = TRUE)
ID：pathway的ID名； GeneRatio：差异基因中富集到该pathway的基因数目/富集到所有pathway的总差异基因数目； BgRatio：该pathway的基因数目/所有pathway的总背景基因数目； pvalue：显著性p值； p.adjust：多重检验矫正后的p值； qvalue：应用FDR校正的q值； Count：富集到该pathway的基因数目； geneID：富集到这个pathway的差异基因（ENTREZ ID）；
我们可以根据GeneRatio与BgRatio信息手动计算出显著性p值：
如果分析物种是人或大鼠，我们修改参数organism为hsa（Homo sapiens，human）或rno（Rattus norvegicus，rat）即可。
enrichKEGG函数适配的物种目录如下：
新的clusterprofiler版本返回结果新增了category信息，更新后，我们就可以很容易的按照通路所属类别展示结果了：
library(ggplot2)
p1 <- dotplot(Enrichment_KEGG,x = "GeneRatio",color = "p.adjust",showCategory = 10) +
facet_grid(rows = vars(category),scales = 'free_y',space = 'free_y') +
scale_color_gradientn(colors = c('#BF1E27','#FEB466','#F9FCCB','#6296C5','#38489D')) +
theme_bw(base_size = 18)
p2 <- barplot(Enrichment_KEGG,x = "GeneRatio",color = "p.adjust",showCategory = 10) +
facet_grid(rows = vars(category),scales = 'free_y',space = 'free_y') +
scale_fill_gradientn(colors = c('#BF1E27','#FEB466','#F9FCCB','#6296C5','#38489D')) +
theme_bw(base_size = 18)
setReadable函数将geneID列的ENTREZ ID重转为Symbol：Enrichment_KEGG2 <- setReadable(Enrichment_KEGG,OrgDb = org.Mm.eg.db,keyType = "ENTREZID")
head(Enrichment_KEGG2)
library(openxlsx)
write.xlsx(Enrichment_KEGG2@result,file = "./Enrichment_KEGG2.xlsx",colNames = TRUE,rowNames = TRUE)
富集因子：差异基因中富集到该通路的基因数（count）/背景基因中富集到该通路的基因数（BgRatio的分子）； 富集倍数：GeneRatio/BgRatio；
library(dplyr)
#计算Rich Factor（富集因子）：
Enrichment_KEGG2 <- mutate(Enrichment_KEGG,
RichFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
#计算Fold Enrichment（富集倍数）：
Enrichment_KEGG2 <- mutate(Enrichment_KEGG2,
FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
head(Enrichment_KEGG2@result)GO富集分析
enrichGO函数可以接受多种ID，包括SYMBOL与ENTREZ ID：#####SYMBOL
Enrichment_GOBP <- enrichGO(gene = rownames(DEG_DESeq2[DEG_DESeq2$Group == "Up",]),
OrgDb = org.Mm.eg.db, #记得修改对应org包
ont = "BP", #"BP"，"MF","CC"三选一,或"ALL"合并
keyType = "SYMBOL",
pvalueCutoff = 0.1,
qvalueCutoff = 0.1,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 10000,
readable = FALSE)
head(Enrichment_GOBP@result)
#####ENTREZ ID
Enrichment_GOBP2 <- enrichGO(gene = unique(ID_entrez$ENTREZID),
OrgDb = org.Mm.eg.db, #记得修改对应org包
ont = "BP", #"BP"，"MF","CC"三选一,或"ALL"合并
keyType = "ENTREZID",
pvalueCutoff = 0.1,
qvalueCutoff = 0.1,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 10000,
readable = TRUE)
head(Enrichment_GOBP2@result)
library(openxlsx)
write.xlsx(Enrichment_GOBP@result,file = "./Enrichment_GOBP.xlsx",colNames = TRUE,rowNames = TRUE)
ID：GO term ID号； Description：GO term的功能描述； GeneRatio：差异基因中富集到该GO term的基因数/差异基因总数； BgRatio：背景基因中富集到该GO term的基因数/背景基因总数； pvalue：显著性p值； p.adjust：多重检验矫正后的p值； qvalue：应用FDR校正的q值； geneID：富集到这个GO term的差异基因（已转换为gene name）； Count：富集到该GO term的基因数量。
org包记录且被注释数据库收录的基因才会被纳入分析（GeneRatio列的总基因数目），因此，这两种方法的分析结果是完全相同的。enrichKEGG函数，enrichGO函数可以直接通过设置参数readable = FALSE将geneID列的ENTREZ ID重转为SYMBOL（如果输入基因列表就是SYMBOL无需转换输出即是SYMBOL）。移除冗余term
#BiocManager::install("topGO")
library(topGO)
plotGOgraph(Enrichment_GOBP)
#goplot(Enrichment_GOBP) #来自enrichplot
详细讲解可以阅读：GO数据库的分类层级说明
clusterProfiler提供了一个simplify函数来消除这种冗余的GO term。Enrichment_GOBP2 <- clusterProfiler::simplify(Enrichment_GOBP, cutoff = 0.7,
by = "p.adjust",select_fun = min)
dim(Enrichment_GOBP@result) #[1] 4929 9
dim(Enrichment_GOBP2@result) #[1] 228 9
三、其他通路数据库
WikiPathways是一个由社区驱动的生物通路数据库，涵盖了不同物种的生物学通路。它是一个开放平台，用户可以提交和编辑通路信息。由于其社区驱动的特点，它可以包含更广泛的或新兴的通路信息。
Pathway Commons是一个汇集了来自众多资源的多个物种的公开获取通路数据，整合了多种通路数据库。包括生化反应，复合物组装、转运和催化事件，以及涉及蛋白质复合物、DNA、RNA、小分子化合物和配合物在内的物理相互作用。
##ORA analysis for WikiPathways
enrichWP(gene, organism, ...)
#enrichWP支持物种
get_wp_organisms()
# [1] "Anopheles gambiae" "Arabidopsis thaliana" "Bos taurus"
# [4] "Caenorhabditis elegans" "Canis familiaris" "Danio rerio"
# [7] "Drosophila melanogaster" "Equus caballus" "Gallus gallus"
#[10] "Homo sapiens" "Mus musculus" "Pan troglodytes"
#[13] "Populus trichocarpa" "Rattus norvegicus" "Saccharomyces cerevisiae"
#[16] "Solanum lycopersicum" "Sus scrofa" "Zea mays"
##ORA analysis for Pathway Commons
enrichPC(gene, source, keyType = "hgnc", ...)
Enrichment_WP <- enrichWP(gene = unique(ID_entrez$ENTREZID),
#OrgDb = org.Mm.eg.db, #记得修改对应org包
organism = "Mus musculus",
pvalueCutoff = 0.1,
qvalueCutoff = 0.1,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 10000)
head(Enrichment_WP@result)
Enrichment_WP2 <- setReadable(Enrichment_WP,OrgDb = org.Mm.eg.db,keyType = "ENTREZID")
head(Enrichment_WP2@result)
save(Enrichment_KEGG,Enrichment_KEGG2,
Enrichment_GOBP,Enrichment_GOBP2,
Enrichment_WP,Enrichment_WP2,
file = "./ORA.Rdata")
四、自定义注释基因集
GeneSet <- read.gmt("D:/data/MSigdb/mouse/msigdb_v2023.2.Mm_GMTs/m2.cp.reactome.v2023.2.Mm.entrez.gmt")
head(GeneSet)
# term gene
#1 REACTOME_ABACAVIR_ADME 75894
#2 REACTOME_ABACAVIR_ADME 11522
#3 REACTOME_ABACAVIR_ADME 76952
#4 REACTOME_ABACAVIR_ADME 20517
#5 REACTOME_ABACAVIR_ADME 20518
#6 REACTOME_ABACAVIR_ADME 20519
基本参数依旧是类似的：
Enrichment_Reactome <- enricher(unique(ID_entrez$ENTREZID), #待富集的基因列表
pvalueCutoff = 0.5, #指定 p 值阈值（可指定 1 以输出全部）
qvalueCutoff = 0.5, #指定 q 值阈值（可指定 1 以输出全部）
pAdjustMethod = "BH", #指定 p 值校正方法
# universe = rownames(DEG_DESeq2), #背景基因
minGSSize = 5, #term最少基因数量
maxGSSize = 10000, #term最大基因数量
TERM2GENE = GeneSet) #annotation of term TO gene mapping
这里有两个个人认为需要注意的点：
Enrichment_Reactome2 <- enricher(unique(ID_entrez$ENTREZID),
pvalueCutoff = 1,
qvalueCutoff = 1,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 10000,
TERM2GENE = GeneSet[GeneSet$term == "REACTOME_DEVELOPMENTAL_BIOLOGY",])
length(unique(GeneSet$gene)) #[1] 9082
length(unique(GeneSet[GeneSet$term == "REACTOME_DEVELOPMENTAL_BIOLOGY","gene"])) #[1] 484
Enrichment_Reactome@result["REACTOME_DEVELOPMENTAL_BIOLOGY",3:7]
# GeneRatio BgRatio pvalue p.adjust qvalue
#REACTOME_DEVELOPMENTAL_BIOLOGY 28/221 484/9082 1.629723e-05 0.002536935 0.002373105
Enrichment_Reactome2@result["REACTOME_DEVELOPMENTAL_BIOLOGY",3:7]
# GeneRatio BgRatio pvalue p.adjust qvalue
#REACTOME_DEVELOPMENTAL_BIOLOGY 28/28 484/484 1 1 NA
五、结果可视化
，有其他需要可以在评论区留言
。参考文章：
1、生信技工：富集分析：（一）概述
2、联川生物：【联川沧海】富集分析ORA（通用富集分析）、FCS（GSEA、ssGSEA、GSVA）一网打尽
3、生信技工：富集分析：（一）概述
4、观科研：Bulk RNA-seq | 第7期. GO富集，我收藏了
5、基迪奥生物：KEGG富集分析及可视化，一把子拿捏！
6、基迪奥生物：GO富集分析及可视化，一把子拿捏！
https://mp.weixin.qq.com/s/in1O4Wte1jz8bqGwuAaAdw