嵌套for循环的data.table替代品？

Question

我有两个矩阵 diff.exp.protein.expr 和 lncRNA.expr 每个矩阵都有 64 列，但行数不同。

>dim(diff.exp.protein.expr)
[1] 14000 64

>dim(lncRNA.expr)
[1] 7600 64

我在两个单独的线性模型中使用这些矩阵作为输入，在其中我将 diff.exp.protein.expr 的每一行与 lncRNA.expr 的所有行进行比较（ =2*1.06 亿次测试）。最后，我想使用 anova 来比较这两个线性模型在统计上是否不同，我为此编写了一个函数，如下所示：

lm.anova <- function(diff.exp.protein.expr,lncRNA.expr,colData) 
{
    tempdff <- data.frame() #create an empty dataframe
    for(i in 1:nrow(diff.exp.protein.expr)) #traverse through 1st matrix
    {
        for(j in 1:nrow(lncRNA.expr)) #traverse through 2nd matrix
        {
            #model 1
            lm1 <- lm(diff.exp.protein.expr[i,]~colData$gender+colData$age+colData$treatment)
            #model 2 (has an additional interaction term at the end)
            lm2 <- lm(diff.exp.protein.expr[i,]~colData$gender+colData$age+colData$treatment+lncRNA.expr[j,]+colData$treatment*lncRNA.expr[j,])
            #get the summary of model 2
            lm.model <- summary(lm2)
            #get the rownames of ith row of matrix 1 and jth row of matrix 2
            #get the pvalue from the anova model
            #get the estimate and std. error for last two terms of model 2
            #add these values as a row to tempdff
            tempdff <- rbind(tempdff,data.frame(rownames(diff.exp.protein.expr)[i],rownames(lncRNA.expr)[j],
                                      anova(lm1,lm2)$"Pr(>F)"[2]),lm.model$coefficients[4,1:2],lm.model$coefficients[6,1:2])
        }
    }
    return(tempdff) #return results
}

lm.anova.res <- lm.anova(diff.exp.protein.expr,lncRNA.expr,colData) #call function

这就是进入函数的数据的样子：

>head(diff.exp.protein.expr)[,1:5] #log transformed values

                     C00060   C00079   C00135   C00150   C00154
ENSG00000000005.5  5.187977 6.323024 6.022986 5.376513 4.810042
ENSG00000000460.12 7.071433 7.302448 7.499133 7.441582 7.439453
ENSG00000000938.8  8.713285 8.584996 8.982816 9.787420 8.823927
ENSG00000001497.12 9.569952 9.658418 9.392670 9.394250 9.087214
ENSG00000001617.7  9.215362 9.417367 8.949943 8.172713 9.198314
ENSG00000001626.10 6.363638 6.038328 6.698733 5.202132 5.336239

>head(lncRNA.expr)[,1:5] #log transformed values
                     C00060   C00079   C00135   C00150   C00154
ENSG00000100181.17 7.477983 7.776463 7.950514 8.098205 7.278615
ENSG00000115934.11 4.104380 4.104380 4.104380 4.104380 4.104380
ENSG00000122043.6  4.104380 4.104380 4.104380 4.104380 4.104380
ENSG00000124915.6  4.104380 4.104380 4.104380 4.104380 4.104380
ENSG00000125514.5  4.104380 4.104380 4.104380 4.104380 4.104380
ENSG00000129816.5  4.104380 4.104380 4.104380 4.104380 4.104380

>head(colData)[,1:5] #sample information for each column of the two input matrices
sample_name status age gender      site treatment
     C00060     NF  48 Female Cleveland      case
     C00079     NF  26 Female Cleveland      case
     C00135     NF  55 Female Cleveland      case
     C00150     NF  61   Male Cleveland      ctrl
     C00154     NF  50   Male Cleveland      case
     C00176     NF  60 Female Cleveland      ctrl

当我几乎没有要执行的测试 (~50) 时，我编写了这段代码。现在，我已经弄清楚了，除了我不知道如何使我的代码高效，因为在这种必须执行 14000*7600 测试的情况下使用 for 循环没有任何意义。它需要永远运行。我想知道我可以真正加快这段代码的其他替代方案。任何建议将不胜感激。

更新 1：我稍微更新了我的线性模型。现在，anova(lm1,lm2)"Pr(>F)" 给出的值与模型 2 中的交互项的值不同。

更新 2：我更新了我的线性模型，使 lm1 嵌套在 lm2 中。

谢谢，

Answer 1

所以有几件事。

首先，调用data.frame(...) 非常昂贵，因此在每次迭代中调用它会大大减慢速度。另一方面，列表非常有效。因此，请尽量将所有内容都保留在列表中，直到最后。

