代码之家 › 专栏 › 技术社区 › El Dude

R函数超慢

bioinformatics r

El Dude · 技术社区 · 9 年前

这是一个具有R的自对准矩阵的非常简单的生物信息学实现。它使用滑动窗口算子在字符串序列上循环两次 frag1 fra2 序列相同。下面的代码非常慢,不知道如何用标准的R语法加速。在python中,这将非常快,但在R中需要1分钟。我已经通过分配 i,j 和 j,i

sequence = 'MNLDIHCEQLSDARWTELLPLLQQYEVVRLDDCGLTEEHCKDIGSALRANPSLTELCLRTNELGDAGVHLVLQGLQSPTCKIQKLSLQNCSLTEAGCGVLPSTLRSLPTLRELHLSDNPLGDAGLRLLCEGLLDPQCHLEKLQLEYCRLTAASCEPLASVLRATRALKELTVSNNDIGEAGARVLGQGLADSACQLETLRLENCGLTPANCKDLCGIVASQASLRELDLGSNGLGDAGIAELCPGLLSPASRLKTLWLWECDITASGCRDL'

if(!exists('BLOSUM50')){
    library(Biostrings)
    data(BLOSUM50)
    #BLOSUM50['A','N']
  }

  windowSize<-24;
  matrixSize<-nchar(sequence) - windowSize;
  defaultValue = -10000000000;
  scoreMatrix <- matrix(defaultValue, nrow = matrixSize, ncol = matrixSize);

  for(i in 1:matrixSize){
    frag1 = substr(sequence,i,i+windowSize);

    for(j in 1:matrixSize){

      frag2 = substr(sequence,j,j+windowSize);

      totalScore = 0;

      if(scoreMatrix[i,j] == defaultValue){
        for(x in 1:windowSize){
          totalScore = totalScore + BLOSUM50[substr(frag1,x,x),substr(frag2,x,x)] / windowSize;
        }

        scoreMatrix[i,j] = totalScore;
        scoreMatrix[j,i] = totalScore;
      }

    }
  }
  return(scoreMatrix);

1 回复 | 直到 9 年前

Gregory Demin 9 年前

你在我不太新的笔记本电脑(2014年的联想瑜伽2,R3.4)上运行的原始代码只需17秒。经过不太繁重的优化后,此时间减少到2秒。我刚转换 sequence 在计算开始时矢量。之后,我在中更改了名称索引 BLOSUM50 通过数字索引进行索引。这导致了0.5秒的执行时间。请参阅下面的代码和基准测试。

fun = function(sequence){
     windowSize<-24
     matrixSize<-nchar(sequence) - windowSize
     defaultValue = -10000000000
     scoreMatrix <- matrix(defaultValue, nrow = matrixSize, ncol = matrixSize)

     for(i in 1:matrixSize){
         frag1 = substr(sequence,i,i+windowSize)

         for(j in 1:matrixSize){

             frag2 = substr(sequence,j,j+windowSize)

             totalScore = 0

             if(scoreMatrix[i,j] == defaultValue){
                 for(x in 1:windowSize){
                     totalScore = totalScore + BLOSUM50[substr(frag1,x,x),substr(frag2,x,x)] / windowSize
                 }

                 scoreMatrix[i,j] = totalScore
                 scoreMatrix[j,i] = totalScore
             }

         }
     }

     scoreMatrix
 }

 fun2 = function(sequence){
     windowSize<-24
     sequence = unlist(strsplit(sequence, split = ""))
     matrixSize<-length(sequence) - windowSize
     defaultValue = -10000000000
     scoreMatrix <- matrix(defaultValue, nrow = matrixSize, ncol = matrixSize)

     for(i in 1:matrixSize){
         frag1 = sequence[i:(i+windowSize)]

         for(j in 1:matrixSize){

             frag2 = sequence[j:(j+windowSize)]

             totalScore = 0

             if(scoreMatrix[i,j] == defaultValue){
                 for(x in 1:windowSize){
                     totalScore = totalScore + BLOSUM50[frag1[x],frag2[x]] / windowSize
                 }

                 scoreMatrix[i,j] = totalScore
                 scoreMatrix[j,i] = totalScore
             }

         }
     }

     scoreMatrix
 }

 fun3 = function(sequence){
     windowSize = 24
     sequence = unlist(strsplit(sequence, split = ""))
     matrixSize = length(sequence) - windowSize
     scoreMatrix = matrix(NA, nrow = matrixSize, ncol = matrixSize)
     sequence_index = match(sequence, colnames(BLOSUM50))
     for(i in seq_len(matrixSize)){
         frag1 = sequence_index[i:(i+windowSize - 1)]

         for(j in seq_len(matrixSize)){

             frag2 = sequence_index[j:(j+windowSize - 1)]

             if(is.na(scoreMatrix[i,j])){
                 totalScore = sum(BLOSUM50[(frag2 - 1)*NROW(BLOSUM50) + frag1])/windowSize
                 scoreMatrix[i,j] = totalScore
                 scoreMatrix[j,i] = totalScore
             }

         }
     }

     scoreMatrix
 }


 if(!exists('BLOSUM50')){
     library(Biostrings)
     data(BLOSUM50)
     #BLOSUM50['A','N']
 }

 sequence = 'MNLDIHCEQLSDARWTELLPLLQQYEVVRLDDCGLTEEHCKDIGSALRANPSLTELCLRTNELGDAGVHLVLQGLQSPTCKIQKLSLQNCSLTEAGCGVLPSTLRSLPTLRELHLSDNPLGDAGLRLLCEGLLDPQCHLEKLQLEYCRLTAASCEPLASVLRATRALKELTVSNNDIGEAGARVLGQGLADSACQLETLRLENCGLTPANCKDLCGIVASQASLRELDLGSNGLGDAGIAELCPGLLSPASRLKTLWLWECDITASGCRDL'
 library(microbenchmark)
 microbenchmark(
     original = fun(sequence),
     fun2 = fun2(sequence),
     fun3 = fun3(sequence),
     times = 5
 )

 # Unit: milliseconds
 # expr          min         lq       mean     median         uq        max neval
 # original 16395.2108 16660.3295 17533.8563 16755.9680 17594.3596 20263.4137     5
 # fun2      1992.7731  2010.4031  2027.7953  2015.9592  2034.9022  2084.9390     5
 # fun3       472.0641   481.9267   496.2656   498.3259   506.6357   522.3755     5