❝Over Representation Analysis❞

超几何分布

这部分内容摘自百度百科。超几何分布是统计学上一种离散概率分布。它描述了从有限N个物件（其中包含M个指定种类的物件）中抽出n个物件，成功抽出该指定种类的物件的次数（不放回）。超几何分布中的参数是N,n,M，上述超几何分布记作X~H(N,n,M)

产品抽样检查中经常遇到一类实际问题，假定在N件产品中有M件不合格品，即不合格率：

在产品中随机抽n件做检查，发现k件不合格品的概率[1]为：

k的取值范围：（k=t，t+1，…，s）

其中s是M与n中的较小者，t在n不大于合格品数（即n≤N-M）时取0，否则t取n减合格品数之差（即t=n-(N-M)） 亦可写作：

（与上式不同的是M可为任意实数，而C表示的组合数M为非负整数）

为古典概型的组合形式，a为下限，b为上限，此时我们称随机变量[2]X服从超几何分布（hypergeometric distribution）

需要注意的是：

• （1）超几何分布的模型是不放回抽样。
• （2）超几何分布中的参数是N,n,M，上述超几何分布记作X~H(N,n,M)。

ORA过表达分析

富集分析的算法有很多，最常用是Over Representation Analysis，ORA过表达分析，其次是gene set enrichment analysis, GSEA基因集富集分析。

具体区别详见这篇知乎文章：https://zhuanlan.zhihu.com/p/534016487[3]

本文讲更简单易懂的ORA过表达分析

过表达分析(ORA)的输入需要数据：

• 背景基因集：一个物种的某个数据库，包含多个子集。如GO数据库中,有三个ontology（本体），分别描述基因的细胞组分（cellular component，CC）、分子功能（molecular function，MF）、参与的生物过程（biological process，BP）。在R包GO.db，version=3.18中，BP、MF、CC分别包含了15709、1977、5055个通路。
• 特定基因集：常常是差异表达基因集(differentially expressed genes, DEGs)，也可以使用WGCNA识别到的模块中的基因等其他来源

ORA过表达分析实质上就是Fisher's 精确检验：

我们可以形象的理解富集分析：

例如GO数据库的1w8k个基因中，1000个基因在一个通路中，17000个基因在这个通路外，我们有2000个差异表达基因，其中在通路中为600个，不在通路外有1400个。就可以理解为从1w8k个球的桶中不放回的抽球，桶中有1000个黑球和17000个白球，我们进行2000次不放回的抽取，抽到了600个黑球和1400个白球。

富集分析的P值：即计算出现抽中600个球（通路中富集到600个基因）及以上（更极端情况）概率之和。并与显著性系数0.05/0.01/0.001进行比较。

推荐看这篇文章：https://zhuanlan.zhihu.com/p/110932405[4]

应用如下公式计算抽到k个球的概率：

其中n为差异表达基因数，k为差异表达基因中富集到通路（子集）中的基因（就是在通路中的差异表达基因），M是通路（子集）中基因数量，N为数据库中所有基因的基因数（在GO的BP，MF，CC中各约为1w8k个）。

上式计算得到的是p为k个基因富集到通路中的概率，在富集分析中，我们要对k，k+1，k+2.......n或M （当差异表达基因数＜通路中基因数时，取n为最大，代表所有差异表达基因都富集在该通路下；

当差异表达基因数＞通路中基因数时，取m为最大，代表该通路中全部基因都为差异表达基因）的概率进行求和。即代表有k及以上个差异表达基因富集到通路中的概率（极端情况），若这个概率p＜0.05/0.01。则认为其为小概率事件，在一次假设检验中不可能发生，拒绝原假设，接受备择假设，也就差异表达基因与该通路有比较强的联系。

代码实现

这里我们以GO数据库为例

因为我是从写好的整个函数里摘出对应的代码的，所以先需要给一些参数赋值，同时我们取500个感兴趣（这里只是示例）的基因作为输入

