【Bioconductor系列】如何用Bioconductor对基因组注释

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

1. 在mapping-by-sequencing的时候，我找到了一些可能的突变位点，我需要知道这些突变分别是那些基因发生突变，这些突变基因有哪些功能？
2. 差异表达分析之后会得到许多的基因，这些基因有什么样的特征？如果要进行基因富集分析，不可避免就需要知道他们的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


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

• 它的数据库版本是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(&＃8220;Homo sapiens&＃8221;, &＃8220;hg19&＃8221;, &＃8220;GTF&＃8221;))

比如说我想查找拟南芥相关的注释数据库，就可以去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就是&＃8221;AH53758&＃8243;，然后根据这个ID可以进行下载。当然下载方式已经出现过了，

retrieve record with &＃8216;object[[&＃8220;AH53758&＃8221;]]&＃8217;

ath <- ah[['AH53758']]


bioconductor除了AnnotationHub能用来查找生物数据，还有一个库叫做biomaRt，可以用来查找biomart中的数据。不过目前biomart网站正在进行服务器的数据迁移，就不在这里演示。

小结：

• AnnotationHub是生物数据的中转站，方面我们搜索目标数据，另一个相似包是biomaRt;
• 我们通过query,subset等方法（图形界面则是display），逐步从AnnotationHub的metadata筛选到所需数据的ID;
• 使用[]是查看目标数据的metadata, 使用[[]]用于下载数据;

探索注释数据库

找到和下载注释数据库只是第一步，学会如何使用这些数据库更加重要。

AnnotationHub对象的通用方法

之前下载完数据后，在R里面被我指向到了&＃8217;ath&＃8217;，那么我们先简单了解一下这个&＃8217;ath&＃8217;

直接输入对象名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类呢？

mode(ath)
[1] "S4"


用class看下它具体继承什么类（面向对象编程的概念）

class(ath)
[1] "OrgDb"
attr(,"package")
[1] "AnnotationDbi"


好了，我们继续调查什么是&＃8221;AnnotationDbi&＃8221;,了解到他主要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中必然涉及到多个基因，我随便找到了一个基因AT1G07420在TAIR搜索。

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

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

知识星球