当前位置：首页 > article >正文

解析富集分析中的过表达分析（ORA）：原理、应用与优化

article 2025/3/14 15:47:50

TLDR: 在基因富集分析中，超几何分布被用于评估特定功能基因集在差异表达基因列表中的过表达程度，而 Fisher 精确检验（Fisher’s Exact Test）正是基于超几何分布计算 p 值。本文详细解析超几何分布的数学原理、概率质量函数（PMF）、累积分布函数（CDF），并探讨其在富集分析中的应用，包括 R 语言实现、算法优化。

什么是超几何分布

基础概念

超几何分布（Hypergeometric distribution）是统计学上一种离散概率分布。它描述了由有限个对象中抽出n个对象，成功抽出k次指定种类的对象的概率（抽出不放回（without replacement））。例如在有N个样本，其中K个是不及格的。超几何分布描述了在该N个样本中抽出n个，其中k个是不及格的个数：

表示所有在N个样本中抽出n个的方法数目；

表示在K个样本中，抽出k个的方法数目；

表示在N-K个样本中，抽出n-k个的方法数目。

若随机变量X服从参数为n, K, N的超几何分布，则记为

X∼H(n, K, N)。

在公式中，括号里的表达式表示组合数（也叫"二项系数"），它的数学含义是，在a个元素中，不考虑顺序地选择b个元素的所有可能的方式数。计算公式是：

其中，a!表示阶乘，即从1乘到a的积。

例子

我们从一副扑克牌出发（总共 52 张牌），其中有 4 张是 A（A♠, A♥, A♦, A♣），剩下的 48 张是其他牌。

从 52 张牌中随机抽出 4 张，想知道这 4 张全是 A 的概率。

从这 52 张牌中随机抽出 26 张， 26 张有4个A 的概率。

从这 52 张牌中随机抽出 26 张，这 26 张中1、2、3、4张A的概率

从 52 张牌中抽取 26 张牌时，抽到 1、2、3、4 张 A 的概率。可以看到，抽到 2 张 A 的概率最大，其他情况依次递减。

从这 52 张牌中随机抽出 26 张，抽到 1 到 13 张方片的概率

概率质量函数

概率质量函数（Probability Mass Function, PMF） 是用来描述离散型随机变量分布的一个函数，它表示随机变量 X 取某个具体值 k 的概率 P(X=k)。 PMF 是离散随机变量的核心概念，用于描述随机变量的具体取值概率。和连续型随机变量的概率密度函数（PDF）类似，但 PMF 只适用于离散的情况，比如骰子点数、抽到某种牌的数量等。

给定参数 N=52, K=13, n=10 的概率质量函数（PMF），展示了从 52 张牌中抽取 10 张牌时，抽到 0 到 10 张方片的概率分布。

累积分布函数

累积分布函数是用来描述随机变量 X 的分布特性的一种函数，它表示随机变量 X 的取值小于或等于某个值 x 的累计概率。

从图中可以看出：

• P(X≤2) 的累计概率约为 0.4，表示抽到 2 张或更少方片的概率约为 40%。
• P(X≤4)的累计概率约为 0.9，表示抽到 4 张或更少方片的概率约为 90%。

与二项分布、几何分布的联系

几何分布和二项分布内容此处不详述。

• 超几何分布和二项分布是相似的，只是抽样方式不同。
• 几何分布关注的是第一次成功的试验次数，与二项分布在目标变量上不同。

Fisher 精确检验

Fisher 精确检验本质上是基于超几何分布公式计算的。当提到 ORA 的统计检验方法时，通常会直接说使用 Fisher 精确检验。但 Fisher 精确检验的计算逻辑依赖于超几何分布

过表达分析（ORA)

基本原理

在过表达分析（Over-Representation Analysis, ORA）中，统计检验的核心假设和计算方式与Fisher 精确检验和超几何分布密切相关。 ORA的目标是判断在一个感兴趣的基因集（例如差异表达基因）中，某个特定功能类别（如 KEGG通路）是否显著富集。

假设：

