Contents

1 Introduction

Accessing the data available in Ensembl is by far most frequent use of the biomaRt package. With that in mind biomaRt provides a number of functions that are tailored to work specifically with the BioMart instances provided by Ensembl. This vignette details this Ensembl specific functionailty and provides a number of example usecases that can be used as the basis for specifying your own queries.

2 Selecting an Ensembl BioMart database and dataset

Every analysis with biomaRt starts with selecting a BioMart database to use.

library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

If this your first time using biomaRt , you might wonder how to find the two arguments we supplied to the useEnsembl() command. This is a two step process, but once you know the setting you need you can use the version shown above as a single command. These initial steps are outlined below.

2.1 Identifying the database you need

The first step is to find the names of the BioMart services Ensembl is currently providing. We can do this using the function listEnsembl(), which will display all available Ensembl BioMart web services. The first column gives us the name we should provide to the biomart argument in useEnsembl(), and the second gives a more comprehensive title for the dataset along with the Ensembl version.

listEnsembl()
##         biomart                version
## 1         genes      Ensembl Genes 101
## 2 mouse_strains      Mouse strains 101
## 3          snps  Ensembl Variation 101
## 4    regulation Ensembl Regulation 101

The useEnsembl() function can now be used to connect to the desired BioMart database. The biomart argument should be given a valid name from the output of listEnsembl(). In the next example we will select the main Ensembl mart, which provides access to gene annotation information.

ensembl <- useEnsembl(biomart = "genes")

If we print the current ensembl object, we can see that the ENSEMBL_MART_ENSEMBL database 1 this is how Ensembl name the database on their server has been selected, but that no dataset has been choosen.

ensembl
## Object of class 'Mart':
##   Using the ENSEMBL_MART_ENSEMBL BioMart database
##   No dataset selected.

2.2 Choosing a dataset

BioMart databases can contain several datasets. For example, within the Ensembl genes mart every species is a different dataset. In the next step we look at which datasets are available in the selected BioMart by using the function listDatasets(). Note: here we use the function head() to display only the first 5 entries as the complete list has 203 entries.

datasets <- listDatasets(ensembl)
head(datasets)
##                      dataset                                  description                version
## 1   acalliptera_gene_ensembl             Eastern happy genes (fAstCal1.2)             fAstCal1.2
## 2 acarolinensis_gene_ensembl               Anole lizard genes (AnoCar2.0)              AnoCar2.0
## 3  acchrysaetos_gene_ensembl              Golden eagle genes (bAquChr1.2)             bAquChr1.2
## 4  acitrinellus_gene_ensembl               Midas cichlid genes (Midas_v5)               Midas_v5
## 5  amelanoleuca_gene_ensembl                        Panda genes (ailMel1)                ailMel1
## 6    amexicanus_gene_ensembl Mexican tetra genes (Astyanax_mexicanus-2.0) Astyanax_mexicanus-2.0

The listDatasets() function will return every available option, however this can be unwieldy when the list of results is long, involving much scrolling to find the entry you are interested in. biomaRt also provides the functions searchDatasets() which will try to find any entries matching a specific term or pattern. For example, if we want to find the details of any datasets in our ensembl mart that contain the term ‘hsapiens’ we could do the following:

searchDatasets(mart = ensembl, pattern = "hsapiens")
##                  dataset              description    version
## 80 hsapiens_gene_ensembl Human genes (GRCh38.p13) GRCh38.p13

To use a dataset we can update our Mart object using the function useDataset(). In the example below we choose to use the hsapiens dataset.

ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)

As mentioned previously, if the dataset one wants to use is known in advance i.e. you’ve gone through this process before, we can select a both the database and dataset in one step:

ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

2.3 Ensembl mirror sites

To improve performance Ensembl provides several mirrors of their site distributed around the globe. When you use the default settings for useEnsembl() your queries will be directed to your closest mirror geographically. In theory this should give you the best performance, however this is not always the case in practice. For example the if the nearest mirror is experiencing many queries from other users it may perform poorly for you. You can use the mirror argument to useEnsembl() to explicitly request a specific mirror.

ensembl <- useEnsembl(biomart = "ensembl", 
                   dataset = "hsapiens_gene_ensembl", 
                   mirror = "useast")

Values for the mirror argument are: useast, uswest, asia, and www.

2.4 Using Ensembl Genomes

Ensembl Genomes expands the effort to provide annotation from the vertebrate genomes provided by the main Ensembl project across taxonomic space, with separate BioMart interfaces for Protists, Plants, Metazoa and Fungi.

Note: Unfortunately there is no BioMart interface to the Ensembl Bacteria data. The number of bacterial genomes is in the tens of thousands and BioMart does not perform well when providing data on that scale.

You can use the functions listEnsemblGenomes() and useEnsemblGenomes() in similar fashion to the functions shown previously. For example first we can list the available Ensembl Genomes marts:

listEnsemblGenomes()
##               biomart                        version
## 1       protists_mart      Ensembl Protists Genes 48
## 2 protists_variations Ensembl Protists Variations 48
## 3          fungi_mart         Ensembl Fungi Genes 48
## 4    fungi_variations    Ensembl Fungi Variations 48
## 5        metazoa_mart       Ensembl Metazoa Genes 48
## 6  metazoa_variations  Ensembl Metazoa Variations 48
## 7         plants_mart        Ensembl Plants Genes 48
## 8   plants_variations   Ensembl Plants Variations 48

