BED: a Biological Entity Dictionary based on a graph data model

Data model

The BED (Biological Entity Dictionary) system relies on a data model inspired by the central dogma of molecular biology (Crick, 1970) and describing relationships between molecular concepts usually manipulated in the frame of genomics studies (Figure 1). 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 et al., 2018) or the identifier status provided by Uniprot (The UniProt Consortium, 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 micro-arrays, allow the identification of biological entities 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.

Figure 1. The BED graph data model.

The model is shown as an Entity/Relationship (ER) diagram: entities correspond to graph nodes and relationships to graph edges. “ID” and “idx” indicate if the corresponding entity property is unique or indexed respectively. Some redundancies occur in this data model. Indeed some “value” properties are duplicated in upper case (“value_up”) in order to improve the performance of case-insensitive searches. Also, the database of a BEID node is provided as a property to ensure uniqueness of the couples of “database” and “value” properties. The same approach has been applied for the “platform” property of ProbeID nodes.

According to this data model, the scope of an identifier is defined by three features:

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 identifiers 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. Relationships corresponds_to and is_associated_to are used to map identifiers from different databases with different transitivity properties (as explained below), whereas is_replaced_by relationships are used to map deprecated identifiers from former versions of the same database. This set of possible relationships allows the indirect mapping of different identifiers not necessarily provided by any integrated resource.

Implementation

In order to efficiently leverage an indirect path through these different relationships, the data model has been implemented in a Neo4j^® graph database (Neo4j inc, 2017).

Two R (R Core Team, 2017) packages have been developed to feed and query the database. The first one, neo2R, provides low level functions to interact with Neo4j^® . The second R package, BED, provides functions to feed and query the BED Neo4j^® graph database according to the data model described above.

Feeding the database

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:

Available database instance

An instance of the BED database (UCB-Human), built using the script provided in the BED R package, is available in a Docker^® image (Docker inc, 2017) available here: https://hub.docker.com/r/patzaw/bed-ucb-human/ (tag 2018-04-30). This instance is focused on Homo sapiens, Mus musculus and Rattus norvegicus organisms and it has been built from these resources:

The results and the use cases described below were obtained and executed using this instance of the BED database.

The numbers of BEID available in this BED database instance and which can be mapped to each other are shown in Table 1. In total, 3,559,720 BEID are available in this BED instance. This number includes deprecated identifiers without successor and which therefore cannot be mapped to any other identifier. All the genomics platforms included in this BED database instance are shown in Table 2. They provide mapping to BEID from 354,205 ProbeID in total.

Transitivity management

In the context of mapping identifiers, transitivity is the inference of a cross-reference between A and C based on existing cross-references between A and B and between B and C. Depending on how biological entities are defined, transitivity is desirable or not. In BED the transitivity mechanism is managed by the two following relationships: corresponds_to and is_associated_to. On one hand, the corresponds_to relationships make the mapping transitive since two BEIDs which are connected through this kind of relationship are considered to identify the same BE. On the other hand, a BEID which is_associated_to another one does not automatically identify the same BE making this kind of relationship not available for transitive mappings. When the BED database is fed, the user chooses which relationship should be of type corresponds_to or of type is_associated_to. For example, in the instance described above, cross-references provided by Ensembl from Ensembl gene identifiers to Entrez, HGNC and Vega gene identifiers are considered as corresponds_to relationships whereas cross-references to miRbase, Unigene and OMIM are considered as is_associated_to relationships.

Querying the database

The BED R package provides functions to connect to the Neo4j^® database (connectToBed) and to directly query it (bedCall). Beside these two functions, others of higher level are provided to explore available data, to manage identifiers and to convert them from one scope to another.

Different functions can be used to explore the data model by listing the type of BE (listBe) or describing their relationships (firstCommonUpstreamBe). Available organisms, databases and platforms can also be retrieved using listOrganisms, listBeIdSources and listPlatforms functions.

All BEID from a specific scope can be obtained with the getBeIds function. A specific BEID and its relationships with others can be graphically explored with the exploreBe function. The use and the results of these two functions are exemplified in Results. Also, the functions guessIdOrigin and checkBeIds are used to guess and to check the scope of any list of identifiers. This set of functions is completed by the getBeIdSymbols, getBeIdNames, getGeneDescription and getBeIdDescription functions which provide different ways to annotate identifiers taking advantage of information related to connected BEID. Other functions, searchId and getRelevantIds, are also provided to seek relevant identifiers for a specific BE. These functions are used by a shiny (Chang et al., 2017) gadget (findBe) providing an interactive dictionary of BEID which is also made available as an Rstudio add-in (Allaire et al., 2017; Cheng, 2016).

As described above, the BED data model has been built to fulfill molecular biology processes in order to ensure the biological relevance of identifier mappings. The is_expressed_as and is_translated_in relationships correspond to the transcription and translation processes respectively whereas codes_for is a fuzzy relationship allowing the mapping of genes on objects not necessarily corresponding to the same kind of biological molecules. These processes are described in different databases with different levels of granularity. For example, Ensembl (Zerbino et al., 2018) provides possible transcripts for each gene specifying which of them is canonical. The getDirectProduct and getDirectOrigin functions allow the user to retrieve direct products or direct origins of such molecular biology processes.

The automatic conversion of identifiers from one scope to another is handled by the convBeIds function. This conversion process can be applied directly on lists of identifiers (convBeIdLists) or on data frames with an identifier column (convDfBeIds). Because such conversion can be intricate, the exploreConvPath function is provided to display the shortest relevant paths between two identifiers. These functions are also exemplified in Results and Use case.

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 Neo4j^® 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 Neo4j^® 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. If needed, it can also be manually flushed by using the lsBedCache and clearBedCache functions.

Operation

Minimal system requirements for running BED and neo2R R packages:

The graph database has been implemented with Neo4j^® version 3 (Neo4j inc, 2017). The BED R package depends on the following packages available in the Comprehensive R Archive Network (CRAN):

Management of identifiers

Identifiers (ID) in BED can identify a biological entity (BEID) or a probe (ProbeID). The getBeIds function returns all ID from a specific scope. As described in Methods a scope is defined by its type (the type of BE or “Probe”), its source (a database for BEID or a platform for ProbeID) and the organism to which it refers. 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)

