Сумма всех попарных продуктов ряда как двухсторонняя матрица

Question

Я ищу данные о высокопроизводительных генах и делаю тип корреляционного анализа на основе байесовской статистики. Одна из вещей, которую мне нужно сделать, - найти каждую попарную комбинацию продуктов в наборе данных и найти сумму каждой результирующей строки.

Так, например, для набора данных матрицы высокой пропускной Dataset

(Dataset <- structure(list(`Condition 1` = c(1L, 3L, 2L, 2L), `Condition 2` = c(2L, 1L, 7L, 2L), `Condition 3` = c(4L, 1L, 2L, 5L)), .Names = c("Condition 1", "Condition 2", "Condition 3"), class = "data.frame", row.names = c("Gene A", "Gene B", "Gene C", "Gene D"))) 
     Condition 1 Condition 2 Condition 3 
Gene A   1   2    4 
Gene B   3   1    1 
Gene C   2   7    2 
Gene D   2   2    5 

Прежде всего я хочу, чтобы умножить каждую возможную пару строк вместе, чтобы получить следующую матрицу под названием Comb:

   Condition 1 Condition 2 Condition 3 
Gene A Gene A   1   4   9 
Gene A Gene B   3   2   4 
Gene A Gene C   2   14   8 
Gene A Gene D   2   4   20 
Gene B Gene B   9   1   1 
Gene B Gene C   6   7   2 
Gene B Gene D   6   2   5 
Gene C Gene C   4   49   4 
Gene C Gene D   4   14   10 
Gene D Gene D   4   4   25 

После I хотите найти суммы строк для каждого продукта и получить суммы в виде матрицы (которую я назову CombSums):

  Gene A  Gene B  Gene C  Gene D 
Gene A   NA   10   24   26 
Gene B   10   NA   15   13 
Gene C   24   15   NA   28 
Gene D   26   13   28   NA 

Когда я попытался сделать это, лучшее, что я мог придумать было

combs <- combn(seq_len(nrow(Dataset)), 2) 
Comb <- Dataset[combs[1,], ] * Dataset[combs[2,], ] 
rownames(Comb) <- apply(combn(rownames(Comb), 2), 2, paste, collapse = " ") 
CombSums <- rowSums(Comb) 

Какой бы дает мне суммы в виде списка, например, как показано ниже:

    [1,] 
Gene A Gene B  10 
Gene A Gene C  24 
Gene A Gene D  26 
Gene B Gene C  15 
Gene B Gene D  13 
Gene C Gene D  28 

К сожалению, я хотите, чтобы это как двухсторонняя матрица, а не список, так что это не работает, поэтому, если кто-нибудь может предложить способ получить суммы в виде матрицы, это будет большой помощью.

Answer 1

Если скорость является важным фактором (например, если вы обрабатываете огромную матрицу), вы могли бы найти реализацию Rcpp полезной. Это только заполняет верхнюю треугольную часть матрицы.

library(Rcpp) 
cppFunction(
"NumericMatrix josilberRcpp(NumericMatrix x) { 
    const int nr = x.nrow(); 
    const int nc = x.ncol(); 
    NumericMatrix y(nr, nr); 
    for (int col=0; col < nc; ++col) { 
    for (int i=0; i < nr; ++i) { 
     for (int j=i; j < nr; ++j) { 
     y(i, j) += x(i, col) * x(j, col); 
     } 
    } 
    } 
    return y; 
}") 
josilberRcpp(as.matrix(Dataset)) 
#  [,1] [,2] [,3] [,4] 
# [1,] 21 9 24 26 
# [2,] 0 11 15 13 
# [3,] 0 0 57 28 
# [4,] 0 0 0 33 

Бенчмаркинг предоставляется в моем другом ответе. Обратите внимание, что в бенчмаркинг не входит время компиляции с использованием cppFunction, что может быть весьма значительным. Поэтому эта реализация, вероятно, полезна только для очень больших входов или когда вам нужно использовать эту функцию много раз.

