|
biomaRt這個(gè)包很久以前我就給它寫過教程(點(diǎn)擊閱讀)，但是排版不好，可讀性很差，所以我用R Markdown重新來一個(gè)。 當(dāng)然了，它本身有官方的英文版教程(點(diǎn)擊閱讀)，我在翻譯的基礎(chǔ)上面，加入了自己的理解， 下面是正文： biomaRt是一個(gè)超級網(wǎng)絡(luò)資源庫，里面的信息非常之多，就是網(wǎng)頁版的biomaRt的R語言接口。谷歌搜索 the biomart user’s guide filetype:pdf 這個(gè)關(guān)鍵詞，就看到關(guān)于這個(gè)包的詳細(xì)介紹以及例子，我這里簡單總結(jié)一下它的用法。
包的安裝Bioconductor系列包的安裝方法都一樣 source("http:///biocLite.R")
biocLite(“biomaRt”)
install.packages('DT')
選擇數(shù)據(jù)庫然后我們就可以加載這個(gè)包來看看它的一些性質(zhì)啦，大家也可以根據(jù)這個(gè)包的說明書里面的代碼一行行代碼敲進(jìn)去看看效果，這樣學(xué)習(xí)最快也最容易記住。
library("biomaRt")
library(DT)
listMarts() #看看有多少數(shù)據(jù)庫資源
## biomart version
## 1 ENSEMBL_MART_ENSEMBL Ensembl Genes 91
## 2 ENSEMBL_MART_MOUSE Mouse strains 91
## 3 ENSEMBL_MART_SNP Ensembl Variation 91
## 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 91
ensembl=useMart("ENSEMBL_MART_ENSEMBL")
all_datasets <- listDatasets(ensembl) #看看選擇的數(shù)據(jù)庫里面有多少數(shù)據(jù)表，這個(gè)跟物種相關(guān)
library(DT)
datatable(all_datasets, options = list(
searching = FALSE,
pageLength = 5,
lengthMenu = c(5, 10, 15, 20)
))
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
# 這是一步法選擇人類的ensembl數(shù)據(jù)庫代碼.
ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")
三個(gè)主要函數(shù)getBM,getSequence,getLDS我們選擇好了數(shù)據(jù)庫就要開始干活啦，這個(gè)數(shù)據(jù)庫的檢索主要是三個(gè)函數(shù)getBM,getSequence,getLDS, 其中g(shù)etBM這個(gè)函數(shù)可以部分用select語句替代。
首先我們講講getBM函數(shù)，它就四個(gè)參數(shù)。 filter來控制根據(jù)什么東西來過濾，可是不同數(shù)據(jù)庫的ID，也可以是染色體定位系統(tǒng)坐標(biāo) Attributes來控制我們想獲得什么，一般是不同數(shù)據(jù)庫的ID Values是我們用來檢索的關(guān)鍵詞向量 Mart是我們前面選擇好的數(shù)據(jù)庫，一般都是ensembl = useMart(“ensembl”,dataset=“hsapiens_gene_ensembl”)
getBM函數(shù)唯一的用處，就是做各種ID轉(zhuǎn)換。 可以查看filter和attribute有哪些東西。 filters = listFilters(ensembl)
#filters[1:5,]
library(DT)
datatable(filters, options = list(
searching = FALSE,
pageLength = 5,
lengthMenu = c(5, 10, 15, 20)
))
attributes = listAttributes(ensembl)
#attributes[1:5,]
datatable(attributes, options = list(
searching = FALSE,
pageLength = 5,
lengthMenu = c(5, 10, 15, 20)
))
getSequence函數(shù)大同小異，建議大家活用help函數(shù)去查看每個(gè)函數(shù)的幫助文檔。 簡單講幾個(gè)例子咯： Ps：這些都是在線注釋，所以都是要網(wǎng)絡(luò)的，網(wǎng)速慢的會(huì)非?？?/strong> 幾個(gè)實(shí)用的例子一．對幾個(gè)芯片探針的ID號，注釋它所捕獲的基因的entrezID# ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
affyids=c("202763_at","209310_s_at","207500_at")
getBM(attributes=c('affy_hg_u133_plus_2', 'entrezgene'),
filters = 'affy_hg_u133_plus_2',
values = affyids,
mart = ensembl
)
可以看到結(jié)果里面已經(jīng)成功的把a(bǔ)ffymetrix的芯片探針I(yè)D，轉(zhuǎn)為了對應(yīng)的基因的entrez ID
二．對剛才的那三個(gè)探針I(yè)D號進(jìn)行多個(gè)內(nèi)容注釋，每個(gè)探針都對應(yīng)著基因名已經(jīng)染色體及起始終止坐標(biāo)。affyids=c("202763_at","209310_s_at","207500_at")
getBM(attributes=c('affy_hg_u133_plus_2', 'hgnc_symbol',
'chromosome_name','start_position','end_position', 'band'),
filters = 'affy_hg_u133_plus_2',
values = affyids,
mart = ensembl
)
三．對給定的基因ID號進(jìn)行GO注釋# library("biomaRt")
# ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
entrez=c("673","837")
goids = getBM(attributes = c('entrezgene', 'go_id'),
filters = 'entrezgene',
values = entrez,
mart = ensembl)
head(goids)
四．通過染色體及起始終止坐標(biāo)來挑選基因getBM(c('affy_hg_u133_plus_2','ensembl_gene_id'),
filters = c('chromosome_name','start','end'),
values=list(16,1100000,1250000),
mart=ensembl
)
五．對特定的GO ID號來查詢該go通路上面的基因是哪些。getBM(c('entrezgene','hgnc_symbol'), filters='go', values='GO:0004707', mart=ensembl)
六．根據(jù)refseq數(shù)據(jù)庫的NM系列ID號來獲取信息refseqids = c("NM_005359","NM_000546")
ipro = getBM(attributes=c("refseq_dna","interpro","interpro_description"), filters="refseq_dna",values=refseqids, mart=ensembl)
這個(gè)例子的代碼有錯(cuò)誤，因?yàn)閞efseq的信息沒有 refseq_dna Error in getBM(attributes = c(“refseq_dna”, “interpro”, “interpro_description”), : Invalid attribute(s): refseq_dna Please use the function 'listAttributes’ to get valid attribute names 我簡單檢查了一下，發(fā)現(xiàn)需要更正為 refseq_mrna
tmp=listAttributes(ensembl)
grep("refseq",tmp[,1])
## [1] 82 83 84 85 86 87
# [1] 81 82 83 84 85 86
tmp[grep("refseq",tmp[,1]),]
## name description page
## 82 refseq_mrna RefSeq mRNA ID feature_page
## 83 refseq_mrna_predicted RefSeq mRNA predicted ID feature_page
## 84 refseq_ncrna RefSeq ncRNA ID feature_page
## 85 refseq_ncrna_predicted RefSeq ncRNA predicted ID feature_page
## 86 refseq_peptide RefSeq peptide ID feature_page
## 87 refseq_peptide_predicted RefSeq peptide predicted ID feature_page
很明顯，是沒有 refseq_dna 的，只有 refseq_mrna 然后我用了新的代碼 refseqids = c("NM_005359","NM_000546")
ipro = getBM(attributes=c("refseq_mrna","interpro","interpro_description"), filters="refseq_mrna",values=refseqids, mart=ensembl)
ipro
## refseq_mrna interpro
## 1 NM_000546 IPR002117
## 2 NM_000546 IPR008967
## 3 NM_000546 IPR010991
## 4 NM_000546 IPR011615
## 5 NM_000546 IPR012346
## 6 NM_000546 IPR013872
## 7 NM_005359 IPR001132
## 8 NM_005359 IPR003619
## 9 NM_005359 IPR008984
## 10 NM_005359 IPR013019
## 11 NM_005359 IPR013790
## 12 NM_005359 IPR017855
## interpro_description
## 1 p53 tumour suppressor family
## 2 p53-like transcription factor, DNA-binding
## 3 p53, tetramerisation domain
## 4 p53, DNA-binding domain
## 5 p53/RUNT-type transcription factor, DNA-binding domain
## 6 p53 transactivation domain
## 7 SMAD domain, Dwarfin-type
## 8 MAD homology 1, Dwarfin-type
## 9 SMAD/FHA domain
## 10 MAD homology, MH1
## 11 Dwarfin
## 12 SMAD domain-like
就成功的完成了轉(zhuǎn)換 七．根據(jù)基因的entrez ID號來挑選該基因的指定上下游區(qū)域信息或者蛋白序列entrez=c("673","7157","837")
getSequence(id = entrez, type="entrezgene",seqType="coding_gene_flank",upstream=100, mart=ensembl)
## coding_gene_flank
## 1 CCTCCGCCTCCGCCTCCGCCTCCGCCTCCCCCAGCTCTCCGCCTCCCTTCCCCCTCCCCGCCCGACAGCGGCCGCTCGGGCCCCGGCTCTCGGTTATAAG
## 2 CACGTTTCCGCCCTTTGCAATAAGGAAATACATAGTTTACTTTCATTTTTGACTCTGAGGCTCTTTCCAACGCTGTAAAAAAGGACAGAGGCTGTTCCCT
## 3 TCCTTCTCTGCAGGCCCAGGTGACCCAGGGTTGGAAGTGTCTCATGCTGGATCCCCACTTTTCCTCTTGCAGCAGCCAGACTGCCTTCCGGGTCACTGCC
## entrezgene
## 1 673
## 2 837
## 3 7157
這樣就得到了三條序列，是給定基因的上游100bp序列 coding_gene_flank entrezgene 其實(shí)這個(gè)getSequence函數(shù)還有非常多的用法，當(dāng)然主要的變化在其readme上面可以看到，主要是seqType可以有多個(gè)選擇。 utr5 = getSequence(chromosome=3, start=185514033, end=185535839,
type="entrezgene",seqType="5utr", mart=ensembl)
utr5
## 5utr
## 1 Sequence unavailable
## 2 AGTCCCTAGGGAACTTCCTGTTGTCACCACACCTCTGAGTCGTCTGAGCTCACTGTGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 3 TGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 4 ATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## entrezgene
## 1 200879
## 2 200879
## 3 200879
## 4 200879
#這樣就查詢到了 chromosome=3, start=185514033, end=185535839,這個(gè)定位里面的所有utr5信息。
# 返回的urt5這個(gè)對象是一個(gè)data.frame，可以用exportFASTA()函數(shù)導(dǎo)出成fasta文件。
protein = getSequence(id=c(100, 5728),type="entrezgene",
seqType="peptide", mart=ensembl)
protein
## peptide
## 1 MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIARL*
## 2 Sequence unavailable
## 3 ALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVS*
## 4 MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 5 MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEAQK*
## 6 Sequence unavailable
## 7 MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 8 MTAIIKEIVSRNKRRYQEDGFDLDLTYIYPNIIAMGFPAERLEGVYRNNIDDVVRFLDSKHKNHYKIYNLCAERHYDTAKFNCRVAQYPFEDHNPPQLELIKPFCEDLDQWLSEDDNHVAAIHCKAGKGRTGVMICAYLLHRGKFLKAQEALDFYGEVRTRDKKGVTIPSQRRYVYYYSYLLKNHLDYRPVALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVKLYFTKTVEEPSNPEASSSTSVTPDVSDNEPDHYRYSDTTDSDPENEPFDEDQHTQITKV*
## entrezgene
## 1 100
## 2 100
## 3 5728
## 4 100
## 5 100
## 6 5728
## 7 100
## 8 5728
#這樣就查到了這兩個(gè)基因轉(zhuǎn)錄本的所有蛋白序列，100這個(gè)基因有5條序列，5728這個(gè)基因有四條序列。
getSequence函數(shù)返回的對象是一個(gè)data.frame，可以用exportFASTA()函數(shù)導(dǎo)出成fasta文件。
八，選擇其它數(shù)據(jù)庫來進(jìn)行查詢，比如snp數(shù)據(jù)庫當(dāng)然還有一些數(shù)據(jù)庫的小技巧，第一個(gè)是參數(shù) archive = TRUE，設(shè)置只用能獲取的數(shù)據(jù)庫 然后是設(shè)置特定選取hg19對應(yīng)的信息。 ensembl = useMart("ensembl_mart_46", dataset="hsapiens_gene_ensembl", archive = TRUE)
listMarts(host='may2009.archive.ensembl.org')
ensembl54=useMart(host='may2009.archive.ensembl.org', biomart='ENSEMBL_MART_ENSEMBL')
ensembl54=useMart(host='may2009.archive.ensembl.org', biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl')
還可以選取其它物種，比如秀麗隱桿線蟲 下面這個(gè)代碼是過時(shí)了的，你們需要自己去看說明書找到現(xiàn)在最新的代碼wormbase=useMart("WS220",dataset="wormbase_gene")
listFilters(wormbase)
listAttributes(wormbase)
getBM(attributes = c("public_name","rnai","rnai_phenotype_phenotype_label"),
filters="gene_name", values=c("unc-26","his-33"),
mart=wormbase)
寫在最后最后我簡單提一下select函數(shù)是如何部分替代getBM函數(shù)的，因?yàn)閎iomaRt是在線數(shù)據(jù)庫，本來只能用它自己的getBM系列函數(shù)，但是為了對接其它bioconductor系列包，也可以用select函數(shù)來操作這個(gè)在線數(shù)據(jù)庫。
mart<-useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")
head(keytypes(mart), n=3)
## [1] "affy_hc_g110" "affy_hg_focus" "affy_hg_u133_plus_2"
length(keytypes(mart))
## [1] 334
length(columns(mart))
## [1] 1918
head(columns(mart), n=3)
## [1] "3_utr_end" "3_utr_end" "3_utr_start"
k = keys(mart, keytype="chromosome_name")
head(k, n=3)
## [1] "1" "2" "3"
k = keys(mart, keytype="chromosome_name", pattern="LRG")
head(k, n=3)
## character(0)
affy=c("202763_at","209310_s_at","207500_at")
select(mart, keys=affy, columns=c('affy_hg_u133_plus_2','entrezgene'), keytype='affy_hg_u133_plus_2')
## affy_hg_u133_plus_2 entrezgene
## 1 202763_at 836
## 2 209310_s_at 837
## 3 207500_at 838
而且這個(gè)包更新也比較頻繁，所以大家如果看到比較老舊的教程，可能無法模仿，或者多年以后，你看到我這個(gè)教程，也會(huì)發(fā)現(xiàn)，重復(fù)不出來咯。
|