其次，可能只是运行 2 * 1.06 亿次回归占用了大部分时间。

我会倾向于这样尝试：

## not tested...
df1 <- t(diff.exp.protein.expr)
df2 <- t(lncRNA.expr)

df  <- cbind(df1,df2,colData)
names <- expand.grid(x=colnames(df1),y=colnames(df2),stringsAsFactors=FALSE)
get.anova <- function(n){
  form.1 <- as.formula(paste0(n[1],"~gender+age+treatment+",n[2]))
  form.2 <- as.formula(paste0(n[1],"~gender+age+treatment+",n[2],"+treatment:",n[2]))
  lm1    <- lm(form.1,data=df)
  lm2    <- lm(form.2,data=df)
  coef <- summary(lm2)$coefficients
  list(n[1],n[2],anova(lm1,lm2)$"Pr(>F)"[2],coef[5,1],coef[5,2],coef[6,1],coef[6,2])
}
result <- do.call(rbind,apply(names,1,get.anova))
result <- data.frame(result)
colnames(result) <- c("protein","lncRNA","p.value","est.1","se.1","est.2","se.2")

这没有经过测试，因为您提供的数据集不够大：模型的误差

编辑（对 OP 的评论和数据提供的回应）。

使用您评论中提供的数据集，我们可以对原始代码（基于您更新的帖子）和新版本（基于上面的代码，更新以反映您的新模型公式）进行基准测试。在这两种情况下，我只使用每个数据集中的 10 行（100 种组合）。

f.original <- function() lm.anova(diff.exp.protein.expr.sub[1:10,],lncRNA.expr.sub[1:10,],colData)
f.new      <- function() do.call(rbind,apply(names,1,get.anova))
library(microbenchmark)
microbenchmark(f.original(), f.new(), times=10)
# Unit: seconds
#          expr      min       lq   median       uq      max neval
#  f.original() 2.622461 2.712527 2.824419 2.869193 2.914099    10
#       f.new() 2.049756 2.074909 2.144422 2.191546 2.224831    10

所以我们可以看到，返回列表而不是数据帧的新版本快了大约 25%。

使用Rprof 对这两种方法进行分析表明，原始版本在lm(...) 中花费了大约50% 的时间，而新版本在lm(...) 中花费了大约60%（更短的）时间。这是有道理的，并且表明您可以做的最好的事情是比新版本提高约 30%。换句话说：lm(...) 的调用是瓶颈：lm(...) 的 2 亿次调用只需要很长时间。

您可以考虑使用并行计算方法，例如使用foreach package，但在走这条路之前，IMO 您应该考虑其他策略来获得最终结果。

Answer 2

我的包裹MatrixEQTL似乎解决了你的问题。

http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/

http://cran.r-project.org/web/packages/MatrixEQTL/index.html

您需要使用modelLINEAR_CROSS 模型来测试交互项的重要性。

请随时询问有关包裹的任何问题。

Answer 3

这里的主要时间消费者不是嵌套循环，而是lm。以下是您可以优化的几件事（但另请参阅 Andrey Shabalin 的回答 - 他可能会为您的情况提供更简单的解决方案）。

您有两个 lm，简化形式：

        lm1 <- lm(Y~ X1 + X2 + X3 + X4)
        lm2 <- lm(Y~ X1 + X2 + X3 + X4 + X3:X4)

然后您将lm2 中的summary 进行比较，并将lm1 与anova 进行比较。因为您仅从 anova 中提取 p 值，但 p 值与 lm2 中交互项的 p 值相同，所以执行 anova 是多余的。因此lm1 也是多余的。然后使用rbind 增加数据框是浪费时间，所以我建议使用 list 并在每次迭代时添加一个新元素。因此循环中的代码可以简化为：

# lm1 <- ... # skip this 
#model 2 (has an additional interaction term at the end)
lm2 <- lm(diff.exp.protein.expr[i,]~colData$gender+colData$age+colData$treatment+lncRNA.expr[j,]+lncRNA.expr[j,]:colData$treatment)
#get the summary of model 2, * and coefficients from that summary *
lms <- coef(summary(lm2))
tempdff[[length(tempdff)+1]]  <- data.frame(rownames(diff.exp.protein.expr)[i],
     rownames(lncRNA.expr)[j], lms[6,4],lms[5,1:2],lms[6,1:2]) 

此外，您将结果作为 data.frame - 使用 list 也可以节省一些时间。

下一步可能是检查lm 和summary.lm 正在做什么。你不需要它做的所有事情，你只需要使用 b 和它们的标准错误，以及最后一行的 p 值。您也许可以跳过 lm 和 summary.lm 所做的一些计算......也许借助一些矩阵代数。