KEGG、GO富集分析与GSEA并不是二选一
测试开头
测试结尾
今天是生信星球陪你的第733天
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我!
最广为人知的富集分析做法是把上调、下调基因分别或者合并，拿来做GO和KEGG富集分析。经常有一些数据集，拿差异基因做得不到结果，那是因为确实富集不到任何通路，尤其KEGG，收录的基因数量只有8000左右，你富集不到其实是正常的，只是你不愿接受这个残酷的事实。不妨试试gsea，不是拿差异基因，而是拿全部基因作为输入哦。
1.R包和数据
输入数据仍然是差异分析结果表格，需要用到logFC值。示例数据在生信星球公众号回复“富集输入”即可拿到。
library(stringr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
library(enrichplot)
load("step4output.Rdata")
head(deg,3)
## logFC AveExpr t P.Value adj.P.Val B probe_id
## 1 5.780170 7.370282 82.94833 3.495205e-12 1.163798e-07 16.32898 8133876
## 2 -4.212683 9.106625 -68.40113 1.437468e-11 2.393169e-07 15.71739 7965335
## 3 5.633027 8.763220 57.61985 5.053466e-11 4.431880e-07 15.04752 7972259
## symbol change ENTREZID
## 1 CD36 up 948
## 2 DUSP6 down 1848
## 3 DCT up 1638
2.已有数据整理为要求的输入数据
data(geneList,package = "DOSE")
head(geneList)
## 4312 8318 10874 55143 55388 991
## 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
包的示例数据就是上面展示的geneList的样子，将我们的数据整理成这样:
-
向量
-
内容是logFC
-
名字是entrezid
-
从大到小排序
操作起来，就几句代码，但R语言基础知识要到位哦。
geneList = deg$logFC
names(geneList) = deg$ENTREZID
geneList = sort(geneList,decreasing = T)
head(geneList)
## 948 1638 158471 10610 6947 100133941
## 5.780170 5.633027 4.683610 3.875120 3.357670 3.322533
3.gsea做GO富集分析
gsea的官网数据集里有go，不过，也不用那么麻烦了，因为贴心的Y蜀黍写了函数，kegg也是一样有现成的函数：
ego <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "BP")
head(ego@result,3)
## ID Description
## GO:0048667 GO:0048667 cell morphogenesis involved in neuron differentiation
## GO:1901615 GO:1901615 organic hydroxy compound metabolic process
## GO:0002446 GO:0002446 neutrophil mediated immunity
## setSize enrichmentScore NES pvalue p.adjust qvalues
## GO:0048667 481 0.4025225 1.663780 0.001239157 0.03893605 0.03149766
## GO:1901615 478 0.3677970 1.517653 0.001248439 0.03893605 0.03149766
## GO:0002446 469 0.3420296 1.409709 0.001250000 0.03893605 0.03149766
## rank leading_edge
## GO:0048667 3323 tags=27%, list=18%, signal=23%
## GO:1901615 2003 tags=15%, list=11%, signal=14%
## GO:0002446 3315 tags=24%, list=18%, signal=20%
## core_enrichment
## GO:0048667 288/10512/4208/57556/51062/6696/22906/5803/6792/6092/2044/9839/2549/57698/23043/214/3908/5764/4747/8828/4035/3676/2263/8633/22854/6091/5295/56956/1952/3913/10090/1942/9820/5797/5530/6285/2049/2064/4147/2737/4133/1813/22902/3397/7869/1954/2674/4915/8609/3897/4131/1123/80208/55558/323/2045/2675/54664/1949/1749/26050/23114/1523/79600/4983/10013/4891/55816/59277/1020/1267/65981/220164/23303/9175/55619/23363/84628/347733/51057/2817/3792/120892/22871/5595/10752/659/6709/5663/55715/658/57498/10211/8239/57216/886/23334/6654/2534/8153/10458/22903/7155/7474/5728/79955/10006/4038/6324/5727/393/10500/7436/8660/627/26259/9499/5802/23032/5909/1382/56965/113246/3475/203228/1742/9798/815/2909/5578
## GO:1901615 1638/1545/6591/3293/2494/2246/84812/653/18/9839/8029/221/1734/27122/1593/7108/10404/6121/51477/8291/64788/51196/5950/5071/285440/6609/9227/5360/3158/5333/2903/7042/3635/32/151056/427/120227/1813/5208/403313/114876/66002/19/284161/6256/6799/2629/217/2539/3636/538/23236/6342/9677/3990/54884/7220/84532/8877/4157/114881/57016/1001/7157/1181/2034/9249/113026/6783/37/5092/10005/4435/27120
## GO:0002446 948/6947/57554/25840/10396/221/9961/718/5836/5787/3958/51646/51135/4758/51316/6286/2896/968/57153/8673/11240/175/25797/4126/427/55/10493/7077/5337/5660/4311/6868/2548/10043/138050/83716/1774/55904/4125/932/1992/823/1512/1509/23114/2760/11031/387921/23385/5869/53917/5269/374395/10312/158747/2588/210/5476/2720/2799/410/93100/178/22875/1832/2444/6402/57020/6814/10855/6709/5663/2512/7130/1773/966/29952/5265/5547/84418/728/230/9145/9545/201294/9342/3482/8649/6727/10159/3615/4048/567/1520/2934/8773/57126/3685/58485/2517/11010/5580/1778/1508/55860/25801/23593/9725/127829/2495/535/51071/9798
ekk <- gseKEGG(geneList = geneList,
verbose = T)
head(ekk@result,3)
## ID Description setSize enrichmentScore NES
## hsa04142 hsa04142 Lysosome 121 0.5584081 2.008820
## hsa05204 hsa05204 Chemical carcinogenesis 58 0.6233126 2.014398
## hsa00140 hsa00140 Steroid hormone biosynthesis 45 0.5790892 1.797418
## pvalue p.adjust qvalues rank leading_edge
## hsa04142 0.001582278 0.03560372 0.02863195 3101 tags=44%, list=17%, signal=37%
## hsa05204 0.001680672 0.03560372 0.02863195 1784 tags=24%, list=10%, signal=22%
## hsa00140 0.001686341 0.03560372 0.02863195 2119 tags=18%, list=11%, signal=16%
## core_enrichment
## hsa04142 1513/1519/26503/6609/4758/968/6272/8906/175/4126/427/8722/3073/5660/54/2629/2581/2548/1777/138050/1200/4125/8218/1512/1509/2760/84572/3425/4891/10312/2588/5476/2720/2799/9516/410/53/950/1203/79158/23457/537/3423/3482/23659/20/1520/9114/1497/2517/4669/1508/256471
## hsa05204 1545/7364/221/2952/4258/2941/4259/2052/1571/6799/373156/405/10/2944
## hsa00140 1545/3293/7364/7923/1571/6783/79154/1583
4.可视化
以前讲过GSEA：
有很漂亮的图，以及新版的单条通路、多条通路的展示图
gseaplot2(ego, geneSetID = 1, title = ego$Description[1])
gseaplot2(ego, geneSetID = 1:3)
gseaplot2(ekk, geneSetID = 1:3)
妈妈再也不用担心你的KEGG富集没有结果。
插个小广告！
答疑公告：生信星球答疑公告-2020年全年有效
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