• 零假设（H0）：感兴趣的基因集（差异表达基因）中，某个功能类别的基因数与其在背景基因集中的分布一致，没有富集。
• 备择假设（H1）：感兴趣的基因集中，某个功能类别的基因数显著高于背景中的比例，存在富集。

P(X≥k) 是在零假设成立时，比观察值 k 更少见的情况发生的总概率。

如果这个概率很小，说明观察到的 k 和零假设的随机分布差异显著，可以拒绝零假设。

注意：这里我们要找到的是X≥k的累计概率，直接反映了观察结果及更极端情况的可能性，如果 p-值很小（如 <0.05），说明观察值极难由随机分布产生，因此可以拒绝零假设，认为有富集现象。

具体实现

R中的超几何分布概率计算

R 中对应的函数是 phyper()，其参数含义如下：

phyper(x, m, n, k, lower.tail = FALSE)

其中

• x抽出的白球数量（基因集/通路中显著的基因总数)
• m白球的数量（基因集/通路中基因的总数）
• n黑球的数量（不属于该基因集/通路的基因总数）
• k抽球的数量（显著基因的总数)

在实际计算时，为了包括观察值本身，需将参数 x 调整为 x - 1，最终公式为：

phyper(x - 1, m, n, k, lower.tail = FALSE)

函数理解

假设有一个袋子，里面有 20 个球：

• 白球数量 m=8，代表白球总数（有某个特定属性的球）。
• 黑球数量 n=12，代表非白球总数（没有该属性的球）。
• 从袋子里随机抽取 5 个球（总抽取数 k=5）。

问：抽到 3 个白球（x=3）或更多白球的概率是多少？

假设有一个基因背景库，共包含 10,000 个基因，其中：

• 某个通路包含 500 个基因（m=500），相当于袋子中的白球。
• 不属于该通路的基因有 9,500 个（n=9500），相当于袋子中的黑球。
• 在实验中识别出 100 个显著基因（k=100），相当于抽取的球数。
• 其中 10 个显著基因属于该通路（x=10）。

问：在随机分布下，观察到 10 个或更多显著基因属于该通路的概率是多少？

对应超几何分布公式：

在 R 中计算：

phyper(10 - 1, 500, 9500, 100, lower.tail = FALSE)

以下代码均来自公众号“方圆之处”的文章《提速ORA富集分析20倍》，本文对部分原理进行了解释。

提速ORA富集分析20倍

btw这是另外一个生信大佬的公众号。

后续需要要做大批量的模拟分析，因此ORA的提速真是帮了大忙了。

单个基因集/通路的ORA实现

# 假设 `genes`, `gene_set` 和 `universe` 使用相同的基因 ID 类型
ora_single =function(genes, gene_set, universe){
    n_universe =length(universe)# 背景基因总数

    # 确保显著基因和基因集均属于背景基因集
    genes = intersect(genes, universe)
    gene_set = intersect(gene_set, universe)

    # 参数计算
    x =length(intersect(genes, gene_set))# 显著基因中属于该基因集的基因数
    m =length(gene_set)                    # 基因集中的总基因数
    n = n_universe - m                      # 不属于该基因集的基因数
    k =length(genes)                       # 显著基因的总数

    # 计算 p 值
    phyper(x -1, m, n, k, lower.tail =FALSE)
}

注意：在实际分析中，genes 和 gene_set 不一定完全属于背景基因集，因此需要使用 intersect() 手动去除不属于背景基因的元素，这里的输入均为ID的向量形式。

多个基因集/通路的ORA实现

对于一组基因集或通路，可以使用 sapply() 或 for 循环将 ora_single() 应用于每个基因集或通路，返回对应的 p 值。

ora_v1 = function(genes, gene_sets, universe) {
    # genes: 显著基因向量
    # gene_sets: 基因集或通路列表，每个元素是一个基因集向量
    # universe: 背景基因向量

    # 对每个基因集计算 p 值
    p_values = sapply(gene_sets, function(gene_set) {
        ora_single(genes, gene_set, universe)
    })

    return(p_values)
}

向量化输入参数提升效率

在R中，向量化计算往往是推荐的。注意phyper()同样也支持以向量作为参数。在版本2中，在运行phyper()之前，直接生成参数x，m和n向量。

