数据整理
因为这个是我学生的真实数据，所以就不方便放出来给大家做示例文件了，可以只看代码，或者拿自己的数据练练手。
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
## <chr> <int>
## 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
搞起来啦。