#注意在调试代码的时候不要与以下变量重名
min_gene=10
max_gene=500
org="org.Hs.eg.db"
method="BH"
pvalue_cutoff=0.05
padjust_cutoff=0.05

#这里我们随便取一些基因作为差异表达基因，可以使用自己找的
#我写的代码需要entrezid作为数据输入，所以

library(dplyr)
library(AnnotationDbi)
library(org,character.only = T)
#随便取前500个id赋给向量gene
AnnotationDbi::keys(org.Hs.eg.db,keytype = "ENTREZID") %>% head(.,2000)->gene

现在向量gene里有500个基因的entrezid

因为BP，MF，CC中有注释的总的基因数不同，我们要对这三个ontology分别进行过表达分析，所以需要计算他们下面的所包含基因的数量

#gs是该物种org.db中包含的通路，就是这个物种有的通路
gs<-as.list(get(paste0(strsplit(org,".db")[[1]],"GO2ALLEGS")))
#去重复
gs = lapply(gs, unique)

gs是一个列表，其下面每一个元素包含一个向量，向量内容为该通路下的所有基因

#载入GO.db包
library(GO.db)
as.list(GOTERM)->GO

GO也是一个列表，其每个元素下是一个S4类对象，包含信息有：

• GOID：该通路在GO数据库中的编号
• Term：通路名称
• Ontology：属于哪个“本体”，就是是BP，MF还是CC
• Defination：详细描述
• Synonym：同义词
• Secondary：Alternate term identifier就是备用的通路识别号，就是搜这个编号也能出来这个通路

其实这个列表就等于GO数据库搜索ontology的网页界面

  #获取该物种有的GO_accession
  GO[names(gs)]->GO
  #提取GO列表中的每一个元素下的GOID,Term,Ontology,Definition
  #这里是写了一个函数，用lapply循环列表下的每一个元素，使用@提取其信息，以向量的形式返回赋给go_anno
  go_anno<-lapply(GO,function(element){
    d<-c(element@GOID,element@Term,element@Ontology,element@Definition)
    return(d)
  })
  t(as.data.frame(go_anno))->go_anno
  colnames(go_anno)<-c("GOID","Term","Ontology","Definition")
  gsub("\\.", ":", rownames(go_anno))->rownames(go_anno)
  #这里有个小坑，GO数据库的编号都是GO:XXXXXXXX,而go_anno中的行名为GO.XXXXXXXX，用gsub函数替换一下

go_anno中的信息将用于我们最后输出的表格，毕竟我们也需要知道自己富集到了什么东西（生物学意义）

上文说到BP，CC，MF中包含基因数量不同。所以我们富集分析时所有基因的数量N应该根据这三个ontology发生改变：

  #获取BP下包含的全部基因
  #从go_anno信息中取出第三列为BP的GOid
  rownames(go_anno[go_anno[,3]=="BP",])->GO_BP
  #然后用unlist将gs列表中所有对应BP的元素（内容前文说到是包含基因entrezid的向量）从头到尾连成一个向量
  BP_gene<-unlist(gs[GO_BP],use.names = F)
  #然后进行去重，就得到ontology BP中的所有基因了
  unique(BP_gene)->BP_gene
  
  #如法炮制，获取CC下包含的全部基因
  rownames(go_anno[go_anno[,3]=="CC",])->GO_CC
  CC_gene<-unlist(gs[GO_CC],use.names = F)
  unique(CC_gene)->CC_gene
  #如法炮制，获取MF下包含的全部基因
  rownames(go_anno[go_anno[,3]=="MF",])->GO_MF
  MF_gene<-unlist(gs[GO_MF],use.names = F)
  unique(MF_gene)->MF_gene

富集分析的限速步骤是判断差异表达基因是否在通路中，由上文知道BP+MF+CC总共有2w多通路，R的intersect()函数过于低效，所以这里我们使用顾叔提速富集分析的方法。

  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中每一个向量中重复的元素。

