用Bioconductor对基因组注释-阿里云开发者社区

这一次，我们来聊聊基因组注释。首先问自己一个问题，为什么要进行基因注释。
就我目前而言，它用来解决如下问题：

在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")</pre>

使用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

上述结果告诉了我们以下信息：

它的数据库版本是2017-4-25，目前有39213条记录
你可以用ah$dataprovider的方式查看数据来源，还可以看有哪些物种和数据类型可以用。
你可以用mcols看更多的元信息。
获取数据的方式是object[["AH2"]]

根据这些知识点，我们就可以问第一个问题：

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[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 <- 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                :<externalptr> 
  .. ..@ dbname             : chr "D:\\xuzho\\Documents\\AppData\\.AnnotationHub\\60496"
  .. ..@ loadable.extensions: logi TRUE
  .. ..@ flags              : int 1
  .. ..@ vfs                : chr ""
  .. ..@ ref                :<environment: 0x000000002195a428> 
 $ packageName: chr(0) 
 and 14 methods.

用mode看下它的数据模式(Get or set the type or storage mode of an object.)，发现它是一个S4类。所有bioconductor的包都是S4类，然而什么是S4类呢？

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)

由于随机取样，很有可能找不到任何富集，可也是有有可能的。

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()

但是拟南芥的参考基因组好像一直没有更新，我表示很尴尬。

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中必然涉及到多个基因，我随便找到了一个基因<font size="+1">AT1G07420</font>在TAIR搜索。

最终发现一个页面能搜到path，一个页面不能搜索到path。但是都没有KEGG的说明，我很好奇。

顺便放一下自己的知识星球，如果你觉得我对你有帮助的话。

知识星球

用Bioconductor对基因组注释