We can the select the Ensembl Plants database, and search for the dataset name for Arabidopsis.

ensembl_plants <- useEnsemblGenomes(biomart = "plants_mart")
searchDatasets(ensembl_plants, pattern = "Arabidopsis")
##             dataset                         description version
## 4  ahalleri_eg_gene Arabidopsis halleri genes (Ahal2.2) Ahal2.2
## 5   alyrata_eg_gene    Arabidopsis lyrata genes (v.1.0)   v.1.0
## 7 athaliana_eg_gene Arabidopsis thaliana genes (TAIR10)  TAIR10

We can then use this information to create our Mart object that will access the correct database and dataset.

ensembl_arabidopsis <- useEnsemblGenomes(biomart = "plants_mart", 
                                         dataset = "athaliana_eg_gene")

3 How to build a biomaRt query

The getBM() function has three arguments that need to be introduced: filters, attributes and values. Filters define a restriction on the query. For example you want to restrict the output to all genes located on the human X chromosome then the filter chromosome_name can be used with value ‘X’. The listFilters() function shows you all available filters in the selected dataset.

filters = listFilters(ensembl)
filters[1:5,]
##              name              description
## 1 chromosome_name Chromosome/scaffold name
## 2           start                    Start
## 3             end                      End
## 4      band_start               Band Start
## 5        band_end                 Band End

Attributes define the values we are interested in to retrieve. For example we want to retrieve the gene symbols or chromosomal coordinates. The listAttributes() function displays all available attributes in the selected dataset.

attributes = listAttributes(ensembl)
attributes[1:5,]
##                            name                  description         page
## 1               ensembl_gene_id               Gene stable ID feature_page
## 2       ensembl_gene_id_version       Gene stable ID version feature_page
## 3         ensembl_transcript_id         Transcript stable ID feature_page
## 4 ensembl_transcript_id_version Transcript stable ID version feature_page
## 5            ensembl_peptide_id            Protein stable ID feature_page

The getBM() function is the main query function in biomaRt. It has four main arguments:

Note: for some frequently used queries to Ensembl, wrapper functions are available: getGene() and getSequence(). These functions call the getBM() function with hard coded filter and attribute names.

Now that we selected a BioMart database and dataset, and know about attributes, filters, and the values for filters; we can build a biomaRt query. Let’s make an easy query for the following problem: We have a list of Affymetrix identifiers from the u133plus2 platform and we want to retrieve the corresponding EntrezGene identifiers using the Ensembl mappings.

The u133plus2 platform will be the filter for this query and as values for this filter we use our list of Affymetrix identifiers. As output (attributes) for the query we want to retrieve the EntrezGene and u133plus2 identifiers so we get a mapping of these two identifiers as a result. The exact names that we will have to use to specify the attributes and filters can be retrieved with the listAttributes() and listFilters() function respectively. Let’s now run the query:

affyids <- c("202763_at","209310_s_at","207500_at")
getBM(attributes = c('affy_hg_u133_plus_2', 'entrezgene_id'),
      filters = 'affy_hg_u133_plus_2',
      values = affyids, 
      mart = ensembl)
##   affy_hg_u133_plus_2 entrezgene_id
## 1         209310_s_at           837
## 2           207500_at           838
## 3           202763_at           836

3.1 Searching for filters and attributes

The functions listAttributes() and listFilters() will return every available option for their respective types, which can produce a very long output where it is hard to find the value you are interested in. biomaRt also provides the functions searchAttributes() and searchFilters() which will try to find any entries matching a specific term or pattern, in a similar fashion to searchDatasets() seen previously. You can use these functions to find available attributes and filters that you may be interested in. The example below returns the details for all attributes that contain the pattern ‘hgnc’.

searchAttributes(mart = ensembl, pattern = "hgnc")
##               name        description         page
## 63         hgnc_id            HGNC ID feature_page
## 64     hgnc_symbol        HGNC symbol feature_page
## 94 hgnc_trans_name Transcript name ID feature_page

For advanced use, note that the pattern argument takes a regular expression. This means you can create more complex queries if required. Imagine, for example, that we have the string ENST00000577249.1, which we know is an Ensembl ID of some kind, but we aren’t sure what the appropriate filter term is. The example shown next uses a pattern that will find all filters that contain the terms ‘ensembl’ and ‘id’. This allows us to reduced the list of filters to only those that might be appropriate for our example.

searchFilters(mart = ensembl, pattern = "ensembl.*id")
##                             name                                                    description
## 56               ensembl_gene_id                       Gene stable ID(s) [e.g. ENSG00000000003]
## 57       ensembl_gene_id_version       Gene stable ID(s) with version [e.g. ENSG00000000003.15]
## 58         ensembl_transcript_id                 Transcript stable ID(s) [e.g. ENST00000000233]
## 59 ensembl_transcript_id_version Transcript stable ID(s) with version [e.g. ENST00000000233.10]
## 60            ensembl_peptide_id                    Protein stable ID(s) [e.g. ENSP00000000233]
## 61    ensembl_peptide_id_version     Protein stable ID(s) with version [e.g. ENSP00000000233.5]
## 62               ensembl_exon_id                              Exon ID(s) [e.g. ENSE00000000003]

From this we can compare ENST00000577249.1 with the examples given in the description column, and see it is a Transcript ID with version. Thus the appropriate filter value to use with it is ensembl_transcript_id_version.

3.2 Using predefined filter values

