用Bioconductor对基因组注释
這一次，我們來聊聊基因組注釋。首先問自己一個問題，為什么要進行基因注釋。
就我目前而言，它用來解決如下問題：
如果一個基因沒有注釋信息，那么他就只是一段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.2bitAH54630 | Xiphophorus_maculatus.Xipmac4.4.2.dna_sm.toplevel.2bitAH54631 | 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.sqliteAH52246 | TxDb.Athaliana.BioMart.plantsmart25.sqliteAH52247 | TxDb.Athaliana.BioMart.plantsmart28.sqliteAH53758 | 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 columnsSYMBOL 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 columnsGO 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包
由于隨機取樣，很有可能找不到任何富集，可也是有有可能的。
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，比較一下，有結果放到這里。
函數有以下幾種
# 知由來 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.knownGenestep2: （可選）讀取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对基因组注释的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 京东618：六年历程步步为营，京东商城的
- 下一篇： 青龙峡蹦极攻略