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 matrixcomb2mat <- 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 matrixseq2mat <- 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 allcomb_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 100010    1  NA   NA100  NA   1   NA1000 NA  NA    1$kegg     10 100 100010    1   1    0100   1   1    01000  0   0    1$go             10        100      100010   1.00000000 0.05424063 0.1012658100  0.05424063 1.00000000 0.22087751000 0.10126582 0.22087746 1.0000000Profiling the function for this example is:$by.self                        self.time self.pct total.time total.pctmatrix                    10.34    33.77      10.46     34.16order                      7.98    26.06       7.98     26.06.Call                      2.74     8.95       2.74      8.95make.unique                1.58     5.16       1.76      5.75unique.default             1.14     3.72       1.14      3.72[.data.frame               1.08     3.53      22.10     72.18data.frame                 1.06     3.46       1.54      5.03paste                      1.02     3.33       1.32      4.31FUN                        0.90     2.94      13.86     45.26anyDuplicated.default      0.58     1.89       0.58      1.89edgeMatrix                 0.36     1.18       0.46      1.50unlist                     0.34     1.11       0.46      1.50match                      0.32     1.05       0.32      1.05as.character               0.18     0.59       0.18      0.59standardGeneric            0.16     0.52      30.62    100.00lapply                     0.10     0.33       3.32     10.84rbind                      0.08     0.26       2.68      8.75c                          0.06     0.20       0.06      0.20duplicated.default         0.06     0.20       0.06      0.20is.list                    0.06     0.20       0.06      0.20eval                       0.04     0.13       0.94      3.07.findInheritedMethods      0.04     0.13       0.20      0.65unique                     0.02     0.07      27.32     89.22[                          0.02     0.07      22.10     72.18merge.data.frame           0.02     0.07      12.02     39.26==                         0.02     0.07      11.26     36.77GOGraph                    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.26asNamespace                0.02     0.07       0.04      0.13isUpToDate                 0.02     0.07       0.04      0.13.sigLabel                  0.02     0.07       0.04      0.13<                          0.02     0.07       0.02      0.07check_valid                0.02     0.07       0.02      0.07degree                     0.02     0.07       0.02      0.07el                         0.02     0.07       0.02      0.07.identC                    0.02     0.07       0.02      0.07is.data.frame              0.02     0.07       0.02      0.07is.na                      0.02     0.07       0.02      0.07isNamespace                0.02     0.07       0.02      0.07L2Rchain.Lkeyname          0.02     0.07       0.02      0.07structure                  0.02     0.07       0.02      0.07$by.total                              total.time total.pct self.time self.pctstandardGeneric                  30.62    100.00      0.16     0.52bio.cor2                         30.62    100.00      0.00     0.00unique                           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.07FUN                              13.86     45.26      0.90     2.94react_genes                      12.98     42.39      0.00     0.00merge                            12.04     39.32      0.00     0.00merge.data.frame                 12.02     39.26      0.02     0.07==                               11.26     36.77      0.02     0.07Ops.data.frame                   11.24     36.71      0.00     0.00matrix                           10.46     34.16     10.34    33.77apply                            10.36     33.83      0.00     0.00genes.info                       10.36     33.83      0.00     0.00do.call                           8.92     29.13      0.00     0.00order                             7.98     26.06      7.98    26.06.local                            3.38     11.04      0.00     0.00lapply                            3.32     10.84      0.10     0.33comb2mat                          3.32     10.84      0.00     0.00func                              3.32     10.84      0.00     0.00sapply                            3.32     10.84      0.00     0.00dbGetQuery                        2.78      9.08      0.00     0.00dbQuery                           2.78      9.08      0.00     0.00sqliteGetQuery                    2.78      9.08      0.00     0.00simLL                             2.76      9.01      0.00     0.00.Call                             2.74      8.95      2.74     8.95sqliteFetch                       2.74      8.95      0.00     0.00rbind                             2.68      8.75      0.08     0.26comb_biopath                      2.62      8.56      0.00     0.00GOGraph                           2.54      8.30      0.02     0.07mget                              2.44      7.97      0.00     0.00dbUniqueVals                      2.34      7.64      0.00     0.00Lkeys                             2.18      7.12      0.00     0.00keys                              1.84      6.01      0.00     0.00make.unique                       1.76      5.75      1.58     5.16data.frame                        1.54      5.03      1.06     3.46cbind                             1.54      5.03      0.00     0.00paste                             1.32      4.31      1.02     3.33flatten                           1.32      4.31      0.00     0.00unique.data.frame                 1.32      4.31      0.00     0.00unique.default                    1.14      3.72      1.14     3.72duplicated                        1.02      3.33      0.00     0.00duplicated.data.frame             0.98      3.20      0.00     0.00eval                              0.94      3.07      0.04     0.13toTable                           0.94      3.07      0.00     0.00dbSelectFromL2Rchain              0.78      2.55      0.00     0.00sp.between                        0.64      2.09      0.00     0.00anyDuplicated.default             0.58      1.89      0.58     1.89anyDuplicated                     0.58      1.89      0.00     0.00s.path                            0.56      1.83      0.00     0.00dijkstra.sp                       0.54      1.76      0.00     0.00edgeMatrix                        0.46      1.50      0.36     1.18unlist                            0.46      1.50      0.34     1.11as.list                           0.46      1.50      0.00     0.00data.row.names                    0.42      1.37      0.00     0.00match                             0.32      1.05      0.32     1.05<Anonymous>                       0.28      0.91      0.00     0.00initialize                        0.22      0.72      0.00     0.00new                               0.22      0.72      0.00     0.00.findInheritedMethods             0.20      0.65      0.04     0.13as.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.00Rkeys                             0.16      0.52      0.00     0.00graphNEL_init_edgeL_weights       0.14      0.46      0.00     0.00eg2gofun                          0.12      0.39      0.00     0.00.getWHEC                          0.12      0.39      0.00     0.00getPartialSubmap                  0.10      0.33      0.00     0.00simLP                             0.10      0.33      0.00     0.00.getClassFromCache                0.08      0.26      0.02     0.07edgeData<-                        0.08      0.26      0.00     0.00edgeWeights                       0.08      0.26      0.00     0.00c                                 0.06      0.20      0.06     0.20duplicated.default                0.06      0.20      0.06     0.20is.list                           0.06      0.20      0.06     0.20getClass                          0.06      0.20      0.00     0.00isDirected                        0.06      0.20      0.00     0.00outerLabels                       0.06      0.20      0.00     0.00asNamespace                       0.04      0.13      0.02     0.07isUpToDate                        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.00dbGetInfo                         0.04      0.13      0.00     0.00dbHasCompleted                    0.04      0.13      0.00     0.00edgemode                          0.04      0.13      0.00     0.00elNamed                           0.04      0.13      0.00     0.00.formatList                       0.04      0.13      0.00     0.00get                               0.04      0.13      0.00     0.00.makeSQL                          0.04      0.13      0.00     0.00match.fun                         0.04      0.13      0.00     0.00<                                 0.02      0.07      0.02     0.07check_valid                       0.02      0.07      0.02     0.07degree                            0.02      0.07      0.02     0.07el                                0.02      0.07      0.02     0.07.identC                           0.02      0.07      0.02     0.07is.data.frame                     0.02      0.07      0.02     0.07is.na                             0.02      0.07      0.02     0.07isNamespace                       0.02      0.07      0.02     0.07L2Rchain.Lkeyname                 0.02      0.07      0.02     0.07structure                         0.02      0.07      0.02     0.07anyStrings                        0.02      0.07      0.00     0.00as.factor                         0.02      0.07      0.00     0.00colnames                          0.02      0.07      0.00     0.00deparse                           0.02      0.07      0.00     0.00.deparseOpts                      0.02      0.07      0.00     0.00edgeData                          0.02      0.07      0.00     0.00edgeDataDefaults<-                0.02      0.07      0.00     0.00edgeParser                        0.02      0.07      0.00     0.00extends                           0.02      0.07      0.00     0.00factor                            0.02      0.07      0.00     0.00.getAllEdges                      0.02      0.07      0.00     0.00getClassDef                       0.02      0.07      0.00     0.00getw                              0.02      0.07      0.00     0.00ifelse                            0.02      0.07      0.00     0.00.inheritedArgsExpression          0.02      0.07      0.00     0.00is                                0.02      0.07      0.00     0.00L2Rchain.Rattribnames<-           0.02      0.07      0.00     0.00Lkeyname                          0.02      0.07      0.00     0.00match.arg                         0.02      0.07      0.00     0.00nrow                              0.02      0.07      0.00     0.00Rattribnames<-                    0.02      0.07      0.00     0.00Rattribnames                      0.02      0.07      0.00     0.00sort.list                         0.02      0.07      0.00     0.00split                             0.02      0.07      0.00     0.00split.default                     0.02      0.07      0.00     0.00validGraph                        0.02      0.07      0.00     0.00validityMethod                    0.02      0.07      0.00     0.00validObject                       0.02      0.07      0.00     0.00$sample.interval[1] 0.02$sampling.time[1] 30.62The 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)
Calculating how similar two objects are according to a database
matrix;graph;combinatorics;r;bioinformatics
I'll address the two main bottle necks in your code.First bottle neckTo 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 neckWith 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.