1、安裝，加載所用到到R包
用BiocManager安裝，可同時加載依賴包
source("https://bioconductor.org/biocLite.R")
BiocManager::install("clusterProfiler")
library(clusterProfiler) ##富集分析
library(topGO) ###畫GO圖
library(AnnotationHub) ##獲取數據庫
library(BiocFileCache) ##依賴包
library(dbplyr) ##依賴包
library(pathview) ##看KEGG pathway
2、利用annotataionHub去抓取目標orgDb
ah <- AnnotationHub() ##收索所有orgdb，到ah
unique(ah$dataprovider) ##可查看數據注釋來源
query(ah, "Apis cerana") ##查找目標物種
tar_org <- ah[["AH62635"]] ##下載目標物種到org數據
3、了解org數據庫
主要有5個函數
columns(x)： 顯示當前對象有哪些數據 keytypes(x)： 有哪些keytypes可以用作select或keys的keytypes參數 keys(x, keytype, ...)：返回當前數據對象的keys select(x, keys, columns, keytype, ...)：基於keys, columns和keytype以data.frame數據類型返回數據，可以是一對多的關系 mapIds(x, keys, column, keytype, ..., multiVals)： 類似於select，只不過就返回一個列。3.1 以symbol形式展示
head(keys(tar_org,keytype = "SYMBOL"),30) ##默認為ENTREZID
3.2、可以查看gene類型
keytypes(tar_org) [1] "ARACYC" "ARACYCENZYME" "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" [8] "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PMID" "REFSEQ" [15] "SYMBOL" "TAIR" 3.3、select則是根據你提供的key值去查找注釋數據庫，返回你需要的columns信息
> select(tar_org, keys= "AGO1", columns=c("TAIR","GO"),keytype = "SYMBOL") 'select()' returned 1:many mapping between keys and columns SYMBOL TAIR GO EVIDENCE ONTOLOGY 1 AGO1 AT1G48410 GO:0004521 IDA MF 2 AGO1 AT1G48410 GO:0005515 IPI MF 3 AGO1 AT1G48410 GO:0005634 IDA CC 4 AGO1 AT1G48410 GO:0005737 IDA CC 5 AGO1 AT1G48410 GO:0005737 ISM CC 6 AGO1 AT1G48410 GO:0005737 TAS CC 7 AGO1 AT1G48410 GO:0005829 IDA CC 8 AGO1 AT1G48410 GO:0006306 RCA BP 9 AGO1 AT1G48410 GO:0006342 RCA BP 10 AGO1 AT1G48410 GO:0006346 RCA BP⚠️找到差異基因后，必須得確定你得基因號是對應ENTREZID 或者SYMBOL，，是屬於哪種類型，若沒有符合上述類型，可自行找到NCBI上得數據名稱，進行blast更換名字。
4、進行作圖
統一將差異基因名字改為ENTREZID，防止在做GO分析的時候出現報錯，需要將symbolID轉換成ENTREZID：用mapIds函數就可以轉換ID。
DEG.entrez_id = mapIds(x = tar_org, #### 數據庫 keys = DEG.gene_symbol, #####差異基因名字 keytype = "SYMBOL", ####差異基因類型是SYMBOL column = "ENTREZID") #####轉換為ENTREZID 這時就已經把symbolID轉換成ENTREZID了，但會出現個別的轉換不成功的情況，就是圖中NA的地方，我們進行以下操作即可去掉：
DEG.entrez_id = na.omit(DEG.entrez_id)
4.1、GO分析代碼
BP（Biological process）層面上的富集分析：
erich.go.BP = enrichGO(gene = DEG.entrez_id, ###差異基因ID OrgDb = tar_org, ###數據庫 keyType = "ENTREZID", ##基因ID類型 ont = "BP", ##對BP進行GO分析 pvalueCutoff = 0.5, ###fisher檢驗對p值 qvalueCutoff = 0.5) ###對p值進行校對對q值，一般大於p值
作圖：
dotplot(erich.go.BP)解讀BP層面富集分析圖：
橫坐標是GeneRatio，意思是說輸入進去的基因，它每個term（縱坐標）站整體基因的百分之多少。圓圈的大小代表基因的多少，圖中給出了最大的圓圈代表60個基因，圓圈的顏色代表P-value，也就是說P-value越小gene count圈越大，這事就越可信。
也可以畫柱形圖
barplot(erich.go.CC)
一般GO分析畫這兩個圖就可以了，有時也把GO分析畫成樹形圖，可以更加幫助我們理解。
plotGOgraph(erich.go.BP)樹狀圖很大，所以我們用代碼把它存成pdf，學習下如何用代碼
pdf(file="./enrich.go.bp.tree.pdf",width = 10,height = 15) plotGOgraph(erich.go.BP) dev.off()
至此，GO分析就做完了 ----> over
4.2、KEGG pathway 介紹
KEGG pathway是最常用的功能注釋數據庫之一，可以利用KEGG 的API獲取一個物種所有基因對應的pathway注釋，human對應的API 鏈接如下
http://rest.kegg.jp/link/hsa/pathway 人類
http://rest.kegg.jp/link/soe/pathway 菠菜
path:hsa00010 hsa:10327
path:hsa00010 hsa:124
path:hsa00010 hsa:125
第一列為pathway編號，第二列為基因編號。這里只提供了pathway編號，我們還需要pathway對應的描述信息，同樣也可以通過以下API鏈接得到
http://rest.kegg.jp/list/
通過該鏈接可以獲得如下內容
path:map00010 Glycolysis / Gluconeogenesis
path:map00020 Citrate cycle (TCA cycle)
path:map00030 Pentose phosphate pathway
path:map00040 Pentose and glucuronate interconversions
path:map00051 Fructose and mannose metabolism
第一列為pathway編號，第二列為具體的描述信息。需要注意的是，pathway是一個跨物種的概念，原始的pathway編號為map或者ko加數字，對於特定物種，改成物種對應的三字母縮寫, 比如human對應hsa, 所有擁有pathway信息的物種和對應的三字母縮寫見如下鏈接
https://www.genome.jp/kegg/catalog/org_list.html
clusterProfiler也是通過KEGG API去獲取物種對應的pathway注釋，對於已有pathway注釋的物種，我們只需要知道對應的三字母縮寫， clusterProfiler就會聯網自動獲取該物種的pathway注釋信息。
和GO富集分析類似，對於KEGG的富集分析也包含以下兩種
-
過表征分析 (over representation analysis, ORA) ###先會篩選，並挑選出我們感興趣對基因
-
基因富集分析 (gene set enrichment analysis, GSEA) ###不進行篩選
enrich.KEGG.BP <- enrichKEGG(gene = test_sample, #### 差異基因ID ENTREZID
keyType = "kegg", ####key類型
organism = "soe", ###物種3字母
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
qvalueCutoff = 0.1,)
4.2.1、柱狀圖
barplot(enrich.KEGG.BP, showCategory = 10)
橫軸為該pathway的差異基因個數，縱軸為富集到的pathway的描述信息， showCategory指定展示的pathway的個數，默認展示顯著富集的top10個，即p.adjust最小的10個。注意的顏色對應p.adjust值，從小到大，對應藍色到紅色。
4.2.2、點圖
dotplot(enrich.KEGG.BP, showCategory = 10)
橫軸為GeneRatio, 代表該pathway下的差異基因個數占差異基因總數的比例，縱軸為富集到的pathway的描述信息， showCategory指定展示的pathway的個數，默認展示顯著富集的top10個，即p.adjust最小的10個。圖中點的顏色對應p.adjust的值，從小到大，對應藍色到紅色，大小對應該GO terms下的差異基因個數，個數越多，點越大。
4.2.3、Gene-Concept Network
前面的 兩款神器 兩個函數，都只能展示富集最顯著的 GO term，而函數 cnetplot() 可以將基因與生物學概念 (e.g.* GO terms or KEGG pathways) 的關系繪制成網狀圖。對於基因和富集的pathways之間的對應關系進行展示，如果一個基因位於一個pathway下，則將該基因與pathway連線，用法如下
cnetplot(enrich.KEGG.BP, showCategory = 5)
圖中灰色的點代表基因，黃色的點代表富集到的pathways, 默認畫top5富集到的pathwayss, pathways節點的大小對應富集到的基因個數。數字就是基因ID，如果需要更換，可以更換keytype，或者直接在enrich.KEGG.BP 的結果中進行相同ID更換
cnetplot(enrich.KEGG.BP,circular=T, ###畫為圈圖 colorEdge=T) ##線條用顏色區分
4.2.4、Enrichment Map
Enrichment Map 可以將富集條目和重疊的基因集整合為一個網絡圖，相互重疊的基因集則趨向於成簇，從而易於分辨功能模型。對於富集到的pathways之間的基因重疊關系進行展示，如果兩個pathway的差異基因存在重疊，說明這兩個節點存在overlap關系，在圖中用線條連接起來，用法如下
emapplot(enrich.KEGG.BP)
每個節點是一個富集到的pathway, 默認畫top30個富集到的pathways, 節點大小對應該pathway下富集到的差異基因個數，節點的顏色對應p.adjust的值，從小到大，對應藍色到紅色。
4.2.5 browseKEGG
browseKEGG(enrich.KEGG.BP,"soe00564") ###畫出某一特定pathway的圖
關於KEGG解釋，可以查看鏈接🔗
https://www.jianshu.com/p/f90ed1c52079
4.2.6 pathview 包里的 上帝視角 PATHVIEW!
pathview(gene.data = test_sample, ##是需要提供的基因向量，默認是Entrez_ID。其由gene.idtype決定
pathway.id = "soe00564", ###指的是在KEGG中的ID
species = "soe",
kegg.native = TRUE,###默認是TRUE輸出完整pathway的png格式文件，反之輸出僅是輸入的基因列表的pdf文件。
)
感覺結果和 browseKEGG 差不多
持續整理。。。
關注下方公眾號可獲得更多精彩
參考：https://www.jianshu.com/p/ae94178918bc GO
https://www.jianshu.com/p/47b5ea646932?utm_source=desktop&utm_medium=timeline GO
https://www.cnblogs.com/djx571/p/10271874.html GO
https://blog.csdn.net/weixin_43569478/article/details/83744384 KEGG
https://www.jianshu.com/p/f90ed1c52079 KEGG
https://www.jianshu.com/p/e133ab3169fa KEGG