I want to calculate how similar two objects are according to a database. However the code is very slow; it takes around 2 minutes to analyze just 3 objects. How can I speed it up? I have tried to avoid for loops, and I'm also avoiding recalculating things.
This functions returns the ith comparison of the n elements in r groups, in my case, always 2.
".combinadic" <- function(n, r, i) {
# http://msdn.microsoft.com/en-us/library/aa289166(VS.71).aspx
# http://en.wikipedia.org/wiki/Combinadic
n0 <- length(n)
if (i < 1 | i > choose(n0,r)) {
stop("'i' must be 0 < i <= n0!/(n0-r)!")
}
largestV <- function(n, r, i) {
v <- n # Adjusted for one-based indexing
while (choose(v,r) >= i) { # Adjusted for one-based indexing
v <- v - 1
}
return(v)
}
res <- rep(NA,r)
for (j in 1:r) {
res[j] <- largestV(n0,r,i)
i <- i - choose(res[j],r)
n0 <- res[j]
r <- r - 1
}
res <- res + 1
res <- n[res]
return(res)
}
This functions compares graphs by calculating how many nodes are in both graphs, or if a character vector is given then it compares how many are in both objects.
compare_graphs <- function(g1, g2){
# Function to estimate how much two graphs overlap by looking if the nodes
# are the same
# Check which case are we using
if (is(g1, "graph") & is(g2, "graph")) {
prot1 <- nodes(g1)
prot2 <- nodes(g2)
if (length(prot1) == 0 | length(prot2) == 0) {
return(NA)
}
} else if (is(g1, "graph") & is.character(g2)) {
prot1 <- nodes(g1)
prot2 <- g2
} else if (is(g2, "graph") & is.character(g1)) {
prot2 <- nodes(g2)
prot1 <- g1
} else {
prot1 <- g1
prot2 <- g2
}
score <- (length(intersect(prot1, prot2)))*2/(
length(prot2) + length(prot1))
score
}
This functions calculates the degree of overlap of Gene Ontologies under the class of Biological Process (BP), it is different from compare graph because the structure of this graph or its subgraphs is known, is a DAG, and I use that to compare the paths of both graphs.
# Calculates the degree of overlap of the GO BP ontologies of entrez ids.
go_cor <- function(e_a, e_b, chip = "hgu133plus2.db", mapfun = NULL, ...){
# https://support.bioconductor.org/p/85702/#85732
if (is.na(e_a) | is.na(e_b)) {
return(NA)
}
# Ensure proper format
e_a <- as.character(e_a)
e_b <- as.character(e_b)
if (mapfun) {
mapfunc <- function(x) mget(x, revmap(org.Hs.egGO2EG), ifnotfound = NA)
LP <- simLL(e_a, e_b, "BP", measure = "LP", mapfun = mapfunc, ...)
UI <- simLL(e_a, e_b, "BP", measure = "UI", mapfun = mapfunc, ...)
} else {
LP <- simLL(e_a, e_b, "BP", measure = "LP", chip = chip, ...)
UI <- simLL(e_a, e_b, "BP", measure = "UI", chip = chip, ...)
}
if (length(LP) > 1 | length(UI) > 1) {
if (is.na(LP["sim"]) | is.na(UI["sim"])) {
return(NA)
}
} else if (is.na(LP) | is.na(UI)) {
return(NA)
}
s.path <- function(ig){
# The longest of the shortest path of a graph
lfi <- leaves(ig, "in")
degs <- degree(ig)
root <- names(degs$outDegree)[degs$outDegree == 0]
paths <- sp.between(ig, lfi, root)
plens <- subListExtract(paths, "length", simplify = TRUE)
max(plens)
}
# Calculates the score taking into account the size and the middle path
# Taking advantage of the fact that in GO there is a root and leaves
# UI: Union intersect, is the size of the intersection of the node
# sets divided by the size of the union of the node sets
# LP: longest path, is the longest path in the intersection graph of
# the two supplied graph.
(UI$sim/LP$sim)*max(s.path(LP$g1), s.path(LP$g2))
}
Funcion to perform "efficiently" the conversion from combinations to symmetric matrix
comb2mat <- function(input, func, ...){
# Perform all the combinations of 2 from the input
cobs <- list()
for (i in 1:length(input)) {
cobs[[i]] <- .combinadic(input, 2, i)
}
# cobs <- combn(input, 2)
N <- sapply(cobs, function(x){func(x[1], x[2], ...)})
# Function that performs the calculus
# N <- seq_len(ncol(combs))
out <- matrix(ncol = length(input), nrow = length(input))
out[lower.tri(out)] <- N
out <- t(out)
out[lower.tri(out)] <- N
out <- t(out)
diag(out) <- 1
rownames(out) <- colnames(out) <- input
return(out)
}
Transform a vector to a symmetric matrix, using dat to fill it and x to set names on the matrix
seq2mat <- function(x, dat) {
out <- matrix(ncol = length(x), nrow = length(x))
out[lower.tri(out)] <- unlist(dat)
out <- t(out)
out[lower.tri(out)] <- unlist(dat)
out <- t(out)
diag(out) <- 1
rownames(out) <- colnames(out) <- x
return(out)
}
Extract all the ids of a given column of a df (pathways) for each element on the combination and compare them all
comb_biopath <- function(comb, info, by, biopath){
# react_path <- apply(comb, 2, function(y){
a <- unique(info[info[by] == comb[1], biopath])
a <- a[a != ""]
a <- a[!is.na(a)]
b <- unique(info[info[by] == comb[2], biopath])
b <- b[b != ""]
b <- b[!is.na(b)]
if (all(sapply(a, is.na)) | all(sapply(b, is.na))) {
return(NA)
} else if (length(a) == 0 | length(b) == 0) {
return(NA)
}
expand.grid(a, b)
# })
}
Using data correlates biologically pairwise all elements on x, ids define which type of identification is used, and the other arguments define which comparison is made.
bio.cor2 <- function(x, ids = "Entrez Gene",
go = FALSE, react = TRUE, kegg = FALSE, all = FALSE) {
# Using data correlates biologically two genes or probes
# From the graphite package
# x should be entrez id
# or change the internals from "Entrez Gene" to "Symbols"
if (!ids %in% c("Entrez Gene", "Symbol")) {
stop("Please check the input of genes in Symbol or Entrez Gene format")
}
if (all) {
go <- kegg <- react <- all
}
genes_id <- x
# Obtain data from the annotation packages
gene.symbol <- unique(toTable(org.Hs.egSYMBOL2EG))
gene.symbol <- gene.symbol[gene.symbol$gene_id %in% genes_id, ]
colnames(gene.symbol) <- c("Entrez Gene", "Symbol")
if (kegg) {
# Obtain data
gene.kegg <- unique(toTable(org.Hs.egPATH2EG))
colnames(gene.kegg) <- c("Entrez Gene", "KEGG")
# Merge data
genes <- merge(gene.symbol, gene.kegg, all = TRUE)
kegg.bio <- rep(NA, choose(length(genes_id), 2))
}
if (react) {
gene.reactome <- unique(toTable(reactomePATHID2EXTID))
colnames(gene.reactome) <- c("Entrez Gene", "Reactome")
genes <- unique(merge(genes, gene.reactome, all = TRUE))
react.bio <- rep(NA, choose(length(genes_id), 2))
}
## Calculate for each combination the values
# Calculate the GO correlation
if (go) {
go_mat <- comb2mat(genes_id, go_cor, mapfun = TRUE)
}
if (kegg | react) {
# Calculate the pathways correlation
for (i in 1:choose(length(genes_id), 2)) {
comb <- .combinadic(genes_id, 2, i)
# Kegg calculus
if (kegg) {
kegg.bio[i] <- react_genes(comb, genes, "KEGG", ids)
}
# Reactome calculus
if (react) {
react.bio[i] <- react_genes(comb, genes, "Reactome", ids)
}
}
}
if (react) {
react_mat <- seq2mat(genes_id, react.bio)
}
if (kegg) {
kegg_mat <- seq2mat(genes_id, kegg.bio)
}
if (all) {
cor_mat <- list(reactome = react_mat, kegg = kegg_mat, go = go_mat)
} else if (kegg & react) {
cor_mat <- list(reactome = react_mat, kegg = kegg_mat)
} else if (kegg & go) {
cor_mat <- list(kegg = kegg_mat, go = go_mat)
} else if (go & react) {
cor_mat <- list(reactome = react_mat, go = go_mat)
}
return(cor_mat)
}
Extract which genes are from which pathway on a db.
genes.info <- function(genes, colm, id) {
# Genes is the df, colm is the column you want, id is the id of the pathway
out <- unique(genes[genes[colm] == id, "Symbol"])
out[!is.na(out)]
}
Calculates for all the combinations of pathway the max score shared between genes, on the comb (pair in question)
react_genes <- function(comb, genes, react, id) {
if (!react %in% colnames(genes)) {
stop("Please check which type of reaction you want")
}
react_path <- comb_biopath(comb, genes, id, react)
# Check that we have pathways info for this combination
if (is.null(react_path)) {
return(NA)
} else if (length(react_path) == 2) {
if (nrow(react_path) == 0) {
return(NA)
}
} else if (is.na(react_path)) {
return(NA)
}
react <- apply(react_path, 1, function(x){
a <- genes.info(genes, react, x[1])
b <- genes.info(genes, react, x[2])
out <- compare_graphs(a, b)
out
})
# If NA returns NA
if (length(react) != sum(is.na(react))) {
out <- max(react, na.rm = TRUE)
} else {
out <- NA
}
out
}
The main function is bio.cor2, which can be used as bio.cor2(c("10", "100", "1000"), all = TRUE). Or it can use any other numbers as bio.cor2(c("1000", "135489", "4892346"), all = TRUE).
The result of bio.cor2 is a list of matrix of nxn dimensions where n is the length of the first argument. The numeric values can be from 0 to 1 or if no information is available NA:
> a <- bio.cor2(c("10", "100", "1000"), all = TRUE)
> a
$reactome
10 100 1000
10 1 NA NA
100 NA 1 NA
1000 NA NA 1
$kegg
10 100 1000
10 1 1 0
100 1 1 0
1000 0 0 1
$go
10 100 1000
10 1.00000000 0.05424063 0.1012658
100 0.05424063 1.00000000 0.2208775
1000 0.10126582 0.22087746 1.0000000
Profiling the function for this example is:
$by.self
self.time self.pct total.time total.pct
"matrix" 10.34 33.77 10.46 34.16
"order" 7.98 26.06 7.98 26.06
".Call" 2.74 8.95 2.74 8.95
"make.unique" 1.58 5.16 1.76 5.75
"unique.default" 1.14 3.72 1.14 3.72
"[.data.frame" 1.08 3.53 22.10 72.18
"data.frame" 1.06 3.46 1.54 5.03
"paste" 1.02 3.33 1.32 4.31
"FUN" 0.90 2.94 13.86 45.26
"anyDuplicated.default" 0.58 1.89 0.58 1.89
"edgeMatrix" 0.36 1.18 0.46 1.50
"unlist" 0.34 1.11 0.46 1.50
"match" 0.32 1.05 0.32 1.05
"as.character" 0.18 0.59 0.18 0.59
"standardGeneric" 0.16 0.52 30.62 100.00
"lapply" 0.10 0.33 3.32 10.84
"rbind" 0.08 0.26 2.68 8.75
"c" 0.06 0.20 0.06 0.20
"duplicated.default" 0.06 0.20 0.06 0.20
"is.list" 0.06 0.20 0.06 0.20
"eval" 0.04 0.13 0.94 3.07
".findInheritedMethods" 0.04 0.13 0.20 0.65
"unique" 0.02 0.07 27.32 89.22
"[" 0.02 0.07 22.10 72.18
"merge.data.frame" 0.02 0.07 12.02 39.26
"==" 0.02 0.07 11.26 36.77
"GOGraph" 0.02 0.07 2.54 8.30
".makeSQLchunks" 0.02 0.07 0.18 0.59
".getClassFromCache" 0.02 0.07 0.08 0.26
"asNamespace" 0.02 0.07 0.04 0.13
"isUpToDate" 0.02 0.07 0.04 0.13
".sigLabel" 0.02 0.07 0.04 0.13
"<" 0.02 0.07 0.02 0.07
"check_valid" 0.02 0.07 0.02 0.07
"degree" 0.02 0.07 0.02 0.07
"el" 0.02 0.07 0.02 0.07
".identC" 0.02 0.07 0.02 0.07
"is.data.frame" 0.02 0.07 0.02 0.07
"is.na" 0.02 0.07 0.02 0.07
"isNamespace" 0.02 0.07 0.02 0.07
"L2Rchain.Lkeyname" 0.02 0.07 0.02 0.07
"structure" 0.02 0.07 0.02 0.07
$by.total
total.time total.pct self.time self.pct
"standardGeneric" 30.62 100.00 0.16 0.52
"bio.cor2" 30.62 100.00 0.00 0.00
"unique" 27.32 89.22 0.02 0.07
"[.data.frame" 22.10 72.18 1.08 3.53
"[" 22.10 72.18 0.02 0.07
"FUN" 13.86 45.26 0.90 2.94
"react_genes" 12.98 42.39 0.00 0.00
"merge" 12.04 39.32 0.00 0.00
"merge.data.frame" 12.02 39.26 0.02 0.07
"==" 11.26 36.77 0.02 0.07
"Ops.data.frame" 11.24 36.71 0.00 0.00
"matrix" 10.46 34.16 10.34 33.77
"apply" 10.36 33.83 0.00 0.00
"genes.info" 10.36 33.83 0.00 0.00
"do.call" 8.92 29.13 0.00 0.00
"order" 7.98 26.06 7.98 26.06
".local" 3.38 11.04 0.00 0.00
"lapply" 3.32 10.84 0.10 0.33
"comb2mat" 3.32 10.84 0.00 0.00
"func" 3.32 10.84 0.00 0.00
"sapply" 3.32 10.84 0.00 0.00
"dbGetQuery" 2.78 9.08 0.00 0.00
"dbQuery" 2.78 9.08 0.00 0.00
"sqliteGetQuery" 2.78 9.08 0.00 0.00
"simLL" 2.76 9.01 0.00 0.00
".Call" 2.74 8.95 2.74 8.95
"sqliteFetch" 2.74 8.95 0.00 0.00
"rbind" 2.68 8.75 0.08 0.26
"comb_biopath" 2.62 8.56 0.00 0.00
"GOGraph" 2.54 8.30 0.02 0.07
"mget" 2.44 7.97 0.00 0.00
"dbUniqueVals" 2.34 7.64 0.00 0.00
"Lkeys" 2.18 7.12 0.00 0.00
"keys" 1.84 6.01 0.00 0.00
"make.unique" 1.76 5.75 1.58 5.16
"data.frame" 1.54 5.03 1.06 3.46
"cbind" 1.54 5.03 0.00 0.00
"paste" 1.32 4.31 1.02 3.33
"flatten" 1.32 4.31 0.00 0.00
"unique.data.frame" 1.32 4.31 0.00 0.00
"unique.default" 1.14 3.72 1.14 3.72
"duplicated" 1.02 3.33 0.00 0.00
"duplicated.data.frame" 0.98 3.20 0.00 0.00
"eval" 0.94 3.07 0.04 0.13
"toTable" 0.94 3.07 0.00 0.00
"dbSelectFromL2Rchain" 0.78 2.55 0.00 0.00
"sp.between" 0.64 2.09 0.00 0.00
"anyDuplicated.default" 0.58 1.89 0.58 1.89
"anyDuplicated" 0.58 1.89 0.00 0.00
"s.path" 0.56 1.83 0.00 0.00
"dijkstra.sp" 0.54 1.76 0.00 0.00
"edgeMatrix" 0.46 1.50 0.36 1.18
"unlist" 0.46 1.50 0.34 1.11
"as.list" 0.46 1.50 0.00 0.00
"data.row.names" 0.42 1.37 0.00 0.00
"match" 0.32 1.05 0.32 1.05
"<Anonymous>" 0.28 0.91 0.00 0.00
"initialize" 0.22 0.72 0.00 0.00
"new" 0.22 0.72 0.00 0.00
".findInheritedMethods" 0.20 0.65 0.04 0.13
"as.character" 0.18 0.59 0.18 0.59
".makeSQLchunks" 0.18 0.59 0.02 0.07
"%in%" 0.18 0.59 0.00 0.00
"Rkeys" 0.16 0.52 0.00 0.00
"graphNEL_init_edgeL_weights" 0.14 0.46 0.00 0.00
"eg2gofun" 0.12 0.39 0.00 0.00
".getWHEC" 0.12 0.39 0.00 0.00
"getPartialSubmap" 0.10 0.33 0.00 0.00
"simLP" 0.10 0.33 0.00 0.00
".getClassFromCache" 0.08 0.26 0.02 0.07
"edgeData<-" 0.08 0.26 0.00 0.00
"edgeWeights" 0.08 0.26 0.00 0.00
"c" 0.06 0.20 0.06 0.20
"duplicated.default" 0.06 0.20 0.06 0.20
"is.list" 0.06 0.20 0.06 0.20
"getClass" 0.06 0.20 0.00 0.00
"isDirected" 0.06 0.20 0.00 0.00
"outerLabels" 0.06 0.20 0.00 0.00
"asNamespace" 0.04 0.13 0.02 0.07
"isUpToDate" 0.04 0.13 0.02 0.07
".sigLabel" 0.04 0.13 0.02 0.07
":::" 0.04 0.13 0.00 0.00
".contextualizeColnames" 0.04 0.13 0.00 0.00
"dbGetInfo" 0.04 0.13 0.00 0.00
"dbHasCompleted" 0.04 0.13 0.00 0.00
"edgemode" 0.04 0.13 0.00 0.00
"elNamed" 0.04 0.13 0.00 0.00
".formatList" 0.04 0.13 0.00 0.00
"get" 0.04 0.13 0.00 0.00
".makeSQL" 0.04 0.13 0.00 0.00
"match.fun" 0.04 0.13 0.00 0.00
"<" 0.02 0.07 0.02 0.07
"check_valid" 0.02 0.07 0.02 0.07
"degree" 0.02 0.07 0.02 0.07
"el" 0.02 0.07 0.02 0.07
".identC" 0.02 0.07 0.02 0.07
"is.data.frame" 0.02 0.07 0.02 0.07
"is.na" 0.02 0.07 0.02 0.07
"isNamespace" 0.02 0.07 0.02 0.07
"L2Rchain.Lkeyname" 0.02 0.07 0.02 0.07
"structure" 0.02 0.07 0.02 0.07
"anyStrings" 0.02 0.07 0.00 0.00
"as.factor" 0.02 0.07 0.00 0.00
"colnames" 0.02 0.07 0.00 0.00
"deparse" 0.02 0.07 0.00 0.00
".deparseOpts" 0.02 0.07 0.00 0.00
"edgeData" 0.02 0.07 0.00 0.00
"edgeDataDefaults<-" 0.02 0.07 0.00 0.00
"edgeParser" 0.02 0.07 0.00 0.00
"extends" 0.02 0.07 0.00 0.00
"factor" 0.02 0.07 0.00 0.00
".getAllEdges" 0.02 0.07 0.00 0.00
"getClassDef" 0.02 0.07 0.00 0.00
"getw" 0.02 0.07 0.00 0.00
"ifelse" 0.02 0.07 0.00 0.00
".inheritedArgsExpression" 0.02 0.07 0.00 0.00
"is" 0.02 0.07 0.00 0.00
"L2Rchain.Rattribnames<-" 0.02 0.07 0.00 0.00
"Lkeyname" 0.02 0.07 0.00 0.00
"match.arg" 0.02 0.07 0.00 0.00
"nrow" 0.02 0.07 0.00 0.00
"Rattribnames<-" 0.02 0.07 0.00 0.00
"Rattribnames" 0.02 0.07 0.00 0.00
"sort.list" 0.02 0.07 0.00 0.00
"split" 0.02 0.07 0.00 0.00
"split.default" 0.02 0.07 0.00 0.00
"validGraph" 0.02 0.07 0.00 0.00
"validityMethod" 0.02 0.07 0.00 0.00
"validObject" 0.02 0.07 0.00 0.00
$sample.interval
[1] 0.02
$sampling.time
[1] 30.62
The information is extracted from the following libraries from Bioconductor:
library("biomaRt")
library("hgu133plus2.db")
library("testthat")
library("GOstats")
library("graphite")
library("WGCNA")
library("KEGGgraph")
library("KEGG.db")
library("RBGL")
library("org.Hs.eg.db")
library("graph")
library("Rgraphviz")
library("reactome.db")
