R序列比对函数运行时间过长

Question

所以，我对R相当陌生，我有一个运行时问题。我使用"Biostrings"软件包（biocLite）编写了以下嵌套while循环，以便在两个物种的比对得分>90%的情况下连接它们的蛋白质序列。

基本上，我输入两个蛋白质基因组，将SeqData1中的每个氨基酸序列与SeqData2中的每个核酸序列进行比较，计算比对得分，如果得分>90%，我将匹配的蛋白质名称列表和SeqData2蛋白质的序列连接起来。

该功能确实像它应该的那样工作，唯一的问题是，根据它必须扫描的蛋白质数量，我预计整个功能的运行时间约为1.4个月。T

有人对如何以指数级速度加快此函数的运行时间有什么建议吗？

谢谢！

R代码：

SeqScore <- function()
{
source("http://bioconductor.org/biocLite.R")
biocLite()
require("Biostrings")
data(BLOSUM100)

SeqData1 <- readDNAStringSet("SeqData1.fasta")
proNum1 = 84390     # number of proteins in Seq1
SeqData2 <- readDNAStringSet("SeqData2.fasta")
proNum2 = 15194     # number of proteins in Seq

#Create empty list to fill with percent scores and matching sequences:
DList=NULL
QueSeqList = NULL
TotList = NULL
#initiating the counters:
i=1
j=1
c=0
#Perform alignment and generate percent identity scores:
while(i<=proNum1)
  {
  while(j<=proNum2)
    {
  SeqAlign <- pairwiseAlignment(SeqData1[i], SeqData2[j], substitutionMatrix=BLOSUM100, gapOpening=0, gapExtension=-5)
  PercAlign <- pid(SeqAlign)
  if(PercAlign>=90)
  {
    DList = c(DList, names(SeqData1[i]), names(SeqData2[j]))
    QueSeqList = c(QueSeqList, toString(SeqData2[j]))
    c=c+1
  }
  else{c=c+1; print(c)}
  j=j+1
  }
  i=i+1
  j=1 #to reset the inner while loop
  }
unlist<-t(sapply(DList, unlist));
outputMatrix<-cbind(DList,QueSeqList)
outputMatrix<-as.matrix(outputMatrix, ncol=3)
write.csv(outputMatrix, "outputMatrix.csv")
}

Answer 1

从帮助页面中，我认为当pairwiseAlignment的第一个参数是任何长度的DNAStringSet时，它都能工作，所以外循环可以用代替

SeqAlign <- pairwiseAlignment(SeqData1, SeqData2[j],
    substitutionMatrix=BLOSUM100, gapOpening=0, gapExtension=-5)

我不会使用while循环，而是按照的路线瞄准lapply

results <- lapply(SeqData2, function(elt_j) {
        SeqAlign <- pairwiseAlignment(SeqData1, elt_j,
            substitutionMatrix=BLOSUM100, gapOpening=0, gapExtension=-5)
        pid(SeqAlign)
})
data.frame(toString(SeqData1), rep(toString(SeqData2), each=length(SeqData1)),
    unlist(pid, use.names=FALSE))

这将您从R推断的第二个循环中拯救出来，并且还建议了一个简单的策略，如果计算仍然很慢，可以在非Windows上并行化：

library(parallel)
options(mc.cores=detectCores())
results <- mclapply(seq_along(... ## as before

与其单独键入proNum1 = length(SeqData1)，不如询问数据的长度。此外，我希望您不要每次运行脚本时都使用biocLite()；只需一次即可将软件包安装到任何给定的R安装中。

你会在Bioconductor邮件列表上得到一个更权威的答案——无需订阅。