# 版本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 = 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)
               }

哈希表改进交集操作提升效率

什么是哈希表

哈希表（Hash Table）是一种基于哈希函数实现的高效数据结构，用于快速地查找、插入和删除数据。它的核心思想是：将数据通过哈希函数映射到一个数组中的某个位置（称为“桶”），从而在需要查找一个元素时，直接根据哈希值找到它，而无需遍历整个数据结构。

在传统的 intersect() 实现中，需要对 x 和 universe 两个向量进行完全扫描。假设 x 的长度为 m ，universe 的长度为 n ，复杂度为 O(m×n)。而在哈希表中，查找一个元素的平均复杂度是 O(1)（常数时间复杂度）。因此，将 universe 转换为哈希表后，只需将 x 中的每个元素在哈希表中查找一次，复杂度为 O(m)。

例子： 假设 universe = c("geneA", "geneB", "geneC", "geneD")，我们构建一个哈希表：

哈希函数映射：
"geneA" → 位置 0
"geneB" → 位置 1
"geneC" → 位置 2
"geneD" → 位置 3

当我们检查某个基因（如 geneB）是否在 universe 中时，只需通过哈希函数计算出索引，直接跳转到对应位置检查即可，而不需要遍历整个 universe。

在上述的ORA分析中，每个基因集合都会重新对 universe 进行完整扫描。通过将 universe 转换为哈希表，只需对 universe 进行一次初始化，而不需要对每个基因集合重复扫描，这节省了大量时间。

当 n 和 m 都很大时（如 universe 有 10,000 个基因，x 有 1,000 个基因集合），使用哈希表可以显著减少计算时间。

为什么哈希表的查找复杂度是 O(n+m)？

在传统方法中，判断某个基因是否在背景基因集（universe）中，需要对整个 universe 逐一扫描，其复杂度为 O(n)。但在哈希表中：

• 构建阶段：
- • 将所有的背景基因（universe）插入哈希表。每个插入操作的平均复杂度为 O(1)，因此构建整个哈希表的复杂度为 O(n)（n 是 universe 的大小）。
• 查找阶段：
- • 对每个基因集合（如 x），只需逐一检查它的每个基因是否存在于哈希表中。查找每个基因的平均复杂度为 O(1)，所以对 m 个基因查找的总复杂度为 O(m)。

总复杂度为 O(n+m)，其中 n 是构建哈希表的时间，m 是查找的时间。

哈希表和署名向量、列表的区别

哈希表的功能和 R 中的 named vector 或 list 有相似之处，因为它们都可以通过“键”（名字）快速访问“值”。不过，哈希表和 named vector 之间有一些关键的不同点：

1. 相似点

• 键值对结构：二者都存储键值对，可以通过键快速定位到对应的值。

# Named vector
v <- c(A = 1, B = 2, C = 3)
v["B"]  # 输出: 2

在哈希表中类似：

cpp
std::unordered_map<std::string, int> hashTable;
hashTable["A"] = 1;
hashTable["B"] = 2;
hashTable["C"] = 3;
// hashTable["B"] 会返回 2

• 快速查找：通过键查找值的速度都比直接遍历快。R 中的 named vector 使用内部索引机制加速访问，哈希表则通过哈希函数直接找到位置。

2. 不同点

(1) 底层机制

• named vector：
- • 本质上是普通向量，只是在属性上添加了 names。
- • 查找机制是通过线性搜索或二分查找实现的（具体取决于是否有索引优化）。
- • 如果有很多元素，查找效率通常是 O(log⁡n)（二分查找）或者 O(n)（线性查找）。
• 哈希表：
- • 使用哈希函数直接将键映射到数组索引，实现快速定位。
- • 查找效率是 O(1) 平均复杂度，极端情况下（哈希冲突严重）可能退化到 O(n)，但通常很少发生。

(2) 去重

• named vector：允许重复值，也可以有重复的名字（尽管可能不是好实践）。

v <- c(A = 1, A = 2, B = 3)

访问 v["A"] 时可能会产生问题，因为 R 会返回多个结果。

• 哈希表：键是唯一的，重复插入会覆盖之前的值。