Many filters have a predefined list of values that are known to be in the dataset they are associated with. An common example would be the names of chromosomes when searching a dataset at Ensembl. In this online interface to BioMart these available options are displayed as a list as shown in Figure 1.

The options available to the Chromosome/Scaffold field are limited to a pretermined list based on the values in this dataset.

Figure 1: The options available to the Chromosome/Scaffold field are limited to a pretermined list based on the values in this dataset

You can list this in an R session with the function listFilterValues(), passing it a mart object and the name of the filter. For example, to list the possible chromosome names you could run the following:

listFilterValues(mart = ensembl, filter = "chromosome_name")

It is also possible to search the list of available values via searchFilterValues(). In the two examples below, the first returns all chromosome names starting with “GL”, while the second will find any phenotype descriptions that contain the string “Crohn”.

searchFilterValues(mart = ensembl, filter = "chromosome_name", pattern = "^GL")
searchFilterValues(mart = ensembl, filter = "phenotype_description", pattern = "Crohn")

4 Using archived versions of Ensembl

It is possible to query archived versions of Ensembl through biomaRt, so you can maintain consistent annotation throughout the duration of a project.

biomaRt provides the function listEnsemblArchives() to view the available Ensembl archives. This function takes no arguments, and produces a table containing the name and version number of the available archives, the date they were first released, and the URL where they can be accessed.

listEnsemblArchives()
##              name     date                                url version current_release
## 1  Ensembl GRCh37 Feb 2014          http://grch37.ensembl.org  GRCh37                
## 2     Ensembl 101 Aug 2020 http://aug2020.archive.ensembl.org     101               *
## 3     Ensembl 100 Apr 2020 http://apr2020.archive.ensembl.org     100                
## 4      Ensembl 99 Jan 2020 http://jan2020.archive.ensembl.org      99                
## 5      Ensembl 98 Sep 2019 http://sep2019.archive.ensembl.org      98                
## 6      Ensembl 97 Jul 2019 http://jul2019.archive.ensembl.org      97                
## 7      Ensembl 96 Apr 2019 http://apr2019.archive.ensembl.org      96                
## 8      Ensembl 95 Jan 2019 http://jan2019.archive.ensembl.org      95                
## 9      Ensembl 94 Oct 2018 http://oct2018.archive.ensembl.org      94                
## 10     Ensembl 93 Jul 2018 http://jul2018.archive.ensembl.org      93                
## 11     Ensembl 92 Apr 2018 http://apr2018.archive.ensembl.org      92                
## 12     Ensembl 91 Dec 2017 http://dec2017.archive.ensembl.org      91                
## 13     Ensembl 90 Aug 2017 http://aug2017.archive.ensembl.org      90                
## 14     Ensembl 89 May 2017 http://may2017.archive.ensembl.org      89                
## 15     Ensembl 88 Mar 2017 http://mar2017.archive.ensembl.org      88                
## 16     Ensembl 87 Dec 2016 http://dec2016.archive.ensembl.org      87                
## 17     Ensembl 86 Oct 2016 http://oct2016.archive.ensembl.org      86                
## 18     Ensembl 85 Jul 2016 http://jul2016.archive.ensembl.org      85                
## 19     Ensembl 84 Mar 2016 http://mar2016.archive.ensembl.org      84                
## 20     Ensembl 83 Dec 2015 http://dec2015.archive.ensembl.org      83                
## 21     Ensembl 82 Sep 2015 http://sep2015.archive.ensembl.org      82                
## 22     Ensembl 77 Oct 2014 http://oct2014.archive.ensembl.org      77                
## 23     Ensembl 75 Feb 2014 http://feb2014.archive.ensembl.org      75                
## 24     Ensembl 67 May 2012 http://may2012.archive.ensembl.org      67                
## 25     Ensembl 54 May 2009 http://may2009.archive.ensembl.org      54

Alternatively, one can use the http://www.ensembl.org website to find an archived version. From the main page scroll down the bottom of the page, click on ‘view in Archive’ and select the archive you need.

You will notice that there is an archive URL even for the current release of Ensembl. It can be useful to use this if you wish to ensure that script you write now will return exactly the same results in the future. Using www.ensembl.org will always access the current release, and so the data retrieved may change over time as new releases come out.

Whichever method you use to find the URL of the archive you wish to query, copy the url and use that in the host argument as shown below to connect to the specified BioMart database. The example below shows how to query Ensembl 54.

listEnsembl(version = 54)
##            biomart              version
## 1            genes           Ensembl 54
## 2             snps Ensembl Variation 54
## 3             vega              Vega 35
## 4         REACTOME   Reactome(CSHL US) 
## 5 wormbase_current   WormBase (CSHL US)
## 6            pride       PRIDE (EBI UK)
ensembl54 <- useEnsembl(biomart = 'genes', 
                       dataset = 'hsapiens_gene_ensembl',
                       version = 54)

5 Result Caching

To save time and computing resources biomaRt will attempt to identify when you are re-running a query you have executed before. Each time a new query is run, the results will be saved to a cache on your computer. If a query is identified as having been run previously, rather than submitting the query to the server, the results will be loaded from the cache.

You can get some information on the size and location of the cache using the function biomartCacheInfo():

biomartCacheInfo()
## biomaRt cache
## - Location: /home/biocbuild/.cache/biomaRt
## - No. of files: 702
## - Total size: 88.8 Mb

The cache can be deleted using the command biomartCacheClear(). This will remove all cached files.