虽然我们不会c++，但是我们只需要会用顾叔的代码就好了，这就是一名生（代）信（码）工（裁）程（缝）师应有的修养。

以下是用于将entrezid转换成gene symbol，作为输出使结果可读性更高的代码

  #获取ENTREZID和symbol对应关系
  k<-keys(org.Hs.eg.db,keytype = "ENTREZID")
  e2s<-AnnotationDbi::select(org.Hs.eg.db,keys=k,columns="SYMBOL",keytype = "ENTREZID")
  symbol<-e2s[,2];names(symbol)<-e2s[,1]
  
  #定义一个函数，其输入为entrezid的向量，转化成symbol的向量
  symboltrans<-function(element,symbol) {
    element<-paste(symbol[element],collapse = "/")
    return(element)}

接下来就是计算富集分析的p value，adjust p value，generatio，bgratio，fold enrich，enrich factor了：

• GeneRatio：差异表达基因富集到目的通路基因数占差异表达基因数的比例
• BgRatio：目标通路基因数占通路集总基因数的比例
• pvalue：假设检验结果的显著性
• p.adj：多重假设检验校正后p值，默认使用“BH”法
• Fold_Enrichment：GeneRatio / BgRatioRich_Factor：富集到通路中的基因数/通路中基因数（即GeneRatio的分子除以BgRatio的分母）

我们定义一个函数ora_v3（函数名是沿用顾叔推文的取名），然后写一个循环调用函数进行富集分析

  ora_v3 = function(genes, gene_sets, universe ,symbol ,min_gene , max_gene ,
                    method ,pvalue_cutoff ,padjust_cutoff) {
    
    #m为基因集中每个通路包含的基因数量
    m = sapply(gene_sets, length)
    m >= min_gene & m <= max_gene ->keep
    m[keep]->m
    gene_sets[keep]->gene_sets
    
    #基因集的名字
    gs_names = names(gene_sets)
    #取基因集中在MF/BP/CC中的基因
    genes = intersect(genes, universe)
    #函数intersectToList()可以在R中使用。它接受两个参数。
    #第一个是一个包含了若干向量的列表（lt），第二个参数是一个向量（x）。
    #其中x会和lt中的每一个向量进行intersection。在Cpp代码中，顾叔也去掉了lt中每一个向量中重复的元素
    n_universe = length(universe)
    n_genes = length(genes)
    
    #对基因集列表每一个子列表进行循环，判断基因是否在通路中，返回在通路中的差异表达基因
    genes_in_set<-intersectToList(gene_sets, genes)
    #x为每个通路中富集到基因的数量
    x = sapply(genes_in_set, length)
    
    x!=0 ->keep
    x[keep]->x
    m[keep]->m
    gene_sets[keep]->gene_sets
    gs_names[keep]->gs_names
    genes_in_set[keep]->genes_in_set
    
    #创建一个向量，其每一个元素对应富集到的通路中的所有基因的gene symbol，使用/作为分隔
    genesymbol<-c()
    for(i in 1:length(genes_in_set)){
      genesymbol[i]<-paste(symboltrans(genes_in_set[[i]],symbol),collapse = "/")
    }
    #GeneRatio:x/k,BgRatio:m/n_universe
    #包含在ontology中但不在m对应通路中的基因数量
    n = n_universe - m
    #差异表达基因数量
    k = n_genes
    
    GeneRatio = paste(x,k,sep="/")
    BgRatio = paste(m,n_universe,sep="/")
    Fold_Enrichment=(x*n_universe)/(k*m)
    Rich_Factor=x/k
    
    #可以尝试下不用phyer函数，自己写，phyer函数就是计算：
    #通路中富集到基因的数量大于等于x的极端情况的概率，即富集分析的p值
    p = phyper(x - 1, m, n, k, lower.tail = FALSE)
    #method默认是BH法
    p.adjust=p.adjust(p,method = method)
    
    #将结果输出为数据框
    output<-data.frame(GeneRatio=GeneRatio,BgRatio=BgRatio,
                       gene_symbol=genesymbol,count=x,pvalue=p,p.adjust=p.adjust,
                       Fold_Enrichment=Fold_Enrichment,Rich_Factor=Rich_Factor)
    rownames(output)<-gs_names
    
    #这个判断是使用p value还是p adjust进行结果筛选，默认为p value<0.05
    if(is.null(padjust_cutoff)){
      output[output$pvalue <= pvalue_cutoff,]->output
    }
    else{
      output[output$p.adjust <= padjust_cutoff,]->output
    }
    
    return(output)
  }