(3) 可扩展性

• named vector：扩展或删除元素的成本较高，因为需要重新分配内存和更新索引。
• 哈希表：支持动态扩展，插入或删除的平均时间复杂度是O(1)。

(4) 专用功能

• 哈希表支持复杂的键类型（如字符串、整数组合），并有更丰富的功能，例如统计键的存在性、高效合并等。
• named vector 主要是简单的键值存储，功能相对较少。

R 中更接近哈希表的结构：environment

如果需要在 R 中模仿哈希表的行为，可以使用 environment，它是 R 中类似哈希表的数据结构。示例如下：

# 创建环境
hash_env <- new.env(hash =TRUE)

# 插入键值对
hash_env[["A"]]<-1
hash_env[["B"]]<-2
hash_env[["C"]]<-3

# 查找
hash_env[["B"]]# 输出: 2

# 检查键是否存在
exists("A", envir = hash_env)# 输出: TRUE
exists("D", envir = hash_env)  # 输出: FALSE

与 named vector 不同，environment 使用哈希表实现，所以查找效率很高，适合大规模键值对操作。

注：常用的clusterprofiler，从代码上看也许是在使用此方式(environment)进行interaction的操作，但是有待通过实际数据验证。

具体实现

library(Rcpp)
sourceCpp(code = '
// [[Rcpp::plugins(cpp11)]]

#include <Rcpp.h>
#include <unordered_set>
using namespace Rcpp;

// [[Rcpp::export]]
List intersectToList(List lt, StringVector x) {

    int n = lt.size();
    List out(n);

    std::unordered_set<String> seen;
    seen.insert(x.begin(), x.end());

    for(int i = 0; i < n; i++) {
      
        StringVector v = as<StringVector>(lt[i]);
        LogicalVector l(v.size());

        std::unordered_set<String> seen2;

        for(int j = 0; j < v.size(); j ++) {
            l[j] = seen.find(v[j]) != seen.end() && seen2.insert(v[j]).second;
        }

        out[i] = v[l];
    }

    return out;
}
')

这个新编写的函数intersectToList()可以在R中使用。它接受两个参数。第一个是一个包含了若干向量的列表（lt），第二个参数是一个向量（x），其中x会和lt中的每一个向量进行intersection。在Cpp代码中，去掉了lt中每一个向量中重复的元素。

现在我们将ora_v2()中原来的这行

gene_sets = lapply(gene_sets, function(x) intersect(x, universe))

替换成优化过的代码：

gene_sets = intersectToList(gene_sets, universe)

除了对背景基因的intersection做了优化，在对每个基因集合和DE基因做intersection时，也使用intersectToList()。下面这个新函数ora_v3()基于ora_v2()，其中只修改了两行。

# 版本 3
ora_v3 =function(genes, gene_sets, universe){

  gs_names =names(gene_sets)

  genes = intersect(genes, universe)
# 下面这行优化过了
  gene_sets = intersectToList(gene_sets, universe)

  n_universe =length(universe)
  n_genes =length(genes)

# 下面这行优化过了
  x = sapply(intersectToList(gene_sets, genes),length)
  m = sapply(gene_sets,length)
  n = n_universe - m
  k = n_genes

  p = phyper(x -1, m, n, k, lower.tail =FALSE)
names(p)= gs_names
  p
}

因为在Cpp代码中，基因集合名丢失了，因此在ora_v3()的最后三行，手动将基因集合名添加到变量p中。

参考

https://blog.csdn.net/Luciferchang/article/details/115684092

https://zhuanlan.zhihu.com/p/110932405

https://cloud.tencent.com/developer/article/2376798

https://mp.weixin.qq.com/s?__biz=Mzg4MjUzODM4Ng==&mid=2247484287&idx=1&sn=2c8cc00888c559403f7ecabc298b6ef9&scene=21#wechat_redirect

https://montilab.github.io/BS831/articles/docs/HyperEnrichment.html

看完本文应该能手搓一个ORA分析了

🔥 点赞、收藏、转发，一键三连助你基因富集分析起飞！🚀🧬

查看全文

http://www.kler.cn/a/584317.html