In patzaw/BED: Biological Entity Dictionary (BED)
Introduction
This document shows how to feed the Biological Entity Dictionary (BED).
It can be adapted according to specific needs and DB access.
The BED functions used to feed the DB are not exported to avoid
unintended modifications of the DB. To call them, they are preceded
by BED:::.
In this example several source databases are dumped and their content
is integrated in BED.
Some helper functions are provided to get information from famous databases.
The following chunk is used to configure source versions.
The reDumpThr object is used to define time intervals during which some
data sources should not be re-downloaded.
##
library(knitr)
library(BED)
library(jsonlite)
library(rvest)
##
config <- jsonlite::read_json("build_config.json")
config <- lapply(
config, function(x){
if(!is.character(x)){
return(x)
}else{
sub(pattern="___HOME___", replacement=Sys.getenv("HOME"), x = x)
}
}
)
config <- lapply(
config, function(x){
if(!is.character(x)){
return(x)
}else{
sub(pattern="___ROOT___", replacement=config$ROOT, x = x)
}
}
)
##
wd <- config$BED_WORKING
##
# opts_knit$set(root.dir=wd)
opts_chunk$set(
eval=TRUE,
message=FALSE
# root.dir=wd
)
## Specific config
bedInstance <- config$BED_INSTANCE
bedVersion <- config$BED_VERSION
## Identify the latest ensembl release here: https://www.ensembl.org/
ensembl_org_dir <- grep(
"_core_",
rvest::html_table(
rvest::read_html("https://ftp.ensembl.org/pub/current/mysql/")
)[[1]]$Name,
value=TRUE
)
ensembl_org_table <- do.call(rbind, lapply(
strsplit(ensembl_org_dir, split="_"),
function(x){
gv <- sub("/", "", x[length(x)])
e <- x[length(x)-1]
core <- which(x=="core")
org <- paste(x[1:(core-1)], collapse=" ")
substr(org, 1, 1) <- toupper(substr(org, 1,1))
return(data.frame(
organism=org,
release=e,
genome_version=gv
))
}
))
ensembl_org_table <- ensembl_org_table[
order(as.numeric(ensembl_org_table$genome_version), decreasing=T),
]
ensembl_org_table <- ensembl_org_table[
which(!duplicated(ensembl_org_table$organism)),
]
ensembl_release <- unique(ensembl_org_table$release)
stopifnot(length(ensembl_release)==1)
org <-"Homo sapiens"
ensembl_Hsapiens <- list(
release=ensembl_release,
organism=org,
gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)],
gdbCref=c( # Gene cross-references DBs
"HGNC"="HGNC",
"EntrezGene"="EntrezGene",
"Vega_gene"="Vega_gene",
"Ens_Hs_gene"="Ens_gene"
),
gdbAss=c( # Gene associated IDs (DB)
"miRBase"="miRBase",
"MIM_GENE"="MIM_GENE",
"UniGene"="UniGene"
),
tdbCref=c( # Transcript cross-references DBs
"RefSeq_mRNA"="RefSeq",
"RefSeq_ncRNA"="RefSeq",
"RefSeq_mRNA_predicted"="RefSeq",
"RefSeq_ncRNA_predicted"="RefSeq",
"Vega_transcript"="Vega_transcript",
"Ens_Hs_transcript"="Ens_transcript"
),
pdbCref=c( # Peptide cross-references DBs
"RefSeq_peptide"="RefSeq_peptide",
"RefSeq_peptide_predicted"="RefSeq_peptide",
"Uniprot/SPTREMBL"="Uniprot",
"Uniprot/SWISSPROT"="Uniprot",
"Vega_translation"="Vega_translation",
"Ens_Hs_translation"="Ens_translation"
),
canChromosomes=c(1:22, "X", "Y", "MT")
)
org <- "Mus musculus"
ensembl_Mmusculus <- list(
release=ensembl_release,
organism=org,
gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)],
gdbCref=c( # Gene cross-references DBs
"MGI"="MGI",
"EntrezGene"="EntrezGene",
"Vega_gene"="Vega_gene",
"Ens_Mm_gene"="Ens_gene"
),
gdbAss=c( # Gene associated IDs (DB)
"miRBase"="miRBase",
"UniGene"="UniGene"
),
tdbCref=c( # Transcript cross-references DBs
"RefSeq_mRNA"="RefSeq",
"RefSeq_ncRNA"="RefSeq",
"RefSeq_mRNA_predicted"="RefSeq",
"RefSeq_ncRNA_predicted"="RefSeq",
"Vega_transcript"="Vega_transcript",
"Ens_Mm_transcript"="Ens_transcript"
),
pdbCref=c( # Peptide cross-references DBs
"RefSeq_peptide"="RefSeq_peptide",
"RefSeq_peptide_predicted"="RefSeq_peptide",
"Uniprot/SPTREMBL"="Uniprot",
"Uniprot/SWISSPROT"="Uniprot",
"Vega_translation"="Vega_translation",
"Ens_Mm_translation"="Ens_translation"
),
canChromosomes=c(1:19, "X", "Y", "MT")
)
org <- "Rattus norvegicus"
ensembl_Rnorvegicus <- list(
release=ensembl_release,
organism=org,
gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)],
gdbCref=c( # Gene cross-references DBs
"RGD"="RGD",
"EntrezGene"="EntrezGene",
"Vega_gene"="Vega_gene",
"Ens_Rn_gene"="Ens_gene"
),
gdbAss=c( # Gene associated IDs (DB)
"miRBase"="miRBase",
"UniGene"="UniGene"
),
tdbCref=c( # Transcript cross-references DBs
"RefSeq_mRNA"="RefSeq",
"RefSeq_ncRNA"="RefSeq",
"RefSeq_mRNA_predicted"="RefSeq",
"RefSeq_ncRNA_predicted"="RefSeq",
"Vega_transcript"="Vega_transcript",
"Ens_Rn_transcript"="Ens_transcript"
),
pdbCref=c( # Peptide cross-references DBs
"RefSeq_peptide"="RefSeq_peptide",
"RefSeq_peptide_predicted"="RefSeq_peptide",
"Uniprot/SPTREMBL"="Uniprot",
"Uniprot/SWISSPROT"="Uniprot",
"Vega_translation"="Vega_translation",
"Ens_Rn_translation"="Ens_translation"
),
canChromosomes=c(1:20, "X", "Y", "MT")
)
org <- "Sus scrofa"
ensembl_Sscrofa <- list(
release=ensembl_release,
organism=org,
gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)],
gdbCref=c( # Gene cross-references DBs
"EntrezGene"="EntrezGene",
"Vega_gene"="Vega_gene",
"Ens_Ss_gene"="Ens_gene"
),
gdbAss=c( # Gene associated IDs (DB)
"miRBase"="miRBase",
"UniGene"="UniGene"
),
tdbCref=c( # Transcript cross-references DBs
"RefSeq_mRNA"="RefSeq",
"RefSeq_ncRNA"="RefSeq",
"RefSeq_mRNA_predicted"="RefSeq",
"RefSeq_ncRNA_predicted"="RefSeq",
"Vega_transcript"="Vega_transcript",
"Ens_Ss_transcript"="Ens_transcript"
),
pdbCref=c( # Peptide cross-references DBs
"RefSeq_peptide"="RefSeq_peptide",
"RefSeq_peptide_predicted"="RefSeq_peptide",
"Uniprot/SPTREMBL"="Uniprot",
"Uniprot/SWISSPROT"="Uniprot",
"Vega_translation"="Vega_translation",
"Ens_Ss_translation"="Ens_translation"
),
canChromosomes=c(1:18, "X", "Y", "MT")
)
org <- "Danio rerio"
ensembl_Drerio <- list(
release=ensembl_release,
organism=org,
gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)],
gdbCref=c( # Gene cross-references DBs
"EntrezGene"="EntrezGene",
"ZFIN_ID"="ZFIN_gene",
"Vega_gene"="Vega_gene",
"Ens_Dr_gene"="Ens_gene"
),
gdbAss=c( # Gene associated IDs (DB)
"miRBase"="miRBase",
"UniGene"="UniGene"
),
tdbCref=c( # Transcript cross-references DBs
"RefSeq_mRNA"="RefSeq",
"RefSeq_ncRNA"="RefSeq",
"RefSeq_mRNA_predicted"="RefSeq",
"RefSeq_ncRNA_predicted"="RefSeq",
"Vega_transcript"="Vega_transcript",
"Ens_Dr_transcript"="Ens_transcript"
),
pdbCref=c( # Peptide cross-references DBs
"RefSeq_peptide"="RefSeq_peptide",
"RefSeq_peptide_predicted"="RefSeq_peptide",
"Uniprot/SPTREMBL"="Uniprot",
"Uniprot/SWISSPROT"="Uniprot",
"Vega_translation"="Vega_translation",
"Ens_Dr_translation"="Ens_translation"
),
canChromosomes=c(1:25, "MT")
)
## General config
reDumpThr <- as.difftime(config$SOURCE_REDUMP_THR, units="days")
curDate <- Sys.Date()
BED initialization
BED is based on Neo4j.
The S01-NewBED-Container.sh shows how to run it in
a docker container.
Because the import functions use massively the LOAD CSV Neo4j query,
the feeding of the BED database can only be down from the
computer hosting the Neo4j relevant instance.
The chunk below shows how to connect to BED. In this example,
neo4j authentication is disabled.
connectToBed(
url=sprintf("localhost:%s", config$NJ_HTTP_PORT),
remember=FALSE,
useCache=TRUE,
importPath=config$BED_IMPORT
)
Check empty DB
Do not go further if your BED DB is not empty.
dbSize <- bedCall(cypher, 'MATCH (n) RETURN count(n)')[,1]
if(dbSize!=0){
stop("BED DB is not empty ==> clean it before loading the content below")
}
Set BED instance and version
print(bedInstance)
print(bedVersion)
BED:::setBedVersion(bedInstance=bedInstance, bedVersion=bedVersion)
Load Data model
Start: r Sys.time()
BED:::loadBedModel()
End: r Sys.time()
Loading taxonomy from NCBI
Information is downloaded if older than r reDumpThr days
according to the reDumpThr object.
Start: r Sys.time()
BED:::loadNcbiTax(
reDumpThr=reDumpThr,
ddir=wd,
orgOfInt=c(
"Homo sapiens", "Rattus norvegicus", "Mus musculus",
"Sus scrofa", "Danio rerio"
),
curDate=curDate
)
End: r Sys.time()
Loading data from Ensembl
Register Ensembl DBs
Genes
BED:::registerBEDB(
name="Ens_gene",
description="Ensembl gene",
currentVersion=ensembl_release,
idURL='http://www.ensembl.org/id/%s'
)
Transcripts
BED:::registerBEDB(
name="Ens_transcript",
description="Ensembl transcript",
currentVersion=ensembl_release,
idURL='http://www.ensembl.org/id/%s'
)
Peptides
BED:::registerBEDB(
name="Ens_translation",
description="Ensembl peptides",
currentVersion=ensembl_release,
idURL='http://www.ensembl.org/id/%s'
)
Danio rerio
ensembl <- ensembl_Drerio
print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$gdbCref,
dbAss=ensembl$gdbAss,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$tdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$pdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Homo sapiens
ensembl <- ensembl_Hsapiens
print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$gdbCref,
dbAss=ensembl$gdbAss,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$tdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$pdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Mus musculus
ensembl <- ensembl_Mmusculus
print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$gdbCref,
dbAss=ensembl$gdbAss,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$tdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$pdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Rattus norvegicus
ensembl <- ensembl_Rnorvegicus
print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$gdbCref,
dbAss=ensembl$gdbAss,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$tdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$pdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Sus scrofa
ensembl <- ensembl_Sscrofa
print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$gdbCref,
dbAss=ensembl$gdbAss,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$tdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds(
organism=ensembl$organism,
release=ensembl$release,
gv=ensembl$gv,
ddir=wd,
dbCref=ensembl$pdbCref,
canChromosomes=ensembl$canChromosomes
)
gc()
End: r Sys.time()
Loading data from NCBI
Information is downloaded if older than r reDumpThr days
according to the reDumpThr object.
Register NCBI DBs
BED:::dumpNcbiDb(
taxOfInt = c(), reDumpThr=reDumpThr,
ddir=wd,
toLoad=c(), curDate=curDate
)
Genes
BED:::registerBEDB(
name="EntrezGene",
description="NCBI gene",
currentVersion=format(dumpDate, "%Y%m%d"),
idURL='https://www.ncbi.nlm.nih.gov/gene/%s'
)
Transcripts
BED:::registerBEDB(
name="RefSeq",
description="NCBI nucleotide",
currentVersion=format(dumpDate, "%Y%m%d"),
idURL='https://www.ncbi.nlm.nih.gov/nuccore/%s'
)
Peptides
BED:::registerBEDB(
name="RefSeq_peptide",
description="NCBI protein",
currentVersion=format(dumpDate, "%Y%m%d"),
idURL='https://www.ncbi.nlm.nih.gov/protein/%s'
)
Danio rerio data
Start: r Sys.time()
BED:::getNcbiGeneTransPep(
organism="Danio rerio",
ddir=wd,
curDate=curDate
)
gc()
End: r Sys.time()
Homo sapiens data
Start: r Sys.time()
BED:::getNcbiGeneTransPep(
organism="Homo sapiens",
ddir=wd,
curDate=curDate
)
gc()
End: r Sys.time()
Mus musculus data
Start: r Sys.time()
BED:::getNcbiGeneTransPep(
organism="Mus musculus",
ddir=wd,
curDate=curDate
)
gc()
End: r Sys.time()
Rattus norvegicus data
Start: r Sys.time()
BED:::getNcbiGeneTransPep(
organism="Rattus norvegicus",
ddir=wd,
curDate=curDate
)
gc()
End: r Sys.time()
Sus scrofa data
Start: r Sys.time()
BED:::getNcbiGeneTransPep(
organism="Sus scrofa",
ddir=wd,
curDate=curDate
)
gc()
End: r Sys.time()
Direct cross-references with Uniprot
Start: r Sys.time()
message("Direct cross-references with Uniprot")
BED:::dumpNcbiDb(
taxOfInt="",
reDumpThr=Inf,
ddir=wd,
toLoad="gene_refseq_uniprotkb_collab",
curDate=Sys.Date()
)
gene_refseq_uniprotkb_collab <- gene_refseq_uniprotkb_collab[
gene_refseq_uniprotkb_collab$method %in% c("uniprot", "identical"),
]
gene_refseq_uniprotkb_collab$NCBI_protein_accession <- sub(
"[.].*$", "", gene_refseq_uniprotkb_collab$NCBI_protein_accession
)
for(org in listOrganisms()){
message(" ", org)
curRS <- getBeIds(
be="Peptide", source="RefSeq_peptide", organism=org,
restricted=TRUE
)
toAdd <- gene_refseq_uniprotkb_collab[
which(
gene_refseq_uniprotkb_collab$NCBI_protein_accession %in% curRS$id
),
]
## External DB IDs
toImport <- unique(toAdd[, "UniProtKB_protein_accession", drop=F])
colnames(toImport) <- "id"
BED:::loadBE(
d=toImport, be="Peptide",
dbname="Uniprot",
taxId=NA
)
## The cross references
toImport <- toAdd[
, c("NCBI_protein_accession", "UniProtKB_protein_accession")
]
colnames(toImport) <- c("id1", "id2")
BED:::loadCorrespondsTo(
d=toImport,
db1="RefSeq_peptide",
db2="Uniprot",
be="Peptide"
)
}
End: r Sys.time()
Loading data from Uniprot
Release is defined according to the reldate.txt file on the Uniprot FTP
and data is downloaded only if not already done for the current release.
ftp <- "ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions"
avRel <- readLines(file.path(ftp, "reldate.txt"), n=1)
avRel <- sub(
"^UniProt Knowledgebase Release ", "",
sub(" consists of:$", "", avRel)
)
if(is.na(as.Date(paste0(avRel, "_01"), format="%Y_%m_%d"))){
print(avRel)
stop(sprintf("Check reldate.txt file on %s", ftp))
}
BED:::registerBEDB(
name="Uniprot",
description="Uniprot",
currentVersion=avRel,
idURL='http://www.uniprot.org/uniprot/%s'
)
Danio rerio data
Start: r Sys.time()
BED:::getUniprot(
organism="Danio rerio", taxDiv="vertebrates", release=avRel, ddir=wd
)
gc()
End: r Sys.time()
Homo sapiens data
Start: r Sys.time()
BED:::getUniprot(
organism="Homo sapiens", taxDiv="human", release=avRel, ddir=wd
)
gc()
End: r Sys.time()
Mus musculus data
Start: r Sys.time()
BED:::getUniprot(
organism="Mus musculus", taxDiv="rodents", release=avRel, ddir=wd
)
gc()
End: r Sys.time()
Rattus norvegicus data
Start: r Sys.time()
BED:::getUniprot(
organism="Rattus norvegicus", taxDiv="rodents", release=avRel, ddir=wd
)
gc()
End: r Sys.time()
Sus scrofa data
Start: r Sys.time()
BED:::getUniprot(
organism="Sus scrofa", taxDiv="mammals", release=avRel, ddir=wd
)
gc()
End: r Sys.time()
Indirect cross-references with EntrezGene
Start: r Sys.time()
message("Indirect cross-references with Uniprot")
dumpDir <- file.path(wd, "NCBI-gene-DATA")
f <- "gene2accession.gz"
if(file.exists(dumpDir)){
load(file.path(dumpDir, "dumpDate.rda"))
message("Last download: ", dumpDate)
if(curDate - dumpDate > reDumpThr | !file.exists(file.path(dumpDir, f))){
toDownload <- TRUE
}else{
toDownload <- FALSE
}
}else{
message("Not downloaded yet")
toDownload <- TRUE
}
if(toDownload){
ftp <- "ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/"
dlok <- try(download.file(
url=paste0(ftp, f),
destfile=file.path(dumpDir, f),
method="wget",
quiet=T
), silent=T)
}else{
message("Existing data are going to be used")
}
cn <- readLines(file.path(dumpDir, f), n=1)
cn <- sub("^#", "", cn)
cn <- unlist(strsplit(cn, split="[ \t]"))
for(org in listOrganisms()){
message(" ", org)
tid <- getTaxId(org)
toAdd <- read.table(
text=system(
sprintf("zgrep ^%s %s", tid, file.path(dumpDir, f)),
intern=TRUE
),
sep="\t",
header=F,
stringsAsFactors=F,
quote="", comment.char=""
)
colnames(toAdd) <- cn
toAdd <- toAdd[
which(toAdd$tax_id==tid),
c("tax_id", "GeneID", "protein_accession.version")
]
toAdd$pacc <- sub("[.].*$", "", toAdd$protein_accession.version)
curUP <- getBeIdConvTable(
from="Gene", from.source="BEDTech_gene", organism=org,
to="Peptide", to.source="Uniprot",
restricted=TRUE
)
toAdd <- merge(
toAdd[,c("GeneID", "pacc")],
curUP[,c("from", "to")],
by.x="pacc", by.y="to",
all=FALSE
)
toAdd <- toAdd[,c("from", "GeneID")]
toAdd$from <- as.character(toAdd$from)
toAdd$GeneID <- as.character(toAdd$GeneID)
colnames(toAdd) <- c("id1", "id2")
BED:::loadIsAssociatedTo(
d=toAdd,
db1="BEDTech_gene", db2="EntrezGene",
be="Gene"
)
}
End: r Sys.time()
Loading homologs
Orthologs from biomaRt
Start: r Sys.time()
library(biomaRt)
loadBmHomologs <- function(mart, org2){
#mattr <- listAttributes(mart)
toImport <- getBM(
mart=mart,
attributes=c(
"ensembl_gene_id",
paste0(
org2,
c("_homolog_ensembl_gene", "_homolog_orthology_confidence")
)
)
)
colnames(toImport) <- c("id1", "id2", "cs")
toImport <- unique(toImport[
which(toImport$id1 != "" & toImport$id2 != "" & toImport$cs==1),
c("id1", "id2")
])
BED:::loadIsHomologOf(
d=toImport,
db1="Ens_gene", db2="Ens_gene",
be="Gene"
)
}
#########################################
orgOfInt <- c("hsapiens", "mmusculus", "rnorvegicus", "sscrofa", "drerio")
marts <-listMarts()
bm <- "ENSEMBL_MART_ENSEMBL"
version <- ensembl_release
if(
grep(
sprintf(" %s$", version),
marts[which(marts$biomart==bm), "version"]
)==1
){
version <- NULL
}
for(i in 1:(length(orgOfInt)-1)){
#########################################
## The mart
org1 <- orgOfInt[i]
mart <- useEnsembl(
biomart=bm,
dataset=paste0(org1, "_gene_ensembl"),
version=version
)
for(j in (i+1):length(orgOfInt)){
loadBmHomologs(
mart=mart,
org2=orgOfInt[j]
)
}
}
End: r Sys.time()
Orthologs from NCBI
Start: r Sys.time()
#####################################
gdbname <- "EntrezGene"
taxOfInt <- unlist(lapply(
c(
"Homo sapiens", "Mus musculus", "Rattus norvegicus",
"Sus scrofa", "Danio rerio"
),
getTaxId
))
for(i in 1:length(taxOfInt)){
BED:::dumpNcbiDb(
taxOfInt=taxOfInt[i],
reDumpThr=reDumpThr,
ddir=wd,
toLoad=c("gene_orthologs"),
curDate=curDate
)
toImport <- gene_orthologs[
which(
gene_orthologs$tax_id %in% taxOfInt &
gene_orthologs$Other_tax_id %in% taxOfInt &
gene_orthologs$relationship == "Ortholog"
),
c("GeneID", "Other_GeneID")
]
if(nrow(toImport)>0){
colnames(toImport) <- c("id1", "id2")
toImport <- dplyr::mutate_all(toImport, as.character)
BED:::loadIsHomologOf(
d=toImport,
db1=gdbname, db2=gdbname,
be="Gene"
)
}
}
gc()
End: r Sys.time()
Loading probes
Probes from GEO
library(GEOquery)
geoPath <- file.path(wd, "geo")
dir.create(geoPath, showWarnings=FALSE)
GPL1708: Agilent-012391 Whole Human Genome Oligo Microarray G4112A (Feature Number version)
Start: r Sys.time()
## Import plateform
platName <- "GPL1708"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping with Entrez
d <- Table(gds)
toImport <- d[which(!is.na(d$SPOT_ID)), c("SPOT_ID", "GENE")]
colnames(toImport) <- c("probeID", "id")
toImport$probeID <- as.character(toImport$probeID)
toImport$id <- as.character(toImport$id)
toImport <- toImport[which(!is.na(toImport$id)),]
toImport <- unique(toImport)
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
# ## Import mapping with UniGene
# toImport <- d[which(!is.na(d$SPOT_ID)), c("SPOT_ID", "UNIGENE_ID")]
# colnames(toImport) <- c("probeID", "id")
# toImport$probeID <- as.character(toImport$probeID)
# toImport$id <- as.character(toImport$id)
# toImport <- toImport[which(!is.na(toImport$id) & toImport$id!=""),]
# dbname <- "UniGene"
# ##
# BED:::loadProbes(
# d=toImport,
# be=be,
# platform=platName,
# dbname=dbname
# )
End: r Sys.time()
GPL6480: Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version)
Start: r Sys.time()
## Import plateform
platName <- "GPL6480"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping with Entrez
d <- Table(gds)
toImport <- d[which(!is.na(d$ID)), c("ID", "GENE")]
colnames(toImport) <- c("probeID", "id")
toImport$probeID <- as.character(toImport$probeID)
toImport$id <- as.character(toImport$id)
toImport <- toImport[which(!is.na(toImport$id)),]
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
# ## Import mapping with UniGene
# toImport <- d[which(!is.na(d$ID)), c("ID", "UNIGENE_ID")]
# colnames(toImport) <- c("probeID", "id")
# toImport$probeID <- as.character(toImport$probeID)
# toImport$id <- as.character(toImport$id)
# toImport <- toImport[which(!is.na(toImport$id)),]
# dbname <- "UniGene"
# ##
# BED:::loadProbes(
# d=toImport,
# be=be,
# platform=platName,
# dbname=dbname
# )
End: r Sys.time()
GPL570: Affymetrix Human Genome U133 Plus 2.0 Array
Start: r Sys.time()
## Import plateform
platName <- "GPL570"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL571: Affymetrix Human Genome U133A 2.0 Array
Start: r Sys.time()
## Import plateform
platName <- "GPL571"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL13158: Affymetrix HT HG-U133+ PM Array Plate
Start: r Sys.time()
## Import plateform
platName <- "GPL13158"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL96: Affymetrix Human Genome U133A Array
Start: r Sys.time()
## Import plateform
platName <- "GPL96"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL1261: Affymetrix Mouse Genome 430 2.0 Array
Start: r Sys.time()
## Import plateform
platName <- "GPL1261"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL1355: Affymetrix Rat Genome 230 2.0 Array
Start: r Sys.time()
## Import plateform
platName <- "GPL1355"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- strsplit(
as.character(d$ENTREZ_GENE_ID),
split=" /// "
)
names(toImport) <- d$ID
toImport <- stack(toImport)
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL10558: Illumina HumanHT-12 V4.0 expression beadchip
Start: r Sys.time()
## Import plateform
platName <- "GPL10558"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- d[,c("Entrez_Gene_ID", "ID")]
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),]
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL6947: Illumina HumanHT-12 V3.0 expression beadchip
Start: r Sys.time()
## Import plateform
platName <- "GPL6947"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- d[,c("Entrez_Gene_ID", "ID")]
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),]
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL6885: Illumina MouseRef-8 v2.0 expression beadchip
Start: r Sys.time()
## Import plateform
platName <- "GPL6885"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
d <- Table(gds)
# e <- getBeIds(
# be="Gene", source="EntrezGene", organism="mouse", restricted=FALSE
# )
# sum(d$Entrez_Gene_ID %in% e$id) < sum(sub("[.].*$", "", d$RefSeq_ID) %in% f$id)
be <- "Transcript"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
toImport <- d[,c("RefSeq_ID", "ID")]
toImport[,1] <- as.character(toImport[,1])
toImport[,1] <- sub("[.].*$", "", toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),]
dbname <- "RefSeq"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL6887: Illumina MouseWG-6 v2.0 expression beadchip
Start: r Sys.time()
## Import plateform
platName <- "GPL6887"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
d <- Table(gds)
# e <- getBeIds(
# be="Gene", source="EntrezGene", organism="mouse", restricted=FALSE
# )
# sum(d$Entrez_Gene_ID %in% e$id) > sum(sub("[.].*$", "", d$RefSeq_ID) %in% f$id)
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
toImport <- d[,c("Entrez_Gene_ID", "ID")]
toImport[,1] <- as.character(toImport[,1])
# toImport[,1] <- sub("[.].*$", "", toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),]
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
GPL6101: Illumina ratRef-12 v1.0 expression beadchip
Start: r Sys.time()
## Import plateform
platName <- "GPL6101"
gds <- getGEO(platName, destdir=geoPath)
platDesc <- Meta(gds)$title
be <- "Gene"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping
d <- Table(gds)
toImport <- d[,c("Entrez_Gene_ID", "ID")]
toImport[,1] <- as.character(toImport[,1])
toImport[,2] <- as.character(toImport[,2])
colnames(toImport) <- c("id", "probeID")
toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),]
dbname <- "EntrezGene"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
End: r Sys.time()
Probes from biomaRt
Start: r Sys.time()
library(biomaRt)
bm <- "ENSEMBL_MART_ENSEMBL"
marts <-listMarts()
version <- ensembl_release
if(grep(
sprintf(" %s$", version),
marts[which(marts$biomart==bm), "version"]
) == 1) {
version <- NULL
}
orgOfInt <- c("hsapiens", "mmusculus", "rnorvegicus", "sscrofa", "drerio")
for(org in orgOfInt){
message(org)
mart <- useEnsembl(
biomart = bm,
dataset = paste0(org, "_gene_ensembl"),
version = version
)
at <- listAttributes(mart) %>%
dplyr::filter(
stringr::str_detect(
description, stringr::regex("probe", ignore_case=TRUE)
)
) %>%
dplyr::filter(
!name %in%
c(
"affy_huex_1_0_st_v2",
"affy_moex_1_0_st_v1",
"affy_raex_1_0_st_v1"
)
)
for(i in 1:nrow(at)){
message(" ", i, "/", nrow(at), " platforms")
message(" ", Sys.time())
platName <- at$name[i]
platDesc <- paste(at$description[i], "(Ensembl BioMart mapping)")
be <- "Transcript"
##
BED:::loadPlf(name=platName, description=platDesc, be=be)
## Import mapping with Ens_transcript
toImport <- getBM(
mart=mart,
attributes=c(
"ensembl_transcript_id",
platName
)
) %>%
dplyr::as_tibble() %>%
magrittr::set_colnames(c("id", "probeID")) %>%
dplyr::filter(!is.na(probeID) & probeID!="" & !is.na(id) & id!="") %>%
dplyr::select(probeID, id) %>%
dplyr::distinct()
dbname <- "Ens_transcript"
##
BED:::loadProbes(
d=toImport,
be=be,
platform=platName,
dbname=dbname
)
}
}
End: r Sys.time()
Loading Clarivate Analytics data from TKCat
library(TKCat)
.get_tk_headers <- do.call(function(){
## Choose the relevant credentials
credentials <- readRDS("~/etc/kmt_authorization.rds")
credentials$refresh()
return(function(){
if(!credentials$validate()){
credentials$refresh()
}
list(
"Authorization"=paste("Bearer", credentials$credentials$access_token)
)
})
}, list())
.tkcon <- chTKCat(
"tkcat.ucb.com",
password="",
port=443, https=TRUE,
extended_headers=.get_tk_headers()
)
.db_reconnect <- function(x){
xn <- deparse(substitute(x))
nv <- db_reconnect(x, extended_headers=.get_tk_headers())
assign(xn, nv, envir=parent.frame(n=1))
invisible(nv)
}
MetaBase
Start: r Sys.time()
.db_reconnect(.tkcon)
if("MetaBase" %in% list_MDBs(.tkcon)$name){
source("helpers/loadMBObjects_fromTKCat.R")
}else{
stop("MetaBase is not available in this TKCat instance.")
}
tkmb <- get_MDB(.tkcon, "MetaBase")
Register MetaBase DB
# mbInfo <- mbquery("select * from zzz_System")
BED:::registerBEDB(
name="MetaBase_gene",
description="Clarivate Analytics MetaBase",
currentVersion=tkmb$MetaBase_sourceDatabases$current,
idURL='https://portal.genego.com/cgi/entity_page.cgi?term=20&id=%s'
)
BED:::registerBEDB(
name="MetaBase_object",
description="Clarivate Analytics MetaBase",
currentVersion=tkmb$MetaBase_sourceDatabases$current,
idURL='https://portal.genego.com/cgi/entity_page.cgi?term=100&id=%s'
)
Homo sapiens data
.db_reconnect(tkmb)
loadMBObjects_fromTKCat(
orgOfInt=c("Homo sapiens"),
tkmb=tkmb
)
gc()
Mus musculus data
.db_reconnect(tkmb)
loadMBObjects_fromTKCat(
orgOfInt=c("Mus musculus"),
tkmb=tkmb
)
gc()
Rattus norvegicus data
.db_reconnect(tkmb)
loadMBObjects_fromTKCat(
orgOfInt=c("Rattus norvegicus"),
tkmb=tkmb
)
gc()
End: r Sys.time()
Cortellis
###############################################################################@
## CortellisTAR targets "Cortellis_target" coded by genes ----
# ctar <- get_MDB(.tkcon, "CortellisTAR")
# ctarv <- as.character(ctar$CortellisTAR_APIdump$date)
# targets <- ctar$CortellisTAR_targets %>%
# filter(organism %in% listOrganisms())
# prtargets <- targets %>% filter(type!="Nonspecified") %>%
# select(id, nameMain)
# xref <- ctar$CortellisTAR_externalIdentifiers
# entXref <- xref %>%
# filter(be=="Gene", source=="EntrezGene", targetId %in% prtargets$id)
# stopifnot(all(prtargets$id %in% entXref$targetId))
# allent <- getBeIds(be="Gene", source="EntrezGene", restricted=FALSE)
# entXref <- entXref %>% filter(externalId %in% allent$id)
#
# ### Add Cortellis_target database ----
# ddName <- "Cortellis_target"
# BED:::registerBEDB(
# name=ddName,
# description="Clarivate Analytics Cortellis Targets",
# currentVersion=ctarv,
# idURL=NA
# )
#
# ### Add CortellisTAR targets ----
# toImport <- prtargets[, "id", drop=F]
# colnames(toImport) <- "id"
# BED:::loadBE(
# d=toImport, be="Object",
# dbname=ddName,
# version=NA,
# taxId=NA
# )
#
# ### Add CortellisTAR target names ----
# toImport <- prtargets
# colnames(toImport) <- c("id", "name")
# BED:::loadBENames(d=toImport, be="Object", dbname=ddName)
#
#
# ### Add "codes_for" edges ----
# toImport <- entXref %>%
# select(externalId, targetId)
# colnames(toImport) <- c("gid", "oid")
# BED:::loadCodesFor(
# d=toImport,
# gdb="EntrezGene",
# odb=ddName
# )
###############################################################################@
## CortellisId targets "Cortellis_idtarget" coded by genes ----
cid <- get_MDB(.tkcon, "CortellisID")
cidv <- as.character(cid$CortellisID_APIdump$date)
alltargets <- cid$CortellisID_targets
up <- cid$CortellisID_targetUniprot
### Get uniprot entrez gene identifiers ----
upEntrez <- c()
for(org in listOrganisms()){
message(org)
toAdd <- convBeIds(
unique(up$uniprot),
from="Peptide", from.source="Uniprot", from.org=org,
to="Gene", to.source="EntrezGene",
prefFilter=TRUE, restricted=FALSE
)
upEntrez <- bind_rows(
upEntrez,
toAdd %>%
as_tibble() %>%
filter(!is.na(to)) %>%
select("uniprot"=from, "entrez"=to) %>%
mutate(organism=!!org)
)
}
fup <- left_join(upEntrez, up, by="uniprot", relationship = "many-to-many")
### Add Cortellis_idtarget database ----
ddName <- "Cortellis_idtarget"
BED:::registerBEDB(
name=ddName,
description="Clarivate Analytics Cortellis Investigational Drugs",
currentVersion=cidv,
idURL=NA
)
### Add CortellisID targets ----
targets <- alltargets # %>%
# filter(id %in% fup$target)
toImport <- targets[, "id", drop=F]
colnames(toImport) <- "id"
BED:::loadBE(
d=toImport, be="Object",
dbname=ddName,
version=NA,
taxId=NA
)
### Add CortellisID target names ----
toImport <- targets
colnames(toImport) <- c("id", "name")
BED:::loadBENames(d=toImport, be="Object", dbname=ddName)
### Add "codes_for" edges ----
toImport <- fup %>%
select(entrez, target)
colnames(toImport) <- c("gid", "oid")
BED:::loadCodesFor(
d=toImport,
gdb="EntrezGene",
odb=ddName
)
# ## CortellisTAR ID targets "Cortellis_idtarget" related to "Cortellis_target"
# idtargets <- ctar$CortellisTAR_IDTargets
# idtargets2 <- cid$CortellisID_targetXref %>%
# filter(source=="Cortellis_target") %>%
# select(cortellisTAR=xref, cortellisID=targetId)
# idtargets <- bind_rows(idtargets, idtargets2) %>%
# distinct(cortellisTAR, cortellisID)
# toImport <- idtargets %>%
# filter(cortellisTAR %in% prtargets$id, cortellisID %in% targets$id) %>%
# select(id1=cortellisID, id2=cortellisTAR)
# BED:::loadIsAssociatedTo(
# d=toImport,
# db1="Cortellis_idtarget",
# db2="Cortellis_target",
# be="Object"
# )
## CortellisID targets "Cortellis_idtarget" related to "MetaBase_object"
idObjects <- cid$CortellisID_targetXref %>%
filter(source=="MetaBase_object") %>%
select(metabase=xref, cortellisID=targetId)
mbo <- getBeIds(be="Object", source="MetaBase_object", restricted=FALSE)
toImport <- idObjects %>%
filter(metabase %in% mbo$id, cortellisID %in% targets$id) %>%
select(id1=cortellisID, id2=metabase)
BED:::loadIsAssociatedTo(
d=toImport,
db1="Cortellis_idtarget",
db2="MetaBase_object",
be="Object"
)
Other information
Databases ID URL
otherIdURL <- list(
"HGNC"='http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=%s',
"miRBase"='http://www.mirbase.org/cgi-bin/mirna_entry.pl?acc=%s',
"Vega_gene"='http://vega.sanger.ac.uk/id/%s',
"UniGene"='https://www.ncbi.nlm.nih.gov/unigene?term=%s',
"Vega_transcript"='http://vega.sanger.ac.uk/id/%s',
"MGI"='http://www.informatics.jax.org/marker/MGI:%s',
"Vega_translation"='http://vega.sanger.ac.uk/id/%s',
"RGD"='https://rgd.mcw.edu/rgdweb/report/gene/main.html?id=%s',
"MIM_GENE"='http://www.omim.org/entry/%s',
"ZFIN_gene"='http://zfin.org/%s'
)
for(db in names(otherIdURL)){
BED:::registerBEDB(
name=db,
idURL=otherIdURL[[db]]
)
}
Load Lucene Indexes
Start: r Sys.time()
BED:::loadLuceneIndexes()
End: r Sys.time()
Session info
sessionInfo()
patzaw/BED documentation built on Dec. 13, 2024, 4:29 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.
Introduction
This document shows how to feed the Biological Entity Dictionary (BED).
It can be adapted according to specific needs and DB access.
The BED functions used to feed the DB are not exported to avoid
unintended modifications of the DB. To call them, they are preceded
by BED:::.
In this example several source databases are dumped and their content
is integrated in BED.
Some helper functions are provided to get information from famous databases.
The following chunk is used to configure source versions.
The reDumpThr object is used to define time intervals during which some
data sources should not be re-downloaded.
## library(knitr) library(BED) library(jsonlite) library(rvest) ## config <- jsonlite::read_json("build_config.json") config <- lapply( config, function(x){ if(!is.character(x)){ return(x) }else{ sub(pattern="___HOME___", replacement=Sys.getenv("HOME"), x = x) } } ) config <- lapply( config, function(x){ if(!is.character(x)){ return(x) }else{ sub(pattern="___ROOT___", replacement=config$ROOT, x = x) } } ) ## wd <- config$BED_WORKING ## # opts_knit$set(root.dir=wd) opts_chunk$set( eval=TRUE, message=FALSE # root.dir=wd ) ## Specific config bedInstance <- config$BED_INSTANCE bedVersion <- config$BED_VERSION ## Identify the latest ensembl release here: https://www.ensembl.org/ ensembl_org_dir <- grep( "_core_", rvest::html_table( rvest::read_html("https://ftp.ensembl.org/pub/current/mysql/") )[[1]]$Name, value=TRUE ) ensembl_org_table <- do.call(rbind, lapply( strsplit(ensembl_org_dir, split="_"), function(x){ gv <- sub("/", "", x[length(x)]) e <- x[length(x)-1] core <- which(x=="core") org <- paste(x[1:(core-1)], collapse=" ") substr(org, 1, 1) <- toupper(substr(org, 1,1)) return(data.frame( organism=org, release=e, genome_version=gv )) } )) ensembl_org_table <- ensembl_org_table[ order(as.numeric(ensembl_org_table$genome_version), decreasing=T), ] ensembl_org_table <- ensembl_org_table[ which(!duplicated(ensembl_org_table$organism)), ] ensembl_release <- unique(ensembl_org_table$release) stopifnot(length(ensembl_release)==1) org <-"Homo sapiens" ensembl_Hsapiens <- list( release=ensembl_release, organism=org, gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)], gdbCref=c( # Gene cross-references DBs "HGNC"="HGNC", "EntrezGene"="EntrezGene", "Vega_gene"="Vega_gene", "Ens_Hs_gene"="Ens_gene" ), gdbAss=c( # Gene associated IDs (DB) "miRBase"="miRBase", "MIM_GENE"="MIM_GENE", "UniGene"="UniGene" ), tdbCref=c( # Transcript cross-references DBs "RefSeq_mRNA"="RefSeq", "RefSeq_ncRNA"="RefSeq", "RefSeq_mRNA_predicted"="RefSeq", "RefSeq_ncRNA_predicted"="RefSeq", "Vega_transcript"="Vega_transcript", "Ens_Hs_transcript"="Ens_transcript" ), pdbCref=c( # Peptide cross-references DBs "RefSeq_peptide"="RefSeq_peptide", "RefSeq_peptide_predicted"="RefSeq_peptide", "Uniprot/SPTREMBL"="Uniprot", "Uniprot/SWISSPROT"="Uniprot", "Vega_translation"="Vega_translation", "Ens_Hs_translation"="Ens_translation" ), canChromosomes=c(1:22, "X", "Y", "MT") ) org <- "Mus musculus" ensembl_Mmusculus <- list( release=ensembl_release, organism=org, gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)], gdbCref=c( # Gene cross-references DBs "MGI"="MGI", "EntrezGene"="EntrezGene", "Vega_gene"="Vega_gene", "Ens_Mm_gene"="Ens_gene" ), gdbAss=c( # Gene associated IDs (DB) "miRBase"="miRBase", "UniGene"="UniGene" ), tdbCref=c( # Transcript cross-references DBs "RefSeq_mRNA"="RefSeq", "RefSeq_ncRNA"="RefSeq", "RefSeq_mRNA_predicted"="RefSeq", "RefSeq_ncRNA_predicted"="RefSeq", "Vega_transcript"="Vega_transcript", "Ens_Mm_transcript"="Ens_transcript" ), pdbCref=c( # Peptide cross-references DBs "RefSeq_peptide"="RefSeq_peptide", "RefSeq_peptide_predicted"="RefSeq_peptide", "Uniprot/SPTREMBL"="Uniprot", "Uniprot/SWISSPROT"="Uniprot", "Vega_translation"="Vega_translation", "Ens_Mm_translation"="Ens_translation" ), canChromosomes=c(1:19, "X", "Y", "MT") ) org <- "Rattus norvegicus" ensembl_Rnorvegicus <- list( release=ensembl_release, organism=org, gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)], gdbCref=c( # Gene cross-references DBs "RGD"="RGD", "EntrezGene"="EntrezGene", "Vega_gene"="Vega_gene", "Ens_Rn_gene"="Ens_gene" ), gdbAss=c( # Gene associated IDs (DB) "miRBase"="miRBase", "UniGene"="UniGene" ), tdbCref=c( # Transcript cross-references DBs "RefSeq_mRNA"="RefSeq", "RefSeq_ncRNA"="RefSeq", "RefSeq_mRNA_predicted"="RefSeq", "RefSeq_ncRNA_predicted"="RefSeq", "Vega_transcript"="Vega_transcript", "Ens_Rn_transcript"="Ens_transcript" ), pdbCref=c( # Peptide cross-references DBs "RefSeq_peptide"="RefSeq_peptide", "RefSeq_peptide_predicted"="RefSeq_peptide", "Uniprot/SPTREMBL"="Uniprot", "Uniprot/SWISSPROT"="Uniprot", "Vega_translation"="Vega_translation", "Ens_Rn_translation"="Ens_translation" ), canChromosomes=c(1:20, "X", "Y", "MT") ) org <- "Sus scrofa" ensembl_Sscrofa <- list( release=ensembl_release, organism=org, gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)], gdbCref=c( # Gene cross-references DBs "EntrezGene"="EntrezGene", "Vega_gene"="Vega_gene", "Ens_Ss_gene"="Ens_gene" ), gdbAss=c( # Gene associated IDs (DB) "miRBase"="miRBase", "UniGene"="UniGene" ), tdbCref=c( # Transcript cross-references DBs "RefSeq_mRNA"="RefSeq", "RefSeq_ncRNA"="RefSeq", "RefSeq_mRNA_predicted"="RefSeq", "RefSeq_ncRNA_predicted"="RefSeq", "Vega_transcript"="Vega_transcript", "Ens_Ss_transcript"="Ens_transcript" ), pdbCref=c( # Peptide cross-references DBs "RefSeq_peptide"="RefSeq_peptide", "RefSeq_peptide_predicted"="RefSeq_peptide", "Uniprot/SPTREMBL"="Uniprot", "Uniprot/SWISSPROT"="Uniprot", "Vega_translation"="Vega_translation", "Ens_Ss_translation"="Ens_translation" ), canChromosomes=c(1:18, "X", "Y", "MT") ) org <- "Danio rerio" ensembl_Drerio <- list( release=ensembl_release, organism=org, gv=ensembl_org_table$genome_version[which(ensembl_org_table$organism==org)], gdbCref=c( # Gene cross-references DBs "EntrezGene"="EntrezGene", "ZFIN_ID"="ZFIN_gene", "Vega_gene"="Vega_gene", "Ens_Dr_gene"="Ens_gene" ), gdbAss=c( # Gene associated IDs (DB) "miRBase"="miRBase", "UniGene"="UniGene" ), tdbCref=c( # Transcript cross-references DBs "RefSeq_mRNA"="RefSeq", "RefSeq_ncRNA"="RefSeq", "RefSeq_mRNA_predicted"="RefSeq", "RefSeq_ncRNA_predicted"="RefSeq", "Vega_transcript"="Vega_transcript", "Ens_Dr_transcript"="Ens_transcript" ), pdbCref=c( # Peptide cross-references DBs "RefSeq_peptide"="RefSeq_peptide", "RefSeq_peptide_predicted"="RefSeq_peptide", "Uniprot/SPTREMBL"="Uniprot", "Uniprot/SWISSPROT"="Uniprot", "Vega_translation"="Vega_translation", "Ens_Dr_translation"="Ens_translation" ), canChromosomes=c(1:25, "MT") ) ## General config reDumpThr <- as.difftime(config$SOURCE_REDUMP_THR, units="days") curDate <- Sys.Date()
BED initialization
BED is based on Neo4j.
The S01-NewBED-Container.sh shows how to run it in a docker container.
Because the import functions use massively the LOAD CSV Neo4j query,
the feeding of the BED database can only be down from the
computer hosting the Neo4j relevant instance.
The chunk below shows how to connect to BED. In this example, neo4j authentication is disabled.
connectToBed( url=sprintf("localhost:%s", config$NJ_HTTP_PORT), remember=FALSE, useCache=TRUE, importPath=config$BED_IMPORT )
Check empty DB
Do not go further if your BED DB is not empty.
dbSize <- bedCall(cypher, 'MATCH (n) RETURN count(n)')[,1] if(dbSize!=0){ stop("BED DB is not empty ==> clean it before loading the content below") }
Set BED instance and version
print(bedInstance) print(bedVersion) BED:::setBedVersion(bedInstance=bedInstance, bedVersion=bedVersion)
Load Data model
Start: r Sys.time()
BED:::loadBedModel()
End: r Sys.time()
Loading taxonomy from NCBI
Information is downloaded if older than r reDumpThr days
according to the reDumpThr object.
Start: r Sys.time()
BED:::loadNcbiTax( reDumpThr=reDumpThr, ddir=wd, orgOfInt=c( "Homo sapiens", "Rattus norvegicus", "Mus musculus", "Sus scrofa", "Danio rerio" ), curDate=curDate )
End: r Sys.time()
Loading data from Ensembl
Register Ensembl DBs
Genes
BED:::registerBEDB( name="Ens_gene", description="Ensembl gene", currentVersion=ensembl_release, idURL='http://www.ensembl.org/id/%s' )
Transcripts
BED:::registerBEDB( name="Ens_transcript", description="Ensembl transcript", currentVersion=ensembl_release, idURL='http://www.ensembl.org/id/%s' )
Peptides
BED:::registerBEDB( name="Ens_translation", description="Ensembl peptides", currentVersion=ensembl_release, idURL='http://www.ensembl.org/id/%s' )
Danio rerio
ensembl <- ensembl_Drerio print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$gdbCref, dbAss=ensembl$gdbAss, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$tdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$pdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Homo sapiens
ensembl <- ensembl_Hsapiens print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$gdbCref, dbAss=ensembl$gdbAss, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$tdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$pdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Mus musculus
ensembl <- ensembl_Mmusculus print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$gdbCref, dbAss=ensembl$gdbAss, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$tdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$pdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Rattus norvegicus
ensembl <- ensembl_Rnorvegicus print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$gdbCref, dbAss=ensembl$gdbAss, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$tdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$pdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Sus scrofa
ensembl <- ensembl_Sscrofa print(ensembl)
Genes
Start: r Sys.time()
BED:::getEnsemblGeneIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$gdbCref, dbAss=ensembl$gdbAss, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Transcripts
Start: r Sys.time()
BED:::getEnsemblTranscriptIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$tdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Peptides
Start: r Sys.time()
BED:::getEnsemblPeptideIds( organism=ensembl$organism, release=ensembl$release, gv=ensembl$gv, ddir=wd, dbCref=ensembl$pdbCref, canChromosomes=ensembl$canChromosomes ) gc()
End: r Sys.time()
Loading data from NCBI
Information is downloaded if older than r reDumpThr days
according to the reDumpThr object.
Register NCBI DBs
BED:::dumpNcbiDb( taxOfInt = c(), reDumpThr=reDumpThr, ddir=wd, toLoad=c(), curDate=curDate )
Genes
BED:::registerBEDB( name="EntrezGene", description="NCBI gene", currentVersion=format(dumpDate, "%Y%m%d"), idURL='https://www.ncbi.nlm.nih.gov/gene/%s' )
Transcripts
BED:::registerBEDB( name="RefSeq", description="NCBI nucleotide", currentVersion=format(dumpDate, "%Y%m%d"), idURL='https://www.ncbi.nlm.nih.gov/nuccore/%s' )
Peptides
BED:::registerBEDB( name="RefSeq_peptide", description="NCBI protein", currentVersion=format(dumpDate, "%Y%m%d"), idURL='https://www.ncbi.nlm.nih.gov/protein/%s' )
Danio rerio data
Start: r Sys.time()
BED:::getNcbiGeneTransPep( organism="Danio rerio", ddir=wd, curDate=curDate ) gc()
End: r Sys.time()
Homo sapiens data
Start: r Sys.time()
BED:::getNcbiGeneTransPep( organism="Homo sapiens", ddir=wd, curDate=curDate ) gc()
End: r Sys.time()
Mus musculus data
Start: r Sys.time()
BED:::getNcbiGeneTransPep( organism="Mus musculus", ddir=wd, curDate=curDate ) gc()
End: r Sys.time()
Rattus norvegicus data
Start: r Sys.time()
BED:::getNcbiGeneTransPep( organism="Rattus norvegicus", ddir=wd, curDate=curDate ) gc()
End: r Sys.time()
Sus scrofa data
Start: r Sys.time()
BED:::getNcbiGeneTransPep( organism="Sus scrofa", ddir=wd, curDate=curDate ) gc()
End: r Sys.time()
Direct cross-references with Uniprot
Start: r Sys.time()
message("Direct cross-references with Uniprot") BED:::dumpNcbiDb( taxOfInt="", reDumpThr=Inf, ddir=wd, toLoad="gene_refseq_uniprotkb_collab", curDate=Sys.Date() ) gene_refseq_uniprotkb_collab <- gene_refseq_uniprotkb_collab[ gene_refseq_uniprotkb_collab$method %in% c("uniprot", "identical"), ] gene_refseq_uniprotkb_collab$NCBI_protein_accession <- sub( "[.].*$", "", gene_refseq_uniprotkb_collab$NCBI_protein_accession ) for(org in listOrganisms()){ message(" ", org) curRS <- getBeIds( be="Peptide", source="RefSeq_peptide", organism=org, restricted=TRUE ) toAdd <- gene_refseq_uniprotkb_collab[ which( gene_refseq_uniprotkb_collab$NCBI_protein_accession %in% curRS$id ), ] ## External DB IDs toImport <- unique(toAdd[, "UniProtKB_protein_accession", drop=F]) colnames(toImport) <- "id" BED:::loadBE( d=toImport, be="Peptide", dbname="Uniprot", taxId=NA ) ## The cross references toImport <- toAdd[ , c("NCBI_protein_accession", "UniProtKB_protein_accession") ] colnames(toImport) <- c("id1", "id2") BED:::loadCorrespondsTo( d=toImport, db1="RefSeq_peptide", db2="Uniprot", be="Peptide" ) }
End: r Sys.time()
Loading data from Uniprot
Release is defined according to the reldate.txt file on the Uniprot FTP and data is downloaded only if not already done for the current release.
ftp <- "ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions" avRel <- readLines(file.path(ftp, "reldate.txt"), n=1) avRel <- sub( "^UniProt Knowledgebase Release ", "", sub(" consists of:$", "", avRel) ) if(is.na(as.Date(paste0(avRel, "_01"), format="%Y_%m_%d"))){ print(avRel) stop(sprintf("Check reldate.txt file on %s", ftp)) } BED:::registerBEDB( name="Uniprot", description="Uniprot", currentVersion=avRel, idURL='http://www.uniprot.org/uniprot/%s' )
Danio rerio data
Start: r Sys.time()
BED:::getUniprot( organism="Danio rerio", taxDiv="vertebrates", release=avRel, ddir=wd ) gc()
End: r Sys.time()
Homo sapiens data
Start: r Sys.time()
BED:::getUniprot( organism="Homo sapiens", taxDiv="human", release=avRel, ddir=wd ) gc()
End: r Sys.time()
Mus musculus data
Start: r Sys.time()
BED:::getUniprot( organism="Mus musculus", taxDiv="rodents", release=avRel, ddir=wd ) gc()
End: r Sys.time()
Rattus norvegicus data
Start: r Sys.time()
BED:::getUniprot( organism="Rattus norvegicus", taxDiv="rodents", release=avRel, ddir=wd ) gc()
End: r Sys.time()
Sus scrofa data
Start: r Sys.time()
BED:::getUniprot( organism="Sus scrofa", taxDiv="mammals", release=avRel, ddir=wd ) gc()
End: r Sys.time()
Indirect cross-references with EntrezGene
Start: r Sys.time()
message("Indirect cross-references with Uniprot") dumpDir <- file.path(wd, "NCBI-gene-DATA") f <- "gene2accession.gz" if(file.exists(dumpDir)){ load(file.path(dumpDir, "dumpDate.rda")) message("Last download: ", dumpDate) if(curDate - dumpDate > reDumpThr | !file.exists(file.path(dumpDir, f))){ toDownload <- TRUE }else{ toDownload <- FALSE } }else{ message("Not downloaded yet") toDownload <- TRUE } if(toDownload){ ftp <- "ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/" dlok <- try(download.file( url=paste0(ftp, f), destfile=file.path(dumpDir, f), method="wget", quiet=T ), silent=T) }else{ message("Existing data are going to be used") } cn <- readLines(file.path(dumpDir, f), n=1) cn <- sub("^#", "", cn) cn <- unlist(strsplit(cn, split="[ \t]")) for(org in listOrganisms()){ message(" ", org) tid <- getTaxId(org) toAdd <- read.table( text=system( sprintf("zgrep ^%s %s", tid, file.path(dumpDir, f)), intern=TRUE ), sep="\t", header=F, stringsAsFactors=F, quote="", comment.char="" ) colnames(toAdd) <- cn toAdd <- toAdd[ which(toAdd$tax_id==tid), c("tax_id", "GeneID", "protein_accession.version") ] toAdd$pacc <- sub("[.].*$", "", toAdd$protein_accession.version) curUP <- getBeIdConvTable( from="Gene", from.source="BEDTech_gene", organism=org, to="Peptide", to.source="Uniprot", restricted=TRUE ) toAdd <- merge( toAdd[,c("GeneID", "pacc")], curUP[,c("from", "to")], by.x="pacc", by.y="to", all=FALSE ) toAdd <- toAdd[,c("from", "GeneID")] toAdd$from <- as.character(toAdd$from) toAdd$GeneID <- as.character(toAdd$GeneID) colnames(toAdd) <- c("id1", "id2") BED:::loadIsAssociatedTo( d=toAdd, db1="BEDTech_gene", db2="EntrezGene", be="Gene" ) }
End: r Sys.time()
Loading homologs
Orthologs from biomaRt
Start: r Sys.time()
library(biomaRt) loadBmHomologs <- function(mart, org2){ #mattr <- listAttributes(mart) toImport <- getBM( mart=mart, attributes=c( "ensembl_gene_id", paste0( org2, c("_homolog_ensembl_gene", "_homolog_orthology_confidence") ) ) ) colnames(toImport) <- c("id1", "id2", "cs") toImport <- unique(toImport[ which(toImport$id1 != "" & toImport$id2 != "" & toImport$cs==1), c("id1", "id2") ]) BED:::loadIsHomologOf( d=toImport, db1="Ens_gene", db2="Ens_gene", be="Gene" ) } ######################################### orgOfInt <- c("hsapiens", "mmusculus", "rnorvegicus", "sscrofa", "drerio") marts <-listMarts() bm <- "ENSEMBL_MART_ENSEMBL" version <- ensembl_release if( grep( sprintf(" %s$", version), marts[which(marts$biomart==bm), "version"] )==1 ){ version <- NULL } for(i in 1:(length(orgOfInt)-1)){ ######################################### ## The mart org1 <- orgOfInt[i] mart <- useEnsembl( biomart=bm, dataset=paste0(org1, "_gene_ensembl"), version=version ) for(j in (i+1):length(orgOfInt)){ loadBmHomologs( mart=mart, org2=orgOfInt[j] ) } }
End: r Sys.time()
Orthologs from NCBI
Start: r Sys.time()
##################################### gdbname <- "EntrezGene" taxOfInt <- unlist(lapply( c( "Homo sapiens", "Mus musculus", "Rattus norvegicus", "Sus scrofa", "Danio rerio" ), getTaxId )) for(i in 1:length(taxOfInt)){ BED:::dumpNcbiDb( taxOfInt=taxOfInt[i], reDumpThr=reDumpThr, ddir=wd, toLoad=c("gene_orthologs"), curDate=curDate ) toImport <- gene_orthologs[ which( gene_orthologs$tax_id %in% taxOfInt & gene_orthologs$Other_tax_id %in% taxOfInt & gene_orthologs$relationship == "Ortholog" ), c("GeneID", "Other_GeneID") ] if(nrow(toImport)>0){ colnames(toImport) <- c("id1", "id2") toImport <- dplyr::mutate_all(toImport, as.character) BED:::loadIsHomologOf( d=toImport, db1=gdbname, db2=gdbname, be="Gene" ) } } gc()
End: r Sys.time()
Loading probes
Probes from GEO
library(GEOquery) geoPath <- file.path(wd, "geo") dir.create(geoPath, showWarnings=FALSE)
GPL1708: Agilent-012391 Whole Human Genome Oligo Microarray G4112A (Feature Number version)
Start: r Sys.time()
## Import plateform platName <- "GPL1708" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping with Entrez d <- Table(gds) toImport <- d[which(!is.na(d$SPOT_ID)), c("SPOT_ID", "GENE")] colnames(toImport) <- c("probeID", "id") toImport$probeID <- as.character(toImport$probeID) toImport$id <- as.character(toImport$id) toImport <- toImport[which(!is.na(toImport$id)),] toImport <- unique(toImport) dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname ) # ## Import mapping with UniGene # toImport <- d[which(!is.na(d$SPOT_ID)), c("SPOT_ID", "UNIGENE_ID")] # colnames(toImport) <- c("probeID", "id") # toImport$probeID <- as.character(toImport$probeID) # toImport$id <- as.character(toImport$id) # toImport <- toImport[which(!is.na(toImport$id) & toImport$id!=""),] # dbname <- "UniGene" # ## # BED:::loadProbes( # d=toImport, # be=be, # platform=platName, # dbname=dbname # )
End: r Sys.time()
GPL6480: Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version)
Start: r Sys.time()
## Import plateform platName <- "GPL6480" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping with Entrez d <- Table(gds) toImport <- d[which(!is.na(d$ID)), c("ID", "GENE")] colnames(toImport) <- c("probeID", "id") toImport$probeID <- as.character(toImport$probeID) toImport$id <- as.character(toImport$id) toImport <- toImport[which(!is.na(toImport$id)),] dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname ) # ## Import mapping with UniGene # toImport <- d[which(!is.na(d$ID)), c("ID", "UNIGENE_ID")] # colnames(toImport) <- c("probeID", "id") # toImport$probeID <- as.character(toImport$probeID) # toImport$id <- as.character(toImport$id) # toImport <- toImport[which(!is.na(toImport$id)),] # dbname <- "UniGene" # ## # BED:::loadProbes( # d=toImport, # be=be, # platform=platName, # dbname=dbname # )
End: r Sys.time()
GPL570: Affymetrix Human Genome U133 Plus 2.0 Array
Start: r Sys.time()
## Import plateform platName <- "GPL570" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL571: Affymetrix Human Genome U133A 2.0 Array
Start: r Sys.time()
## Import plateform platName <- "GPL571" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL13158: Affymetrix HT HG-U133+ PM Array Plate
Start: r Sys.time()
## Import plateform platName <- "GPL13158" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL96: Affymetrix Human Genome U133A Array
Start: r Sys.time()
## Import plateform platName <- "GPL96" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL1261: Affymetrix Mouse Genome 430 2.0 Array
Start: r Sys.time()
## Import plateform platName <- "GPL1261" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL1355: Affymetrix Rat Genome 230 2.0 Array
Start: r Sys.time()
## Import plateform platName <- "GPL1355" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- strsplit( as.character(d$ENTREZ_GENE_ID), split=" /// " ) names(toImport) <- d$ID toImport <- stack(toImport) toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL10558: Illumina HumanHT-12 V4.0 expression beadchip
Start: r Sys.time()
## Import plateform platName <- "GPL10558" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- d[,c("Entrez_Gene_ID", "ID")] toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),] dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL6947: Illumina HumanHT-12 V3.0 expression beadchip
Start: r Sys.time()
## Import plateform platName <- "GPL6947" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- d[,c("Entrez_Gene_ID", "ID")] toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),] dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL6885: Illumina MouseRef-8 v2.0 expression beadchip
Start: r Sys.time()
## Import plateform platName <- "GPL6885" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title d <- Table(gds) # e <- getBeIds( # be="Gene", source="EntrezGene", organism="mouse", restricted=FALSE # ) # sum(d$Entrez_Gene_ID %in% e$id) < sum(sub("[.].*$", "", d$RefSeq_ID) %in% f$id) be <- "Transcript" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping toImport <- d[,c("RefSeq_ID", "ID")] toImport[,1] <- as.character(toImport[,1]) toImport[,1] <- sub("[.].*$", "", toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),] dbname <- "RefSeq" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL6887: Illumina MouseWG-6 v2.0 expression beadchip
Start: r Sys.time()
## Import plateform platName <- "GPL6887" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title d <- Table(gds) # e <- getBeIds( # be="Gene", source="EntrezGene", organism="mouse", restricted=FALSE # ) # sum(d$Entrez_Gene_ID %in% e$id) > sum(sub("[.].*$", "", d$RefSeq_ID) %in% f$id) be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping toImport <- d[,c("Entrez_Gene_ID", "ID")] toImport[,1] <- as.character(toImport[,1]) # toImport[,1] <- sub("[.].*$", "", toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),] dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
GPL6101: Illumina ratRef-12 v1.0 expression beadchip
Start: r Sys.time()
## Import plateform platName <- "GPL6101" gds <- getGEO(platName, destdir=geoPath) platDesc <- Meta(gds)$title be <- "Gene" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping d <- Table(gds) toImport <- d[,c("Entrez_Gene_ID", "ID")] toImport[,1] <- as.character(toImport[,1]) toImport[,2] <- as.character(toImport[,2]) colnames(toImport) <- c("id", "probeID") toImport <- toImport[which(!is.na(toImport$id) & toImport$id != ""),] dbname <- "EntrezGene" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname )
End: r Sys.time()
Probes from biomaRt
Start: r Sys.time()
library(biomaRt) bm <- "ENSEMBL_MART_ENSEMBL" marts <-listMarts() version <- ensembl_release if(grep( sprintf(" %s$", version), marts[which(marts$biomart==bm), "version"] ) == 1) { version <- NULL } orgOfInt <- c("hsapiens", "mmusculus", "rnorvegicus", "sscrofa", "drerio") for(org in orgOfInt){ message(org) mart <- useEnsembl( biomart = bm, dataset = paste0(org, "_gene_ensembl"), version = version ) at <- listAttributes(mart) %>% dplyr::filter( stringr::str_detect( description, stringr::regex("probe", ignore_case=TRUE) ) ) %>% dplyr::filter( !name %in% c( "affy_huex_1_0_st_v2", "affy_moex_1_0_st_v1", "affy_raex_1_0_st_v1" ) ) for(i in 1:nrow(at)){ message(" ", i, "/", nrow(at), " platforms") message(" ", Sys.time()) platName <- at$name[i] platDesc <- paste(at$description[i], "(Ensembl BioMart mapping)") be <- "Transcript" ## BED:::loadPlf(name=platName, description=platDesc, be=be) ## Import mapping with Ens_transcript toImport <- getBM( mart=mart, attributes=c( "ensembl_transcript_id", platName ) ) %>% dplyr::as_tibble() %>% magrittr::set_colnames(c("id", "probeID")) %>% dplyr::filter(!is.na(probeID) & probeID!="" & !is.na(id) & id!="") %>% dplyr::select(probeID, id) %>% dplyr::distinct() dbname <- "Ens_transcript" ## BED:::loadProbes( d=toImport, be=be, platform=platName, dbname=dbname ) } }
End: r Sys.time()
Loading Clarivate Analytics data from TKCat
library(TKCat) .get_tk_headers <- do.call(function(){ ## Choose the relevant credentials credentials <- readRDS("~/etc/kmt_authorization.rds") credentials$refresh() return(function(){ if(!credentials$validate()){ credentials$refresh() } list( "Authorization"=paste("Bearer", credentials$credentials$access_token) ) }) }, list()) .tkcon <- chTKCat( "tkcat.ucb.com", password="", port=443, https=TRUE, extended_headers=.get_tk_headers() ) .db_reconnect <- function(x){ xn <- deparse(substitute(x)) nv <- db_reconnect(x, extended_headers=.get_tk_headers()) assign(xn, nv, envir=parent.frame(n=1)) invisible(nv) }
MetaBase
Start: r Sys.time()
.db_reconnect(.tkcon) if("MetaBase" %in% list_MDBs(.tkcon)$name){ source("helpers/loadMBObjects_fromTKCat.R") }else{ stop("MetaBase is not available in this TKCat instance.") } tkmb <- get_MDB(.tkcon, "MetaBase")
Register MetaBase DB
# mbInfo <- mbquery("select * from zzz_System") BED:::registerBEDB( name="MetaBase_gene", description="Clarivate Analytics MetaBase", currentVersion=tkmb$MetaBase_sourceDatabases$current, idURL='https://portal.genego.com/cgi/entity_page.cgi?term=20&id=%s' ) BED:::registerBEDB( name="MetaBase_object", description="Clarivate Analytics MetaBase", currentVersion=tkmb$MetaBase_sourceDatabases$current, idURL='https://portal.genego.com/cgi/entity_page.cgi?term=100&id=%s' )
Homo sapiens data
.db_reconnect(tkmb) loadMBObjects_fromTKCat( orgOfInt=c("Homo sapiens"), tkmb=tkmb ) gc()
Mus musculus data
.db_reconnect(tkmb) loadMBObjects_fromTKCat( orgOfInt=c("Mus musculus"), tkmb=tkmb ) gc()
Rattus norvegicus data
.db_reconnect(tkmb) loadMBObjects_fromTKCat( orgOfInt=c("Rattus norvegicus"), tkmb=tkmb ) gc()
End: r Sys.time()
Cortellis
###############################################################################@ ## CortellisTAR targets "Cortellis_target" coded by genes ---- # ctar <- get_MDB(.tkcon, "CortellisTAR") # ctarv <- as.character(ctar$CortellisTAR_APIdump$date) # targets <- ctar$CortellisTAR_targets %>% # filter(organism %in% listOrganisms()) # prtargets <- targets %>% filter(type!="Nonspecified") %>% # select(id, nameMain) # xref <- ctar$CortellisTAR_externalIdentifiers # entXref <- xref %>% # filter(be=="Gene", source=="EntrezGene", targetId %in% prtargets$id) # stopifnot(all(prtargets$id %in% entXref$targetId)) # allent <- getBeIds(be="Gene", source="EntrezGene", restricted=FALSE) # entXref <- entXref %>% filter(externalId %in% allent$id) # # ### Add Cortellis_target database ---- # ddName <- "Cortellis_target" # BED:::registerBEDB( # name=ddName, # description="Clarivate Analytics Cortellis Targets", # currentVersion=ctarv, # idURL=NA # ) # # ### Add CortellisTAR targets ---- # toImport <- prtargets[, "id", drop=F] # colnames(toImport) <- "id" # BED:::loadBE( # d=toImport, be="Object", # dbname=ddName, # version=NA, # taxId=NA # ) # # ### Add CortellisTAR target names ---- # toImport <- prtargets # colnames(toImport) <- c("id", "name") # BED:::loadBENames(d=toImport, be="Object", dbname=ddName) # # # ### Add "codes_for" edges ---- # toImport <- entXref %>% # select(externalId, targetId) # colnames(toImport) <- c("gid", "oid") # BED:::loadCodesFor( # d=toImport, # gdb="EntrezGene", # odb=ddName # ) ###############################################################################@ ## CortellisId targets "Cortellis_idtarget" coded by genes ---- cid <- get_MDB(.tkcon, "CortellisID") cidv <- as.character(cid$CortellisID_APIdump$date) alltargets <- cid$CortellisID_targets up <- cid$CortellisID_targetUniprot ### Get uniprot entrez gene identifiers ---- upEntrez <- c() for(org in listOrganisms()){ message(org) toAdd <- convBeIds( unique(up$uniprot), from="Peptide", from.source="Uniprot", from.org=org, to="Gene", to.source="EntrezGene", prefFilter=TRUE, restricted=FALSE ) upEntrez <- bind_rows( upEntrez, toAdd %>% as_tibble() %>% filter(!is.na(to)) %>% select("uniprot"=from, "entrez"=to) %>% mutate(organism=!!org) ) } fup <- left_join(upEntrez, up, by="uniprot", relationship = "many-to-many") ### Add Cortellis_idtarget database ---- ddName <- "Cortellis_idtarget" BED:::registerBEDB( name=ddName, description="Clarivate Analytics Cortellis Investigational Drugs", currentVersion=cidv, idURL=NA ) ### Add CortellisID targets ---- targets <- alltargets # %>% # filter(id %in% fup$target) toImport <- targets[, "id", drop=F] colnames(toImport) <- "id" BED:::loadBE( d=toImport, be="Object", dbname=ddName, version=NA, taxId=NA ) ### Add CortellisID target names ---- toImport <- targets colnames(toImport) <- c("id", "name") BED:::loadBENames(d=toImport, be="Object", dbname=ddName) ### Add "codes_for" edges ---- toImport <- fup %>% select(entrez, target) colnames(toImport) <- c("gid", "oid") BED:::loadCodesFor( d=toImport, gdb="EntrezGene", odb=ddName ) # ## CortellisTAR ID targets "Cortellis_idtarget" related to "Cortellis_target" # idtargets <- ctar$CortellisTAR_IDTargets # idtargets2 <- cid$CortellisID_targetXref %>% # filter(source=="Cortellis_target") %>% # select(cortellisTAR=xref, cortellisID=targetId) # idtargets <- bind_rows(idtargets, idtargets2) %>% # distinct(cortellisTAR, cortellisID) # toImport <- idtargets %>% # filter(cortellisTAR %in% prtargets$id, cortellisID %in% targets$id) %>% # select(id1=cortellisID, id2=cortellisTAR) # BED:::loadIsAssociatedTo( # d=toImport, # db1="Cortellis_idtarget", # db2="Cortellis_target", # be="Object" # ) ## CortellisID targets "Cortellis_idtarget" related to "MetaBase_object" idObjects <- cid$CortellisID_targetXref %>% filter(source=="MetaBase_object") %>% select(metabase=xref, cortellisID=targetId) mbo <- getBeIds(be="Object", source="MetaBase_object", restricted=FALSE) toImport <- idObjects %>% filter(metabase %in% mbo$id, cortellisID %in% targets$id) %>% select(id1=cortellisID, id2=metabase) BED:::loadIsAssociatedTo( d=toImport, db1="Cortellis_idtarget", db2="MetaBase_object", be="Object" )
Other information
Databases ID URL
otherIdURL <- list( "HGNC"='http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=%s', "miRBase"='http://www.mirbase.org/cgi-bin/mirna_entry.pl?acc=%s', "Vega_gene"='http://vega.sanger.ac.uk/id/%s', "UniGene"='https://www.ncbi.nlm.nih.gov/unigene?term=%s', "Vega_transcript"='http://vega.sanger.ac.uk/id/%s', "MGI"='http://www.informatics.jax.org/marker/MGI:%s', "Vega_translation"='http://vega.sanger.ac.uk/id/%s', "RGD"='https://rgd.mcw.edu/rgdweb/report/gene/main.html?id=%s', "MIM_GENE"='http://www.omim.org/entry/%s', "ZFIN_gene"='http://zfin.org/%s' ) for(db in names(otherIdURL)){ BED:::registerBEDB( name=db, idURL=otherIdURL[[db]] ) }
Load Lucene Indexes
Start: r Sys.time()
BED:::loadLuceneIndexes()
End: r Sys.time()
Session info
sessionInfo()
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.