只有基因名，能做什么富集分析？

对于我这样的小白来说，芯片和转录组数据最常见的做法的取表达矩阵做差异分析获得logFC,标识change后取DEG再进行富集分析，但是随着学习的深入，还会有相关基因或其他目标基因需要进行富集分析，此时即使只有感兴趣的基因名，也是能做ORA富集分析的。

ORA分析原理

推荐生信技能树文章：

ORA富集分析 (qq.com)

ORA分析通俗来说就是将获得的基因与所选择的数据库中的通路进行比对，计算输入的基因在通路基因中的比例。

通常会得到GenRatio(属于该通路的输入基因个数M/在数据库中的输入基因个数N)和BgRatio(该通路的所有基因a/收集在数据库中的所有通路基因b)。

可以看作在有b个基因的数据库(含有a个x通路基因)的暗箱中进行N次不放回的抽取，获得了M个属于x通路的基因（超几何模型：不放回的抽取）。

然后通过假设检验的方式来确定来判断我们选择的基因和通路是否有关系，假设检验的步骤如下：

设立H0：输入的基因与某通路无关。

在H0成立的前提下，计算获得我们的抽取结果的概率p，若p<0.05(自定义的一类错误，即假阳性控制概率α，通常为0.1,0.05,0.01),则说明该结果为小概率事件。

正常情况（即H0成立时）我们获得小概率事件概率为p，由于p是那么小，和买彩票中五百万差不多大的概率。

我们就可以认为不是我们的运气太好获得了小概率事件，而是H0它就不成立，即从我们获得输入基因的源头开始，就指向了现在我们在数据库里找到的通路。

比如我们找了100个从年龄小到年龄大的实验者进行测序，得到某些变化具有一定规律，于是取其进行ORA分析，找到了一系列与衰老相关的通路，我们能在其中找到一定量的基因导致获得p值小于α，不是因为我们运气导致的，是由于这一系列基因本身就与衰老相关导致的。

不过还有假阳性呢，我们判断H0不成立的出错概率还是存在的，也许我们的运气就是这么好呢？

所以即使p值都很小，但找到了多个通路，相当于进行多重假设检验，那么其错误期望很可能大于1，即一次假设检验都不出错的概率非常小，那么我们通常选择进行p值校正，常用的包enrichGO中使用的是BH算法，我们通常更相信校正后的p值，毕竟其在一定程度上排除了我们由于随机性导致的判断失误。

GO/KEGG数据库ORA分析及可视化

此处使用的数据是相同样本在不同流程的上游分析后得到的差异基因，结果奇怪情有可原

添加ENTREZID列

library(clusterProfiler)
library(org.Mm.eg.db)
library(enrichplot)
library(ggplot2)
library(GOplot)
load("gene_list")
# 样本为小鼠
go_organism <- "org.Mm.eg.db" # org.Hs.eg.db / org.Mm.eg.db
kegg_organism <- 'mmu' # hsa / mmu
anno_organism <- 'mouse' # human / mouse 
​
# 转换基因名称格式
ids <- bitr(gene_list,fromType="SYMBOL", toType="ENTREZID", OrgDb = go_organism)
head(ids)
ids.unique <- ids %>% dplyr::distinct(ENTREZID, .keep_all = TRUE)
rich<-ids

GO数据库分析及可视化

#GO富集分析#
go_enrich = enrichGO(gene = rich[,2], #表示前景基因，即待富集的基因列表;[,1]表示对entrezid_all数据集的第1列进行处理
                     OrgDb = org.Mm.eg.db, 
                     keyType = "ENTREZID", #输入数据的类型
                     ont = "ALL", #可以输入CC\MF\BP\ALL
                     pAdjustMethod = "fdr", # 指定多重假设检验矫正的方法，选项包含 "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
                     pvalueCutoff = 0.3,  #指定 p 值阈值（指定 1 以输出全部），默认为0.05
                     qvalueCutoff = 0.3, #指定 q 值阈值（指定 1 以输出全部），默认0.2
                     readable = T) # 是否将gene ID转换到 gene symbol]
​
##可视化
##条形图
pdf(file="GO-barplot.pdf",width = 14,height = 5)
barplot(go_enrich, drop = TRUE, showCategory =10,label_format=100,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
​
##气泡图
pdf(file="GO-bubble.pdf",width = 14,height = 5)
dotplot(go_enrich,showCategory = 10,label_format=100,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()

KEGG数据库分析及可视化

#KEGG富集分析#
KEGG_enrich = enrichKEGG(gene = rich[,2], #即待富集的基因列表
                         keyType = "kegg",
                         organism= "mouse",  #hsa，可根据你自己要研究的物种更改，可在https://www.kegg.jp/brite/br08611中寻找
                         pAdjustMethod = "fdr", # 指定多重假设检验矫正的方法，选项包含 "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
                         pvalueCutoff = 0.3,  #指定 p 值阈值（指定 1 以输出全部），默认为0.05
                         qvalueCutoff =0.3)  #指定 q 值阈值（指定 1 以输出全部），默认0.2
​
save(ids,go_enrich,KEGG_enrich,file="GO.Rdata")
​
##可视化
##条形图
pdf(file="KEGG-barplot.pdf",width = 14,height = 5)
barplot(KEGG_enrich, drop = TRUE, showCategory = 15,label_format=100)
dev.off()
​
##气泡图
pdf(file="KEGG-bubble.pdf",width = 14,height = 5)
dotplot(KEGG_enrich, showCategory = 15,label_format=100)
dev.off()

可以看到一些奇怪的结果，p.adjust相同是由于BH算法可能产生连续相同的p值

GO:

KEGG: