Citing BED

This package and the underlying research has been published in this peer reviewed article:

Patrice Godard and Jonathan van Eyll (2018). BED: a Biological Entity Dictionary based on a graph data model (version 3; peer review: 2 approved). F1000Research, 7:195.

Installation

devtools::install_github("patzaw/BED")

Error: package or namespace load failed for ‘BED’:
 .onLoad failed in loadNamespace() for 'BED', details:
  call: connections[[connection]][["cache"]]
  error: subscript out of bounds

file.exists(file.path(Sys.getenv("HOME"), "R", "BED"))

Connection

Before using BED, the connection needs to be established with the underlying Neo4j DB. url, username and password should be adapted.

library(BED)

## [1] FALSE

connectToBed(url="localhost:5454", remember=FALSE, useCache=FALSE)

The remember parameter can be set to TRUE in order to save connection information that will be automatically used the next time the connectToBed() function is called. By default, this parameter is set to FALSE to comply with CRAN policies. Saved connection can be managed with the lsBedConnections() and the forgetBedConnection() functions.

The useCache parameter is by default set to FALSE to comply with CRAN policies. However, it is recommended to set it to TRUE to improve the speed of recurrent queries: the results of some large queries are saved locally in a file.

The connection can be checked the following way.

checkBedConn(verbose=TRUE)

If the verbose parameter is set to TRUE, the URL and the content version are displayed as messages.

The following function list saved connections.

lsBedConnections()

The connection param of the connectToBed function can be used to connect to a saved connection other than the last one.

Data model

The BED underlying data model can be shown at any time using the following command.

showBedDataModel()

Direct calls

Cypher queries can be run directly on the Neo4j database using the cypher function from the neo2R package through the bedCall function.

results <- bedCall(
    cypher,
    query=prepCql(
       'MATCH (n:BEID)',
       'WHERE n.value IN $values',
       'RETURN DISTINCT n.value AS value, labels(n), n.database'
    ),
    parameters=list(values=c("10", "100"))
)
results

Feeding the database

Many functions are provided within the package to build your own BED database instance. These functions are not exported in order to avoid their use when interacting with BED normally. Information about how to get an instance of the BED neo4j database is provided here:

• https://github.com/patzaw/BED#bed-database-instance-available-as-a-docker-image
• https://github.com/patzaw/BED#build-a-bed-database-instance

It can be adapted to user needs.

Caching

This part is relevant if the useCache parameter is set to TRUE when calling connectToBed().

Functions of the BED package used to retrieve thousands of identifiers can take some time (generally a few seconds) before returning a result. Thus for this kind of query, the query is run for all the relevant ID in the DB and thanks to a cache system implemented in the package same queries with different filters should be much faster the following times.

By default the cache is flushed when the system detect inconsistencies with the BED database. However, it can also be manualy flushed if needed using the clearBedCache() function.

Queries already in cache can be listed using the lsBedCache() function which also return the occupied disk space.

Biological entities

BED is organized around the central concept of Biological Entity (BE). All supported types of BE can be listed.

listBe()

These BE are organized according to how they are related to each other. For example a Gene is_expressed_as a Transcript. This organization allows to find the first upstream BE common to a set of BE.

firstCommonUpstreamBe(c("Object", "Transcript"))
firstCommonUpstreamBe(c("Peptide", "Transcript"))

Organisms

Several organims can be supported by the BED underlying database. They can be listed the following way.

listOrganisms()

Common names are also supported and the corresponding taxonomic identifiers can be retrieved. Conversely the organism names corresponding to a taxonomic ID can be listed.

getOrgNames(getTaxId("human"))

Identifiers of biological entities

The main aim of BED is to allow the mapping of identifiers from different sources such as Ensembl or Entrez. Supported sources can be listed the following way for each supported organism.

listBeIdSources(be="Transcript", organism="human")

The database gathering the largest number of BE of specific type can also be identified.

largestBeSource(be="Transcript", organism="human", restricted=TRUE)

