In patzaw/BED: Biological Entity Dictionary (BED)
reComp <- FALSE
library(knitr)
library(RColorBrewer)
library(VennDiagram)
library(gridExtra)
flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")
opts_chunk$set(
include=FALSE,
warning=FALSE, echo=FALSE, message=FALSE,
concordance=TRUE
)
bigFormat <- function(x) {format(x, big.mark=",")}
library(BED)
connectToBed(url="localhost:5454", username="neo4j", password="1234")
\newcommand{\tm}{\textsuperscript{\textregistered}}
\newcommand{\neo}{Neo4j\tm{}}
\newcommand{\cypher}{Cypher\tm{}}
\newcommand{\docker}{Docker\tm{}}
\newcommand{\metabase}{MetaBase\tm{}}
Introduction
Since the advent of genome sequencing projects, many technologies have been
developped to get access to different molecular information at a large scale
and with high throughput. DNA microarrays are probably the archetype of such
technology because of their historical impact on gathering data related
to nucleic acids: genomic DNA and RNA. They triggered the emergence of
"omics" fields of research such as genomics, epigenomics or transcriptomics.
Lately massive parallel sequencing
further increased the throughput of data generation related to nucleic acids
by several orders of magnitude.
In a different way mass spectrometry-related technologies allow the
identification and the quantification of many kinds of molecular entities
such as metabolites and proteins.
Many information systems have been developped to manage
the exploding amount of data and knowledge related to biological
molecular entities.
These resources manage different aspects of
the knowledge.
For example some are genome or proteome centered whereas other
are focused on molecular interactions and pathways.
Thus all these resources rely on different identifier systems to organise
the concepts of interest.
The value of all the experimental data and all the knowledge
collected in public or private resources is very high as such
but is also often synegistically leveraged by their cross comparison
in a dedicated manner. Indeed many datasets can be relevant when
adressing the understanding of a specific biological system, a phenotypic trait
or a disease for example. These datasets can focus on different biological
entities such as transcripts or proteins in different tissues, conditions
or organisms. Comparing all these data and integrating them with
available knowledge require the ability to map the identifiers on which
each resource relies.
To achieve this task public and proprietary information systems
provide mapping tables between their own identifiers and those
from other resources.
Furthermore many tools have been developped to facilitate the access
to this information.
Ensembl BioMarts [@kinsella_ensembl_2011],
mygene [@wu_biogps_2013],
and g:Profiler [@reimand_g:profiler-web_2016]
are popular examples among many others.
However, as pointed by @van_iersel_bridgedb_2010, these tools
are generally dedicated to a particular domain not necessarly
relevant or complete for all research project and keeping them
up-to-date can also be an issue.
Recognizing these challenges @van_iersel_bridgedb_2010 proposed
the BridgeDb framework providing to bioinformatics developers
a standard interface between tools and
mapping services and also allowing the easy integration of custom data
by a transitivity mechanism.
Here we present BED: a biological entity dictionary.
BED has been developped to adress three main challenges.
The first one is related to the completeness of identifier mappings.
Indeed direct mapping information provided by the different
systems are not always complete and can be enriched by mappings provided
by other resources. More interestingly direct mappings not identified by
any of these resources can be indireclty infered by using mappings to a
third reference. For example many human Ensembl gene identifiers not
directly mapped to any Entrez gene identifiers by neither Ensembl nor the NCBI
can be indirectly mapped using mappings to HGNC identifiers.
The second challenge is related to the mapping of deprecated identifiers.
Indeed entity identifiers can change from one resource release to another.
The identifier history is provided by some resources, such as Ensembl or
the NCBI, but it is generally not used by mapping tools.
The third challenge is related to the automation of the mapping process
according to the relationships between the biological entities of interest.
Indeed mapping between gene and protein identifier scopes should not be done
the same way than two scopes of gene identifiers.
Also converting identifiers from different organism should be possible
using gene ortholog information.
To meet these challenges we designed a graph data model describing possible
relationships between different biological entities and their identifiers.
This data model has been implemented with the \neo{} graph
database [@neo4j_inc_neo4j_2017]
and conversion rules have been defined and coded in
an R [@r_core_team_r:_2017] package.
We provide an instance of the BED database focused on human, mouse and
rat organism but many functions are available to construct other instances
tailored to other needs.
Methods
Data model
include_graphics('img/BED-Data-Model.pdf')
The BED (Biological Entity Dictionary) system relies on a data model
inspired by the central dogma of
molecular biology [@crick_central_1970] and describing
relationships between molecular concepts
usually manipulated in the frame of genomics studies
(Figure \@ref(fig:Data-Model)).
A biological entity identifier (BEID) can identify either a Gene (GeneID),
a Transcript (TranscriptID), a Peptide (PeptideID) or
an Object (ObjectID). Object entities can correspond to complex concepts
coded by any number of genes (i.e. a protein complex or a molecular function).
BEID are extracted from public or private databases (BEDB).
BEDB can provide an Attribute related to each BEID.
For example it can be the
sequencing region provided by the Ensembl database [@zerbino_ensembl_2017]
or the identifier
status provided by Uniprot [@the_uniprot_consortium_uniprot:_2017].
BEID can have one or several associated names (BENames) and symbols
(BESymbol). GeneID can have one or several homologs in other organisms
belonging to the same GeneIDFamily.
Many genomics platforms, such as microarray, allow the identification of
biological entity by using probes identified by ProbeID. In general
BEID can be targeted by several probes belonging to a Platform which
is focused on one and only one type of entity (BEType) among those described
above: Gene, Transcript, Peptide or Object. A BEType can have several
BEType products but can be the product of at most one BEType.
This constraint allows the unambiguous identification of the most relevant path
to convert identfiers from one scope to another and
is fulfilled by the current data model:
peptides are only produced from transcripts which are only produced
from genes which can also code for objects.
BEID identifying the same biological entity are related through three
different kinds of relationship according to the information available
in the source databases and to the decision made by the database administrator
about how to use them. Two BEID which corresponds_to each other both
identify the same biological entity. A BEID which is_associated_to or
which is_replaced_by another BEID does not directly
identify any biological entity: the link is always indirect through
one or several other BEID. Therefore by design a BEID
which is_associated_to or which is_replaced_by another BEID can be
related to several different biological entities. It is not the case for
other BEID which identify one and only one biological entity.
This set of possible relationship allows the indirect mapping of
different identifiers not necessarily provided by any integrated resource.
In order to efficiently leverage indirect path through these different
relationships the data model has been implemented in
a \neo{} graph database [@neo4j_inc_neo4j_2017].
Feeding the database
Two R [@r_core_team_r:_2017] packages have been developed to feed and query
the database.
The first one, neo2R, provides low level functions to interact with
\neo{}. The second R package, BED,
provides functions to feed and query the BED
\neo{} graph database according to the
data model described above.
Many functions are provided within the package to build a tailored BED database
instance. These functions are not exported in order not to mislead
the user when querying the database (which is the expected most frequent
usage of the system). An R markdown document showing how to build a BED database
instance for human, mouse and rat organisms is provided within the
package. It can be adapted to other organisms or needs.
Briefly these functions can be divided according to three main levels:
- The lowest level function is the
bedImport function
which load a table in the \neo{} database according to a \cypher{} query.
- Functions of the second level allow loading identifiers and relationships
tables ensuring the integrity of the data model.
- Highest level functions are helpers for loading information provided by
some public resources in different specific format.
Querying the database
The BED R package provides several functions to retrieve identifiers
from different resources and also to convert identifiers from one
reference to another.
These functions generate and call \cypher{} queries on the \neo{} database.
Converting thousands of identifiers can take some time
(generally a few seconds). Also such conversions are often recurrent and
redundant. In order to improve the performance for such recurrent and redundant
queries, a cache system has been implemented. The first time, the query is run
on \neo{} for all the relevant ID related to user input and the result is saved
in a local file. Next time similar queries are requested, the system does
not call \neo{} but loads the cached results and filters it according
to user input.
By default the cache is flushed when the system detects inconsistencies
with the BED database. It can also be manually flushed if needed.
Operation
The graph database has been implemented
with \neo{} version 3 [@neo4j_inc_neo4j_2017].
The BED R package depends on the following packages available in the
Comprehensive R Archive Network [@cran_comprehensive_nodate]:
- visNetwork [@almende_b.v._visnetwork:_2017]
- dplyr [@wickham_dplyr:_2017]
- htmltools [@rstudio_inc_htmltools:_2017]
- DT [@xie_dt:_2016]
- shiny [@chang_shiny:_2017]
- miniUI [@cheng_miniui:_2016]
- rstudioapi [@allaire_rstudioapi:_2017]
Use Cases
Available database instance
An instance of the BED database (UCB-Human)
has been built using the script provided
in the BED R package and made available in a \docker{}
image [@docker_inc_docker_2017] available here:
https://hub.docker.com/r/patzaw/bed-ucb-human/
This instance used to exemplify the following use cases
is focused on Homo sapiens, Mus musculus and Rattus norvegicus organisms
and it has been built from the following resources:
- Ensembl [@zerbino_ensembl_2017]
- NCBI [@ncbi_resource_coordinators_database_2017]
- Uniprot [@the_uniprot_consortium_uniprot:_2017]
- biomaRt [@durinck_mapping_2009]
- GEOquery [@davis_geoquery:_2007]
- Clarivate Analytics \metabase{} [@clarivate_analytics_metacore_2017]
toShow <- c()
for(be in listBe()){
for(org in BED::listOrganisms()){
toAdd <- listBeIdSources(
be=be,
organism=org,
exclude=c("BEDTech_gene", "BEDTech_transcript"),
verbose=FALSE
)
colnames(toAdd) <- c("Database", "nbBE", "BEID", "BE")
toAdd$Organism <- org
toShow <- rbind(
toShow,
toAdd[,c("BE", "Organism", "Database", "BEID")]
)
}
}
urls <- c()
for(i in 1:nrow(toShow)){
urls <- c(urls, getBeIdURL("", toShow[i, "Database"]))
}
urls <- substr(
urls, start=1, stop=regexpr("[^/][/][^/].*$", urls)
)
# urls <- sub("^https*[:][/][/]", "", urls)
toShow$URL <- urls
kable(
toShow,
row.names=FALSE,
format='pandoc',
format.args=list(big.mark=","),
caption='Numbers of BEID available in the BED UCB-Human database instance.
Numbers have been split according to the BE type
and the organism.
Only BEID which can be mapped to each other are taken into account.'
)
nbBeIds <- 3519181
nbProbeIds <- 354205
## Total number of BEID
nbBeIds <- bedCall(
cypher,
query <- prepCql(c(
'MATCH (n:BEID)',
#'-[:is_replaced_by|is_associated_to*0..]->(:BEID)-[:identifies]->()',
'WHERE NOT n.database IN ["BEDTech_gene", "BEDTech_transcript"]',
'RETURN count(distinct n);'
))
)
nbProbeIds <- bedCall(
cypher,
query <- prepCql(c(
'MATCH (n:ProbeID)',
'RETURN count(distinct n);'
))
)
The numbers of biological entity (BE) identifiers (BEID) available in this
BED database instance and which can be mapped to each other are shown
in table \@ref(tab:ID-Numbers).
In total, r bigFormat(nbBeIds) BEID are available in this
BED instance but many
of them cannot be mapped to other BEID because corresponding to deprecated
identifiers not corresponding to any identifier currently in use.
Also all the genomics platforms included in this BED database instance are
shown in table \@ref(tab:Platforms) . They provide mapping to BEID from
r bigFormat(nbProbeIds) ProbeID in total.
toShow <- listPlatforms()
colnames(toShow) <- c("Name", "Description", "BE")
toShow$Description <- sub(" [(].*$", "", toShow$Description)
kable(
toShow,
row.names=FALSE,
format='pandoc',
caption='Genomics platforms available in the BED UCB-Human database instance.'
)
Exploring identifiers of biological entities
The getBeIds function returns all BE identifiers from a specific scope.
A scopes is defined
by the type of BE or probe, the source of the identifiers (database or platform)
and the organism.
For example the following code
returns all the Ensembl identifiers of human genes.
beids <- getBeIds(
be="Gene", source="Ens_gene", organism="human",
restricted=FALSE
)
head(beids)
redUnRestIds <- table(table(beids$Gene))
redUnRestIds <- redUnRestIds[order(as.numeric(names(redUnRestIds)))]
ensVersion <- bedCall(
cypher,
query=prepCql(c(
'MATCH (n:BEDB {name:"Ens_gene"}) RETURN n.currentVersion'
))
)
layout(matrix(c(1,1,1,1,2,2,2), 1, 7, byrow = TRUE))
par(mar=c(5.1, 4.5, 1.1, 0.5))
barplot(
redUnRestIds,
xlab="Number of BEID",
ylab="Number of BE", log="y",
las=2,
cex.axis=0.7,
cex.names=0.7
)
mtext(text="(a)", side=2, line=3, at=10^par()$usr[4], las=2, font=2)
##
rbeids <- getBeIds(
be="Gene", source="Ens_gene", organism="human",
restricted = TRUE
)
redIds <- table(table(rbeids$Gene))
redIds <- redIds[order(as.numeric(names(redIds)))]
barplot(
redIds,
xlab="Number of BEID",
ylab="Number of BE", log="y",
las=2,
cex.axis=0.7,
cex.names=0.7
)
mtext(text="(b)", side=2, line=3, at=10^par()$usr[4], las=2, font=2)
The id column corresponds to the BEID from the source of interest.
The column named according to the BE type (in this case Gene)
corresponds to the internal identifiers of the related BE.
This internal identifier is not a stable reference that can be used as such.
Nevertheless it is useful to identify BEID identifying the
same BE.
In the example above even if most of Gene BE are identified by only
one Ensembl gene BEID, many of them are identified by two or more
(r bigFormat(sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=2)]))
/ r bigFormat(sum(redUnRestIds))
= r round(100*sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=2)])/sum(redUnRestIds)) %);
r sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=10)]) BE
are even identified by more than 10 Ensembl BEID
(Figure \@ref(fig:beidTables).a).
In this case, most of these redundancies come from deprecated ID from former
versions of the Ensembl database (version in used here: r ensVersion)
and can be excluded by setting the restricted parameter to TRUE when calling
the getBeIds function (Figure \@ref(fig:beidTables).b).
However many BE are still identified by two or more current Ensembl BEID
(r bigFormat(sum(redIds[which(as.numeric(names(redIds))>=2)]))
/ r bigFormat(sum(redIds))
= r round(100*sum(redIds[which(as.numeric(names(redIds))>=2)])/sum(redIds))~%).
This result comes from the way the BED database is constructed:
When two identifiers from the same resource correspond to the same identifier
in another resource (correspond_to relationship in the data model),
all these BEID are considered to identify the same BE.
A complex example of such mapping is shown in figure \@ref(fig:TAS2R8)
mapping all the BEID of the human TAS2R8 gene which codes for a protein
of the family of candidate taste receptors. There are three identifiers
corresponding to this gene symbol in Ensembl. All these three identifiers
correspond to the same Entrez gene and the same HGNC identifiers.
All these BEID are thus considered to identify the same gene. It turns out
that the three Ensembl BEID correspond to the same gene mapped on different
sequence version of the chromosome 12: the canonical (ENSG00000121314),
CHR_HSCHR12_2_CTG2 (ENSG00000272712)
and CHR_HSCHR12_3_CTG2 (ENSG00000277316).
This information provided by Ensembl is encoded in the seq_region
attribute for each Ensembl BEID (see data model)
and is used to define preferred BEID which are mapped on canonical version
of chromosome sequences.
The ENSG00000272712 identifier shows also a complex history in former
Ensembl versions.
```r identifiers. Solid arrows represent \emph{correspond\_to} and \emph{is\_known\_as} relationships. Dotted arrows represent \emph{is\_replaced\_by} and \emph{is\_associated\_to} relationships. This graph has been drawn with the \texttt{exploreBe} function.'}
include_graphics('img/TAS2R8-Identifiers.png')
## Converting identifiers
```r
## Save/load results from commmands below to not recompute everything
load("~/Tmp/ConvTables-3rdTools.rda")
The main goal of BED is to convert identifiers from one scope to another
easily, rapidly and with high completeness.
It has been thought in order to allow recurring comparisons to each other of
many lists of biological entities from various origins.
The function guessIdOrigin can be used to guess the scope
of any list of identifiers.
A simple example regarding the conversion of human Ensembl
gene to human Entrez gene identifiers is shown below and discussed hereafter.
By setting the restricted parameter to TRUE the converted BEID are
restricted to current - non-deprecated - version of Entrez gene identifiers.
Nevertheless all the input BEID are taken into account,
current and deprecated ones.
runTimes <- list()
##
clearBedCache()
s <- Sys.time()
bedConv <- convBeIds(
ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human",
to.source="EntrezGene", restricted=TRUE
)
e <- Sys.time()
runTimes[["BED (Not cached)"]] <- e-s
##
s <- Sys.time()
bedConv <- convBeIds(
ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human",
to.source="EntrezGene", restricted=TRUE
)
e <- Sys.time()
runTimes[["BED (Cached)"]] <- e-s
bedConv <- convBeIds(
ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human",
to.source="EntrezGene", restricted=TRUE
)
curBeid <- unique(beids$id[which(beids$db.deprecated=="FALSE")])
depBeid <- unique(beids$id[which(beids$db.deprecated!="FALSE")])
Among all the r bigFormat(length(curBeid)+ length(depBeid))
human Ensembl gene identifiers available in the
database, r bigFormat(sum(is.na(bedConv$to)))
(r round(100*sum(is.na(bedConv$to))/length(unique(bedConv$from))) %)
were not converted to any human Entrez gene
identifier: r bigFormat(sum(is.na(bedConv$to[which(bedConv$from %in% curBeid)])))
(r round(100*sum(is.na(bedConv$to[which(bedConv$from %in% curBeid)]))/length(curBeid)) %)
of the r bigFormat(length(curBeid)) non-deprecated
and r bigFormat(sum(is.na(bedConv$to[which(bedConv$from %in% depBeid)])))
(r round(100*sum(is.na(bedConv$to[which(bedConv$from %in% depBeid)]))/length(depBeid)) %)
of the r bigFormat(length(depBeid)) deprecated identifiers.
## biomaRt
library(biomaRt)
bmEnsembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
s <- Sys.time()
bmConv <- getBM(
values = unique(beids$id),
filters = 'ensembl_gene_id',
attributes=c('ensembl_gene_id', 'entrezgene'),
mart = bmEnsembl
)
e <- Sys.time()
runTimes[["biomaRt"]] <- e-s
## mygene
library(mygene)
s <- Sys.time()
mgConv <- queryMany(
qterm=unique(beids$id),
scopes = "ensembl.gene",
fields="entrezgene",
species="human"
)
e <- Sys.time()
runTimes[["mygene"]] <- e-s
## gProfileR
library(gProfileR)
s <- Sys.time()
gpConv <- gconvert(
query=unique(beids$id),
target="ENTREZGENE_ACC",
organism="hsapiens",
filter_na=FALSE,
df=FALSE
)
gpConv <- do.call(
rbind,
lapply(
gpConv,
function(x){
data.frame(
alias=as.character(x[,"alias"]),
target=as.character(x[,"target"]),
stringsAsFactors=FALSE
)
}
)
)
gpConv$target <- sub("ENTREZGENE_ACC:", "", gpConv$target)
gpConv$target <- ifelse(gpConv$target=="N/A", NA, gpConv$target)
rownames(gpConv) <- c()
gpConv <- unique(gpConv)
e <- Sys.time()
runTimes[["gProfileR"]] <- e-s
evalDate <- Sys.Date()
save(
bedConv, bmConv, mgConv, gpConv, evalDate, runTimes,
file="~/Tmp/ConvTables-3rdTools.rda"
)
Three other tools were used on r format(evalDate, "%b %d, %Y")
to perform the same conversion task:
biomaRt [@kinsella_ensembl_2011; @durinck_mapping_2009],
mygene [@wu_biogps_2013; @mark_mygene:_2014],
and gProfileR [@reimand_g:profiler-web_2016; @reimand_gprofiler:_2016].
At that time, biomaRt and mygene were based on the Ensembl 91 release
whereas gProfileR was based on release 90.
gq.all <- unique(beids$id)
toComp.all <- list(
"BED"=unique(bedConv$from[which(!is.na(bedConv$to))]),
"biomaRt"=unique(bmConv$ensembl_gene_id[which(!is.na(bmConv$entrezgene))]),
"mygene"=unique(mgConv$query[which(!is.na(mgConv$entrezgene))]),
"gProfileR"=unique(gpConv$alias[which(!is.na(gpConv$target))])
)
gq.cur <- unique(beids$id[which(beids$db.deprecated=="FALSE")])
toComp.cur <- lapply(
toComp.all,
intersect,
gq.cur
)
venn.plot.all <- venn.diagram(
toComp.all, NULL,
fill=brewer.pal(4, "Accent"),
alpha=rep(0.5, 4),
main="(a)", main.pos=c(0.05,0.95),
main.fontface=2, main.cex=1.8
)
venn.plot.cur <- venn.diagram(
toComp.cur, NULL,
fill=brewer.pal(4, "Accent"),
alpha=rep(0.5, 4),
margin=0.1,
main="(b)", main.pos=c(0.05,0.95),
main.fontface=2, main.cex=1.8
)
grid.newpage()
grid.arrange(
gTree(children=venn.plot.all, vp=viewport(width=0.8,height=0.9)),
gTree(children=venn.plot.cur, vp=viewport(width=0.8,height=0.9)),
ncol=2
)
onlyInBed <- bedConv[
which(
!bedConv$from %in% unlist(toComp.all[-1]) &
!is.na(bedConv$to)
),
]
onlyInBedCur <- onlyInBed[
which(
onlyInBed$from %in% beids$id[which(beids$db.deprecated=="FALSE")]
),
]
onlyInBedDep <- onlyInBed[
which(onlyInBed$from %in% setdiff(onlyInBed$from, onlyInBedCur$from)),
]
head(onlyInBedDep[order(as.numeric(onlyInBedDep$to)),])
head(onlyInBedCur[order(as.numeric(onlyInBedCur$to)),])
sum(unique(onlyInBedCur$from) %in% beids$id[which(beids$preferred)])
fromNcbi <- bedCall(
cypher,
query=prepCql(c(
'MATCH (f:GeneID {database:"Ens_gene"}) WHERE f.value IN $from',
'MATCH (t:GeneID {database:"EntrezGene"})',
'MATCH (f)-[:corresponds_to]-(t)',
'RETURN f.value as from, t.value as to'
)),
parameters=list(from=as.list(unique(onlyInBedCur$from)))
)
all(
paste(fromNcbi$from, fromNcbi$to) %in%
paste(onlyInBedCur$from, onlyInBedCur$to)
)
##
fromHgnc <- bedCall(
cypher,
query=prepCql(c(
'MATCH (f:GeneID {database:"Ens_gene"}) WHERE f.value IN $from',
'MATCH (i:GeneID {database:"HGNC"})',
'MATCH (t:GeneID {database:"EntrezGene"})',
'MATCH (f)-[:corresponds_to]-(i)-[:corresponds_to]-(t)',
'RETURN f.value as from, t.value as to'
)),
parameters=list(from=as.list(setdiff(onlyInBedCur$from, fromNcbi$from)))
)
##
remaining <- setdiff(
onlyInBedCur$from, c(fromHgnc$from, fromNcbi$from)
)
The numbers of human Ensembl gene identifiers successfully converted by each
method are compared in figure \@ref(fig:vennTools).
Five identifiers were only converted by gProfileR.
They were provided by former versions of Ensembl or NCBI
but are now deprecated in the current releases of these two resources.
All the other gene identifiers converted by the different methods
were also converted by BED. However BED was able to map at least
r bigFormat(length(unique(onlyInBed$from))) more identifiers than all the
other tools (figure \@ref(fig:vennTools).a). A few of these mappings
(r bigFormat(length(unique(onlyInBedDep$from))))
are explained by the fact that BED is the only tool mapping
deprecated identifiers to current versions.
Nevertheless, even when focusing on the mapping of current versions of
Ensembl identifiers BED was able to
map r bigFormat(length(unique(onlyInBedCur$from)))
more identifiers than all the
other tools (figure \@ref(fig:vennTools).b).
A few of these mappings (r bigFormat(length(unique(fromNcbi$from))))
are directly provided by the NCBI. But most of them
(r bigFormat(length(unique(fromHgnc$from))))
are inferred from a mapping of the Ensembl and Entrez gene identifiers
to the same HGNC [@gray_genenames.org:_2015] identifier.
A rough approximation of running times of the different methods is provided
in table \@ref(tab:runTimes). The aim of this table is to show that BED,
as a dedicated and localy available tool,
is a very efficient option to convert large lists of identifiers on the fly and
recurrently.
The aim of BED is to improve the efficiency of identifier conversion
in a well defined context (organism, information resources of interest...)
and not to replace biomaRt, mygene, gProfileR
or other tools which provide many more features
for many organisms and which should not be narrowed to this task for a
complete comparison.
# runTimes <- c(
# "BED (Not cached)"="~10 seconds",
# "BED (Cached)"="~3 seconds",
# "biomaRt"="~1 minute",
# "mygene"="~4 minutes",
# "gProfileR"="~1 minute"
# )
runTimes <- data.frame(
"Method"=names(runTimes),
"Running time"=paste0(
"~",
unlist(lapply(runTimes, function(x) format(signif(x, 2))))
),
stringsAsFactors=FALSE,
check.names=FALSE
)
kable(
runTimes,
row.names=FALSE,
format='pandoc',
format.args=list(big.mark=","),
caption='Rough approximation of running time of different methods
to convert human Ensembl gene identifiers in human Entrez gene identifiers.'
)
The BED convBeIds function can be used to convert identifiers from
any available scope to any other one. It automaticaly find the most
relevant path according to the considered biological entities.
It allows elaborate mapping such as the conversion between
probe identifiers from a platform focused on mouse transcripts into
human protein identifiers.
Because such mappings can be intricated,
BED also provides a function to show the shortest relevant path between
two different identifiers (Figure \@ref(fig:explConv)).
```r relationship. This graph has been drawn with the \texttt{exploreConvPath} function.'}
include_graphics('img/ILMN_1220595-Conversion.png')
## Additional features
Some additional use cases and examples are provided in
the BED R package vignette.
Several functions are available for annotating BEID with symbols and names,
again taking advantage of information related to connected identifiers.
Other functions are also provided to seek relevant identifier of a specific
biological entity. These functions are used by a shiny [@chang_shiny:_2017]
gadget (Figure \@ref(fig:findBe))
providing an interactive dictionary of BEID which is also made
available as an Rstudio addin [@cheng_miniui:_2016; @allaire_rstudioapi:_2017].
```r Shiny gadget to seek relevant identifier of a specific biological entity. In this example the user is looking after human Ensembl transcript identifiers corresponding to \`\`il6\'\'.'}
include_graphics('img/findBe.png')
Summary
BED is a system dedicated to the mapping between identifiers of molecular
biological entities. It relies on a graph data model implemented with
\neo{} and on rules coded in an R package.
BED leverages mapping information provided by different resources in order
to increase the mapping efficiency between each of them.
It also allows the mapping of deprecated identifiers.
Rules are used to automatically convert identifiers from one scope to another
using the most appropriate path.
The intend of BED is to be tailored to specific needs and
beside functions for querying the system the BED R package provides functions
to build custom instances of the database.
Database instances can be localy installed or shared accross a community.
This design combined with a cache system makes BED performant for converting
large lists of identifiers from and to a large variety of scopes.
Because of our research field we provide an instance focused on human,
mouse and rat organisms. This database instance can be directly used
in relevant projects but it can also be enriched depending on user or
community needs.
Software availability
This section will be generated by the Editorial Office before publication. Authors are asked to provide some initial information to assist the Editorial Office, as detailed below.
- URL link to where the software can be downloaded from or used by a non-coder (AUTHOR TO PROVIDE; optional): NO
- URL link to the author's version control system repository containing the source code (AUTHOR TO PROVIDE; required):
- https://github.com/patzaw/BED
- https://github.com/patzaw/neo2R
- Link to source code as at time of publication (F1000Research TO GENERATE)
- Link to archived source code as at time of publication (F1000Research TO GENERATE)
- Software license (AUTHOR TO PROVIDE; required): GPL-3
Author contributions
PG designed and implemented the BED graph database,
the neo2R and BED R packages.
PG and JVE defined use cases, tested the whole system and
wrote the manuscript.
Competing interests
No competing interests were disclosed.
Grant information
This work was entirely supported by UCB Pharma.
The authors declared that no grants were involved in supporting this work.
Acknowledgments
We are grateful to Frédéric Vanclef, Malte Lucken, Liesbeth François,
Matthew Page,
Massimo de Francesco, and Marina Bessarabova for fruitful discussions
and constructive criticisms.
patzaw/BED documentation built on April 30, 2024, 5:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
reComp <- FALSE library(knitr) library(RColorBrewer) library(VennDiagram) library(gridExtra) flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger") opts_chunk$set( include=FALSE, warning=FALSE, echo=FALSE, message=FALSE, concordance=TRUE ) bigFormat <- function(x) {format(x, big.mark=",")} library(BED) connectToBed(url="localhost:5454", username="neo4j", password="1234")
\newcommand{\tm}{\textsuperscript{\textregistered}} \newcommand{\neo}{Neo4j\tm{}} \newcommand{\cypher}{Cypher\tm{}} \newcommand{\docker}{Docker\tm{}} \newcommand{\metabase}{MetaBase\tm{}}
Introduction
Since the advent of genome sequencing projects, many technologies have been developped to get access to different molecular information at a large scale and with high throughput. DNA microarrays are probably the archetype of such technology because of their historical impact on gathering data related to nucleic acids: genomic DNA and RNA. They triggered the emergence of "omics" fields of research such as genomics, epigenomics or transcriptomics. Lately massive parallel sequencing further increased the throughput of data generation related to nucleic acids by several orders of magnitude. In a different way mass spectrometry-related technologies allow the identification and the quantification of many kinds of molecular entities such as metabolites and proteins. Many information systems have been developped to manage the exploding amount of data and knowledge related to biological molecular entities. These resources manage different aspects of the knowledge. For example some are genome or proteome centered whereas other are focused on molecular interactions and pathways. Thus all these resources rely on different identifier systems to organise the concepts of interest. The value of all the experimental data and all the knowledge collected in public or private resources is very high as such but is also often synegistically leveraged by their cross comparison in a dedicated manner. Indeed many datasets can be relevant when adressing the understanding of a specific biological system, a phenotypic trait or a disease for example. These datasets can focus on different biological entities such as transcripts or proteins in different tissues, conditions or organisms. Comparing all these data and integrating them with available knowledge require the ability to map the identifiers on which each resource relies.
To achieve this task public and proprietary information systems provide mapping tables between their own identifiers and those from other resources. Furthermore many tools have been developped to facilitate the access to this information. Ensembl BioMarts [@kinsella_ensembl_2011], mygene [@wu_biogps_2013], and g:Profiler [@reimand_g:profiler-web_2016] are popular examples among many others. However, as pointed by @van_iersel_bridgedb_2010, these tools are generally dedicated to a particular domain not necessarly relevant or complete for all research project and keeping them up-to-date can also be an issue. Recognizing these challenges @van_iersel_bridgedb_2010 proposed the BridgeDb framework providing to bioinformatics developers a standard interface between tools and mapping services and also allowing the easy integration of custom data by a transitivity mechanism.
Here we present BED: a biological entity dictionary. BED has been developped to adress three main challenges. The first one is related to the completeness of identifier mappings. Indeed direct mapping information provided by the different systems are not always complete and can be enriched by mappings provided by other resources. More interestingly direct mappings not identified by any of these resources can be indireclty infered by using mappings to a third reference. For example many human Ensembl gene identifiers not directly mapped to any Entrez gene identifiers by neither Ensembl nor the NCBI can be indirectly mapped using mappings to HGNC identifiers. The second challenge is related to the mapping of deprecated identifiers. Indeed entity identifiers can change from one resource release to another. The identifier history is provided by some resources, such as Ensembl or the NCBI, but it is generally not used by mapping tools. The third challenge is related to the automation of the mapping process according to the relationships between the biological entities of interest. Indeed mapping between gene and protein identifier scopes should not be done the same way than two scopes of gene identifiers. Also converting identifiers from different organism should be possible using gene ortholog information.
To meet these challenges we designed a graph data model describing possible relationships between different biological entities and their identifiers. This data model has been implemented with the \neo{} graph database [@neo4j_inc_neo4j_2017] and conversion rules have been defined and coded in an R [@r_core_team_r:_2017] package. We provide an instance of the BED database focused on human, mouse and rat organism but many functions are available to construct other instances tailored to other needs.
Methods
Data model
include_graphics('img/BED-Data-Model.pdf')
The BED (Biological Entity Dictionary) system relies on a data model inspired by the central dogma of molecular biology [@crick_central_1970] and describing relationships between molecular concepts usually manipulated in the frame of genomics studies (Figure \@ref(fig:Data-Model)). A biological entity identifier (BEID) can identify either a Gene (GeneID), a Transcript (TranscriptID), a Peptide (PeptideID) or an Object (ObjectID). Object entities can correspond to complex concepts coded by any number of genes (i.e. a protein complex or a molecular function). BEID are extracted from public or private databases (BEDB). BEDB can provide an Attribute related to each BEID. For example it can be the sequencing region provided by the Ensembl database [@zerbino_ensembl_2017] or the identifier status provided by Uniprot [@the_uniprot_consortium_uniprot:_2017]. BEID can have one or several associated names (BENames) and symbols (BESymbol). GeneID can have one or several homologs in other organisms belonging to the same GeneIDFamily. Many genomics platforms, such as microarray, allow the identification of biological entity by using probes identified by ProbeID. In general BEID can be targeted by several probes belonging to a Platform which is focused on one and only one type of entity (BEType) among those described above: Gene, Transcript, Peptide or Object. A BEType can have several BEType products but can be the product of at most one BEType. This constraint allows the unambiguous identification of the most relevant path to convert identfiers from one scope to another and is fulfilled by the current data model: peptides are only produced from transcripts which are only produced from genes which can also code for objects.
BEID identifying the same biological entity are related through three different kinds of relationship according to the information available in the source databases and to the decision made by the database administrator about how to use them. Two BEID which corresponds_to each other both identify the same biological entity. A BEID which is_associated_to or which is_replaced_by another BEID does not directly identify any biological entity: the link is always indirect through one or several other BEID. Therefore by design a BEID which is_associated_to or which is_replaced_by another BEID can be related to several different biological entities. It is not the case for other BEID which identify one and only one biological entity. This set of possible relationship allows the indirect mapping of different identifiers not necessarily provided by any integrated resource.
In order to efficiently leverage indirect path through these different relationships the data model has been implemented in a \neo{} graph database [@neo4j_inc_neo4j_2017].
Feeding the database
Two R [@r_core_team_r:_2017] packages have been developed to feed and query the database. The first one, neo2R, provides low level functions to interact with \neo{}. The second R package, BED, provides functions to feed and query the BED \neo{} graph database according to the data model described above.
Many functions are provided within the package to build a tailored BED database instance. These functions are not exported in order not to mislead the user when querying the database (which is the expected most frequent usage of the system). An R markdown document showing how to build a BED database instance for human, mouse and rat organisms is provided within the package. It can be adapted to other organisms or needs.
Briefly these functions can be divided according to three main levels:
- The lowest level function is the
bedImportfunction which load a table in the \neo{} database according to a \cypher{} query. - Functions of the second level allow loading identifiers and relationships tables ensuring the integrity of the data model.
- Highest level functions are helpers for loading information provided by some public resources in different specific format.
Querying the database
The BED R package provides several functions to retrieve identifiers from different resources and also to convert identifiers from one reference to another. These functions generate and call \cypher{} queries on the \neo{} database. Converting thousands of identifiers can take some time (generally a few seconds). Also such conversions are often recurrent and redundant. In order to improve the performance for such recurrent and redundant queries, a cache system has been implemented. The first time, the query is run on \neo{} for all the relevant ID related to user input and the result is saved in a local file. Next time similar queries are requested, the system does not call \neo{} but loads the cached results and filters it according to user input. By default the cache is flushed when the system detects inconsistencies with the BED database. It can also be manually flushed if needed.
Operation
The graph database has been implemented with \neo{} version 3 [@neo4j_inc_neo4j_2017]. The BED R package depends on the following packages available in the Comprehensive R Archive Network [@cran_comprehensive_nodate]:
- visNetwork [@almende_b.v._visnetwork:_2017]
- dplyr [@wickham_dplyr:_2017]
- htmltools [@rstudio_inc_htmltools:_2017]
- DT [@xie_dt:_2016]
- shiny [@chang_shiny:_2017]
- miniUI [@cheng_miniui:_2016]
- rstudioapi [@allaire_rstudioapi:_2017]
Use Cases
Available database instance
An instance of the BED database (UCB-Human) has been built using the script provided in the BED R package and made available in a \docker{} image [@docker_inc_docker_2017] available here: https://hub.docker.com/r/patzaw/bed-ucb-human/
This instance used to exemplify the following use cases is focused on Homo sapiens, Mus musculus and Rattus norvegicus organisms and it has been built from the following resources:
- Ensembl [@zerbino_ensembl_2017]
- NCBI [@ncbi_resource_coordinators_database_2017]
- Uniprot [@the_uniprot_consortium_uniprot:_2017]
- biomaRt [@durinck_mapping_2009]
- GEOquery [@davis_geoquery:_2007]
- Clarivate Analytics \metabase{} [@clarivate_analytics_metacore_2017]
toShow <- c() for(be in listBe()){ for(org in BED::listOrganisms()){ toAdd <- listBeIdSources( be=be, organism=org, exclude=c("BEDTech_gene", "BEDTech_transcript"), verbose=FALSE ) colnames(toAdd) <- c("Database", "nbBE", "BEID", "BE") toAdd$Organism <- org toShow <- rbind( toShow, toAdd[,c("BE", "Organism", "Database", "BEID")] ) } } urls <- c() for(i in 1:nrow(toShow)){ urls <- c(urls, getBeIdURL("", toShow[i, "Database"])) } urls <- substr( urls, start=1, stop=regexpr("[^/][/][^/].*$", urls) ) # urls <- sub("^https*[:][/][/]", "", urls) toShow$URL <- urls kable( toShow, row.names=FALSE, format='pandoc', format.args=list(big.mark=","), caption='Numbers of BEID available in the BED UCB-Human database instance. Numbers have been split according to the BE type and the organism. Only BEID which can be mapped to each other are taken into account.' )
nbBeIds <- 3519181 nbProbeIds <- 354205
## Total number of BEID nbBeIds <- bedCall( cypher, query <- prepCql(c( 'MATCH (n:BEID)', #'-[:is_replaced_by|is_associated_to*0..]->(:BEID)-[:identifies]->()', 'WHERE NOT n.database IN ["BEDTech_gene", "BEDTech_transcript"]', 'RETURN count(distinct n);' )) ) nbProbeIds <- bedCall( cypher, query <- prepCql(c( 'MATCH (n:ProbeID)', 'RETURN count(distinct n);' )) )
The numbers of biological entity (BE) identifiers (BEID) available in this
BED database instance and which can be mapped to each other are shown
in table \@ref(tab:ID-Numbers).
In total, r bigFormat(nbBeIds) BEID are available in this
BED instance but many
of them cannot be mapped to other BEID because corresponding to deprecated
identifiers not corresponding to any identifier currently in use.
Also all the genomics platforms included in this BED database instance are
shown in table \@ref(tab:Platforms) . They provide mapping to BEID from
r bigFormat(nbProbeIds) ProbeID in total.
toShow <- listPlatforms() colnames(toShow) <- c("Name", "Description", "BE") toShow$Description <- sub(" [(].*$", "", toShow$Description) kable( toShow, row.names=FALSE, format='pandoc', caption='Genomics platforms available in the BED UCB-Human database instance.' )
Exploring identifiers of biological entities
The getBeIds function returns all BE identifiers from a specific scope.
A scopes is defined
by the type of BE or probe, the source of the identifiers (database or platform)
and the organism.
For example the following code
returns all the Ensembl identifiers of human genes.
beids <- getBeIds( be="Gene", source="Ens_gene", organism="human", restricted=FALSE ) head(beids)
redUnRestIds <- table(table(beids$Gene)) redUnRestIds <- redUnRestIds[order(as.numeric(names(redUnRestIds)))] ensVersion <- bedCall( cypher, query=prepCql(c( 'MATCH (n:BEDB {name:"Ens_gene"}) RETURN n.currentVersion' )) )
layout(matrix(c(1,1,1,1,2,2,2), 1, 7, byrow = TRUE)) par(mar=c(5.1, 4.5, 1.1, 0.5)) barplot( redUnRestIds, xlab="Number of BEID", ylab="Number of BE", log="y", las=2, cex.axis=0.7, cex.names=0.7 ) mtext(text="(a)", side=2, line=3, at=10^par()$usr[4], las=2, font=2) ## rbeids <- getBeIds( be="Gene", source="Ens_gene", organism="human", restricted = TRUE ) redIds <- table(table(rbeids$Gene)) redIds <- redIds[order(as.numeric(names(redIds)))] barplot( redIds, xlab="Number of BEID", ylab="Number of BE", log="y", las=2, cex.axis=0.7, cex.names=0.7 ) mtext(text="(b)", side=2, line=3, at=10^par()$usr[4], las=2, font=2)
The id column corresponds to the BEID from the source of interest.
The column named according to the BE type (in this case Gene)
corresponds to the internal identifiers of the related BE.
This internal identifier is not a stable reference that can be used as such.
Nevertheless it is useful to identify BEID identifying the
same BE.
In the example above even if most of Gene BE are identified by only
one Ensembl gene BEID, many of them are identified by two or more
(r bigFormat(sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=2)]))
/ r bigFormat(sum(redUnRestIds))
= r round(100*sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=2)])/sum(redUnRestIds)) %);
r sum(redUnRestIds[which(as.numeric(names(redUnRestIds))>=10)]) BE
are even identified by more than 10 Ensembl BEID
(Figure \@ref(fig:beidTables).a).
In this case, most of these redundancies come from deprecated ID from former
versions of the Ensembl database (version in used here: r ensVersion)
and can be excluded by setting the restricted parameter to TRUE when calling
the getBeIds function (Figure \@ref(fig:beidTables).b).
However many BE are still identified by two or more current Ensembl BEID
(r bigFormat(sum(redIds[which(as.numeric(names(redIds))>=2)]))
/ r bigFormat(sum(redIds))
= r round(100*sum(redIds[which(as.numeric(names(redIds))>=2)])/sum(redIds))~%).
This result comes from the way the BED database is constructed:
When two identifiers from the same resource correspond to the same identifier
in another resource (correspond_to relationship in the data model),
all these BEID are considered to identify the same BE.
A complex example of such mapping is shown in figure \@ref(fig:TAS2R8) mapping all the BEID of the human TAS2R8 gene which codes for a protein of the family of candidate taste receptors. There are three identifiers corresponding to this gene symbol in Ensembl. All these three identifiers correspond to the same Entrez gene and the same HGNC identifiers. All these BEID are thus considered to identify the same gene. It turns out that the three Ensembl BEID correspond to the same gene mapped on different sequence version of the chromosome 12: the canonical (ENSG00000121314), CHR_HSCHR12_2_CTG2 (ENSG00000272712) and CHR_HSCHR12_3_CTG2 (ENSG00000277316). This information provided by Ensembl is encoded in the seq_region attribute for each Ensembl BEID (see data model) and is used to define preferred BEID which are mapped on canonical version of chromosome sequences. The ENSG00000272712 identifier shows also a complex history in former Ensembl versions.
```r identifiers. Solid arrows represent \emph{correspond\_to} and \emph{is\_known\_as} relationships. Dotted arrows represent \emph{is\_replaced\_by} and \emph{is\_associated\_to} relationships. This graph has been drawn with the \texttt{exploreBe} function.'} include_graphics('img/TAS2R8-Identifiers.png')
## Converting identifiers
```r
## Save/load results from commmands below to not recompute everything
load("~/Tmp/ConvTables-3rdTools.rda")
The main goal of BED is to convert identifiers from one scope to another easily, rapidly and with high completeness. It has been thought in order to allow recurring comparisons to each other of many lists of biological entities from various origins.
The function guessIdOrigin can be used to guess the scope
of any list of identifiers.
A simple example regarding the conversion of human Ensembl
gene to human Entrez gene identifiers is shown below and discussed hereafter.
By setting the restricted parameter to TRUE the converted BEID are
restricted to current - non-deprecated - version of Entrez gene identifiers.
Nevertheless all the input BEID are taken into account,
current and deprecated ones.
runTimes <- list() ## clearBedCache() s <- Sys.time() bedConv <- convBeIds( ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human", to.source="EntrezGene", restricted=TRUE ) e <- Sys.time() runTimes[["BED (Not cached)"]] <- e-s ## s <- Sys.time() bedConv <- convBeIds( ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human", to.source="EntrezGene", restricted=TRUE ) e <- Sys.time() runTimes[["BED (Cached)"]] <- e-s
bedConv <- convBeIds( ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human", to.source="EntrezGene", restricted=TRUE )
curBeid <- unique(beids$id[which(beids$db.deprecated=="FALSE")]) depBeid <- unique(beids$id[which(beids$db.deprecated!="FALSE")])
Among all the r bigFormat(length(curBeid)+ length(depBeid))
human Ensembl gene identifiers available in the
database, r bigFormat(sum(is.na(bedConv$to)))
(r round(100*sum(is.na(bedConv$to))/length(unique(bedConv$from))) %)
were not converted to any human Entrez gene
identifier: r bigFormat(sum(is.na(bedConv$to[which(bedConv$from %in% curBeid)])))
(r round(100*sum(is.na(bedConv$to[which(bedConv$from %in% curBeid)]))/length(curBeid)) %)
of the r bigFormat(length(curBeid)) non-deprecated
and r bigFormat(sum(is.na(bedConv$to[which(bedConv$from %in% depBeid)])))
(r round(100*sum(is.na(bedConv$to[which(bedConv$from %in% depBeid)]))/length(depBeid)) %)
of the r bigFormat(length(depBeid)) deprecated identifiers.
## biomaRt library(biomaRt) bmEnsembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") s <- Sys.time() bmConv <- getBM( values = unique(beids$id), filters = 'ensembl_gene_id', attributes=c('ensembl_gene_id', 'entrezgene'), mart = bmEnsembl ) e <- Sys.time() runTimes[["biomaRt"]] <- e-s ## mygene library(mygene) s <- Sys.time() mgConv <- queryMany( qterm=unique(beids$id), scopes = "ensembl.gene", fields="entrezgene", species="human" ) e <- Sys.time() runTimes[["mygene"]] <- e-s ## gProfileR library(gProfileR) s <- Sys.time() gpConv <- gconvert( query=unique(beids$id), target="ENTREZGENE_ACC", organism="hsapiens", filter_na=FALSE, df=FALSE ) gpConv <- do.call( rbind, lapply( gpConv, function(x){ data.frame( alias=as.character(x[,"alias"]), target=as.character(x[,"target"]), stringsAsFactors=FALSE ) } ) ) gpConv$target <- sub("ENTREZGENE_ACC:", "", gpConv$target) gpConv$target <- ifelse(gpConv$target=="N/A", NA, gpConv$target) rownames(gpConv) <- c() gpConv <- unique(gpConv) e <- Sys.time() runTimes[["gProfileR"]] <- e-s
evalDate <- Sys.Date() save( bedConv, bmConv, mgConv, gpConv, evalDate, runTimes, file="~/Tmp/ConvTables-3rdTools.rda" )
Three other tools were used on r format(evalDate, "%b %d, %Y")
to perform the same conversion task:
biomaRt [@kinsella_ensembl_2011; @durinck_mapping_2009],
mygene [@wu_biogps_2013; @mark_mygene:_2014],
and gProfileR [@reimand_g:profiler-web_2016; @reimand_gprofiler:_2016].
At that time, biomaRt and mygene were based on the Ensembl 91 release
whereas gProfileR was based on release 90.
gq.all <- unique(beids$id) toComp.all <- list( "BED"=unique(bedConv$from[which(!is.na(bedConv$to))]), "biomaRt"=unique(bmConv$ensembl_gene_id[which(!is.na(bmConv$entrezgene))]), "mygene"=unique(mgConv$query[which(!is.na(mgConv$entrezgene))]), "gProfileR"=unique(gpConv$alias[which(!is.na(gpConv$target))]) ) gq.cur <- unique(beids$id[which(beids$db.deprecated=="FALSE")]) toComp.cur <- lapply( toComp.all, intersect, gq.cur )
venn.plot.all <- venn.diagram( toComp.all, NULL, fill=brewer.pal(4, "Accent"), alpha=rep(0.5, 4), main="(a)", main.pos=c(0.05,0.95), main.fontface=2, main.cex=1.8 ) venn.plot.cur <- venn.diagram( toComp.cur, NULL, fill=brewer.pal(4, "Accent"), alpha=rep(0.5, 4), margin=0.1, main="(b)", main.pos=c(0.05,0.95), main.fontface=2, main.cex=1.8 ) grid.newpage() grid.arrange( gTree(children=venn.plot.all, vp=viewport(width=0.8,height=0.9)), gTree(children=venn.plot.cur, vp=viewport(width=0.8,height=0.9)), ncol=2 )
onlyInBed <- bedConv[ which( !bedConv$from %in% unlist(toComp.all[-1]) & !is.na(bedConv$to) ), ] onlyInBedCur <- onlyInBed[ which( onlyInBed$from %in% beids$id[which(beids$db.deprecated=="FALSE")] ), ] onlyInBedDep <- onlyInBed[ which(onlyInBed$from %in% setdiff(onlyInBed$from, onlyInBedCur$from)), ] head(onlyInBedDep[order(as.numeric(onlyInBedDep$to)),]) head(onlyInBedCur[order(as.numeric(onlyInBedCur$to)),]) sum(unique(onlyInBedCur$from) %in% beids$id[which(beids$preferred)])
fromNcbi <- bedCall( cypher, query=prepCql(c( 'MATCH (f:GeneID {database:"Ens_gene"}) WHERE f.value IN $from', 'MATCH (t:GeneID {database:"EntrezGene"})', 'MATCH (f)-[:corresponds_to]-(t)', 'RETURN f.value as from, t.value as to' )), parameters=list(from=as.list(unique(onlyInBedCur$from))) ) all( paste(fromNcbi$from, fromNcbi$to) %in% paste(onlyInBedCur$from, onlyInBedCur$to) ) ## fromHgnc <- bedCall( cypher, query=prepCql(c( 'MATCH (f:GeneID {database:"Ens_gene"}) WHERE f.value IN $from', 'MATCH (i:GeneID {database:"HGNC"})', 'MATCH (t:GeneID {database:"EntrezGene"})', 'MATCH (f)-[:corresponds_to]-(i)-[:corresponds_to]-(t)', 'RETURN f.value as from, t.value as to' )), parameters=list(from=as.list(setdiff(onlyInBedCur$from, fromNcbi$from))) ) ## remaining <- setdiff( onlyInBedCur$from, c(fromHgnc$from, fromNcbi$from) )
The numbers of human Ensembl gene identifiers successfully converted by each
method are compared in figure \@ref(fig:vennTools).
Five identifiers were only converted by gProfileR.
They were provided by former versions of Ensembl or NCBI
but are now deprecated in the current releases of these two resources.
All the other gene identifiers converted by the different methods
were also converted by BED. However BED was able to map at least
r bigFormat(length(unique(onlyInBed$from))) more identifiers than all the
other tools (figure \@ref(fig:vennTools).a). A few of these mappings
(r bigFormat(length(unique(onlyInBedDep$from))))
are explained by the fact that BED is the only tool mapping
deprecated identifiers to current versions.
Nevertheless, even when focusing on the mapping of current versions of
Ensembl identifiers BED was able to
map r bigFormat(length(unique(onlyInBedCur$from)))
more identifiers than all the
other tools (figure \@ref(fig:vennTools).b).
A few of these mappings (r bigFormat(length(unique(fromNcbi$from))))
are directly provided by the NCBI. But most of them
(r bigFormat(length(unique(fromHgnc$from))))
are inferred from a mapping of the Ensembl and Entrez gene identifiers
to the same HGNC [@gray_genenames.org:_2015] identifier.
A rough approximation of running times of the different methods is provided in table \@ref(tab:runTimes). The aim of this table is to show that BED, as a dedicated and localy available tool, is a very efficient option to convert large lists of identifiers on the fly and recurrently. The aim of BED is to improve the efficiency of identifier conversion in a well defined context (organism, information resources of interest...) and not to replace biomaRt, mygene, gProfileR or other tools which provide many more features for many organisms and which should not be narrowed to this task for a complete comparison.
# runTimes <- c( # "BED (Not cached)"="~10 seconds", # "BED (Cached)"="~3 seconds", # "biomaRt"="~1 minute", # "mygene"="~4 minutes", # "gProfileR"="~1 minute" # ) runTimes <- data.frame( "Method"=names(runTimes), "Running time"=paste0( "~", unlist(lapply(runTimes, function(x) format(signif(x, 2)))) ), stringsAsFactors=FALSE, check.names=FALSE ) kable( runTimes, row.names=FALSE, format='pandoc', format.args=list(big.mark=","), caption='Rough approximation of running time of different methods to convert human Ensembl gene identifiers in human Entrez gene identifiers.' )
The BED convBeIds function can be used to convert identifiers from
any available scope to any other one. It automaticaly find the most
relevant path according to the considered biological entities.
It allows elaborate mapping such as the conversion between
probe identifiers from a platform focused on mouse transcripts into
human protein identifiers.
Because such mappings can be intricated,
BED also provides a function to show the shortest relevant path between
two different identifiers (Figure \@ref(fig:explConv)).
```r relationship. This graph has been drawn with the \texttt{exploreConvPath} function.'} include_graphics('img/ILMN_1220595-Conversion.png')
## Additional features Some additional use cases and examples are provided in the BED R package vignette. Several functions are available for annotating BEID with symbols and names, again taking advantage of information related to connected identifiers. Other functions are also provided to seek relevant identifier of a specific biological entity. These functions are used by a shiny [@chang_shiny:_2017] gadget (Figure \@ref(fig:findBe)) providing an interactive dictionary of BEID which is also made available as an Rstudio addin [@cheng_miniui:_2016; @allaire_rstudioapi:_2017]. ```r Shiny gadget to seek relevant identifier of a specific biological entity. In this example the user is looking after human Ensembl transcript identifiers corresponding to \`\`il6\'\'.'} include_graphics('img/findBe.png')
Summary
BED is a system dedicated to the mapping between identifiers of molecular biological entities. It relies on a graph data model implemented with \neo{} and on rules coded in an R package. BED leverages mapping information provided by different resources in order to increase the mapping efficiency between each of them. It also allows the mapping of deprecated identifiers. Rules are used to automatically convert identifiers from one scope to another using the most appropriate path.
The intend of BED is to be tailored to specific needs and beside functions for querying the system the BED R package provides functions to build custom instances of the database. Database instances can be localy installed or shared accross a community. This design combined with a cache system makes BED performant for converting large lists of identifiers from and to a large variety of scopes.
Because of our research field we provide an instance focused on human, mouse and rat organisms. This database instance can be directly used in relevant projects but it can also be enriched depending on user or community needs.
Software availability
This section will be generated by the Editorial Office before publication. Authors are asked to provide some initial information to assist the Editorial Office, as detailed below.
- URL link to where the software can be downloaded from or used by a non-coder (AUTHOR TO PROVIDE; optional): NO
- URL link to the author's version control system repository containing the source code (AUTHOR TO PROVIDE; required):
- https://github.com/patzaw/BED
- https://github.com/patzaw/neo2R
- Link to source code as at time of publication (F1000Research TO GENERATE)
- Link to archived source code as at time of publication (F1000Research TO GENERATE)
- Software license (AUTHOR TO PROVIDE; required): GPL-3
Author contributions
PG designed and implemented the BED graph database, the neo2R and BED R packages. PG and JVE defined use cases, tested the whole system and wrote the manuscript.
Competing interests
No competing interests were disclosed.
Grant information
This work was entirely supported by UCB Pharma. The authors declared that no grants were involved in supporting this work.
Acknowledgments
We are grateful to Frédéric Vanclef, Malte Lucken, Liesbeth François, Matthew Page, Massimo de Francesco, and Marina Bessarabova for fruitful discussions and constructive criticisms.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.