可能是個生物資訊學資料超市吧
阿新 • • 發佈：2022-05-02
biomaRt這個包很久以前我就給它寫過教程(點選閱讀)，但是排版不好，可讀性很差，所以我用R Markdown重新來一個。 當然了，它本身有官方的英文版教程(點選閱讀)，我在翻譯的基礎上面，加入了自己的理解， 下面是正文：
biomaRt是一個超級網路資源庫，裡面的資訊非常之多，就是網頁版的biomaRt的R語言介面。谷歌搜尋 the biomart user’s guide filetype:pdf 這個關鍵詞，就看到關於這個包的詳細介紹以及例子，我這裡簡單總結一下它的用法。
包的安裝
Bioconductor系列包的安裝方法都一樣
source("http://bioconductor.org/biocLite.R") biocLite(“biomaRt”) install.packages('DT')
選擇資料庫
然後我們就可以載入這個包來看看它的一些性質啦，大家也可以根據這個包的說明書裡面的程式碼一行行程式碼敲進去看看效果，這樣學習最快也最容易記住。
library("biomaRt")
library(DT)
listMarts() #看看有多少資料庫資源## 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) #看看選擇的資料庫裡面有多少資料表，這個跟物種相關
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資料庫程式碼. ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")
三個主要函式getBM,getSequence,getLDS
我們選擇好了資料庫就要開始幹活啦，這個資料庫的檢索主要是三個函式getBM,getSequence,getLDS, 其中getBM這個函式可以部分用select語句替代。
首先我們講講getBM函式，它就四個引數。
- filter來控制根據什麼東西來過濾，可是不同資料庫的ID，也可以是染色體定位系統座標
- Attributes來控制我們想獲得什麼，一般是不同資料庫的ID
- Values是我們用來檢索的關鍵詞向量
- Mart是我們前面選擇好的資料庫，一般都是ensembl = useMart(“ensembl”,dataset=“hsapiens_gene_ensembl”)
getBM函式唯一的用處，就是做各種ID轉換。
可以檢視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函式大同小異，建議大家活用help函式去檢視每個函式的幫助文件。
簡單講幾個例子咯： Ps：這些都是線上註釋，所以都是要網路的，網速慢的會非常坑
幾個實用的例子
一．對幾個晶片探針的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
)可以看到結果裡面已經成功的把affymetrix的晶片探針ID，轉為了對應的基因的entrez ID
二．對剛才的那三個探針ID號進行多個內容註釋，每個探針都對應著基因名已經染色體及起始終止座標。
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號進行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)四．通過染色體及起始終止座標來挑選基因
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)六．根據refseq資料庫的NM系列ID號來獲取資訊
refseqids = c("NM_005359","NM_000546")
ipro = getBM(attributes=c("refseq_dna","interpro","interpro_description"), filters="refseq_dna",values=refseqids, mart=ensembl)這個例子的程式碼有錯誤，因為refseq的資訊沒有 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 我簡單檢查了一下，發現需要更正為 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就成功的完成了轉換
七．根據基因的entrez ID號來挑選該基因的指定上下游區域資訊或者蛋白序列
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
其實這個getSequence函式還有非常多的用法，當然主要的變化在其readme上面可以看到，主要是seqType可以有多個選擇。
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,這個定位裡面的所有utr5資訊。
# 返回的urt5這個物件是一個data.frame，可以用exportFASTA()函式匯出成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#這樣就查到了這兩個基因轉錄本的所有蛋白序列，100這個基因有5條序列，5728這個基因有四條序列。getSequence函式返回的物件是一個data.frame，可以用exportFASTA()函式匯出成fasta檔案。
八，選擇其它資料庫來進行查詢，比如snp資料庫
當然還有一些資料庫的小技巧，第一個是引數 archive = TRUE，設定只用能獲取的資料庫
然後是設定特定選取hg19對應的資訊。
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')還可以選取其它物種，比如秀麗隱杆線蟲
下面這個程式碼是過時了的，你們需要自己去看說明書找到現在最新的程式碼
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函式是如何部分替代getBM函式的，因為biomaRt是線上資料庫，本來只能用它自己的getBM系列函式，但是為了對接其它bioconductor系列包，也可以用select函式來操作這個線上資料庫。
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] 334length(columns(mart))## [1] 1918head(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而且這個包更新也比較頻繁，所以大家如果看到比較老舊的教程，可能無法模仿，或者多年以後，你看到我這個教程，也會發現，重複不出來咯。