拯救你那些富集不到结果的基因

数据整理

因为这个是我学生的真实数据，所以就不方便放出来给大家做示例文件了，可以只看代码，或者拿自己的数据练练手。

library(org.Hs.eg.db) 
library(clusterProfiler)
library(dplyr) 
library(ggplot2)
deg=read.csv("deg.csv",row.names = 1)
head(deg)
##           ensembl      logFC     logCPM      PValue        FDR change ENTREZID
## 1 ENSG00000000003  0.1655532  5.0415800 0.467563867 0.79040699    not     7105
## 2 ENSG00000000419  0.3873324  5.8257517 0.074773161 0.30181701    not     8813
## 3 ENSG00000000457 -0.3431328  0.8672989 0.517593456 0.80314607    not    57147
## 4 ENSG00000000460  0.5544268  2.1875783 0.107585303 0.37378191    not    55732
## 5 ENSG00000000971  0.3512430 -0.2066577 0.679314330 0.91718616    not     3075
## 6 ENSG00000001036  0.6269545  5.5618964 0.004645405 0.05314606    not     2519
gene_diff = deg$ENTREZID[deg$change != 'not']
length(gene_diff)
## [1] 975

975个基因做不出GO富集，我的第一反应是孩子学艺不精，然后看了他的代码，排除了可能的代码错误和物种错误，我发现他做的富集分析是没有错的！

ego <- enrichGO(gene = gene_diff,
                OrgDb= org.Hs.eg.db,
                ont = "ALL",
                readable = TRUE)
dim(ego)
## [1] 0 9

如果你是kegg富集不到，那我就会发这个链接(简书放不了链接，感兴趣的话去生信星球公号看吧)

但是这可是975个基因,GO的富集分析居然做不出来，实属诡异啊。让我有了想拯救一下的心思。

拯救计划1

首先是不卡p值，pvalueCutoff = 1。就算不显著，你也得让我看见结果，端详一下嘛。

ego1 <- enrichGO(gene = gene_diff,
                OrgDb= org.Hs.eg.db,
                ont = "ALL",
                pvalueCutoff = 1,
                readable = TRUE)
dim(ego1)
## [1] 73 10
table(ego1@result$pvalue<0.05)
## 
## TRUE 
##   73
table(ego1@result$p.adjust<0.05)
## 
## FALSE 
##    73
par(mfrow = c(1,2))
plot(sort(ego1@result$pvalue))
plot(sort(ego1@result$p.adjust))

adjust之后，就全都不显著了，而clusterProfiler默认是按照p.adjust来做各项分析的。

一个可选的操作，用p值筛选条目，然后自己画后面的条带图和气泡图交差，写文章的时候要写明白你用的是原始p值！会不会被怼就要看运气啦。

拯救计划2

经神奇豆豆指点，发现了两个enrichGO的参数，minGSSize = 10,maxGSSize = 500,这是限制每个term里面基因数量的参数，默认只计算那些10-500个基因的term，就避免了父节点子节点一起出现在你的结果里。但目前我们遇到的是子节点富集不到，那调大一点也是可以的。只不过要结合生物学意义，看这些大的term能不能和自己的实验设计联系起来。

ego2 <- enrichGO(gene = gene_diff,
                 OrgDb= org.Hs.eg.db,
                 ont = "ALL",
                 pvalueCutoff = 1,
                 maxGSSize = 5000,
                 minGSSize = 5,
                 readable = TRUE)
dim(ego2)
## [1] 230  10
table(ego2@result$p.adjust<0.05)
## 
## FALSE  TRUE 
##   210    20

跌跌撞撞有了结果，可以挑一挑啦啊。

拯救计划3

既然超几何分布检验不太行，那不如用GSEA来试试咯。 GO数据库的GSEA分析有两个做法，一个是简单的专用函数gseGO：

#rm(list = ls())
library(GSVA)
library(GSEABase)
library(msigdbr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(limma)
#install.packages("tinyarray")
library(tinyarray)
ge = deg$logFC
names(ge) = deg$ENTREZID
ge = sort(ge,decreasing = T)
head(ge)
##     6751    79585   400566    57615     6638    57540 
## 7.931213 7.096336 6.659524 6.659524 6.553707 6.553707
em <- gseGO(geneList     = ge,
              OrgDb        = org.Hs.eg.db,
              ont          = "ALL")
dim(em)
## [1]  7 12

一个是使用msigdbr里的数据，虽然麻烦一些，但是支持更多的基因集。

GO_df = msigdbr(species = "Homo sapiens",category = "C5") %>% 
  dplyr::select(entrez_gene,gs_exact_source,gs_subcat)
GO_df = GO_df[GO_df$gs_subcat!="HPO",]
table(GO_df$gs_subcat)
## 
##  GO:BP  GO:CC  GO:MF 
## 722231 112249 116388
GO_df = GO_df[,c(2,1)]
head(GO_df)
## # A tibble: 6 x 2
##   gs_exact_source entrez_gene
##                    
## 1 GO:0009256            10840
## 2 GO:0009256           160428
## 3 GO:0009256             4522
## 4 GO:0009256            25902
## 5 GO:0009256           441024
## 6 GO:0006103            51166
em2 <- GSEA(ge, TERM2GENE = GO_df)
dim(em2)
## [1]  5 11

这下没有缺憾啦，岁月静好，可以交差。

gseaplot2(em, geneSetID = 1, title = em$Description[1])

拯救计划4

神奇豆豆最近埋头开发他的R包，他让我用他的包试试，因为支持的geneid比orgDb多一些。

#install.packages("genekitr")
library(genekitr)
gene_diff = as.character(gene_diff)
ego3 <- genGO(gene_diff,
  org = "human",
  ont = "all"
)
## Error in genGO(gene_diff, org = "human", ont = "all"): No GO terms enriched ...
dim(ego3)
## Error in eval(expr, envir, enclos): object 'ego3' not found
# 也没富集到，但是它也有那两个参数
ego4 <- genGO(gene_diff,
  org = "human",
  ont = "all",
  maxGSSize = 5000,
  minGSSize = 5
)
dim(ego4)
## [1] 20 13

搞起来啦。