从biomaRt到clusterProfiler
1.安装
source("http://bioconductor.org/biocLite.R")
biocLite('biomaRt')
biocLite('clusterProfiler')
R 版本为3.5以上的需要以下命令
BiocManager::install('biomaRt')
BiocManager::install('clusterProfiler')
如果提示缺哪个包，则用以下命令安装：
install.packages('')#''内为缺失的包名
2.生物学ID转换
——做KEGG和GO分析之前需要对非NCBI ID的转换为NCBI ID，如山羊的Hoxc13基因的NCBI ID为102178245，人的Hoxc13基因的NCBI ID为3229。如果你在NCBI https://www.ncbi.nlm.nih.gov/的gene这一栏中输入102178245，那么它反馈回来的为山羊的Hoxc13基因，如果你输入3229，那么它反馈回来的为人的Hoxc13基因。
2.1 选择数据库
library(biomaRt)
mart <- useMart("ensembl","hsapiens_gene_ensembl")
——useMart第二个双引号内为需要使用的数据库，根据物种的不同，选择不同的数据库。如，人的选择hsapiens_gene_ensembl，鼠的选择mmusculus_gene_ensembl
2.2 选择输入ID类型
——然后，你需要知道输入的是什么类型的ID，即转换前的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类型的输入，可以通过以下命令查看，共380种ID。
listFilters(mart)#输入的ID类型
2.3 选择输出ID类型
——要想知道biomaRt支持哪些ID类型的输出，可以通过以下命令查看，共支持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 导入文件和选定向量
gene <- read.csv(“data.csv”,sep=',',header=T)[,4]
因为需要转换的ID在数据集data.csv的第四列。
题外话
“data.csv”只是我自己的一个例子，如果你的为xlxs格式，可以另存为.csv格式，然后在通过read.csv函数导入；如果你的为txt格式，通过read.table函数导入，如：
gene <- read.table(“data.txt”,sep='\t',header=T)
其中双引号内的“data.txt”，有两种书写方法：
1，输入文件的完整路径（注意反斜线），如：“G:/Download/data.csv”
2，通过R菜单将工作目录切换到data.csv所在目录，“”内只写文件名即可。
gene # 查看结果为向量
#[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 开始转换
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从矩阵变为向量，因为clusterProfiler包的enrichKEGG函数需要输入为向量。
这次结果中，输入49个"external_gene_name"，输出48个"entrezgene"。
2.6 clusterProfiler使用
2.6.1 KEGG分析
library('clusterProfiler')
kegg <- enrichKEGG(geneID, organism = "hsa", keyType = "kegg", use_internal_data = FALSE)
如果是使用clusterProfiler包的bitr函数进行ID转换，则可以这样做：
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个 "SYMBOL",输出45个"ENTREZID"。所以有4.35% of input gene IDs are fail to map。"SYMBOL"和biomaRt包的"external_gene_name"，"ENTREZID"和biomaRt包的"entrezgene"其实是同一种类型，只不过两个包名字不同而已。
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,显示条目的数目
barplot(GO_BP,showCategory=20,drop=T)
barplot(GO_CC,showCategory=20,drop=T)
barplot(GO_MF,showCategory=20,drop=T)
