應用BiomaRt獲取實驗數據與註釋信息

應用BiomaRt獲取實驗數據與註釋信息

獲取數據是生物信息分析的第一步,經常使用的一個示例數據是Bioconductor 中的慢性淋巴細胞白血病(Chronic Lymphocytic Leukemia, CLL)實驗數據包,其安裝辦法是在啟動R主程序並運行:

#安裝CLL數據包。

BiocManager::install("CLL")

#載入CLL數據包。

library(CLL)

#載入數據(庫文件中附帶的示例數據)。

data(CLLbatch)

#查看數據內容與結構。

phenoData(CLLbatch)
應用BiomaRt獲取實驗數據與註釋信息

從上圖,我們可以看到CLL數據集一共有24個樣本,在實際應用中, 需要讀取所有樣本對應的CEL文件。

把一些實驗數據製作成包,便於用戶下載和使用,也便於開發者開發和測試算法。但是對於比較大的數據或者經常動態更新的數據或註釋,包就顯得力不從心。

一個簡單的解決辦法:

Bioconductor只需要提供一個包就能實現Bioconductor與所有在線生物數據庫的接口,用戶就可以根據需要通過操作這個包獲取數據和註釋信息。

這樣大大簡化了數據和註釋信息的獲取,但是不同的數據庫對外提供數據的標準和方式都不同,如何統一生物信息數據庫的數據訪問就成為了一個關鍵問題。

因此,有必要介紹一下由歐洲生物信息研究所(EBI)和冷泉港實驗室(ESHL)共同開發的BioMart數據管理系統。BioMart系統可以管理任意格式的數據庫,也可以按照不同的需求安裝不同的查詢工具和界面,由於其內部採用關係型的數據組織模式,所以更易於進行復雜的數據挖掘研究,讀者可以通過官方網站(htp://www.biomart.org) 進一步瞭解BioMart。

當前各主要生物信息數據庫(如EBI維護的Ensembl數據庫)都提供了基於BioMart管理系統的批量數據訪問服務,這些數據庫也可以統稱為BioMart數據庫。

下載數據有兩種方式:

  • ①登陸官方網站htp//ww.ensembl.org/下載
  • ②R語言程序下載

我們今天學習的是R語言程序下載數據。

biomaRt包不僅可以輕鬆獲取數據庫中的數據和註釋,還能夠獲取相關數據或註釋在不同數據庫間的關聯信息,為生物數據和註釋的獲取提供了極大的便利,這是一般的 R擴展包無法實現的。

#安裝biomaRt包。

BiocManager::install("biomaRt")

#載入biomaRt包。

library(biomaRt)

#獲取當前可用的數據源,一個數據源叫做一個mart。

marts 

#只查看前幾個。

head(marts)
應用BiomaRt獲取實驗數據與註釋信息

#使用ensembl數據源,如果知道用這個,前面沒必要查看所有數據源。

ensembl_mart 

#獲取ensembl_mart 中可用數據集。

datasets 

#查看前10個。

datasets[1:10, ]
應用BiomaRt獲取實驗數據與註釋信息

#使用豬基因組數據集。

dataset_pig 

#獲取dataset_pig數據集上可用的篩選器。

filters 

#只查看前幾個,後面沒用到任何篩選器。

head(filters)
應用BiomaRt獲取實驗數據與註釋信息

attributes 

#只查看前幾個

head(attributes)
 
應用BiomaRt獲取實驗數據與註釋信息

#從dataset_pig數據集中提取ensembl_transcript_id和description信息。

idlist 

#從dataset_pig數據集中根據ensembl_transcript_id提取序列。

seqs=getSequence(id=idlist["ensembl_transcript_id"],
type="ensembl_transcript_id",
seqType="3utr", mart = dataset_pig)

#去除沒有序列內容的數據記錄。

seqs = seqs[!seqs[,1]== "Sequence unavailable", ]

#去除沒有UTR註釋的數據記錄。

seqs = seqs[!seqs[,1]=="No UTR is annotated for this transcript", ]

#提取序列的內容。

x=seqs[ ,1]

#提取序列的ID。

names(x)=seqs[ ,2]

#結果存入文件"UTR3seqs-1.fa",格式為fasta。

writeXStringSet(DNAStringSet(x, use.names=TRUE),"UTR3seqs-1.fa")

#同時提取對應3'UTR序列的cDNA序列。

cDNAseqs = getSequence(id=idlist"ensembl_transcript_id"],type="ensembl_transcript_id",seqType="cdna", mart = dataset_pig)
x=cDNAseqs[ ,1]
names(x)=cDNAseqs[ ,2]

#結果存入文件"UTR3seqs-1.fa",格式為fasta。

writeXStringSet(DNAStringSet(x, use.names=TRUE), "transcriptom.fasta")

讀者運行並對比兩次下載的結果會發現,編程下載的數據沒有添加註釋信息,其實該信息已經包含在變量"idlist" 中,只需稍稍改動一點程序,即可以得到含有註釋信息的數據。

使用biomaRt 包應特別注意兩點:

  • 其一是該包很多函數運行速度嚴重依賴網絡速度,有些在線數據庫對每次下載的數據量的大小作出了限制,在批處理數據的過程中,還可能發生不可意料的錯誤
  • 其二是biomaRt包中可用的數據庫和數據集合的版本與該數據庫官方網站的版本不保證同步更新,如biomaRt包中Ensembl 數據庫版本已經更新到"ENSEMBL GENES 73",而官方網站的BioMart數據庫則可能還是"ENSEMBL GENES 70",因此使用數據時,一定要記錄好數據庫版本,否則將來無法重複計算結果。


分享到:


相關文章: