富集分析，Over-representation analysis（ORA）是用来分析一个基因列表（例如差异基因，或者任意形式的基因列表，本文中我们以差异基因，也就是DE基因为例）是否富集在一个基因集合中的方法。我们一般使用超几何分布来计算p值，对应的函数是phyper()，参数为：
phyper(x, m, n, k, lower.tail = FALSE)
其中
x是基因集合中DE基因的个数，m是基因集合中基因的总数，n是不属于此基因集合的基因总数，k是DE基因的总数。
参数lower.tail = FALSE计算的是$\mathrm{Pr}(X > x)$。在计算p值时，我们常也包含$X = x$的情形，此时最终的p值应当为$\mathrm{Pr}(X \ge x)$。那么在使用phyper()时，我们需要将x调整为x - 1。
phyper(x - 1, m, n, k, lower.tail = FALSE)
使用for/sapply循环
已知一个DE基因列表和一个基因集合，很容易能够计算出phyper()所需的四个参数的值。首先我们编写一个ora_single()函数用来计算对一个单独基因集合做富集分析的p值。如下代码所示，ora_single()需要三个参数：一个DE基因向量genes，一个基因集合基因向量gene_set和一个背景基因向量universe。在很多情况下，并不能保证背景基因能够完全包含DE基因或者基因集合，因此，我们需要手动去除不包含在背景基因中的基因（使用intersect()函数）。
在下面的例子中，我们假设所有的基因向量中的基因都具有相同类型的ID类型。
ora_single = function(genes, gene_set, universe) {
n_universe = length(universe)
# 将genes和gene_set都限制在universe中
genes = intersect(genes, universe)
gene_set = intersect(gene_set, universe)
x = length(intersect(genes, gene_set)) # DE genes in the gene set
m = length(gene_set) # total genes in the gene set
n = n_universe - m # total genes not in the gene set
k = length(genes) # total DE genes
phyper(x - 1, m, n, k, lower.tail = FALSE)
}
然后，对于一组基因集合，我们只需要使用sapply()或者for循环将ora_single()应用到每一个基因集合上。在ora_v1()中，我们假设基因集合的格式是一个list，其中每一个成员向量表示一个基因集合。ora_v1()返回一个p值向量。
# version 1
ora_v1 = function(genes, gene_sets, universe) {
sapply(gene_sets, function(x) ora_single(genes, x, universe))
}
可见在R中编写一个支持ORA分析的函数非常简单。为了执行ora_v1()，我们使用GO基因集合，使用蛋白编码基因作为背景基因，从其中随机抽取1000个基因作为DE基因。
library(org.Hs.eg.db)
gs = as.list(org.Hs.egGO2ALLEGS)
gs = lapply(gs, unique) # 基因有可能重复
library(GO.db)
gs = gs[Ontology(names(gs)) == "BP"] # 取GO BP基因集合
pc_genes = select(org.Hs.eg.db, key = "protein-coding",
keytype = "GENETYPE", column = c("ENTREZID"))[, 2] # 获取所有的protein-coding基因
genes = sample(pc_genes, 1000)
注意在上述的例子中pc_genes作为背景基因，但是并不是所有gs中的基因均能被包含在pc_genes中。因为在GO注释中，也包含了非蛋白编码基因，如microRNA。
现在让我们运行ora_v1()：
system.time(p1 <- ora_v1(genes, gs, pc_genes))
## 用户 系统 流逝
## 10.946 3.357 14.314
使用向量化的hyper()
在版本1中，我们使用循环的方式进行代码编写。在R中，我们往往推荐向量化计算。注意phyper()同样也支持以向量作为参数，同时计算多个p值。在版本2中，在运行phyper()之前，我们直接生成参数x，m和n向量。
# version 2
ora_v2 = function(genes, gene_sets, universe) {
genes = intersect(genes, universe)
gene_sets = lapply(gene_sets, function(x) intersect(x, universe))
n_universe = length(universe)
n_genes = length(genes)
# 注意现在x, m, n, k均为向量
x = sapply(gene_sets, function(x) length(intersect(x, genes)))
m = sapply(gene_sets, length)
n = n_universe - m
k = n_genes
phyper(x - 1, m, n, k, lower.tail = FALSE)
}
使用相同的genes，gs和pc_genes变量，运行ora_v2()。
system.time(p2 <- ora_v2(genes, gs, pc_genes))
## 用户 系统 流逝
## 4.954 1.090 6.049
可见和ora_v1()相比，ora_v2()提速了一倍。
再看intersect()函数
下面我们对ora_v2()进行代码profiling，看看其中最耗时的是哪一部分。
Rprof()
p2 = ora_v2(genes, gs, pc_genes)
Rprof(NULL)
summaryRprof()$by.self
## self.time self.pct total.time total.pct
## "base::intersect" 5.04 93.33 5.32 98.52
## "duplicated.default" 0.18 3.33 0.18 3.33
## "FUN" 0.04 0.74 5.40 100.00
## "intersect" 0.04 0.74 5.36 99.26
## "duplicated" 0.04 0.74 0.22 4.07
## "numeric" 0.04 0.74 0.04 0.74
## "c" 0.02 0.37 0.02 0.37
可见intersect()是代码中最耗时的一步。
数学中，intersect(A, B)和intersect(B, A)是等价的，但是在R中，intersect()内部是基于哈希表，其第一个参数是输入向量，而第二个参数会被视为“字典”用来内部构建一个哈希表。
让我们做如下实验：
system.time(lapply(gs, function(x) intersect(x, pc_genes)))
## 用户 系统 流逝
## 4.402 0.866 5.272
system.time(lapply(gs, function(x) intersect(pc_genes, x)))
## 用户 系统 流逝
## 11.896 2.114 14.021
在上述代码的第一行中，因为pc_genes只需要在第一次运行intersect()时生成一份临时的哈希表，而在以后的运行时可以重复使用。而在第二行中，由于intersect()的第二个参数为x，在每次循环中都会发生变化，因此每次都会重新生成关于x的新的哈希表。
那么到这里的讨论告诉我们，如果以后我们需要对一组向量对同一个向量做intersect，永远把那个单独的向量作为intersect()的第二个参数。
fastmatch包
如果你查看intersect()的源代码，会发现其内部使用了match()函数。在R中，fastmatch提供了一个更快的match操作，名为fmatch()，或者与%in%类似的%fin%。我们来把ora_v2()修改一下：
# version 3
library(fastmatch)
ora_v3 = function(genes, gene_sets, universe) {
genes = genes[genes %fin% universe]
gene_sets = lapply(gene_sets, function(x) x[x %fin% universe])
n_universe = length(universe)
n_genes = length(genes)
# 注意现在x, m, n, k均为向量
x = sapply(gene_sets, function(x) sum(x %fin% genes))
m = sapply(gene_sets, length)
n = n_universe - m
k = n_genes
phyper(x - 1, m, n, k, lower.tail = FALSE)
}
system.time(p3 <- ora_v3(genes, gs, pc_genes))
## 用户 系统 流逝
## 0.134 0.002 0.136
速度得到了超大幅度的提升！
三个方法得到的p值是相同的。
identical(p1, p2)
## [1] TRUE
identical(p1, p3)
## [1] TRUE
结语
下面综合比较一下：
library(microbenchmark)
microbenchmark(
"loop(v1)" = ora_v1(genes, gs, pc_genes),
"vectorized(v2)" = ora_v2(genes, gs, pc_genes),
"fastmatch(v3)" = ora_v3(genes, gs, pc_genes),
times = 10
)
## Unit: milliseconds
## expr min lq mean median uq max neval
## loop(v1) 15239.0657 16168.2848 16441.5549 16347.8581 16638.9054 17934.5915 10
## vectorized(v2) 6624.4581 6650.8479 6988.7350 6782.4591 7039.8733 8211.7624 10
## fastmatch(v3) 107.9137 110.6585 120.8098 116.0022 126.0807 147.1155 10
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS 26.3.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] microbenchmark_1.5.0 fastmatch_1.1-4 GO.db_3.17.0 org.Hs.eg.db_3.17.0
## [5] AnnotationDbi_1.62.2 IRanges_2.36.0 S4Vectors_0.40.2 Biobase_2.60.0
## [9] BiocGenerics_0.48.1 knitr_1.45
##
## loaded via a namespace (and not attached):
## [1] bit_4.0.5 jsonlite_1.8.8 compiler_4.3.3 crayon_1.5.2
## [5] blob_1.2.4 bitops_1.0-7 Biostrings_2.68.1 jquerylib_0.1.4
## [9] png_0.1-8 yaml_2.3.8 fastmap_1.1.1 R6_2.5.1
## [13] XVector_0.40.0 GenomeInfoDb_1.36.4 bookdown_0.39 GenomeInfoDbData_1.2.10
## [17] DBI_1.2.2 bslib_0.7.0 rlang_1.1.3 KEGGREST_1.40.1
## [21] cachem_1.0.8 xfun_0.43 sass_0.4.9 bit64_4.0.5
## [25] RSQLite_2.3.6 memoise_2.0.1 cli_3.6.2 zlibbioc_1.46.0
## [29] digest_0.6.35 lifecycle_1.0.4 vctrs_0.6.5 evaluate_0.23
## [33] blogdown_1.19 RCurl_1.98-1.14 rmarkdown_2.26 httr_1.4.7
## [37] pkgconfig_2.0.3 tools_4.3.3 htmltools_0.5.8.1