#这个for循环是对BP,MF,CC三个ontology分别进行过表达分析，结果赋给BP_Res,MF_Res,CC_Res
    for (k in c("BP","MF","CC")) {
      assign(paste0(k,"_Res"),ora_v3(gene,gs[get(paste0("GO_",k))],
                                    get(paste0(k,"_gene")),
                                    symbol=symbol ,min_gene=min_gene ,max_gene=max_gene, method=method,
                                    pvalue_cutoff=pvalue_cutoff,padjust_cutoff=padjust_cutoff))
}

最后合并结果输出为一个表格

  rbind(BP_Res,MF_Res,CC_Res)->Res
  cbind(go_anno[match(rownames(Res),rownames(go_anno)),],Res)->Res
 write.csv(Res,file = "ORA_res.csv")

我将上文所有内容封装成了一个函数，最简单的用法就是输入一个是基因entrezid的向量，返回结果就是过表达分析的表格，函数代码如下：

• 参数gene就是输入的entrezid向量，
• min_gene是默认过滤掉富集到基因数量＜10的通路，
• org就是研究物种的org.db包，注意，只有在bioconductor上有才能做过表达分析，否则无法使用我的代码，
• method是p值校正方法，默认BH法（其他可用方法详见stats包的p.adjust函数），
• pvalue_cutoff是过滤掉p值＜0.05的结果，
• padjust_cutoff是过滤掉p值校正值＜0.05的结果。（默认使用pvalue_cutoff过滤，当设置padjust_cutoff时，优先使用padjust_cutoff进行过滤）

#GO数据库ORA富集分析函数 enrichment function written by Xiao Chen

