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

96
徐洲更hoptop 22d8d123 271c 4d80 9c59 6990844a9e37
2.0 2017.06.18 12:53* 字数 2902

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

  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")</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的说明,我很好奇。

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


知识星球
R语言学习笔记