trjhtr

I am evaluating if a list of pathways and its genes has pathways included in another pathways. So the way I thought about it is to create a function that does just a single comparison.

all_in <- function(x, y) 
 if (length(x) > length(y)) 
 0
 else 
 ifelse(all(x %in% y), 1, 0)
 

And then Vectorize to be able to use outer

all_in_vec <- Vectorize(all_in, vectorize.args = c("x", "y"))

And use it as follows:

nested <- outer(paths2genes, paths2genes, all_in_vec)

example:

paths2genes <- structure(list(`1430728` = c("10", "9"), `156580` = c("10", "3", 
"9"), `156582` = c("10", "9"), `194840` = c("2", "3"), `211859` = c("10", 
"9")), .Names = c("1430728", "156580", "156582", "194840", "211859"
))

Tests:

library(testthat)
expect_true(all(diag(nested) == 1L))
expect_equal(nested[1, 2], 1L)

However I have recently found that Vectorize should be avoided in favor of a built-in vectorization. Is there any built-in method I am missing? How can I vectorize this code ?

If you were to add a cat("hello") at the top of your all_in function, you would find that your function is called 25 times, once for each combination (pair) of pathways. So yes, despite having used Vectorized, it is still essentially a big old loop you have under the hood... Here is how I would write a vectorized function so the heavy-lifting function (%in%, or in my case, match) is only called once or twice:

all_in_outer <- function(list_x, list_y) 
 uniq_x <- unique(unlist(list_x, use.names = FALSE))
 len_x <- vapply(list_x, length, integer(1L))
 as_mat <- function(list_a, ids = uniq_x) 
 vec <- unlist(list_a, use.names = FALSE)
 len <- vapply(list_a, length, integer(1L))
 idx <- rep(seq_along(list_a), len)
 mat <- matrix(0L, nrow = length(list_a), ncol = length(ids),
 dimnames = list(names(list_a), ids))
 mat[cbind(idx, match(vec, ids))] <- 1L
 mat
 
 (as_mat(list_x) %*% t(as_mat(list_y)) == len_x) * 1


all_in_outer(paths2genes, paths2genes)
# 1430728 156580 156582 194840 211859
# 1430728 1 1 1 0 1
# 156580 0 1 0 0 0
# 156582 1 1 1 0 1
# 194840 0 0 0 1 0
# 211859 1 1 1 0 1

Some explanation: After we have found the list of unique genes across all pathways in the first argument (10, 9, 3, 2), we turn both arguments into a matrix A where A[i,j] = 1 if the pathway i contains the gene j, or 0 otherwise. For paths2genes, this matrix is as_mat(paths2genes):

# 10 9 3 2
# 1430728 1 1 0 0
# 156580 1 1 1 0
# 156582 1 1 0 0
# 194840 0 0 1 1
# 211859 1 1 0 0

Then, using matrix multiplication between two such matrices, you get the number of gene matches for each possible combination of pathways. You then just have to compare that number of matches with the length of the pathway: (as_mat(list_x) %*% t(as_mat(list_y)) == len_x).

In the conversion of each input into a matrix of zero and ones, see this particular line of code: mat[cbind(idx, match(vec, ids))] <- 1L. This is what fills the ones in the matrix. It is completely vectorized, via a single call to match.