#输入为ENTREZID，通路中默认最少要包含10个基因，物种默认为人类
#p值校正方法BH法，结果显示阈值：p值小于0.05，padjust结果显示阈值为空
ORA_GO<-function(gene,min_gene=10,max_gene=500,org="org.Hs.eg.db",method="BH",
            pvalue_cutoff=0.05,padjust_cutoff=NULL){
  library(AnnotationDbi)
  library(org,character.only = T)
  
  #gs是该物种org.db中包含的通路，就是这个物种有的通路
  gs<-as.list(get(paste0(strsplit(org,".db")[[1]],"GO2ALLEGS")))
  #去重复
  gs = lapply(gs, unique)
  
  #方法2.选择BP,MF,CC各自包含的基因分别作为背景基因
  library(GO.db)
  as.list(GOTERM)->GO
  
  #取该物种有的GO_accession
  GO[names(gs)]->GO
  #提取GO列表中的每一个元素下的GOID,Term,Ontology,Definition
  #这里是写了一个函数，用lapply循环列表下的每一个元素，使用@提取其信息，以向量的形式返回赋给go_anno
  go_anno<-lapply(GO,function(element){
    d<-c(element@GOID,element@Term,element@Ontology,element@Definition)
    return(d)
  })
  t(as.data.frame(go_anno))->go_anno
  colnames(go_anno)<-c("GOID","Term","Ontology","Definition")
  #这里有个小坑，GO数据库的编号都是GO:XXXXXXXX,而go_anno中的行名为GO.XXXXXXXX，用gsub函数替换一下
  gsub("\\.", ":", rownames(go_anno))->rownames(go_anno)
  
  #获取BP下包含的全部基因
  #从go_anno信息中取出第三列为BP的GOid
  rownames(go_anno[go_anno[,3]=="BP",])->GO_BP
  #然后用unlist将gs列表中所有对应BP的元素（内容前文说到是包含基因entrezid的向量）从头到尾连成一个向量
  BP_gene<-unlist(gs[GO_BP],use.names = F)
  #然后进行去重，就得到ontology BP中的所有基因了
  unique(BP_gene)->BP_gene
  
  #如法炮制，获取CC下包含的全部基因
  rownames(go_anno[go_anno[,3]=="CC",])->GO_CC
  CC_gene<-unlist(gs[GO_CC],use.names = F)
  unique(CC_gene)->CC_gene
  #如法炮制，获取MF下包含的全部基因
  rownames(go_anno[go_anno[,3]=="MF",])->GO_MF
  MF_gene<-unlist(gs[GO_MF],use.names = F)
  unique(MF_gene)->MF_gene
  
  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;
}
')
  
  #获取ENTREZID和symbol对应关系
  k<-keys(org.Hs.eg.db,keytype = "ENTREZID")
  e2s<-AnnotationDbi::select(org.Hs.eg.db,keys=k,columns="SYMBOL",keytype = "ENTREZID")
  symbol<-e2s[,2];names(symbol)<-e2s[,1]
  
  #定义一个函数，其输入为entrezid的向量，转化成symbol的向量
  symboltrans<-function(element,symbol) {
    element<-paste(symbol[element],collapse = "/")
    return(element)}
  
  ora_v3 = function(genes, gene_sets, universe ,symbol ,min_gene , max_gene ,
                    method ,pvalue_cutoff ,padjust_cutoff) {
    
    #m为基因集中每个通路包含的基因数量
    m = sapply(gene_sets, length)
    m >= min_gene & m <= max_gene ->keep
    m[keep]->m
    gene_sets[keep]->gene_sets
    
    #基因集的名字
    gs_names = names(gene_sets)
    #取差异表达基因中在MF/BP/CC中的基因
    genes = intersect(genes, universe)
    #函数intersectToList()可以在R中使用。它接受两个参数。
    #第一个是一个包含了若干向量的列表（lt），第二个参数是一个向量（x）。
    #其中x会和lt中的每一个向量进行intersection。在Cpp代码中，顾叔也去掉了lt中每一个向量中重复的元素
    n_universe = length(universe)
    n_genes = length(genes)
    
    #对基因集列表每一个子列表进行循环，判断基因是否在通路中，返回在通路中的差异表达基因
    genes_in_set<-intersectToList(gene_sets, genes)
    #x为每个通路中富集到基因的数量
    x = sapply(genes_in_set, length)
    
    x!=0 ->keep
    x[keep]->x
    m[keep]->m
    gene_sets[keep]->gene_sets
    gs_names[keep]->gs_names
    genes_in_set[keep]->genes_in_set
    
    #创建一个向量，其每一个元素对应富集到的通路中的所有基因的gene symbol，使用/作为分隔
    genesymbol<-c()
    for(i in 1:length(genes_in_set)){
      genesymbol[i]<-paste(symboltrans(genes_in_set[[i]],symbol),collapse = "/")
    }
    #GeneRatio:x/k,BgRatio:m/n_universe
    #包含在ontology中但不在m对应通路中的基因数量
    n = n_universe - m
    #差异表达基因数量
    k = n_genes
    
    GeneRatio = paste(x,k,sep="/")
    BgRatio = paste(m,n_universe,sep="/")
    Fold_Enrichment=(x*n_universe)/(k*m)
    Rich_Factor=x/k
    
    #可以尝试下不用phyer函数，自己写
    p = phyper(x - 1, m, n, k, lower.tail = FALSE)
    p.adjust=p.adjust(p,method = method)
    
    #将结果输出为数据框
    output<-data.frame(GeneRatio=GeneRatio,BgRatio=BgRatio,
                       gene_symbol=genesymbol,count=x,pvalue=p,p.adjust=p.adjust,
                       Fold_Enrichment=Fold_Enrichment,Rich_Factor=Rich_Factor)
    rownames(output)<-gs_names
    
    
    #这个判断是使用p value还是p adjust进行结果筛选，默认为p value<0.05
    
    if(is.null(padjust_cutoff)){
      output[output$pvalue <= pvalue_cutoff,]->output
    }
    else{
      output[output$p.adjust <= padjust_cutoff,]->output
    }
    
    
    
    return(output)
  }
  #这个for循环是对BP,MF,CC三个ontology分别进行过表达分析，结果赋给BP_Res,MF_Res,CC_Res
    for (k in c("BP","MF","CC")) {
      assign(paste0(k,"_Res"),ora_v3(gene,gs[get(paste0("GO_",k))],
                                    get(paste0(k,"_gene")),
                                    symbol=symbol ,min_gene=min_gene ,
                                    max_gene=max_gene,method=method,
                                    pvalue_cutoff=pvalue_cutoff,padjust_cutoff=padjust_cutoff))
}
  
  rbind(BP_Res,MF_Res,CC_Res)->Res
  cbind(go_anno[match(rownames(Res),rownames(go_anno)),],Res)->Res
  
  return(Res)
}

