# Calculating how similar two objects are according to a database

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
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")

• I suggest that you tell us in detail what the intended calculation is, so that we don't have to reverse-engineer your code to review it. – 200_success Sep 26 '16 at 14:52
• It's more code than most of us are willing to review, and there is no data to make it a reproducible example... I'd suggest you look into profiling your code so you can identify the bottle neck. – flodel Sep 26 '16 at 23:21
• @200_success I tried to explain more in detail what I want to do. (BTW thanks for the edit), flodel I did the profiling but I don't know how to reduce my bottle neck, which seems to be the matrix manipulation with [ and unique – llrs Sep 27 '16 at 8:04

I'll address the two main bottle necks in your code.

# First bottle neck

To help understand the issue, let's first remind ourselves the difference between the [ and [[ operators:

• when applied to a list, [ returns a sub-list, while [[ returns a list element.
• when applied to a data.frame (which is a form of a list), [ returns a data.frame, while [[ returns a vector (the data in a column).

Inside genes.info, where you do:

out <- unique(genes[genes[colm] == id, "Symbol"])


genes is a data.frame (i.e. a list), so genes[colm] is also a data.frame (a sub-list). When you then do genes[colm] == id, the == operator has to convert your one-column data.frame into a matrix before it can compare it to id, which is very expensive. This is where the matrix item at the top of your profile comes from. Instead, you meant to do:

out <- unique(genes[genes[[colm]] == id, "Symbol"])


where genes[[colm]] is a vector, so == does not have to do any conversion.

Note that you have a similar issue twice inside comb_biopath where you meant to use info[[by]] instead of info[by].

# Second bottle neck

With the iterative merge calls, you end up with pretty large data. What comes as pretty costly in that these merge calls, by default, also sort your data. That's where the second item in your profile (order) comes from. To get rid of it, which should not affect your results, add sort = FALSE to all your merge() calls.

On my machine, these two changes cut the computation times by roughly two thirds. I hope this puts you on the right track.

• Thank you very much Flodel for your review! I wasn't aware of the difference between [ and [[. This will definitely help me a lot, not just in this case. – llrs Sep 29 '16 at 6:54