Finally, the getAllBeIdSources() function returns all the source databases of BE identifiers whatever the BE type.

Experimental platforms and probes

BED also supports experimental platforms and provides mapping betweens probes and BE identifiers (BEID).

The supported platforms can be listed the following way. The getTargetedBe() function returns the type of BE on which a specific platform focus.

head(listPlatforms())
getTargetedBe("GPL570")

Retrieving all identifiers from a source

All identifiers of an organism BEs from one source can be retrieved.

beids <- getBeIds(
    be="Gene", source="EntrezGene", organism="human",
    restricted=FALSE
)
dim(beids)
head(beids)

The first column, id, corresponds to the identifiers of the BE in the source. The column named according to the BE type (in this case Gene) corresponds to the internal identifier of the related BE. BE CAREFUL, THIS INTERNAL ID IS NOT STABLE AND CANNOT BE USED AS A REFERENCE. This internal identifier is useful to identify BEIDS corresponding to the same BE. The following code can be used to have an overview of such redundancy.

sort(table(table(beids$Gene)), decreasing = TRUE)
ambId <- sum(table(table(beids$Gene)[which(table(beids$Gene)>=10)]))

In the example above we can see that most of Gene BE are identified by only one EntrezGene ID. However many of them are identified by two or more ID; XXX BE are even identified by 10 or more EntrezGeneID. In this case, most of these redundancies come from ID history extracted from Entrez. Legacy ID can be excluded from the retrieved ID by setting the restricted parameter to TRUE.

beids <- getBeIds(
    be="Gene", source="EntrezGene", organism="human",
    restricted = TRUE
)
dim(beids)

The same code as above can be used to identify remaining redundancies.

sort(table(table(beids$Gene)), decreasing = TRUE)

In the example above we can see that allmost all Gene BE are identified by only one EntrezGene ID. However some of them are identified by two or more ID. This result comes from how the BED database is constructed according to the ID mapping provided by the different source databases. The graph below shows how the mapping was done for such a BE with redundant EntrezGene IDs.

This issue has been mainly solved by not taking into account ambigous mappings between NCBI Entrez gene identifiers and Ensembl gene identifier provided by Ensembl. It has been achieved using the cleanDubiousXRef() function from the 2019.10.11 version of the BED-UCB-Human database.

eid <- beids$id[which(beids$Gene %in% names(which(table(beids$Gene)>=3)))][1]
print(eid)
exploreBe(id=eid, source="EntrezGene", be="Gene") %>%
   visPhysics(solver="repulsion")

The way the ID correspondances are reported in the different source databases leads to this mapping ambiguity which has to be taken into account when comparing identifiers from different databases.

The getBeIds() returns other columns providing additional information about the id. The same function can be used to retrieved symbols or probe identifiers.

mapt <- convBeIds(
   "MAPT", from="Gene", from.source="Symbol", from.org="human",
   to.source="Ens_gene", restricted=TRUE
)
exploreBe(
   mapt[1, "to"],
   source="Ens_gene",
   be="Gene"
)
getBeIds(
   be="Gene", source="Ens_gene", organism="human",
   restricted=TRUE,
   attributes=listDBAttributes("Ens_gene"),
   filter=mapt$to
)

Checking identifiers

The origin of identifiers can be guessed as following.

oriId <- c(
    "17237", "105886298", "76429", "80985", "230514", "66459",
    "93696", "72514", "20352", "13347", "100462961", "100043346",
    "12400", "106582", "19062", "245607", "79196", "16878", "320727",
    "230649", "66880", "66245", "103742", "320145", "140795"
)
idOrigin <- guessIdScope(oriId)
print(idOrigin$be)
print(idOrigin$source)
print(idOrigin$organism)

The best guess is returned as a list but other possible origins are listed in the details attribute.

print(attr(idOrigin, "details"))

If the origin of identifiers is already known, it can also be tested.

checkBeIds(ids=oriId, be="Gene", source="EntrezGene", organism="mouse")