6 biomaRt helper functions

This section describes a set of biomaRt helper functions that can be used to export FASTA format sequences, retrieve values for certain filters and exploring the available filters and attributes in a more systematic manner.

6.1 exportFASTA

The data.frames obtained by the getSequence function can be exported to FASTA files using the exportFASTA() function. One has to specify the data.frame to export and the filename using the file argument.

6.2 Finding out more information on filters

6.2.1 filterType

Boolean filters need a value TRUE or FALSE in biomaRt. Setting the value TRUE will include all information that fulfill the filter requirement. Setting FALSE will exclude the information that fulfills the filter requirement and will return all values that don’t fulfill the filter. For most of the filters, their name indicates if the type is a boolean or not and they will usually start with “with”. However this is not a rule and to make sure you got the type right you can use the function filterType() to investigate the type of the filter you want to use.

filterType("with_affy_hg_u133_plus_2",ensembl)
## [1] "boolean_list"

6.2.2 filterOptions

Some filters have a limited set of values that can be given to them. To know which values these are one can use the filterOptions() function to retrieve the predetermed values of the respective filter.

filterOptions("biotype", ensembl)
## [1] "[IG_C_gene,IG_C_pseudogene,IG_D_gene,IG_J_gene,IG_J_pseudogene,IG_pseudogene,IG_V_gene,IG_V_pseudogene,lncRNA,miRNA,misc_RNA,Mt_rRNA,Mt_tRNA,polymorphic_pseudogene,processed_pseudogene,protein_coding,pseudogene,ribozyme,rRNA,rRNA_pseudogene,scaRNA,scRNA,snoRNA,snRNA,sRNA,TEC,transcribed_processed_pseudogene,transcribed_unitary_pseudogene,transcribed_unprocessed_pseudogene,translated_processed_pseudogene,translated_unprocessed_pseudogene,TR_C_gene,TR_D_gene,TR_J_gene,TR_J_pseudogene,TR_V_gene,TR_V_pseudogene,unitary_pseudogene,unprocessed_pseudogene,vault_RNA]"

If there are no predetermed values e.g. for the entrezgene filter, then filterOptions() will return the type of filter it is. And most of the times the filter name or it’s description will suggest what values one case use for the respective filter (e.g. entrezgene filter will work with enterzgene identifiers as values)

6.3 Attribute Pages

For large BioMart databases such as Ensembl, the number of attributes displayed by the listAttributes() function can be very large. In BioMart databases, attributes are put together in pages, such as sequences, features, homologs for Ensembl. An overview of the attributes pages present in the respective BioMart dataset can be obtained with the attributePages() function.

pages = attributePages(ensembl)
pages
## [1] "feature_page" "structure"    "homologs"     "snp"          "snp_somatic"  "sequences"

To show us a smaller list of attributes which belong to a specific page, we can now specify this in the listAttributes() function. The set of attributes is still quite long, so we use head() to show only the first few items here.

head(listAttributes(ensembl, page="feature_page"))
##                            name                  description         page
## 1               ensembl_gene_id               Gene stable ID feature_page
## 2       ensembl_gene_id_version       Gene stable ID version feature_page
## 3         ensembl_transcript_id         Transcript stable ID feature_page
## 4 ensembl_transcript_id_version Transcript stable ID version feature_page
## 5            ensembl_peptide_id            Protein stable ID feature_page
## 6    ensembl_peptide_id_version    Protein stable ID version feature_page

We now get a short list of attributes related to the region where the genes are located.

7 Using select()

In order to provide a more consistent interface to all annotations in Bioconductor the select(), columns(), keytypes() and keys() have been implemented to wrap some of the existing functionality above. These methods can be called in the same manner that they are used in other parts of the project except that instead of taking a AnnotationDb derived class they take instead a Mart derived class as their 1st argument. Otherwise usage should be essentially the same. You still use columns() to discover things that can be extracted from a Mart, and keytypes() to discover which things can be used as keys with select().

mart <- useMart(dataset="hsapiens_gene_ensembl",biomart='ensembl')
head(keytypes(mart), n=3)
## [1] "affy_hc_g110"        "affy_hg_focus"       "affy_hg_u133_plus_2"
head(columns(mart), n=3)
## [1] "3_utr_end"   "3_utr_end"   "3_utr_start"

And you still can use keys() to extract potential keys, for a particular key type.

k = keys(mart, keytype="chromosome_name")
head(k, n=3)
## [1] "1" "2" "3"

When using keys(), you can even take advantage of the extra arguments that are available for others keys methods.

k = keys(mart, keytype="chromosome_name", pattern="LRG")
head(k, n=3)
## character(0)

Unfortunately the keys() method will not work with all key types because they are not all supported.

But you can still use select() here to extract columns of data that match a particular set of keys (this is basically a wrapper for getBM()).

affy=c("202763_at","209310_s_at","207500_at")
select(mart, keys=affy, columns=c('affy_hg_u133_plus_2','entrezgene_id'),
  keytype='affy_hg_u133_plus_2')
##   affy_hg_u133_plus_2 entrezgene_id
## 1         209310_s_at           837
## 2           207500_at           838
## 3           202763_at           836

So why would we want to do this when we already have functions like getBM()? For two reasons: 1) for people who are familiar with select and it’s helper methods, they can now proceed to use biomaRt making the same kinds of calls that are already familiar to them and 2) because the select method is implemented in many places elsewhere, the fact that these methods are shared allows for more convenient programmatic access of all these resources. An example of a package that takes advantage of this is the OrganismDbi package. Where several packages can be accessed as if they were one resource.