##                    id preferred  Gene db.version db.deprecated
## 83452 ENSG00000276626      TRUE 65397         92         FALSE
## 83453 ENSG00000199595      TRUE 65401         92         FALSE
## 83454 ENSG00000201381     FALSE 65401         54      20090519
## 83468 ENSG00000200605      TRUE 65406         92         FALSE
## 83469 ENSG00000206757     FALSE 65406         54      20090519
## 83455 ENSG00000207395      TRUE 65407         92         FALSE

The id column corresponds to the BEID from the source of interest. The column named according to the type of BE (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 Gene BE are identified by only one Ensembl gene BEID, many of them are identified by two or more (6,031/59,901 ≈ 10%); 277 BE are even identified by more than 10 Ensembl BEID (Figure 2a). In this case, most of these redundancies come from deprecated BEID from former versions of the Ensembl database (version in use here: 92) and can be excluded by setting the restricted parameter to TRUE when calling the getBeIds function (Figure 2b). However, many BE are still identified by two or more current Ensembl BEID (2,944/59,901 ≈ 5%). This result comes from the way the BED database is constructed: When two BEID from the same resource correspond to the same BEID in another resource (correspond_to relationship in the data model), all these BEID are considered to identify the same BE.

Figure 2. Barplots showing the number (log scale) of gene BE identified by one or more Ensembl gene BEID.

a) All Ensembl gene ID. b) Current Ensembl gene ID (version 92).

A complex example of such mapping is shown in Figure 3 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 BEID 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 versions 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 versions of chromosome sequences. The ENSG00000272712 identifier shows also a complex history in former Ensembl versions.

Figure 3. BED relationships between all the different identifiers of the human TAS2R8 gene recorded in the database.

BEID are shown as circle and gene symbol in the rounded box. The color legend is shown to the left of the figure. BEID surrounded in bold correspond to preferred identifiers. Solid arrows represent correspond_to (between two nodes of identical shape) and is_known_as (between two nodes of different shapes) relationships. Dotted arrows represent is_replaced_by (between two nodes of identical color) and is_associated_to (between two nodes of different colors) relationships. This graph has been drawn with the exploreBe function.

Converting gene identifiers

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 BEID from various origins.

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.