Answer 2

Вы можете сделать это путем вычисления попарных произведений для каждого столбца в первоначальном кадре данных с lapply и outer, а затем вы можете добавить все те попарные произведения вместе с Reduce и +.

Reduce("+", lapply(dat, function(x) outer(x, x))) 
#  [,1] [,2] [,3] [,4] 
# [1,] 21 9 24 26 
# [2,] 9 11 15 13 
# [3,] 24 15 57 28 
# [4,] 26 13 28 33 

вариация на эту тему, которая менее интенсивно использующих память (потому что не нужно хранить матрицу каждого столбца в то же время), но более типирование будет:

ret <- outer(dat[,1], dat[,1]) 
for (i in 2:ncol(dat)) ret <- ret + outer(dat[,i], dat[,i]) 
ret 
#  [,1] [,2] [,3] [,4] 
# [1,] 21 9 24 26 
# [2,] 9 11 15 13 
# [3,] 24 15 57 28 
# [4,] 26 13 28 33 

Вот тест из предложенных подходов до сих пор на раме 100 х 100 данных:

# Larger dataset 
set.seed(144) 
dat <- as.data.frame(matrix(rnorm(10000), nrow=100)) 

josilber <- function(dat) Reduce("+", lapply(dat, function(x) outer(x, x))) 
josilber2 <- function(dat) { 
    ret <- outer(dat[,1], dat[,1]) 
    for (i in 2:ncol(dat)) ret <- ret + outer(dat[,i], dat[,i]) 
    ret 
} 
frank <- function(DF) { 
    mat <- as.matrix(DF) 
    pairs <- combn(1:nrow(DF),2) 
    vals <- rowSums(mat[pairs[1,],]*mat[pairs[2,],]) 
    res <- matrix(,nrow(DF),nrow(DF),dimnames=list(rownames(DF),rownames(DF))) 
    res[lower.tri(res)] <- vals 
    res 
} 

library(microbenchmark) 
microbenchmark(josilber(dat), josilber2(dat), josilberRcpp(as.matrix(dat)), frank(dat)) 
# Unit: microseconds 
#       expr  min  lq  mean median  uq  max neval 
#     josilber(dat) 6867.499 45437.277 45506.731 46372.576 47549.834 85494.063 100 
#    josilber2(dat) 6831.692 7982.539 10245.459 9109.023 10883.965 50612.600 100 
# josilberRcpp(as.matrix(dat)) 989.592 1112.316 1290.617 1204.388 1483.638 2384.348 100 
#     frank(dat) 13043.912 53369.804 52488.997 53921.402 54855.583 62566.730 100 

Answer 3

Использование combn, вы можете не делать лишние вычисления:

mat <- as.matrix(DF) 

pairs <- combn(1:nrow(DF),2) 

vals <- rowSums(mat[pairs[1,],]*mat[pairs[2,],]) 
res <- matrix(,nrow(DF),nrow(DF),dimnames=list(rownames(DF),rownames(DF))) 
res[lower.tri(res)] <- vals 

#  GeneA GeneB GeneC GeneD 
# GeneA NA NA NA NA 
# GeneB  9 NA NA NA 
# GeneC 24 15 NA NA 
# GeneD 26 13 28 NA 

Ваша матрица Comb является промежуточным результатом mat[pairs[1,],]*mat[pairs[2,],].

---

Весь расчет можно сделать внутри combn, поочередно:

vals <- combn(rownames(DF),2,FUN=function(x)sum(apply(DF[x,],2,prod))) 

Как @josilber отметил в комментарии ниже, это невероятно медленно, однако.

---

данные:

DF <- read.table(header=TRUE,text="Condition1 Condition2 Condition3 
GeneA   1   2    4 
GeneB   3   1    1 
GeneC   2   7    2 
GeneD   2   2    5")