8 Examples of biomaRt queries

In the sections below a variety of example queries are described. Every example is written as a task, and we have to come up with a biomaRt solution to the problem.

8.1 Annotate a set of Affymetrix identifiers with HUGO symbol and chromosomal locations of corresponding genes

We have a list of Affymetrix hgu133plus2 identifiers and we would like to retrieve the HUGO gene symbols, chromosome names, start and end positions and the bands of the corresponding genes. The listAttributes() and the listFilters() functions give us an overview of the available attributes and filters and we look in those lists to find the corresponding attribute and filter names we need. For this query we’ll need the following attributes: hgnc_symbol, chromsome_name, start_position, end_position, band and affy_hg_u133_plus_2 (as we want these in the output to provide a mapping with our original Affymetrix input identifiers. There is one filter in this query which is the affy_hg_u133_plus_2 filter as we use a list of Affymetrix identifiers as input. Putting this all together in the getBM() and performing the query gives:

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)
##   affy_hg_u133_plus_2 hgnc_symbol chromosome_name start_position end_position  band
## 1         209310_s_at       CASP4              11      104942866    104969366 q22.3
## 2           207500_at       CASP5              11      104994235    105023168 q22.3
## 3           202763_at       CASP3               4      184627696    184649509 q35.1

8.2 Annotate a set of EntrezGene identifiers with GO annotation

In this task we start out with a list of EntrezGene identiers and we want to retrieve GO identifiers related to biological processes that are associated with these entrezgene identifiers. Again we look at the output of listAttributes() and listFilters() to find the filter and attributes we need. Then we construct the following query:

entrez=c("673","837")
goids = getBM(attributes = c('entrezgene_id', 'go_id'), 
              filters = 'entrezgene_id', 
              values = entrez, 
              mart = ensembl)
head(goids)
##   entrezgene_id      go_id
## 1           673 GO:0005524
## 2           673 GO:0007165
## 3           673 GO:0006468
## 4           673 GO:0035556
## 5           673 GO:0004672
## 6           673 GO:0043231

8.3 Retrieve all HUGO gene symbols of genes that are located on chromosomes 17,20 or Y, and are associated with specific GO terms

The GO terms we are interested in are: GO:0051330, GO:0000080, GO:0000114, GO:0000082. The key to performing this query is to understand that the getBM() function enables you to use more than one filter at the same time. In order to do this, the filter argument should be a vector with the filter names. The values should be a list, where the first element of the list corresponds to the first filter and the second list element to the second filter and so on. The elements of this list are vectors containing the possible values for the corresponding filters.

 go=c("GO:0051330","GO:0000080","GO:0000114","GO:0000082")
 chrom=c(17,20,"Y")
 getBM(attributes= "hgnc_symbol",
        filters=c("go","chromosome_name"),
        values=list(go, chrom), mart=ensembl)
##   hgnc_symbol
## 1        RPA1
## 2       CRLF3
## 3        CDC6
## 4     RPS6KB1
## 5        CDK3
## 6        MCM8

8.4 Annotate set of idenfiers with INTERPRO protein domain identifiers

In this example we want to annotate the following two RefSeq identifiers: NM_005359 and NM_000546 with INTERPRO protein domain identifiers and a description of the protein domains.

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                                               interpro_description
## 1    NM_000546 IPR002117                                       p53 tumour suppressor family
## 2    NM_000546 IPR008967                         p53-like transcription factor, DNA-binding
## 3    NM_000546 IPR010991                                        p53, tetramerisation domain
## 4    NM_000546 IPR011615                                            p53, DNA-binding domain
## 5    NM_000546 IPR012346 p53/RUNT-type transcription factor, DNA-binding domain superfamily
## 6    NM_000546 IPR013872                                         p53 transactivation domain
## 7    NM_000546 IPR036674                        p53-like tetramerisation domain superfamily
## 8    NM_000546 IPR040926               Cellular tumor antigen p53, transactivation domain 2
## 9    NM_005359 IPR001132                                          SMAD domain, Dwarfin-type
## 10   NM_005359 IPR003619                                       MAD homology 1, Dwarfin-type
## 11   NM_005359 IPR008984                                        SMAD/FHA domain superfamily
## 12   NM_005359 IPR013019                                                  MAD homology, MH1
## 13   NM_005359 IPR013790                                                            Dwarfin
## 14   NM_005359 IPR017855                                       SMAD-like domain superfamily
## 15   NM_005359 IPR036578                                        SMAD MH1 domain superfamily

8.5 Select all Affymetrix identifiers on the hgu133plus2 chip and Ensembl gene identifiers for genes located on chromosome 16 between basepair 1100000 and 1250000.

In this example we will again use multiple filters: chromosome_name, start, and end as we filter on these three conditions. Note that when a chromosome name, a start position and an end position are jointly used as filters, the BioMart webservice interprets this as return everything from the given chromosome between the given start and end positions.

getBM(attributes = c('affy_hg_u133_plus_2','ensembl_gene_id'), 
      filters = c('chromosome_name','start','end'),
      values = list(16,1100000,1250000), 
      mart = ensembl)
