R语言|clusterprofile超几何分布富集分析 GO,KEGG富集分析，循环Fisher‘s test

• ID转换
• GO
• KEGG
• • 把KEGG里的geneid转回名称（readable table）
• 超几何分布以及生成data frame
• • 生成用来循环fisher的table
• 循环Fisher's test

这几天真的被这个超几何富集分析做的头大。原来单纯用找到的显著基因分析出来的pathway不一定是准确的，或者说不一定是具有显著性的。详细的原理推荐下面两个链接
https://www.jianshu.com/p/3d01a66e235b
https://blog.csdn.net/linkequa/article/details/88189582

这篇小记专注于使用clusterprofile包进行富集分析，再进行fisher’s 检验。btw，也吐槽一下大部分网站做出来的富集分析多多少少会有一些缺失或者滞后性，还是用最新的包跑一下比较严谨

先导入需要用的包

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) #可以转换的ID格式，如uniport，ENTERZID..

导入数据

#先导入全局蛋白list（234个蛋白）
olink<- read.csv('~/PD/string_protein_annotations.tsv', sep = "t")
symbo<- olink$X.node #（提取蛋白名称，其实不单独提取成一个向量也没关系）

表格格式无所谓，要的只是那一列蛋白的名字

ids<- bitr(symbo, #使用olink$X.node 也可以
 fromType = "SYMBOL", #输入的格式
 toType = c("GENENAME","ENTREZID"), #想得到的格式
           OrgDb = 'org.Hs.eg.db')

接下来导入显著蛋白的名称数据，方法同上

olink.sig<- read.csv('~/PD/PD分期/olink_sig_original.csv', check.names = F)
#这个数据格式里列名是每个显著蛋白，大家根据自己的df来提取名称
symbo.sig<- colnames(olink.sig[6:36]) 
ids.sig<- bitr(symbo.sig, fromType = "SYMBOL", toType = c("GENENAME","ENTREZID"),
           OrgDb = 'org.Hs.eg.db')

其实GO和KEGG都是差不多的过程

#全局蛋白GO
go.mf<- enrichGO(ids$ENTREZID, #使用转换后的ENTREZID
			   OrgDb = org.Hs.eg.db,
               keyType ='ENTREZID', ont = "MF",#MF/CC/BP
               pAdjustMethod = 'BH',
               pvalueCutoff = 0.01,
               qvalueCutoff = 0.05,
               readable = T #富集结果更好展示)
#显著蛋白GO
go.mf.sig<- enrichGO(ids.sig$ENTREZID, OrgDb = org.Hs.eg.db,
                 keyType ='ENTREZID', ont = "MF",
                 pAdjustMethod = 'BH',
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05,
                 readable = T)             
         
go.mf@result #查看GO富集的data frame

其实还可以点圈圈这里，也是富集结果的表格
富集分析的结果就是这个样子

KEGG同理，但我看网上很多人会先进行这一步代码，不然会报错，但我没有跑也顺利得出结果了，不知道有没有聚聚知道原因

data(ids, package = "DOSE") #报错的朋友跑这行代码先试一下
kk<- enrichKEGG(ids$ENTREZID, organism = "hsa",#hsa = homo sapine
                keyType = "kegg",
                pAdjustMethod = "BH",
                pvalueCutoff = 0.01,
)

kk.sig<- enrichKEGG(ids.sig$ENTREZID, organism = "hsa",
                keyType = "kegg",
                pAdjustMethod = "BH",
                pvalueCutoff = 0.01,
)

kk@result

KEGG富集结果

可以看到KEGG的和GO不一样的是geneIDs那一列。GO显示的是基因名称，而KEGG是基因ID，这时候不方便查看或者额外的分析。一行代码可以连带整张表将id转换为名称

y <- setReadable(kk, OrgDb = org.Hs.eg.db, keyType="ENTREZID")

因为我找的是我的全局蛋白（234个）和我的差异蛋白（31个）之间的通路差异显著性，所以有需要的朋友将全局蛋白换成你所需要的（比如GO目前有注释的蛋白是18352个

总而言之，这个分析是在检验差异蛋白的通路是否真的具有显著性。

如果事先了解过超几何分布（话说这个名字真的翻译的太吓人），应该知道这大概是个venn diagram、交集并集的关系

这里我用KEGG的通路来当例子，就不做GO的了。我先用下图红框的通路举个例子
这一行Generatio.x是全局基因的比例，输入169个gene，在这条通路上的有32个； Gene ratio.y是显著基因的比例，输入31个gene，在通路上的有8个。

（可能有人注意到我输入的是200多个gene怎么少了那么多，其实KEGG已有注释的gene并不多，和GO库比起来少了一半还多，再加上很多gene的通路研究还处于未知，所以缺失的这个比例也是可以理解的 ，个人的理解，我也去KEGG网站单独搜了没注释到的gene id，的确是没有pathway机制）

N<- 169 #总共的基因
M<- 32 #在通路上的基因
n<- 20 #显著的基因
k<- 6 #在通路上且显著的基因

a <- data.frame(gene.not.interest=c(M-k, N-M-n+k), gene.in.interest=c(k, n-k))
row.names(a) <- c("In_category", "not_in_category")
a

计算一下可以知道，这四个数字相加 = 169 （总共的基因），这样应该可以理解之前说的venn diagram的道理了

如果只有几条通路我们可以一个一个算，但很多通路的情况下怎么形成那么多table并且进行fisher或者卡方检验？
果然人类的进步原因之一是懒，当时一想到我要做那么多费舍人都晕掉了
方法来源于这位大神，但我跑起来有bug，进行了修改
https://zhuanlan.zhihu.com/p/98704640

还是用我的KEGG 结果

kkall<- kk@result #全局的富集结果表
kksig<- kk.sig@result #显著的富集结果表
kkresult<- merge(kkall, kksig, by = 'ID') #根据通路的ID合并两个表，取交集
kkresult<- kkresult[,-c(4:8,13:16)]
kkresult<- kkresult[,-c(5,7)] #删掉一些不用的列

这是我的kkresult table，需要用的其实只有count.x 和count.y两列

N<- 169
n<- 20
allcount<- kkresult$Count.x
sigcount<- kkresult$Count.y
#根据上面刚刚生成table的公式，获得每个通路超结合分布table的数字
A = allcount-sigcount
B = sigcount
C = 169-20-allcount+sigcount
D = 20 - sigcount

data<- data.frame(rbind(A,B,C,D)) #把每个向量rbind成一个大table
colnames(data)<- kkresult$ID#列名是通路的名字

这样就会形成一个 # of pathway x 4的表格，每一列是一个通路的超几何分布数据

整理好了表格格式，接下来就是循环检验了

qq<- list()#建一个空的list用来放结果
for (i in 1:ncol(data)){
  newdat=data[,i]
  test=matrix(newdat,nrow=2,ncol=2,byrow=TRUE)#可以单独跑这行代码看表data里ABCD的顺序对不对，顺序不会会影响检验结果
  qq[[i]]=fisher.test(test)   
}

f.pvalue<- data.frame(lapply(qq, 
	function(x) x$p.value)) #提取每个检验的p value
colnames(f.pvalue)<- colnames(data) #改列名

列数 = # pathway，
这样我们就完成了每条通路的超几何分布富集分析并进行差异检验了

其实不是很复杂，大约就四点：

1. id的转换
2. 富集分析，得到通路和基因数量
3. 生成用于Fisher‘s test的表
4. fisher.test ()

Done