非模式基因GO富集分析：以玉米为例+使用OrgDb
模式生物做什么都简单，非模式生物则很多缺少注释，没有注释你就没法做，只能是借助于各种软件比如blastgo，自己跑电子注释。但今天要讲的不是这种情况，很多物种还是有注释的，只是你有时候不知道该去那里下载，或者你有数据，却不知道该怎么用！很多的软件都是针对模式生物的，或者针对某一些类型的非模式生物，能够支持多种非模式生物，能够支持用户自己的注释文件的软件相对来讲，就非常少有了，然而clusterProfiler就是这类少有的软件之一。
获得OrgDb
今天要讲的是通过OrgDb来做GO分析，这是clusterProfiler的enrichGO函数所支持的背景注释，Bioconductor自带20个OrgDb可供使用，多半是模式生物，难道我们要做的物种不在这20个里面就不行了吗？显然不是的，clusterProfiler能支持的物种我自己都数不过来。
我们可以通过AnnotationHub在线检索并抓取OrgDb，比如这里以玉米为例：
> require(AnnotationHub)
> hub <- AnnotationHub()
> query(hub, "zea")
AnnotationHub with 2 records
# snapshotDate(): 2017-04-25
# $dataprovider: Inparanoid8, ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# $species: Gibberella zeae, Zea mays
# $rdataclass: Inparanoid8Db, OrgDb
# additional mcols(): taxonomyid, genome, description,
# coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
# rdatapath, sourceurl, sourcetype
# retrieve records with, e.g., 'object[["AH10514"]]'
title
AH10514 | hom.Gibberella_zeae.inp8.sqlite
AH55736 | org.Zea_mays.eg.sqlite
通过检索，org.Zea_mays.eg.sqlite就是我们所要的OrgDb，可以通过相应的accession number, AH55736抓取文件，并存入了maize对象中，它包含了51097个基因的注释：
> maize <- hub[['AH55736']]
> length(keys(maize))
[1] 51097
这个OrgDb，包含有以下一些注释信息：
> columns(maize)
[1] "ACCNUM" "ALIAS" "CHR" "ENTREZID" "EVIDENCE"
[6] "EVIDENCEALL" "GENENAME" "GID" "GO" "GOALL"
[11] "ONTOLOGY" "ONTOLOGYALL" "PMID" "REFSEQ" "SYMBOL"
[16] "UNIGENE"
转换ID
我们可以使用bitr来转换ID，甚至于直接检索GO注释：
> require(clusterProfiler)
> bitr(keys(maize)[1], 'ENTREZID', c("REFSEQ", "GO", "ONTOLOGY"), maize)
ENTREZID REFSEQ GO ONTOLOGY
1 541612 XP_008648268.1 GO:0009507 CC
2 541612 XP_008648268.1 GO:0051537 MF
3 541612 XP_008648268.1 GO:0009055 MF
4 541612 XP_008648268.1 GO:0046872 MF
5 541612 XP_008648268.1 GO:0022900 BP
6 541612 NP_001104837.2 GO:0009507 CC
7 541612 NP_001104837.2 GO:0051537 MF
8 541612 NP_001104837.2 GO:0009055 MF
9 541612 NP_001104837.2 GO:0046872 MF
10 541612 NP_001104837.2 GO:0022900 BP
11 541612 XM_008650046.2 GO:0009507 CC
12 541612 XM_008650046.2 GO:0051537 MF
13 541612 XM_008650046.2 GO:0009055 MF
14 541612 XM_008650046.2 GO:0046872 MF
15 541612 XM_008650046.2 GO:0022900 BP
16 541612 NM_001111367.2 GO:0009507 CC
17 541612 NM_001111367.2 GO:0051537 MF
18 541612 NM_001111367.2 GO:0009055 MF
19 541612 NM_001111367.2 GO:0046872 MF
20 541612 NM_001111367.2 GO:0022900 BP
GO富集分析
> sample_genes <- keys(maize)[1:100]
> head(sample_genes)
[1] "541612" "541613" "541614" "541615" "541617" "541618"
这里我只是简单地使用ID列表中前100个ENTREZ基因ID，也可以使用其它的ID，通过借助于bitr进行转换，或者通过给enrichGO指定ID类型(keyType参数）。
有了OrgDb，使用起来，就跟文档中使用人类基因做为例子一样，用法一致，并且也可以通过clusterProfiler所提供的各种可视化函数对结果进行展示：
> require(clusterProfiler)
> res = enrichGO(sample_genes, OrgDb=maize, pvalueCutoff=1, qvalueCutoff=1)
> res
#
# over-representation test
#
#...@organism Zea mays
#...@ontology MF
#...@keytype ENTREZID
#...@gene chr [1:100] "541612" "541613" "541614" "541615" "541617" "541618" "541619" ...
#...pvalues adjusted by 'BH' with cutoff <1
#...114 enriched terms found
'data.frame': 114 obs. of 9 variables:
$ ID : chr "GO:0004871" "GO:0000155" "GO:0004673" "GO:0016775" ...
$ Description: chr "signal transducer activity" "phosphorelay sensor kinase activity" "protein histidine kinase activity" "phosphotransferase activity, nitrogenous group as acceptor" ...
$ GeneRatio : chr "9/80" "5/80" "5/80" "5/80" ...
$ BgRatio : chr "81/14230" "22/14230" "23/14230" "23/14230" ...
$ pvalue : num 6.65e-10 1.21e-07 1.54e-07 1.54e-07 1.90e-07 ...
$ p.adjust : num 7.58e-08 3.60e-06 3.60e-06 3.60e-06 3.60e-06 ...
$ qvalue : num 6.37e-08 3.03e-06 3.03e-06 3.03e-06 3.03e-06 ...
$ geneID : chr "541618/541625/541627/541634/541636/541638/541641/541642/541663" "541627/541634/541641/541642/541663" "541627/541634/541641/541642/541663" "541627/541634/541641/541642/541663" ...
$ Count : int 9 5 5 5 4 4 4 5 5 5 ...
#...Citation
Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
clusterProfiler: an R package for comparing biological themes among
gene clusters. OMICS: A Journal of Integrative Biology
2012, 16(5):284-287
如果你有表达谱数据，你也可以使用gseGO进行GSEA分析，这里我懒得上网找数据来演示了，用法反正跟文档里的一样，只不过换成了你自己的数据和相应物种的OrgDb对象而已。
如果没有OrgDb怎么办？
必须也是可以做的，这个将在以后讲解！