##    affy_hg_u133_plus_2 ensembl_gene_id
## 1                      ENSG00000260702
## 2            215502_at ENSG00000260532
## 3                      ENSG00000273551
## 4            205845_at ENSG00000196557
## 5                      ENSG00000196557
## 6                      ENSG00000260403
## 7                      ENSG00000259910
## 8                      ENSG00000261294
## 9          220339_s_at ENSG00000116176
## 10                     ENSG00000277010
## 11         205683_x_at ENSG00000197253
## 12         207134_x_at ENSG00000197253
## 13         217023_x_at ENSG00000197253
## 14         210084_x_at ENSG00000197253
## 15         215382_x_at ENSG00000197253
## 16         216474_x_at ENSG00000197253
## 17         205683_x_at ENSG00000172236
## 18         207134_x_at ENSG00000172236
## 19         217023_x_at ENSG00000172236
## 20         210084_x_at ENSG00000172236
## 21         215382_x_at ENSG00000172236
## 22         216474_x_at ENSG00000172236

8.6 Retrieve all EntrezGene identifiers and HUGO gene symbols of genes which have a “MAP kinase activity” GO term associated with it.

The GO identifier for MAP kinase activity is GO:0004707. In our query we will use go_id as our filter, and entrezgene_id and hgnc_symbol as attributes. Here’s the query:

getBM(attributes = c('entrezgene_id','hgnc_symbol'), 
      filters = 'go', 
      values = 'GO:0004707', 
      mart = ensembl)
##    entrezgene_id hgnc_symbol
## 1           5599       MAPK8
## 2           5891         MOK
## 3           5597       MAPK6
## 4           5595       MAPK3
## 5           5598       MAPK7
## 6          51701         NLK
## 7           5596       MAPK4
## 8           5594       MAPK1
## 9           6300      MAPK12
## 10          5600      MAPK11
## 11          5602      MAPK10
## 12          5601       MAPK9
## 13          4117         MAK
## 14          6732       SRPK1
## 15          1432      MAPK14
## 16          5603      MAPK13
## 17         22858       CILK1
## 18          6733       SRPK2
## 19        225689      MAPK15
## 20         26576       SRPK3

8.7 Given a set of EntrezGene identifiers, retrieve 100bp upstream promoter sequences

All sequence related queries to Ensembl are available through the getSequence() wrapper function. getBM() can also be used directly to retrieve sequences but this can get complicated so using getSequence is recommended.

Sequences can be retrieved using the getSequence() function either starting from chromosomal coordinates or identifiers.
The chromosome name can be specified using the chromosome argument. The start and end arguments are used to specify start and end positions on the chromosome. The type of sequence returned can be specified by the seqType argument which takes the following values:

  • cdna
  • peptide for protein sequences
  • 3utr for 3’ UTR sequences
  • 5utr for 5’ UTR sequences
  • gene_exon for exon sequences only
  • transcript_exon for transcript specific exonic sequences only
  • transcript_exon_intron gives the full unspliced transcript, that is exons + introns
  • gene_exon_intron gives the exons + introns of a gene
  • coding gives the coding sequence only
  • coding_transcript_flank gives the flanking region of the transcript including the UTRs, this must be accompanied with a given value for the upstream or downstream attribute
  • coding_gene_flank gives the flanking region of the gene including the UTRs, this must be accompanied with a given value for the upstream or downstream attribute
  • transcript_flank gives the flanking region of the transcript exculding the UTRs, this must be accompanied with a given value for the upstream or downstream attribute
  • gene_flank gives the flanking region of the gene excluding the UTRs, this must be accompanied with a given value for the upstream or downstream attribute

In MySQL mode the getSequence() function is more limited and the sequence that is returned is the 5’ to 3’+ strand of the genomic sequence, given a chromosome, as start and an end position.

This task requires us to retrieve 100bp upstream promoter sequences from a set of EntrezGene identifiers. The type argument in getSequence() can be thought of as the filter in this query and uses the same input names given by listFilters(). In our query we use entrezgene_id for the type argument. Next we have to specify which type of sequences we want to retrieve, here we are interested in the sequences of the promoter region, starting right next to the coding start of the gene. Setting the seqType to coding_gene_flank will give us what we need. The upstream argument is used to specify how many bp of upstream sequence we want to retrieve, here we’ll retrieve a rather short sequence of 100bp. Putting this all together in getSequence() gives:

entrez=c("673","7157","837")
getSequence(id = entrez, 
            type="entrezgene_id",
            seqType="coding_gene_flank",
            upstream=100, 
            mart=ensembl) 
##                                                                                      coding_gene_flank
## 1 TCCTTCTCTGCAGGCCCAGGTGACCCAGGGTTGGAAGTGTCTCATGCTGGATCCCCACTTTTCCTCTTGCAGCAGCCAGACTGCCTTCCGGGTCACTGCC
## 2 CACGTTTCCGCCCTTTGCAATAAGGAAATACATAGTTTACTTTCATTTTTGACTCTGAGGCTCTTTCCAACGCTGTAAAAAAGGACAGAGGCTGTTCCCT
## 3 CCTCCGCCTCCGCCTCCGCCTCCGCCTCCCCCAGCTCTCCGCCTCCCTTCCCCCTCCCCGCCCGACAGCGGCCGCTCGGGCCCCGGCTCTCGGTTATAAG
##   entrezgene_id
## 1          7157
## 2           837
## 3           673

8.8 Retrieve all 5’ UTR sequences of all genes that are located on chromosome 3 between the positions 185,514,033 and 185,535,839

