正文
这一次，我们来聊聊基因组注释。首先问自己一个问题，为什么要进行基因注释。
就我目前而言，它用来解决如下问题：
在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()
ah
AnnotationHub with 39213 records
# snapshotDate(): 2017-04-25
# $dataprovider: BroadInstitute, Ensembl, UCSC, Haemcode, Inparanoid8, Pazar, Gencode, ftp://ftp.ncbi.nlm...
# $species: Homo sapiens, Mus musculus, Bos taurus, Pan troglodytes, Danio rerio, Rattus norvegicus, Dros...
# $rdataclass: GRanges, BigWigFile, FaFile, TwoBitFile, ChainFile, Rle, Inparanoid8Db, EnsDb, TxDb, data....
# additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer,
# rdatadateadded, preparerclass, tags, rdatapath, sourceurl, sourcetype
# retrieve records with, e.g., 'object[["AH2"]]'
title
AH2 | Ailuropoda_melanoleuca.ailMel1.69.dna.toplevel.fa
AH3 | Ailuropoda_melanoleuca.ailMel1.69.dna_rm.toplevel.fa
AH4 | Ailuropoda_melanoleuca.ailMel1.69.dna_sm.toplevel.fa
AH5 | Ailuropoda_melanoleuca.ailMel1.69.ncrna.fa
AH6 | Ailuropoda_melanoleuca.ailMel1.69.pep.all.fa
... ...
AH54627 | Xiphophorus_maculatus.Xipmac4.4.2.cdna.all.2bit
AH54628 | Xiphophorus_maculatus.Xipmac4.4.2.dna.toplevel.2bit
AH54629 | Xiphophorus_maculatus.Xipmac4.4.2.dna_rm.toplevel.2bit
AH54630 | Xiphophorus_maculatus.Xipmac4.4.2.dna_sm.toplevel.2bit
AH54631 | Xiphophorus_maculatus.Xipmac4.4.2.ncrna.2bit
上述结果告诉了我们以下信息：
根据这些知识点，我们就可以问第一个问题：
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 <- query(ah, "Arabidopsis thaliana")
grs
...
title
AH10456 | hom.Arabidopsis_thaliana.inp8.sqlite
AH52245 | TxDb.Athaliana.BioMart.plantsmart22.sqlite
AH52246 | TxDb.Athaliana.BioMart.plantsmart25.sqlite
AH52247 | TxDb.Athaliana.BioMart.plantsmart28.sqlite
AH53758 | org.At.tair.db.sqlite
当然我们还可以用R本身的筛选功能
> ah[ahrdataclass == '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 <- ah[['AH53758']]
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返回的结果。因为并不是所有列都能当做查找对象。
keytypes告诉我们可以当做哪些列是keytype类型，那么keys则列出这个keytype下有哪些关键字。
head(keys(ath,keytype = "SYMBOL"))
select则是根据你提供的key值去查找注释数据库，返回你需要的columns信息。
> select(ath, keys= "AGO1", columns=c("TAIR","GO"),keytype = "SYMBOL")
'select()' returned 1:many mapping between keys and columns
SYMBOL TAIR GO EVIDENCE ONTOLOGY
1 AGO1 AT1G48410 GO:0004521 IDA MF
2 AGO1 AT1G48410 GO:0005515 IPI MF
3 AGO1 AT1G48410 GO:0005634 IDA CC
4 AGO1 AT1G48410 GO:0005737 IDA CC
5 AGO1 AT1G48410 GO:0005737 ISM CC
6 AGO1 AT1G48410 GO:0005737 TAS CC
7 AGO1 AT1G48410 GO:0005829 IDA CC
8 AGO1 AT1G48410 GO:0006306 RCA BP
9 AGO1 AT1G48410 GO:0006342 RCA BP
10 AGO1 AT1G48410 GO:0006346 RCA BP
...
比如说一个AGO1就能返回那么多信息.因为一个基因可以有多个功能（GO注释），当然一个GO注释下也有可以多个基因。
> select(ath, keys= "GO:0004521", columns=c("TAIR"),keytype = "GO")
'select()' returned 1:many mapping between keys and columns
GO EVIDENCE ONTOLOGY TAIR
1 GO:0004521 ISS MF AT1G14210
2 GO:0004521 ISS MF AT1G14220
3 GO:0004521 ISS MF AT1G26820
4 GO:0004521 IDA MF AT1G30460
5 GO:0004521 IDA MF AT1G48410
6 GO:0004521 ISS MF AT2G02990
7 GO:0004521 IDA MF AT2G04270
8 GO:0004521 IMP MF AT2G17520
9 GO:0004521 IDA MF AT2G39780
10 GO:0004521 ISS MF AT2G39780
11 GO:0004521 IDA MF AT2G40410
12 GO:0004521 ISS MF AT3G04480
13 GO:0004521 ISS MF AT3G20390
14 GO:0004521 IDA MF AT5G24360
15 GO:0004521 IDA MF AT5G41190
富集分析就是看不同GO，KEGG注释下，你提供的基因集的分布情况。比如说我随机从拟南芥中抽样200个基因，然后观察这些基因的富集情况。
注：这里用的Y叔的clusterProfiler包
library("clusterProfiler")
tair.sample <- sample(keys(ath,keytype = "ENTREZID"), 100)
library("clusterProfiler")
test <- enrichGO(gene = tair.sample,
OrgDb = ath,
keytype = "ENTREZID",
pAdjustMethod = "none",
pvalueCutoff = 0.1,
qvalueCutoff = 0.2)
summary(test)
由于随机取样，很有可能找不到任何富集，可也是有有可能的。
PS: AGO1没有KEGG注释，毕竟他不属于代谢途径。
mapIds功能和select类似，只不过他返回的是一组向量，而不是数据库。
mapIds(ath, keys = tair.sample, column = c("TAIR"), keytype = "ENTREZID")
这部分的小结就是记住5个函数：
columns
keytypes
keys
select
mapIds
注意，基本上所有的注释数据库，都能先用这些函数看个大概。
TxDb对象的专门方法
因为TxDb注释对象包含转录组数据，所以有一些特殊方法。以拟南芥的TxDb为例。
txdb <- ah[['AH52247']]txdb
第一类方法：根据转录本类型提取内容，transcrpits(), exons(), cds(), genes() and promoters()
transcripts(txdb)
cds(txdb)
genes(txdb)
exons(txdb)
第二类方法： 根据某一类特征转录组进行分类，如transcriptsBy(), exonsBy(), cdsBy(), intronsByTranscript(), fiveUTRsByTranscript() and threeUTRsByTranscript().
transcriptsBy(txdb, by="gene")
exonsBy(txdb, by="gene")
第三类方法：染色体相关函数： seqinfo(), seqlevels(), seqlevelsStyle(), isActivateSeq()
# 染色体命名
seqinfo(txdb)
# 染色体水平
seqlevels(txdb)
# 染色体命名方法
seqlevelsStyle(txdb)
# 决定处理那些染色体
isActiveSeq(txdb)
BSgenome对象专门函数
BSgenome存放的是基因组序列数据，好像无法被AnnotationHub找到（但是共享通用函数）,所以需要加载BSgenome对象进行搜索。
librara(BSgenome)
bs <- available.genomes()
但是拟南芥的参考基因组好像一直没有更新，最新的也是TAIR9，我表示很尴尬。
require(stringr)
str_subset(bs, "TAIR")
#[1] "BSgenome.Athaliana.TAIR.04232008" #"BSgenome.Athaliana.TAIR.TAIR9"
然后我觉得TAIR10是2012年左右发布的，所以肯定有人提问了。于是我熟练的打开搜索引擎去找相关信息，于是我发现了如下内容
TAIR9 and TAIR10 correspond to the same genome assembly so there is no need for a BSgenome pkg for TAIR10 :-)
TAIR9的染色体命名是ChX, 也就是NCBI的风格，注意这一点。
于是问题就这样轻松解决了，下载TAIR9就行了。
标记： 我需要下载TAIR9和TAIR10，比较一下，有结果放到这里。
BiocInstaller::biocLite("BSgenome.Athaliana.TAIR.TAIR9")
函数有以下几种
# 知由来
sourceUrl(Athaliana)
# 看的不是基因组大小
length(Athaliana)
# 序列的命名
seqnames(Athaliana)
# 各序列长度
seqlengths(Athaliana)
# 提取第一条染色体
ch1 <- Athaliana[['1']]
根据GenmoicRanges信息提取序列
对GenomicRanges进行注释
我们可以使用variantAnnotation包中的locateVariants()或predictCoding()对GenmoicRanges数据进行注释。
locateVariants主要是判断GenomicRanges在基因组哪里，所以只需要TxDb数据库就行了
predictCodings要根据突变位点和位置判断位点是无义突变还是同义突变，氨基酸又是如何变化。所以还需要BSgenome数据库. 同时如果提供的是GenomicRanges, 你还需要额外提供突变后的碱基信息，DNAStringSet格式。
step1: 加载数据库
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
step2: （可选）读取VCF文件
fl <- system.file("extdata", "gl_chr1.vcf",
package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
step3 : 提供GenmoicRanges， 或直接用vcf进行找到位点所在基因组位置
region <- IntergenicVariants(upstream=70000, downstream=70000)
loc_int <- locateVariants(vcf, txdb, region)
# 或者用genomicRanges
gr <- rowRanges(vcf)
loc_int <- locateVariants(gr, txdb, region)
mcols(loc_int)[c("LOCATION", "PRECEDEID", "FOLLOWID")]
step4: （可选）： 预测碱基突变的影响
library(BSgenome.Hsapiens.UCSC.hg19)
coding <- predictCoding(vcf, txdb, seqSource=Hsapiens)
如果提供的是GenomicRange,那么你还得需要一个varAllele
varallele <- unlist(mcols(vcf)$ALT)
coding <- predictCoding(gr, txdb, seqSource=Hsapiens,
varAllele=varallele)
PS: VariantAnnotation包还提供了一些列函数，如readVcf()，header(), info(), geno(), rowRanges() 用于从加载VCF文件，然后提供特定的数据。类似bcftools
练习：寻找AGO1，然后对其进行GO注释，KEGG注释等
首先加载拟南芥的物种数据库， 因为存放了GENENAME等信息。
tair.org <- ah[['AH53758']]
然后提GENENAME和GENEID
columns(tair.org)
keytypes(tair.org)
genename <- select(tair.org, keys=keys(tair.org,keytype = "TAIR"), columns=c("TAIR","GENENAME"), keytype="TAIR")
然后通过模式匹配，找到含有AGO1的那一行。但是由于不是每个基因都有基因俗名，NA在模式匹配的时候会返回NA，而不是FALSE，所以需要先处理掉NA。要么是替换成空字符串，要么直接删除改行。
genename <- na.omit(genename)
筛选AGO1：
require(stringr)
ago1 <- genename[str_detect(genename$GENENAME, "AGO1"),]
通过肉眼检查，发现“AT1G48410“是要找到，和直接在TAIR搜索一致。
对该基因进行GO和KEGG注释
ago1 <- select(tair.org, keys="AT1G48410", columns = c("TAIR","PATH"), keytype = "TAIR" )
dim(ago1)
[1] 61 5
发现GO注释挺多，但是KEGG没有注释。想想也是吧，毕竟AGO1是参加miRNA加工的，不算是代谢通路的。
不过为了折腾，我觉得找到有KEGG注释的基因，然后去TAIR上看下。
kegg <- select(tair.org, keys = keys(tair.org,"PATH"), columns=c("TAIR","PATH"), keytype = "PATH")
也是一种一对多的关系，同一条PATH中必然涉及到多个基因，我随便找到了一个基因AT1G07420在TAIR搜索。
最终发现一个页面能搜到path，一个页面不能搜索到path。但是都没有KEGG的说明，我很好奇。