checkBeIds(ids=oriId, be="Gene", source="HGNC", organism="human")

Identifier annotation

Identifiers can be annotated with symbols and names according to available information. The following code returns the most relevant symbol and the most relevant name for each ID. Source URL can also be generated with the getBeIdURL() function.

toShow <- getBeIdDescription(
    ids=oriId, be="Gene", source="EntrezGene", organism="mouse"
)
toShow$id <- paste0(
    sprintf(
        '<a href="%s" target="_blank">',
        getBeIdURL(toShow$id, "EntrezGene")
    ),
    toShow$id,
    '<a>'
)
kable(toShow, escape=FALSE, row.names=FALSE)

All possible symbols and all possible names for each ID can also be retrieved using the following functions.

res <- getBeIdSymbols(
    ids=oriId, be="Gene", source="EntrezGene", organism="mouse",
    restricted=FALSE
)
head(res)

res <- getBeIdNames(
    ids=oriId, be="Gene", source="EntrezGene", organism="mouse",
    restricted=FALSE
)
head(res)

Also probes and some biological entities do not have directly associated symbols or names. These elements can also be annotated according to information related to relevant genes.

someProbes <- c(
    "238834_at", "1569297_at", "213021_at", "225480_at",
    "216016_at", "35685_at", "217969_at", "211359_s_at"
)
toShow <- getGeneDescription(
    ids=someProbes, be="Probe", source="GPL570", organism="human"
)
kable(toShow, escape=FALSE, row.names=FALSE)

Products of molecular biology processes

The BED data model has beeing built to fulfill molecular biology processes:

• is_expressed_as relationships correspond to the transcription process.
• is_translated_in relationships correspond to the translation process.
• codes_for is a fuzzy relationship allowing the mapping of genes on object not necessary corresonpding to the same kind of biological molecule.

These processes are described in different databases with different level of granularity. For exemple, Ensembl provides possible transcripts for each gene specifying which one of them is canonical.

The following functions are used to retrieve direct products or direct origins of molecular biology processes.

getDirectProduct("ENSG00000145335", process="is_expressed_as")
getDirectProduct("ENST00000336904", process="is_translated_in")
getDirectOrigin("NM_001146055", process="is_expressed_as")

Same entity and same organism: from one source to another

res <- convBeIds(
    ids=oriId,
    from="Gene",
    from.source="EntrezGene",
    from.org="mouse",
    to.source="Ens_gene",
    restricted=TRUE,
    prefFilter=TRUE
)
head(res)

Same organism: from one entity to another

res <- convBeIds(
    ids=oriId,
    from="Gene",
    from.source="EntrezGene",
    from.org="mouse",
    to="Peptide",
    to.source="Ens_translation",
    restricted=TRUE,
    prefFilter=TRUE
)
head(res)

From one organism to another

res <- convBeIds(
    ids=oriId,
    from="Gene",
    from.source="EntrezGene",
    from.org="mouse",
    to="Peptide",
    to.source="Ens_translation",
    to.org="human",
    restricted=TRUE,
    prefFilter=TRUE
)
head(res)

Converting lists of identifiers

List of identifiers can be converted the following way. Only converted IDs are returned in this case.

humanEnsPeptides <- convBeIdLists(
    idList=list(a=oriId[1:5], b=oriId[-c(1:5)]),
    from="Gene",
    from.source="EntrezGene",
    from.org="mouse",
    to="Peptide",
    to.source="Ens_translation",
    to.org="human",
    restricted=TRUE,
    prefFilter=TRUE
)
unlist(lapply(humanEnsPeptides, length))
lapply(humanEnsPeptides, head)

entrezGenes <- BEIDList(
   list(a=oriId[1:5], b=oriId[-c(1:5)]),
   scope=list(be="Gene", source="EntrezGene", organism="Mus musculus"),
   metadata=data.frame(
      .lname=c("a", "b"),
      description=c("Identifiers in a", "Identifiers in b"),
      stringsAsFactors=FALSE
   )
)
entrezGenes
entrezGenes$a
ensemblGenes <- focusOnScope(entrezGenes, source="Ens_gene")
ensemblGenes$a