现在调用我们的函数，查看富集分析结果，并使用system.time()记录运行时间

#除了gene参数必须输入，其他参数都是可选或可缺省的
assign("Res1",ORA_GO(gene,padjust_cutoff=0.05)) %>% system.time()
write.csv(Res1,file = "ORA_res1.csv")

y叔的clusterProfiler表现

library(clusterProfiler)
#因为我们的代码没有计算qvalue且没有设置maxGSSize，为了方便比对结果，设为最大不影响输出
assign("Res2",enrichGO(gene,
                       OrgDb = org.Hs.eg.db,
                       ont = "ALL",
                       qvalueCutoff = 1,
                       maxGSSize = 99999,
                       readable = T)) %>% system.time()
write.csv(Res2,file = "ORA_res2.csv")

image.png

#检查两个函数富集到结果是否一致
Res1[order(Res1$p.adjust,Res1$GOID,decreasing = F),]->Res1
Res2[order(Res2$p.adjust,Res2$ID,decreasing = F),]->Res2
identical(Res1$GOID,Res2$ID)
#返回TRUE则一致

hh，之前还寻思着能不能不能做的比y叔快，结果越写越慢。相较于clusterprofiler我们的结果没有计算qvalue，但是提供了Fold_enrichment和Rich_Factor。

ORA_res1.csv的内容如下：

有了这样的富集分析结果，就可以使用ggplot2进行可视化了。这部分我们之前和中文互联网上内容比较多，篇幅问题不再详述。

总结：

1. 富集分析可以简单的理解为不放回的抽球问题，其概率符合超几何分布。
2. 富集分析p值为X个基因富集到通路中及大于X个基因富集到通路中（更极端的情况）的概率之和，p.adjust是对多重假设检验的校正，目的是减少假阳性率。
3. R语言中intersect函数在进行一万多次富集分析时明显较慢，我们可以使用顾叔用c++写的相同功能的函数可以提升一点富集分析的效率。

Reference

[1] 概率: https://baike.baidu.com/item/%E6%A6%82%E7%8E%87/828845?fromModule=lemma_inlink

[2] 随机变量: https://baike.baidu.com/item/%E9%9A%8F%E6%9C%BA%E5%8F%98%E9%87%8F/828980?fromModule=lemma_inlink

[3] https://zhuanlan.zhihu.com/p/534016487: https://zhuanlan.zhihu.com/p/534016487

[4] https://zhuanlan.zhihu.com/p/110932405: https://zhuanlan.zhihu.com/p/110932405