2014-11-24 67 views
0

我正在使用函數從多個HWE卡方檢驗中獲取p值。我正在遍歷一個名爲geno.data的大矩陣,它是(313行x 355232列)來執行此操作。我基本上是逐行循環兩列矩陣。它運行非常緩慢。我怎樣才能讓它更快?由於如何讓我的循環在R中運行得更快?

library(genetics) 
geno.data<-matrix(c("a","c"), nrow=313,ncol=355232) 
Num_of_SNPs<-ncol(geno.data) /2 
alleles<- vector(length = nrow(geno.data)) 
HWE_pvalues<-vector(length = Num_of_SNPs) 
j<- 1 

for (count in 1:Num_of_SNPs){ 
    for (i in 1:nrow(geno.data)){ 
     alleles[i]<- levels(genotype(paste(geno.data[i,c(2*j -1, 2*j)], collapse = "/"))) 
    } 
    g2 <- genotype(alleles) 
    HWE_pvalues[count]<-HWE.chisq(g2)[3] 
    j = j + 2 
} 
+2

請參閱http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – csgillespie 2014-11-24 15:47:37

+0

所以你正在做'選擇(355232,2)'chisq測試?你碰巧認識fortran嗎? – rawr 2014-11-24 15:53:30

+0

@rawr我不知道fortran。我正在使用R包中的函數來進行卡方檢驗。這是特定於我的問題。 – cooldood3490 2014-11-24 16:03:42

回答

3

首先,注意張貼的代碼將導致索引出界外的錯誤,因爲畢竟主循環迭代Num_of_SNPsj值將ncol(geno.data)-1和你所訪問的列2*j-12*j。我假設你想要刪除列2*count-12*countj

矢量化對於編寫快速R代碼非常重要。在你的代碼中,你調用paste函數313次,每次傳遞長度爲1的向量。一旦傳遞了長度爲313的向量,R中的paste就會快得多。下面是main for循環的原始矢量化內部:

# Original 
get.pval1 <- function(count) { 
    for (i in 1:nrow(geno.data)){ 
    alleles[i]<- levels(genotype(paste(geno.data[i,c(2*count -1, 2*count)], collapse = "/"))) 
    } 
    g2 <- genotype(alleles) 
    HWE.chisq(g2)[3] 
} 

# Vectorized 
get.pval2 <- function(count) { 
    g2 <- genotype(paste0(geno.data[,2*count-1], "/", geno.data[,2*count])) 
    HWE.chisq(g2)[3] 
} 

我們得到關於從量化20倍速度提升:

library(microbenchmark) 
all.equal(get.pval1(1), get.pval2(1)) 
# [1] TRUE 
microbenchmark(get.pval1(1), get.pval2(1)) 
# Unit: milliseconds 
#   expr  min  lq  mean median  uq  max neval 
# get.pval1(1) 299.24079 304.37386 323.28321 307.78947 313.97311 482.32384 100 
# get.pval2(1) 14.23288 14.64717 15.80856 15.11013 16.38012 36.04724 100 

隨着量化代碼,你的代碼應該在大約177616 * 0.01580856 =2807.853秒,或約45分鐘完成(相比於16小時爲原始代碼)。如果對你來說這還不夠快,那麼我建議你看一下R中的parallel包。mcmapply應該給你一個很好的加速,因爲外部for循環的每次迭代都是獨立的。