1.安裝
source("http://bioconductor.org/biocLite.R")
biocLite('biomaRt')
biocLite('clusterProfiler')
R 版本為3.5以上的需要以下命令
BiocManager::install('biomaRt')
BiocManager::install('clusterProfiler')
如果提示缺哪個(gè)包，則用以下命令安裝：
install.packages('')#''內(nèi)為缺失的包名
2.生物學(xué)ID轉(zhuǎn)換
——做KEGG和GO分析之前需要對(duì)非NCBI ID的轉(zhuǎn)換為NCBI ID，如山羊的Hoxc13基因的NCBI ID為102178245，人的Hoxc13基因的NCBI ID為3229。如果你在NCBI https://www.ncbi.nlm.nih.gov/的gene這一欄中輸入102178245，那么它反饋回來(lái)的為山羊的Hoxc13基因，如果你輸入3229，那么它反饋回來(lái)的為人的Hoxc13基因。
2.1 選擇數(shù)據(jù)庫(kù)
library(biomaRt)
mart <- useMart("ensembl","hsapiens_gene_ensembl")
——useMart第二個(gè)雙引號(hào)內(nèi)為需要使用的數(shù)據(jù)庫(kù)，根據(jù)物種的不同，選擇不同的數(shù)據(jù)庫(kù)。如，人的選擇hsapiens_gene_ensembl，鼠的選擇mmusculus_gene_ensembl
2.2 選擇輸入ID類型
——然后，你需要知道輸入的是什么類型的ID，即轉(zhuǎn)換前的ID類型，如：
ENSG00000000003或ENMUSG000000003，屬于類型為ensembl_gene_id；
ENST00000000233或ENMUST00000000233，屬于類型為ensembl_transcript_id；
102178245，屬于類型為entrezgene；
Hoxc13，屬于類型為external_gene_name；
NM_000014，屬于類型為refseq_mrna；
hsa-let-7a-1，屬于類型為mirbase_id；
要想知道biomaRt支持哪些ID類型的輸入，可以通過(guò)以下命令查看，共380種ID。
listFilters(mart)#輸入的ID類型
2.3 選擇輸出ID類型
——要想知道biomaRt支持哪些ID類型的輸出，可以通過(guò)以下命令查看，共支持2722種ID輸出。
listAttributes(mart)
常用的：
Gene description：
Gene name：
Chromosome/scaffold:
Gene start (bp):
Gene end (bp):
KEGG Pathway and Enzyme:
NCBI gene:
PDB:
entrezgene_id,即NCBI gene ID
external_gene_name
2.4 導(dǎo)入文件和選定向量
gene <- read.csv(“data.csv”,sep=',',header=T)[,4]
因?yàn)樾枰D(zhuǎn)換的ID在數(shù)據(jù)集data.csv的第四列。
題外話
“data.csv”只是我自己的一個(gè)例子，如果你的為xlxs格式，可以另存為.csv格式，然后在通過(guò)read.csv函數(shù)導(dǎo)入；如果你的為txt格式，通過(guò)read.table函數(shù)導(dǎo)入，如：
gene <- read.table(“data.txt”,sep='\t',header=T)
其中雙引號(hào)內(nèi)的“data.txt”，有兩種書寫方法：
1，輸入文件的完整路徑（注意反斜線），如：“G:/Download/data.csv”
2，通過(guò)R菜單將工作目錄切換到data.csv所在目錄，“”內(nèi)只寫文件名即可。
gene # 查看結(jié)果為向量
#[1] NFKB1 CDKN1A MMP16 KIT Card10 Scube2 TRAF6 IRAK1 TLR4 UHRF1 PDGFRA S100A12 ZNF117 MYO6
#[15] CCDC6 IL1RAP IL1RL2 HNRNPD ZNRF3 IL6 PAX8 NOVA1 RARB EGFR ERBB4 SLC5A5 MALAT1 BSG
#[29] YWHAZ WEE1 TRPS1 MBNL2 KLF13 CDKN1A LATS2 ERBB4 NR4A2 VEGFA TGFBR2 RHOC NFIB CDK2
#[43] CCNA1 TNFAIP1 DKK1 BTG1 LEFTY1 TXNIP ATAD2
#47 Levels: ATAD2 BSG BTG1 Card10 CCDC6 CCNA1 CDK2 CDKN1A DKK1 EGFR ERBB4 HNRNPD IL1RAP IL1RL2 IL6 IRAK1 KIT ... ZNRF3
2.5 開始轉(zhuǎn)換
geneID <- getBM(attributes=c("entrezgene_id"),filters = "external_gene_name",values = gene, mart = mart)[,1]
geneID # 查看
#[1] 7126 8808 8900 6528 7048 2066 109504726 51351 26524 84133 7189
#[12] 1026 29775 29128 22943 57758 3569 4325 7227 378938 7422 5156
#[23] 29028 51621 4790 10628 389 694 3815 7849 7465 5915 3556
#[34] 8030 3184 4646 1956 4857 4781 7099 3654 682 10637 4929
#[45] 6283 10150 7534 1017
filters，指定輸入ID的類型，
attributes指定輸出ID的類型，
[,1]，取第一列，是為了把輸出的ID從矩陣變?yōu)橄蛄浚驗(yàn)閏lusterProfiler包的enrichKEGG函數(shù)需要輸入為向量。
這次結(jié)果中，輸入49個(gè)"external_gene_name"，輸出48個(gè)"entrezgene"。
2.6 clusterProfiler使用
2.6.1 KEGG分析
library('clusterProfiler')
kegg <- enrichKEGG(geneID, organism = "hsa", keyType = "kegg", use_internal_data = FALSE)
如果是使用clusterProfiler包的bitr函數(shù)進(jìn)行ID轉(zhuǎn)換，則可以這樣做：
library('org.Hs.eg.db')
library('org.Mm.eg.db')#物種為鼠
ENTREZID<- bitr(gene, fromType = "SYMBOL", toType="ENTREZID",OrgDb = org.Hs.eg.db)
#Warning message:
#In bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) :
#4.35% of input gene IDs are fail to map...
此外還有"ENSEMBL"。
geneID<-ENTREZID[,2]#取第二列
geneID
#[1] "4790" "1026" "4325" "3815" "7189" "3654" "7099" "29128" "5156" "6283" "51351" "4646" "8030"
#[14] "3556" "8808" "3184" "84133" "3569" "7849" "4857" "5915" "1956" "2066" "6528" "378938" "682"
#[27] "7534" "7465" "7227" "10150" "51621" "26524" "4929" "7422" "7048" "389" "4781" "1017" "8900"
#[40] "7126" "22943" "694" "10637" "10628" "29028"
kegg <- enrichKEGG(geneID, organism = "hsa", keyType = "kegg", use_internal_data = TRUE)
可以看到clusterProfiler包輸入49個(gè) "SYMBOL",輸出45個(gè)"ENTREZID"。所以有4.35% of input gene IDs are fail to map。"SYMBOL"和biomaRt包的"external_gene_name"，"ENTREZID"和biomaRt包的"entrezgene"其實(shí)是同一種類型，只不過(guò)兩個(gè)包名字不同而已。
https://www.genome.jp/kegg/catalog/org_list.html查看支持的物種（organisms）。
人：hsa；
鼠：mmu。
2.5.2 GO分析
GO_BP <- enrichGO(geneID, OrgDb = org.Hs.eg.db, ont='BP',pAdjustMethod = 'BH', qvalueCutoff = 0.2, pvalueCutoff = 0.05,keyType = 'ENTREZID')
GO_CC <- enrichGO(geneID, OrgDb = org.Hs.eg.db, ont='CC',pAdjustMethod = 'BH', qvalueCutoff = 0.2, pvalueCutoff = 0.05,keyType = 'ENTREZID')
GO_MF <- enrichGO(geneID, OrgDb = org.Hs.eg.db, ont='MF',pAdjustMethod = 'BH', qvalueCutoff = 0.2, pvalueCutoff = 0.05,keyType = 'ENTREZID')
2.6 作圖
dotplot(kegg,showCategory=20) #氣泡圖
barplot(kegg,showCategory=20,drop=T) #showCategory,顯示條目的數(shù)目
barplot(GO_BP,showCategory=20,drop=T)
barplot(GO_CC,showCategory=20,drop=T)
barplot(GO_MF,showCategory=20,drop=T)