Converting data frames

IDs in data frames can also be converted.

toConv <- data.frame(a=1:25, b=runif(25))
rownames(toConv) <- oriId
res <- convDfBeIds(
    df=toConv,
    from="Gene",
    from.source="EntrezGene",
    from.org="mouse",
    to.source="Ens_gene",
    restricted=TRUE,
    prefFilter=TRUE
)
head(res)

Explore convertion shortest path between two identifiers

Because the conversion process takes into account several resources, it might be useful to explore the path between two identifiers which have been mapped. This can be achieved by the exploreConvPath function.

from.id <- "ILMN_1220595"
res <- convBeIds(
   ids=from.id, from="Probe", from.source="GPL6885", from.org="mouse",
   to="Peptide", to.source="Uniprot", to.org="human",
   prefFilter=TRUE
)
res
exploreConvPath(
   from.id=from.id, from="Probe", from.source="GPL6885",
   to.id=res$to[1], to="Peptide", to.source="Uniprot"
)

The figure above shows how the XXX ProbeID, targeting the mouse NM_010552 transcript, can be associated to the XXX human protein ID in Uniprot.

Notes about converting from and to gene symbols

Canonical and non-canonical symbols are associated to genes. In some cases the same symbol (canonical or not) can be associated to several genes. This can lead to ambiguous mapping. The strategy to apply for such mapping depends on the aim of the user and his knowledge about the origin of the symbols to consider.

The complete mapping between Ensembl gene identifiers and symbols is retrieved by using the getBeIDSymbolTable function.

compMap <- getBeIdSymbolTable(
   be="Gene", source="Ens_gene", organism="rat",
   restricted=FALSE
)
dim(compMap)
head(compMap)

The canonical field indicates if the symbol is canonical for the identifier. The direct field indicates if the symbol is directly associated to the identifier or indirectly through a relationship with another identifier.

As an example, let’s consider the “Snca” symbol in rat. As shown below, this symbol is associated to 2 genes; it is canonical for one gene and not for another. These 2 genes are also associated to other symbols.

sncaEid <- compMap[which(compMap$symbol=="Snca"),]
sncaEid
compMap[which(compMap$id %in% sncaEid$id),]

The getBeIdDescription function described before, reports only one symbol for each identifier. Canonical and direct symbols are prioritized.

getBeIdDescription(
   sncaEid$id,
   be="Gene", source="Ens_gene", organism="rat"
)

The convBeIds works differently in order to provide a mapping as exhaustive as possible. If a symbol is associated to several input identifiers, non-canonical associations with this symbol are removed if a canonical association exists for any other identifier. This can lead to inconsistent results, depending on the user input, as show below.

convBeIds(
   sncaEid$id[1],
   from="Gene", from.source="Ens_gene", from.org="rat",
   to.source="Symbol"
)
convBeIds(
   sncaEid$id[2],
   from="Gene", from.source="Ens_gene", from.org="rat",
   to.source="Symbol"
)
convBeIds(
   sncaEid$id,
   from="Gene", from.source="Ens_gene", from.org="rat",
   to.source="Symbol"
)

In the example above, when the query is run for each identifier independently, the association to the “Snca” symbol is reported for both. However, when running the same query with the 2 identifiers at the same time, the “Snca” symbol is reported only for one gene corresponding to the canonical association. An additional filter can be used to only keep canonical symbols:

convBeIds(
   sncaEid$id,
   from="Gene", from.source="Ens_gene", from.org="rat",
   to.source="Symbol",
   canonical=TRUE
)

Finally, as shown below, when running the query the other way, “Snca” is only associated to the gene for which it is the canonical symbol.

convBeIds(
   "Snca",
   from="Gene", from.source="Symbol", from.org="rat",
   to.source="Ens_gene"
)

Therefore, the user should chose the function to use with care when needing to convert from or to gene symbol.