As described in the provious task getSequence can also use chromosomal coordinates to retrieve sequences of all genes that lie in the given region. We also have to specify which type of identifier we want to retrieve together with the sequences, here we choose for entrezgene identifiers.

utr5 = getSequence(chromosome=3, start=185514033, end=185535839,
                   type="entrezgene_id",
                   seqType="5utr", 
                   mart=ensembl)
utr5
##                                                                                                                    5utr
## 1                                                                                                  Sequence unavailable
## 2 ACCACACCTCTGAGTCGTCTGAGCTCACTGTGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 3                               TGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 4                                                                               ATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
##   entrezgene_id
## 1        200879
## 2        200879
## 3        200879
## 4        200879

8.9 Retrieve protein sequences for a given list of EntrezGene identifiers

In this task the type argument specifies which type of identifiers we are using. To get an overview of other valid identifier types we refer to the listFilters() function.

protein = getSequence(id=c(100, 5728),
                      type="entrezgene_id",
                      seqType="peptide", 
                      mart=ensembl)
protein
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             peptide
## 1 LERGGEAAAAAAAAAAAPGRGSESPVTISRAGNAGELVSPLLLPPTRRRRRRHIQGPGPVLNLPSAAAAPPVARAPEAAGGGSRSEDYSSSPHSAAAAARPLAAEEKQAQSLQPSSSRRSSHYPAAVQSQAAAERGASATAKSRAISILQKKPRHQQLLPSLSSFFFSHRLPDMTAIIKEIVSRNKRRYQEDGFDLDLTYIYPNIIAMGFPAERLEGVYRNNIDDVVRFLDSKHKNHYKIYNLCAERHYDTAKFNCRVAQYPFEDHNPPQLELIKPFCEDLDQWLSEDDNHVAAIHCKAGKGRTGVMICAYLLHRGKFLKAQEALDFYGEVRTRDKKGVTIPSQRRYVYYYSYLLKNHLDYRPVALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVKLYFTKTVEEPSNPEASSSTSVTPDVSDNEPDHYRYSDTTDSDPENEPFDEDQHTQITKV*
## 2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Sequence unavailable
## 3                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIARL*
## 4                                                                                                                                                                                                                                                                                                                                                                                                                                        ALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVS*
## 5                                                                                                                                                                              MTAIIKEIVSRNKRRYQEDGFDLDLTYIYPNIIAMGFPAERLEGVYRNNIDDVVRFLDSKHKNHYKIYNLCAERHYDTAKFNCRVAQYPFEDHNPPQLELIKPFCEDLDQWLSEDDNHVAAIHCKAGKGRTGVMICAYLLHRGKFLKAQEALDFYGEVRTRDKKGVTIPSQRRYVYYYSYLLKNHLDYRPVALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVKLYFTKTVEEPSNPEASSSTSVTPDVSDNEPDHYRYSDTTDSDPENEPFDEDQHTQITKV*
## 6                                                                                                                                                                                                                      MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 7                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Sequence unavailable
## 8                                                                                                                                                                                                                                              MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 9                                                                                                                                                                                                                                                                                                                          MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEAQK*
##   entrezgene_id
## 1          5728
## 2           100
## 3           100
## 4          5728
## 5          5728
## 6           100
## 7          5728
## 8           100
## 9           100

8.10 Retrieve known SNPs located on the human chromosome 8 between positions 148350 and 148612

For this example we’ll first have to connect to a different BioMart database, namely snp.

snpmart = useEnsembl(biomart = "snp", dataset="hsapiens_snp")

The listAttributes() and listFilters() functions give us an overview of the available attributes and filters.
From these we need: refsnp_id, allele, chrom_start and chrom_strand as attributes; and as filters we’ll use: chrom_start, chrom_end and chr_name.
Note that when a chromosome name, a start position and an end position are jointly used as filters, the BioMart webservice interprets this as return everything from the given chromosome between the given start and end positions. Putting our selected attributes and filters into getBM gives:

getBM(attributes = c('refsnp_id','allele','chrom_start','chrom_strand'), 
      filters = c('chr_name','start','end'), 
      values = list(8,148350,148612), 
      mart = snpmart)