bedConv <- convBeIds(
   ids=beids$id, from="Gene", from.source="Ens_gene", from.org="human",
   to.source="EntrezGene", restricted=TRUE
)

Among all the 69,056 human Ensembl gene identifiers available in the database, 22,056 (32%) were not converted to any human Entrez gene identifier: 21,416 (33%) of the 65,256 non-deprecated and 640 (17%) of the 3,800 deprecated identifiers.

In order to assess the improvement of completeness achieved by BED, we compared it to three other tools: biomaRt (Durinck et al., 2009; Kinsella et al., 2011), mygene (Mark et al., 2014; Wu et al., 2013), and gProfileR (Reimand et al., 2016a; Reimand et al., 2016b). All these tools were used on May 03, 2018 to perform the same conversion task. At that time, biomaRt was based on the Ensembl 92 release (as the BED database instance), mygene on release 91 and gProfileR on release 90.

The numbers of human Ensembl gene identifiers successfully converted by each method are compared in Figure 4. Mappings returned only by gProfileR or by mygene (33 + 94 + 76) are available in releases 90 and 91 of Ensembl respectively but not in release 92. They probably correspond to deprecated cross-references. Conversely, mappings returned by both BED and biomaRt but neither by gProfileR nor mygene (319) are available in release 92 of Ensembl but not in releases 91 or 90. All the gene identifiers successfully converted by biomaRt were also converted by BED. However, BED was able to map at least 18,080 more identifiers than all the other tools (Figure 4a). A few of these mappings (3,160) 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 14,920 more identifiers than all the other tools (Figure 4b). A few of these mappings (683) are directly provided by the NCBI. But most of them (14,237) are inferred from a mapping of the Ensembl and Entrez gene identifiers to the same HGNC (Gray et al., 2015) identifier.

Figure 4.

Venn diagrams showing the number of human Ensembl gene identifiers mapped to at least one human Entrez gene identifier by the different tested tools when focusing (a) on all 69,056 or (b) on current 65,256 BEID (Ensembl 92 release).

We assessed the validity of the mappings by comparing the location on the human genome (GRCh38) of Ensembl gene identifiers as reported by Ensembl and the location of the corresponding Entrez gene identifiers as reported by the NCBI. If the mapping between two identifiers is correct their location on the genome should be identical or highly similar (only gene identifiers located on canonical versions of chromosome sequences were considered for this comparison). We compared the following mapping results:

As shown in Figure 5a, the proportion of former cross-referenced identifiers (“Former”) located on different chromosomes is much higher than for the reference (“Reference”). The proportion of BED cross-referenced identifiers (“BED”) located on different chromosomes is higher than for “Reference” but still low. Figures 5.a and b show that inconsistencies in gene start location and in gene size are much more important in “Former” mappings than in “Reference” mappings. Such difference is not observed between the “Reference” mappings and the additional mappings provided by BED. These results show that the general quality of “Former” mappings is lower than the “Reference”, which is in agreement with their deprecation. More interestingly, additional mappings provided by BED are of similar quality to the “Reference” making their use as valuable.

Figure 5. Comparison of gene locations provided by Ensembl for Ensembl identifiers and by the NCBI for the cross-referenced Entrez identifiers.

“Reference” corresponds to mappings provided by biomaRt for identifiers successfully converted by all the tools. “Former” corresponds to mappings provided by mygene for identifiers successfully converted by mygene but neither by biomaRt nor BED. “BED” corresponds to mappings provided by BED for identifiers successfully converted by BED but not by any of the three other tools. Numbers of considered mappings (N) in each group are indicated (see text for details). a) Proportion of cross-referenced identifiers located on different chromosomes. b) Absolute distance (base pairs) between gene start positions for cross-referenced identifiers located on the same chromosome. c) Absolute difference (base pairs) between gene sizes for cross-referenced identifiers located on the same chromosome. Median and third quartile values are indicated. Outliers are not shown. “BED” and “Former” results are compared to Reference on different boxplots because of scale shift.

