這一次，我們來聊聊基因組注釋。首先問自己一個問題，為什麼要進行基因注釋。
就我目前而言，它用來解決如下問題：
- 在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;
- 使用
[]
[[]]
探索注釋資料庫
找到和下載下傳注釋資料庫隻是第一步，學會如何使用這些資料庫更加重要。
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的說明，我很好奇。
順便放一下自己的知識星球，如果你覺得我對你有幫助的話。
知識星球
![【MySQL資料庫】資料庫索引事務1.索引2.事務[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)