前言
關於clusterProfiler這個R包就不介紹了，網紅教授宣傳得很成功，功能也比較強大，主要是做GO和KEGG的功能富集及其可視化。簡單總結下用法，以后用時可直接找來用。
首先考慮一個問題：clusterProfiler做GO和KEGG富集分析的注釋信息來自哪里？
GO的注釋信息來自Bioconductor，提供了19個物種的org類型的GO注釋信息，如下表所示。Bioconductor中更多的注釋包可參考http://www.bioconductor.org/packages/release/data/annotation/，很亂，大多數我都不知道干啥用的。
| packages | organism |
|---|---|
| org.Ag.eg.db | Anopheles |
| org.At.tair.db | Arabidopsis |
| org.Bt.eg.db | Bovine |
| org.Ce.eg.db | Worm |
| org.Cf.eg.db | Canine |
| org.Dm.eg.db | Fly |
| org.Dr.eg.db | Zebrafish |
| org.EcK12.eg.db | E coli strain K12 |
| org.EcSakai.eg.db | E coli strain Sakai |
| org.Gg.eg.db | Chicken |
| org.Hs.eg.db | Human |
| org.Mm.eg.db | Mouse |
| org.Mmu.eg.db | Rhesus |
| org.Pf.plasmo.db | Malaria |
| org.Pt.eg.db | Chimp |
| org.Rn.eg.db | Rat |
| org.Sc.sgd.db | Yeast |
| org.Ss.eg.db | Pig |
| org.Xl.eg.db | Xenopus |
KEGG的注釋信息clusterProfiler通過KEGG 數據庫的API來獲取，https://www.kegg.jp/kegg/rest/keggapi.html。
首先是一個物種所有基因對應的pathway注釋文件，比如人的：http://rest.kegg.jp/link/hsa/pathway。
其次還需要pathway對應的描述信息，比如人的：
http://rest.kegg.jp/list/pathway/hsa。
關於KEGG數據庫全部的物種及其簡寫（三個字母）如下列表：
https://www.genome.jp/kegg/catalog/org_list.html。
因此對於以上已有pathway注釋的物種，只需要將物種簡寫輸入給clusterProfiler， 它會通過聯網自動獲取該物種的pathway注釋信息。
以上都是有物種信息的情況，那么對於無物種信息的項目怎么辦？
GO可以通過讀取外部的GO注釋文件進行分析。關於基因的GO注釋，interproscan、eggnog-mapper和blas2go等軟件都可以做，不過輸出格式有些不同。clusterProfiler需要導入的GO注釋文件的格式如下：
| GeneID | GO | GO_Description |
|---|---|---|
| 1 | GO:0005819 | spindle |
| 2 | GO:0072686 | mitotic spindle |
| 3 | GO:0000776 | kinetochore |
需要包含以上三列信息，這3列信息任意順序都可。
clusterProfiler包只針對含有OrgDb對象，如果是公共數據庫中有該物種注釋信息，只是未制作成org.db數據庫（標准注釋庫），則可以不需要從頭注釋，只需手動制作org.db數據庫類型，完成后直接使用即可，代碼如下：
source("https://bioconductor.org/biocLite.R")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("AnnotationHub") # 一個包含大量注釋信息的數據庫，里面有很多物種及來源於很多數據庫的注釋信息。
BiocManager::install("biomaRt")
library(AnnotationHub)
library(biomaRt)
hub <- AnnotationHub() #建立AnnotationHub對象（視人品，網不行加載不了）
# unique(hub$species) #查看AnonotationHub里面物種
hub$species[which(hub$species=="Solanum")] #看AnonotationHub里是否包含想要的物種
# Solanum是番茄的拉丁名
query(hub, "Solanum") #查看該物種信息
hub[hub$species=="Solanum" & hub$rdataclass == "OrgDb"] #OrgDb屬於rdataclass中，因此查看下該物種有沒有OrgDb
Solanum.OrgDb <- hub[["AH59087"]]#AH59087是番茄對應的編號
#制作為標准注釋庫，就可和模式生物一樣使用了
同樣地，對於pathway數據庫中沒有的物種，也支持讀取基因的pathway注釋文件，然后進行分析，注釋文件的格式如下：
| GeneID | Pathway | Path_Description |
|---|---|---|
| 1 | ko:00001 | spindle |
| 2 | ko:00002 | mitotic spindle |
| 3 | ko:00003 | kinetochore |
以上三列信息的順序也是任意的。
富集分析
通常用的富集分析有ORA、FCS和拓撲三種方法。ORA簡單來說就是超幾何檢驗或Fisher精確檢驗，大同小異，都符合超幾何檢驗，這也是目前用的最多的方法，優劣不談。FCS的代表就是GSEA，即基因集富集分析，優劣亦不談。clusterProfiler提供了這兩種富集分析方法。
1. ORA（Over-Representation Analysis）
GO富集參考代碼：
#標准富集分析
ego <- enrichGO(
gene = gene$entrzID,
keyType = "ENTREZID",
universe = names(geneList), #背景基因集，可省
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
#通過導入外部注釋文件富集分析
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)
gene差異基因對應的向量;
keyType指定基因ID的類型，默認為ENTREZID, 可參考keytypes(org.Hs.eg.db)類型 ；
OrgDb指定該物種對應的org包的名字；
ont代表GO的3大類別，BP, CC, MF，也可是全部ALL；
pAdjustMethod指定多重假設檢驗矯正的方法，有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一種；
cufoff指定對應的閾值；
readable=TRUE代表將基因ID轉換為gene symbol。
KEGG Pathway富集參考代碼：
#標准富集分析
ego <- enrichKEGG(
gene = gene,
keyType = "kegg",
organism = 'hsa',
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
qvalueCutoff = 0.05
)
#通過外部導入注釋文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)
默認基因ID為kegg gene id，也可以是ncbi-geneid, ncbi-proteinid, uniprot等。
organism物種對應的三字母縮寫，其他參數同GO富集。ID轉換函數：
library(clusterProfiler)
bitr_kegg("1",fromType = "kegg",toType = 'ncbi-proteinid',organism='hsa')
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) #支持的ID類型
bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)
#以上看出ID轉換輸入時，可以向量的形式，也可以單列基因名list導入，也可以是內置數據
gene <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")
gene <- data$V1 #字符串
data(geneList, package="DOSE") #富集分析的背景基因集
gene <- names(geneList)[abs(geneList) > 2]
2. GSEA（Gene Set Enrichment Analysis）
GO富集參考代碼：
#標准富集分析
ego <- gseGO(
geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "CC",
nPerm = 1000, #置換檢驗的置換次數
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)
#通過導入外部注釋文件富集分析參考代碼：
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)
KEGG Pathway富集參考代碼：
#標准富集分析
kk <- gseKEGG(
geneList = gene,
keyType = 'kegg',
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.05,
pAdjustMethod = "BH"
)
#通過外部導入注釋文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)
可視化
1.GO富集分析結果可視化
#barplot
barplot(ego, showCategory = 10) #默認展示顯著富集的top10個，即p.adjust最小的10個
#dotplot
dotplot(ego, showCategory = 10)
#DAG有向無環圖
plotGOgraph(ego) #矩形代表富集到的top10個GO terms, 顏色從黃色過濾到紅色，對應p值從大到小。
#igraph布局的DAG
goplot(ego)
#GO terms關系網絡圖（通過差異基因關聯）
emapplot(ego, showCategory = 30)
#GO term與差異基因關系網絡圖
cnetplot(ego, showCategory = 5)
2.Pathway富集分析結果可視化
#barplot
barplot(kk, showCategory = 10)
#dotplot
dotplot(kk, showCategory = 10)
#pathway關系網絡圖（通過差異基因關聯）
emapplot(kk, showCategory = 30)
#pathway與差異基因關系網絡圖
cnetplot(kk, showCategory = 5)
#pathway映射
browseKEGG(kk, "hsa04934") #在pathway通路圖上標記富集到的基因，會鏈接到KEGG官網
Ref：
https://blog.csdn.net/weixin_43569478/article/details/83744242
https://blog.csdn.net/weixin_43569478/article/details/83744384
https://www.jianshu.com/p/065d38c28e2d
https://www.jianshu.com/p/47b5ea646932
https://www.cnblogs.com/yatouhetademao/p/8018252.html
https://zhuanlan.zhihu.com/p/35510434