我有一个命名列表,它代表了一组生物路径,其中的名称是路径名,列表中的载体是属于该通路的蛋白质。一个小例子是:
ann <- structure(list(`GO:0000010` = c("Q33DR2", "Q9CZQ1", "D6RHT8",
"F6ZCX7", "B8JJX0", "Q33DR3", "F6T4Z4", "E0CYM9"), `GO:0000016` = c("Q5XLR9",
"Q3TZ78", "F8VPT3"), `GO:0000026` = c("Q8BTP0", "Q3TZM9", "A0A077K846",
"F6R220", "A0A077K9W9"), `GO:0000032` = c("Q924M7", "Q3V100",
"F6Q3K8", "Q921Z9"), `GO:0000033` = c("Q9DBE8", "F6RBY3", "Q8BMZ4",
"Q8K2A8", "F6XUH0", "D6RCW8", "Q6P8H8", "Q3URN2")), .Names = c("GO:0000010",
"GO:0000016", "GO:0000026", "GO:0000032", "GO:0000033"))
我对成对的路径感兴趣:
pairs <- t(combn(names(ann), 2))
对于每一对路径,我都希望得到所有可能的蛋白质组合,其中蛋白质#1在路径#1,蛋白质#2在路径#2中。理想的输出是一个两列矩阵的列表,其中第1列包含路径#1中的蛋白质,第2列包含路径#2中的蛋白质。到目前为止,我有这样的结论:
protein_pairs <- purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid(ann[[.x]], ann[[.y]])))
但是,由于我感兴趣的对的总数相当大(通常>1,000),因此对所有可能的对进行映射expand.grid
需要很长的时间--按小时计算。
是否有更快的方法从这份清单中获得每对生物途径中所有可能的蛋白质组合?
发布于 2018-06-25 16:06:28
如果您正在寻找速度,那么您可以很容易地找到一个Rcpp
版本:
// [[Rcpp::export]]
CharacterMatrix fast2Expand(CharacterVector x, CharacterVector y) {
unsigned long int lenX = x.size(), lenY = y.size();
CharacterMatrix result = no_init_matrix(lenX * lenY, 2);
for (std::size_t i = 0, count = 0; i < lenY; ++i) {
for (std::size_t j = 0; j < lenX; ++j, ++count){
result(count, 0) = x[j];
result(count, 1) = y[i];
}
}
return result;
}
它是关于10x
比原始版本快,20%
比rep.int
版本快(对于本例来说):
microbenchmark(OP = purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid(ann[[.x]], ann[[.y]]))),
Rcpp = purrr::map2(pairs[, 1], pairs[, 2], ~ fast2Expand(ann[[.x]], ann[[.y]])),
repInt = purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid.jc(ann[[.x]], ann[[.y]]))))
Unit: microseconds
expr min lq mean median uq max neval
OP 1104.700 1136.4370 1536.4048 1188.9990 1481.4940 6730.960 100
Rcpp 105.505 126.9975 149.9009 138.1195 150.2015 663.146 100
repInt 133.044 151.0175 223.9815 165.5435 203.5335 1269.194 100
下面是一个基于OP示例的人为示例,其目的完全是为了比较效率:
annBig <- lapply(1:5, function(x) rep(ann[[x]], 100))
names(annBig) <- names(ann)
microbenchmark(OP = purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid(annBig[[.x]], annBig[[.y]]))),
Rcpp = purrr::map2(pairs[, 1], pairs[, 2], ~ fast2Expand(annBig[[.x]], annBig[[.y]])),
repInt = purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid.jc(annBig[[.x]], annBig[[.y]]))), times = 20)
Unit: milliseconds
expr min lq mean median uq max neval
OP 522.56536 533.39393 562.60750 555.45345 588.4514 640.8584 20
Rcpp 48.12683 56.17155 92.30095 92.23838 125.8065 142.2949 20
repInt 80.28625 107.32329 140.32793 152.13732 160.9656 193.1310 20
发布于 2018-06-25 15:54:19
我认为rep.int()
的工作速度要快得多,就像在另一个question:中所说的那样
尝试以下几点:
expand.grid.jc <- function(seq1,seq2) {
cbind(Var1 = rep.int(seq1, length(seq2)),
Var2 = rep.int(seq2, rep.int(length(seq1),length(seq2))))
}
protein_pairs <- purrr::map2(pairs[, 1], pairs[, 2], ~ as.matrix(expand.grid.jc(ann[[.x]], ann[[.y]])))
干杯!
https://stackoverflow.com/questions/51027078
复制相似问题