R: amélioration de la vitesse du test t pour chaque ligne d'une trame de données

Question

J'essaie d'améliorer la vitesse d'une boucle de test t pour une trame de données.

J'ai un grand cadre de données (~ 15000 lignes et 205 colonnes). Chaque colonne est une cellule et chaque ligne est un gène. Je peux regrouper les colonnes en 2 groupes en fonction de leur identité fournie dans une autre table de référence.

Voici la boucle que je l'ai écrit:

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() 
} 

Cette boucle me oblige à diviser la trame de données d'origine en 2 trames de données basées sur les identités des colonnes. Cependant, cette boucle prend plus de 45 minutes à compléter.

Je me demande si vous avez tous des suggestions sur ce que je peux faire différemment? Comment puis-je utiliser des fonctions d'application pour améliorer la vitesse?

Merci beaucoup!

Answer 1

limma est la solution.

> #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 

avis qu'il imprime également le p-valeur ajustée (vous pouvez choisir la méthode), qui est l'un des critères pour filtrer les gènes d'intérêt.

Answer 2

nc <- 205 
nr <- 15000 

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

moyens de ligne Calculer et vars avant la boucle (ce fut votre plus grande erreur, de mettre cette opération en boucle):

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

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

En utilisant des moyens de ligne précalculées et la ligne vars nous pouvons calculer p. valeur beaucoup plus rapide sans t.test fonction (vous pouvez regarder t.test code pour extraire les pièces dont vous avez besoin):

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) 

permet de calculer avec bout t.test.p.value fonction et par défaut 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 

On voit que les résultats sont égaux

pour ces données Timing en utilisant seulement mon approche:

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