这一次，我们来聊聊基因组注释。首先问自己一个问题，为什么要进行基因注释。
就我目前而言，它用来解决如下问题：
-
在mapping-by-sequencing的时候，我找到了一些可能的突变位点，我需要知道这些突变分别是那些基因发生突变，这些突变基因有哪些功能？
-
差异表达分析之后会得到许多的基因，这些基因有什么样的特征？如果要进行基因富集分析，不可避免就需要知道他们的GO,KEGG等注释信息。
如果一个基因没有注释信息，那么他就只是一段DNA序列，只是一个符号。你可能会很开心，因为你研究的功能并没有被大多数人所发现，说不定这就是一篇CNS级别的文章；你或许会悲伤，因为没有注释，意味着你的工作从全新的工作，也就是说你的工作量会很大。但是不管如何，你看到一个基因后，都会本能的想知道它到底有哪些功能，如同你看到一个漂亮妹子的照片，你也可能想去知道更多有关于她的信息。
对于一个或几个基因而言，NCBI，EBI,TAIR等网站够用了，但是对于高通量数据分析的结果，你还要一个一个查的话，那就是有点费劲了。（尴尬的是，我第一次寻找突变位点就是靠我手工注释结果）。
因此，本文就是介绍如何在R语言中对高通量分析结果中基因信息进行注释。
找到注释信息
目前存在大量的注释信息的数据库，我们需要一个方便的搜索工具，用于找到我们所需要的信息。Biconductor建立在R语言上的一个开源项目，旨在未高通量数据分析提供可靠的工具。项目的一个重要部分就是组织网络上大量的注释信息，方便科研人员使用。
目前最新的工具包叫做
AnnotationHub
，顾名思义，就是注释信息的中装站。通过它，能找到了几乎所有的注释资源。如果没有，你还可以根据已有的数据用它提供的函数进行构建。
安装方式很简单(首先你得装了R）：
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("AnnotationHub")
使用AnnotationHub，我们需要创建AnnotationHub对象（加载AnnotationHub这一步就不多说了).
library(AnnotationHub)
ah
上述结果告诉了我们以下信息：
根据这些知识点，我们就可以问
第一个问题
：
AnnotationHub的数据来源有哪些？
unique(ah$dataprovider)
[1] "Ensembl" "UCSC"
[3] "RefNet" "Inparanoid8"
[5] "NHLBI" "ChEA"
[7] "Pazar" "NIH Pathway Interaction Database"
[9] "Haemcode" "BroadInstitute"
[11] "PRIDE" "Gencode"
[13] "dbSNP" "CRIBI"
[15] "Genoscope" "MISO, VAST-TOOLS, UCSC"
[17] "ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/"
第二个问题
：
AnnotationHub目前支持哪些物种？我想找的物种在这里面么？
unique(ah$species)
由于结果有391个，不方便查询。但是可以通过筛选，找找目标物种是否存在。
ah$species[which(ah$species == "Arabidopsis thaliana")]
[1] "Arabidopsis thaliana" "Arabidopsis thaliana" "Arabidopsis thaliana" "Arabidopsis thaliana"
[5] "Arabidopsis thaliana"
通过它提供的
query
函数，去搜索ah对象，就能判断目标物种是否被AnnotationHub收录。
query(x, pattern, ignore.case=TRUE, pattern.op=
&
)
Return an AnnotationHub subset containing only those elements whose metadata matches pattern. Matching uses pattern as in grepl to search the as.character representation of each column, performing a logical
&
across columns. e.g., query(x, c(“Homo sapiens”, “hg19”, “GTF”))
比如说我想查找拟南芥相关的注释数据库，就可以去query去查找在metadata里面想关信息。
grs
当然我们还可以用R本身的筛选功能
> ah[ah$species == 'Arabidopsis thaliana' & ah$rdataclass == 'OrgDb']
> subset(ah, species == 'Arabidopsis thaliana' & rdataclass == 'OrgDb')
搜索到的记录就只有如下几个了。
AnnotationHub with 1 record
# snapshotDate(): 2017-04-25
# names(): AH53758
# $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# $species: Arabidopsis thaliana
# $rdataclass: OrgDb
# $rdatadateadded: 2017-04-10
# $title: org.At.tair.db.sqlite
# $description: NCBI gene ID based annotations about Arabidopsis thaliana
# $taxonomyid: 3702
# $genome: NCBI genomes
# $sourcetype: NCBI/ensembl
# $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.ensembl.org/pub/current...
# $sourcesize: NA
# $tags: c("NCBI", "Gene", "Annotation")
# retrieve record with 'object[["AH53758"]]'
如果我们酷爱图形界面（GUI），类似于网页搜索那样的操作，可以使用的是
display
函数了。
display(ah)
第三个问题
：
AnnotationHub的注释信息的数据存放格式是什么？
unique(ah$rdataclass)
[1] "FaFile" "GRanges" "data.frame" "Inparanoid8Db" "TwoBitFile"
[6] "ChainFile" "SQLiteConnection" "biopax" "BigWigFile" "AAStringSet"
[11] "MSnSet" "mzRpwiz" "mzRident" "VcfFile" "list"
[16] "TxDb" "Rle" "EnsDb" "OrgDb"
比如说fasta文件是FaFile, 转录组数据库是TxDb, 提供内含子、外显子、UTR区的信息。有物种数据库，OrgDb,用于基因ID，基因名,GO,KEGG ID之间的相互映射。
第四个问题
：
我们如何去下载所需信息
第二个问题后，你会得到一个ID，比如说拟南芥的OrgDb的注释数据库的ID就是”AH53758”，然后根据这个ID可以进行下载。当然下载方式已经出现过了，
retrieve record with ‘object[[“AH53758”]]’
ath
bioconductor除了
AnnotationHub
能用来查找生物数据，还有一个库叫做
biomaRt
，可以用来查找biomart中的数据。不过目前biomart网站正在进行服务器的数据迁移，就不在这里演示。
小结
：
-
AnnotationHub是生物数据的中转站，方面我们搜索目标数据，另一个相似包是
biomaRt
;
-
我们通过query,subset等方法（图形界面则是display），逐步从AnnotationHub的metadata筛选到所需数据的ID;
-
使用
[]
是查看目标数据的metadata, 使用
[[]]
用于下载数据;
探索注释数据库
找到和下载注释数据库只是第一步，学会如何使用这些数据库更加重要。
AnnotationHub对象的通用方法
之前下载完数据后，在R里面被我指向到了’ath’，那么我们先简单了解一下这个’ath’
直接输入对象名
ath
，显示的就是元数据信息，太长不放。
用
str
了解一下它的数据结构.好吧，我承认我自己看不出名堂。只知道他是AnnotationDbi的OrgDb类。
> str(ath)
Reference class 'OrgDb' [package "AnnotationDbi"] with 2 fields
$ conn :Formal class 'SQLiteConnection' [package "RSQLite"] with 6 slots
.. ..@ ptr :
.. ..@ dbname : chr "D:\\xuzho\\Documents\\AppData\\.AnnotationHub\\60496"
.. ..@ loadable.extensions: logi TRUE
.. ..@ flags : int 1
.. ..@ vfs : chr ""
.. ..@ ref :
$ packageName: chr(0)
and 14 methods.
用
mode
看下它的数据模式(Get or set the type or storage mode of an object.)，发现它是一个S4类。大部分bioconductor的包都是S4类，然而什么是S4类呢？在R语言编程艺术，我看到过这个概念，主要和S3类区别，据说更加安全。
mode(ath)
[1] "S4"
用
class
看下它具体继承什么类（面向对象编程的概念）
class(ath)
[1] "OrgDb"
attr(,"package")
[1] "AnnotationDbi"
好了，我们继续调查什么是”AnnotationDbi”,了解到他主要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，只不过就返回一个列。
返回这个数据包都有哪些列:
> columns(ath)
[1] "ARACYC" "ARACYCENZYME" "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[8] "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PMID" "REFSEQ"
[15] "SYMBOL" "TAIR"
返回这个数据包可以当做关键字(key)进行查找的列:
> keytypes(ath)
[1] "ARACYC" "ARACYCENZYME" "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[8] "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PMID" "REFSEQ"
[15] "SYMBOL" "TAIR"
基本上
keytypes
返回的结果是等于或者少于columns返回的结果。因为并不是所有列都能当做查找对象。