As shown above, additional mappings between Ensembl and Entrez gene identifiers inferred thanks to the use of HGNC cross-references by BED are many and of good quality. However, this transitivity mechanism is not always desirable depending on how different resources define biological entities. It is especially true for gene which is an unstable concept as described by Gerstein et al. (2007). For example, in Entrez the Hs.103110 UniGene identifier is mapped to 5465 and to 150383 Entrez gene identifiers which correspond to two different genes in Entrez (PPARA and CDPF1) but also in Ensembl (ENSG00000186951 and ENSG00000205643) and in HGNC (9232 and 33710). These two genes are located closely on the same chromosome but at different positions and on different directions. The same feature has been observed for many UniGene identifiers which makes UniGene unsuitable for transitivity mappings between gene identifiers. Therefore, in this instance of the BED database, as described in Methods, cross-references between Ensembl and Unigene gene identifiers are considered as is_associated_to relationships avoiding transitivity.

Mapping rules

Beside cross-referencing identifiers of identical type of BE, BED uses the biological relationship between genes, transcript and peptides to convert identifiers across different type of BE. For example, when converting peptides identifiers from the same species it uses only mapping done at the peptide level and does not use mapping to transcripts and genes. This strategy seems to be applied by biomaRt but not by mygene nor by gProfileR which map, for example, one Uniprot identifier to all the Ensembl peptide identifier coded by the same gene. For example, the A6NI28 Uniprot identifier is unambiguously mapped to the ENSP00000298815 Ensembl peptide identifier by BED and biomaRt but is wrongly mapped to three additional Ensembl peptide identifiers by mygene and gProfileR (ENSP00000431776, ENSP00000434304 and ENSP00000435961) which are encoded by the same gene (ENSG00000165895).

Furthermore, in biomaRt, mygene and gProfileR, mapping of BEIDs that are not genes from two different organisms using orthologs information requires at least two steps: one to find the ortholog genes and the other to find the relevant BEID. These two steps are integrated and transparent in BED.

In general, thanks to the BED data model and to mapping rules, identifier conversions from one scope to another are biologically relevant and automatic as exemplified in the use case described below.

Because all these indirect mappings can be intricate, BED provides a function to show the shortest relevant paths between two different identifiers (Figure 6).

Figure 6. BED conversion shortest path between the ILMN_1220595 probe identifier targeting a transcript of the mouse Il17a gene and the Uniprot Q16552 identifier of the human IL17 protein.

The legend is shown to the left of the figure. The red edge represents the is_homolog_of relationship. BEID surrounded in bold correspond to preferred identifiers. This graph has been drawn with the exploreConvPath function.

Performance

The time taken by BED to convert different sets of identifiers from one scope to another has been compared to the time taken by biomaRt, mygene and gProfileR. Three different queries have been executed starting from 100 or 20,000 identifiers.

As shown in Figure 7, BED outperforms mygene and gProfileR in all the tested cases even when not using its cache system. It also outperforms biomaRt for converting Affymetrix probe ID into Ensembl mouse peptide ID (Figures 7.e and 7.f). Without using its cache system, BED performs equally well as biomaRt for converting Ensembl human gene ID into Entrez gene ID (Figures 7.a and 7.b). It outperforms biomaRt for converting 100 Uniprot mouse peptide ID into Ensembl transcript ID (Figures 7.c) but is outperformed by biomaRt when 20,000 identifiers are considered (Figures 7.d). Nevertheless, when cache is used, BED outperforms all the tools whatever the conversion scenario (Figures 7.b, 7.d and 7.f). These results show that BED, as a dedicated and locally available tool, is a very efficient option to convert large lists of identifiers on the fly and recurrently.

Figure 7. Comparison of the time taken to convert a set of randomly selected BE identifiers.

Random selections of 100 or 20,000 identifiers have been executed 40 times before conversion by each considered tools. BED performance was assessed without using its cache system for all the queries and also using its cache system for queries regarding 20,000 identifiers. a) Conversion of 100 Ensembl human gene ID into Entrez gene ID. b) Conversion of 20,000 Ensembl human gene ID into Entrez gene ID. c) Conversion of 100 Uniprot mouse peptide ID into Ensembl transcript ID. d) Conversion of 20,000 Uniprot mouse peptide ID into Ensembl transcript ID. e) Conversion of 100 Affymetrix probe ID into Ensembl mouse peptide ID. f) Conversion of 20,000 Affymetrix probe ID into Ensembl mouse peptide ID.