##       refsnp_id  allele chrom_start chrom_strand
## 1  rs1450830176     G/C      148350            1
## 2  rs1360310185     C/T      148352            1
## 3  rs1434776028     A/T      148353            1
## 4  rs1413161474     C/T      148356            1
## 5  rs1410590268     A/G      148365            1
## 6  rs1193735780   T/A/C      148368            1
## 7  rs1409139861     C/T      148371            1
## 8   rs868546642     A/G      148372            1
## 9   rs547420070     A/C      148373            1
## 10 rs1236874674     C/T      148375            1
## 11 rs1207902742     C/T      148376            1
## 12 rs1437239557     T/C      148377            1
## 13 rs1160135941     T/G      148379            1
## 14 rs1229249227     A/T      148380            1
## 15 rs1328678285     C/G      148390            1
## 16   rs77274555   G/A/C      148391            1
## 17  rs567299969   T/A/C      148394            1
## 18 rs1457776094   A/C/G      148395            1
## 19 rs1292078334     C/T      148404            1
## 20 rs1456392146     A/T      148405            1
## 21  rs368076569     G/A      148407            1
## 22 rs1166604298     A/G      148408            1
## 23 rs1393545673   A/G/T      148410            1
## 24 rs1180076939     A/T      148413            1
## 25 rs1476313471     A/G      148414            1
## 26 rs1248424402     T/C      148417            1
## 27 rs1207939741     A/T      148419            1
## 28 rs1442232389     A/T      148421            1
## 29 rs1266518345     A/C      148424            1
## 30 rs1224022305     A/G      148426            1
## 31 rs1322431371     T/C      148427            1
## 32 rs1289196156     T/C      148429            1
## 33 rs1227114952     C/G      148433            1
## 34 rs1367904435     A/C      148436            1
## 35 rs1298979710     A/C      148440            1
## 36 rs1434874738     T/G      148441            1
## 37 rs1363439435     A/C      148443            1
## 38 rs1294944280     T/C      148453            1
## 39 rs1309134024     C/T      148454            1
## 40 rs1348796979     C/T      148460            1
## 41 rs1374829532     T/C      148461            1
## 42 rs1431088173     G/A      148462            1
## 43 rs1390780034     C/T      148468            1
## 44 rs1187817390     C/G      148470            1
## 45 rs1242396084     T/C      148471            1
## 46 rs1447487229   TTT/T      148471            1
## 47 rs1193176841     A/G      148474            1
## 48 rs1488893482     C/T      148478            1
## 49 rs1264956600     T/A      148483            1
## 50 rs1219544655     A/G      148487            1
## 51 rs1344406365 CTCT/CT      148495            1
## 52  rs745318437     C/G      148497            1
## 53 rs1321280897     G/T      148499            1
## 54 rs1220571972     C/T      148504            1
## 55 rs1228245996     G/T      148508            1
## 56 rs1340516169     C/G      148512            1
## 57 rs1263402165     C/T      148513            1
## 58 rs1303775291     C/G      148517            1
## 59 rs1387005982   C/A/T      148520            1
## 60 rs1291424168     T/A      148526            1
## 61 rs1411168240     T/C      148529            1
## 62 rs1370615124     C/T      148536            1
## 63 rs1469788925     T/C      148559            1
## 64 rs1214284488     T/C      148561            1
## 65 rs1237469667     G/T      148563            1
## 66  rs190721891     C/G      148576            1
## 67 rs1427698597     C/G      148594            1
## 68 rs1179966097     T/C      148595            1

8.11 Given the human gene TP53, retrieve the human chromosomal location of this gene and also retrieve the chromosomal location and RefSeq id of its homolog in mouse.

The getLDS() (Get Linked Dataset) function provides functionality to link 2 BioMart datasets which each other and construct a query over the two datasets. In Ensembl, linking two datasets translates to retrieving homology data across species. The usage of getLDS is very similar to getBM(). The linked dataset is provided by a separate Mart object and one has to specify filters and attributes for the linked dataset. Filters can either be applied to both datasets or to one of the datasets. Use the listFilters and listAttributes functions on both Mart objects to find the filters and attributes for each dataset (species in Ensembl). The attributes and filters of the linked dataset can be specified with the attributesL and filtersL arguments. Entering all this information into getLDS() gives:

human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = c("hgnc_symbol","chromosome_name", "start_position"),
       filters = "hgnc_symbol", values = "TP53",mart = human,
      attributesL = c("refseq_mrna","chromosome_name","start_position"), martL = mouse)
##   HGNC.symbol Chromosome.scaffold.name Gene.start..bp. RefSeq.mRNA.ID Chromosome.scaffold.name.1
## 1        TP53                       17         7661779                                        11
## 2        TP53                       17         7661779   NM_001127233                         11
## 3        TP53                       17         7661779      NM_011640                         11
##   Gene.start..bp..1
## 1          69580359
## 2          69580359
## 3          69580359

9 Session Info

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.11-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] biomaRt_2.44.4   BiocStyle_2.16.1
## 
## loaded via a namespace (and not attached):
##  [1] progress_1.2.2       tidyselect_1.1.0     xfun_0.18            purrr_0.3.4         
##  [5] vctrs_0.3.4          generics_0.0.2       htmltools_0.5.0      stats4_4.0.3        
##  [9] BiocFileCache_1.12.1 yaml_2.2.1           blob_1.2.1           XML_3.99-0.5        
## [13] rlang_0.4.8          pillar_1.4.6         glue_1.4.2           DBI_1.1.0           
## [17] rappdirs_0.3.1       BiocGenerics_0.34.0  bit64_4.0.5          dbplyr_1.4.4        
## [21] lifecycle_0.2.0      stringr_1.4.0        codetools_0.2-16     memoise_1.1.0       
## [25] evaluate_0.14        Biobase_2.48.0       knitr_1.30           IRanges_2.22.2      
## [29] parallel_4.0.3       curl_4.3             AnnotationDbi_1.50.3 highr_0.8           
## [33] Rcpp_1.0.5           openssl_1.4.3        BiocManager_1.30.10  S4Vectors_0.26.1    
## [37] bit_4.0.4            hms_0.5.3            askpass_1.1          png_0.1-7           
## [41] digest_0.6.25        stringi_1.5.3        bookdown_0.21        dplyr_1.0.2         
## [45] tools_4.0.3          magrittr_1.5         RSQLite_2.2.1        tibble_3.0.4        
## [49] crayon_1.3.4         pkgconfig_2.0.3      ellipsis_0.3.1       xml2_1.3.2          
## [53] prettyunits_1.1.1    assertthat_0.2.1     rmarkdown_2.4        httr_1.4.2          
## [57] R6_2.4.1             compiler_4.0.3
warnings()