优化数据框中每行t检验速度

Question

我正在尝试提高数据框中t检验循环的速度。

我有一个大型数据框（约15000行和205列）。每列都是一个单元格，每行都是一个基因。我可以根据不同参考表中提供的标识将列分成两组。

以下是我编写的循环：

for (i in 1:nrow(EC)){
  ttest_result[i,2] <- rowMeans(EC)[i]
  ttest_result[i,3] <- rowMeans(CP)[i]
  ttest_result[i,4] <- (ttest_result[i,2] - ttest_result[i,3])
  ttest_result[i,5] <- (ttest_result[i,2] / ttest_result[i,3])
  ttest_result[i,6]<- t.test(EC[i,],CP[i,], var.equal = TRUE)$p.value
  pb$tick()
}

这个循环要求我根据列标识将原始数据框拆分为2个数据框。然而，这个循环需要超过45分钟才能完成。

我想知道你们是否有建议，可以让我做出不同的选择？我如何潜在地使用apply函数来提高速度？

非常感谢！

- Son nguyen

---

如果所有的列都是数字，将其转换为矩阵。在数据框上进行行操作非常缓慢。如果你转置矩阵并在列上进行操作而不是行上，效果会更好一些。 - lmo

Answer 1

limma 是解决方案。

> #library(BiocInstaller)
> #biocLite("limma")
> 
> #create a dataset
> library(limma)
> data <- matrix(rnorm(15000*205),15000,205)
> dim(data)
[1] 15000   205
> rownames(data) <- paste("Gene",1:15000)
> str(data)
 num [1:15000, 1:205] -0.55603 -0.45478 -1.76432 0.05198 0.00844 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:15000] "Gene 1" "Gene 2" "Gene 3" "Gene 4" ...
  ..$ : NULL
> # save the grouping in a factor
> f<-sample(c("ctrl","treat"),size = 205,replace = T)
> 
> # perform the comparison gene per grouping
> t<-Sys.time()
> design <- model.matrix(~0+f)
> colnames(design) <- c("ctrl","treat")
> fit2 <- lmFit(data,design)
> contrast.matrix <- makeContrasts("treat-ctrl", levels=design)
> fit2 <- contrasts.fit(fit2, contrast.matrix)
> fit2 <- eBayes(fit2)
> top_table<-topTable(fit2, adjust="BH",coef=1,number =15000)
> dim(top_table)
[1] 15000     6
> head(top_table)
                logFC     AveExpr         t      P.Value adj.P.Val          B
Gene 12434 -0.6238005  0.07603032 -4.454575 8.459040e-06 0.1268856 -0.5006572
Gene 11609 -0.5827713  0.11178709 -4.156804 3.242956e-05 0.2174629 -1.0670677
Gene 5924   0.5729590 -0.02980352  4.089151 4.349258e-05 0.2174629 -1.1903102
Gene 10460 -0.5274251 -0.07930193 -3.770822 1.632559e-04 0.5747294 -1.7431451
Gene 5950   0.5216678 -0.03304759  3.730682 1.915765e-04 0.5747294 -1.8096840
Gene 14518  0.5053476 -0.05750282  3.612195 3.044821e-04 0.6298752 -2.0019558
> Sys.time()-t
Time difference of 0.3026412 secs

请注意，它还会打印出调整后的p值（您可以选择方法），这是筛选感兴趣基因的标准之一。

- edoardo pedrini

---

你似乎是一个计算生物学的人。 - PesKchan

Answer 2

nc <- 205
nr <- 15000

set.seed(30)
EC <- matrix(rnorm(nr * nc), nr, nc)
CP <- matrix(rnorm(nr * nc), nr, nc)

在循环之前计算行的均值和方差（这是您犯下的最大错误，将此操作放入循环中）：

meansEC <- rowMeans(EC)
meansCP <- rowMeans(CP)

require(matrixStats)
varsEC <- rowVars(EC)
varsCP <- rowVars(CP)

使用预先计算的行均值和行方差，我们可以在不使用t.test函数的情况下更快地计算p.value（您可以查看 .test 代码以提取所需部分）：

t.test.p.value <- function(i, j, nx, ny){
  mu <- 0
  mx <- meansEC[i]
  vx <- varsEC[i]
  my <- meansCP[j]
  vy <- varsCP[j]
  df <- nx + ny - 2
  v <- 0
  if (nx > 1) v <- v + (nx - 1)*vx
  if (ny > 1) v <- v + (ny - 1)*vy
  v <- v/df
  stderr <- sqrt(v*(1/nx + 1/ny))
  tstat <- (mx - my - mu)/stderr
  pval <- 2 * pt(-abs(tstat), df)
  pval
}

# create resulting matrix
ttest_result <- matrix(NA, nrow(EC), 7)
t <- Sys.time()

nx <- ncol(EC)
ny <- ncol(CP)

让我们使用bout t.test.p.value函数和默认的t.test进行计算:

for (i in 1:nrow(EC)) {
  ttest_result[i, 2] <- meansEC[i]
  ttest_result[i, 3] <- meansCP[i]
  ttest_result[i, 4] <- (meansEC[i] - meansCP[i])
  ttest_result[i, 5] <- (meansEC[i] / meansCP[i])
  ttest_result[i, 6] <- t.test(EC[i, ], CP[i, ], var.equal = TRUE)$p.value
  ttest_result[i, 7] <- t.test.p.value(i, i, nx, ny)
}
t <- Sys.time() - t
t

ttest_result[1:5, 5:7]
#            [,1]       [,2]       [,3]
# [1,] -0.3248217 0.35084307 0.35084307
# [2,] -2.3455622 0.11621785 0.11621785
# [3,] -2.1586716 0.01294876 0.01294876
# [4,]  1.1556035 0.98065576 0.98065576
# [5,]  1.9875296 0.92340948 0.92340948
all.equal(ttest_result[,6], ttest_result[, 7])
# [1] TRUE

我是一位有用的助手，可以翻译文本。

我们看到结果是相等的。

仅使用我的方法计算这些数据的时间：

t <- Sys.time()
meansEC <- rowMeans(EC)
meansCP <- rowMeans(CP)
require(matrixStats)
varsEC <- rowVars(EC)
varsCP <- rowVars(CP)
ttest_result <- matrix(NA, nrow(EC), 7)
for (i in 1:nrow(EC)) {
  ttest_result[i, 2] <- meansEC[i]
  ttest_result[i, 3] <- meansCP[i]
  ttest_result[i, 4] <- (meansEC[i] - meansCP[i])
  ttest_result[i, 5] <- (meansEC[i] / meansCP[i])
  ttest_result[i, 6] <- t.test.p.value(i, i, nx, ny)
}
t <- Sys.time() - t
t #Time difference of 0.145169 secs

- minem