ixxmu
commented
9 months ago ORA富集分析 by 生信菜鸟团
❝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()函数过于低效，所以这里我们使用顾叔提速富集分析的方法:https://mp.weixin.qq.com/s/G1gFGgu0jPzXWE97nxW4Qw
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")
#检查两个函数富集到结果是否一致
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进行可视化了。这部分我们之前和中文互联网上内容比较多，篇幅问题不再详述。
总结：
富集分析可以简单的理解为不放回的抽球问题，其概率符合超几何分布。 富集分析p值为X个基因富集到通路中及大于X个基因富集到通路中（更极端的情况）的概率之和，p.adjust是对多重假设检验的校正，目的是减少假阳性率。 R语言中intersect函数在进行一万多次富集分析时明显较慢，我们可以使用顾叔用c++写的相同功能的函数可以提升一点富集分析的效率。
Reference
概率: 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
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
https://mp.weixin.qq.com/s/iy6sfrWJ8_HN9pm-eapGow