R：提高数据帧中每一行的t检验速度

我正在尝试提高数据帧的t检验循环的速度。我有一个大的数据框（约15000行和205列）。每列是一个细胞，每一行都是一个基因。根据不同参考表中提供的身份，我可以将这些列分组为2个组。R：提高数据帧中每一行的t检验速度

这是我写的循环：

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分钟才能完成。

我想知道你们是否都对我可以做不同的事情有什么建议？我如何使用应用函数来提高速度？

非常感谢！

来源

2017-10-07 Son nguyen

如果所有列是数字，将其转换为一个矩阵。数据帧的行操作非常慢。如果转置矩阵并在列上而不是在行上操作，则会更好一些。 – lmo

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值（您可以选择方法），这是筛选感兴趣基因的标准之一。

来源

2017-10-07 18:02:47

nc <- 205 
nr <- 15000 

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

计算行手段和VAR循环之前（这是你最大的错误，把这种操作成环）：

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

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

使用预先计算行手段和行瓦尔我们可以计算页。值得更快，而不t.test功能（你可以看看t.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)

让我们从计算回合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

来源

2017-10-09 07:56:51 minem

R：提高数据帧中每一行的t检验速度

回答

相关问题