R: быстрое сопоставление записей на основе настраиваемой функции расстояния и нескольких критериев

Я сделал некоторые функции в R для сопоставления химических масс-спектров (матрица с двумя столбцами с целыми массами и интенсивностями) с библиотеками таких спектров на основе настраиваемой функции спектрального сходства и сопоставления так называемого индекса удерживания соединений (т. Е. время элюирования) (см. здесь пример, http://webbook.nist.gov/cgi/cbook.cgi?ID=C630035&Mask=200). Для этого элемент списка «RI» каждой записи должен быть сравнен с элементом списка в библиотеке, и когда абсолютное отклонение меньше заданного допуска, он должен добавить лучшее соответствие спектральной библиотеки к моим записям. Ниже приведен код, который я написал для этого, но проблема в том, что он слишком медленный для моих целей (у меня обычно около 1000 образцов спектров и 200 000 библиотечных спектров). Я попытался распараллелить его, но это тоже не помогло. Возможно, есть какие-нибудь мысли о том, как я могу сделать приведенный ниже код более эффективным, например. используя больше векторизации, или используя встроенный код C, или некоторые другие трюки R? Я знаю общий совет по этому поводу, но не совсем понимаю, как его легко реализовать в этом случае (а я, к сожалению, еще не так хорошо разбираюсь в C) ... Есть мысли или советы? Ах да, а как добавить индикатор выполнения при использовании sfLapply? Помогло бы, возможно, сначала поместить мои спектры в одну большую (разреженную, так как много нулей) матрицу, чтобы избежать шага merge в функции спектрального подобия, или использовать дополнительные критерии, такие как рассмотрение спектров только тогда, когда наибольшее / самый интенсивный пик в спектре запроса имеет ту же массу, что и спектр библиотеки (или содержится в наборе, скажем, 5 самых больших пиков в спектре библиотеки)? В любом случае, любые мысли о том, как ускорить эту задачу, были бы очень признательны!

РЕДАКТИРОВАТЬ: один остаточный запрос, который у меня все еще есть, заключается в том, как я мог бы избежать создания полной копии образцов записей recs в функции addbestlibmatches1, а вместо этого изменить только те записи, для которых есть соответствие библиотеки? Кроме того, передача выборки записей библиотеки, для которых есть соответствие индекса хранения, вероятно, неэффективна (в функции addbestlibmatch). Есть мысли, как я могу этого избежать?

# EXAMPLE DATA

rec1=list(RI=1100,spectrum=as.matrix(cbind(mz=c(57,43,71,41,85,56,55,70,42,84,98,99,39,69,58,113,156),int=c(999,684,396,281,249,173,122,107,94,73,51,48,47,47,37,33,32))))
randrec=function() list(RI=runif(1,1000,1200),spectrum=as.matrix(cbind(mz=seq(30,600,1),int=round(runif(600-30+1,0,999)))))

# spectral library
libsize=2000 # my real lib has 200 000 recs
lib=lapply(1:libsize,function(i) randrec())
lib=append(list(rec1),lib) 

# sample spectra
ssize=100 # I usually have around 1000
s=lapply(1:ssize,function(i) randrec())
s=append(s,list(rec1)) # we add the first library record as the last sample spectrum, so this should match



# SPECTRAL SIMILARITY FUNCTION

SpecSim=function (ms1,ms2,log=F) { 
  alignment = merge(ms1,ms2,by=1,all=T)
  alignment[is.na(alignment)]=0
  if (nrow(alignment)!=0) {
    alignment[,2]=100*alignment[,2]/max(alignment[,2]) # normalize base peak intensities to 100
    alignment[,3]=100*alignment[,3]/max(alignment[,3])
    if (log==T) {alignment[,2]=log2(alignment[,2]+1);alignment[,3]=log2(alignment[,3]+1)} # work on log2 intensity scale if requested
    u = alignment[,2]; v = alignment[,3]
    similarity_score = as.vector((u %*% v) / (sqrt(sum(u^2)) * sqrt(sum(v^2))))
    similarity_score[is.na(similarity_score)]=-1
    return(similarity_score)} else return(-1) }



# FUNCTION TO CALCULATE SIMILARITY VECTOR OF SPECTRUM WITH LIBRARY SPECTRA

SpecSimLib=function(spec,lib,log=F) {
  sapply(1:length(lib), function(i) SpecSim(spec,lib[[i]]$spectrum,log=log)) }



# FUNCTION TO ADD BEST MATCH OF SAMPLE RECORD rec IN SPECTRAL LIBRARY lib TO ORIGINAL RECORD
# we only compare spectra when list element RI in the sample record is within tol of that in the library
# when there is a spectral match > specsimcut within a RI devation less than tol,
# we add the record nrs in the library with the best spectral matches, the spectral similarity and the RI deviation to recs

addbestlibmatch=function(rec,lib,xvar="RI",tol=10,log=F,specsimcut=0.8) {
    nohit=list(record=-1,specmatch=NA,RIdev=NA)
    selected=abs(sapply(lib, "[[", xvar)-rec[[xvar]])<tol
    if (sum(selected)!=0) {
    specsims=SpecSimLib(rec$spectrum,lib[selected],log) # HOW CAN I AVOID PASSING THE RIGHT LIBRARY SUBSET EACH TIME?
    maxspecsim=max(specsims)
    if (maxspecsim>specsimcut) {besthsel=which(specsims==maxspecsim)[[1]] # nr of best hit among selected elements, in case of ties we just take the 1st hit
                                idbesth=which(selected)[[besthsel]] # record nr of best hit in library lib
                                return(modifyList(rec,list(record=idbesth,specsim=specsims[[besthsel]],RIdev=rec[[xvar]]-lib[[idbesth]][[xvar]])))}
                                else {return(rec)} } else {return(rec)}
}



# FUNCTION TO ADD BEST LIBRARY MATCHES TO RECORDS RECS

library(pbapply)
addbestlibmatches1=function(recs,lib,xvar="RI",tol=10,log=F,specsimcut=0.8) {
  pblapply(1:length(recs), function(i) addbestlibmatch(recs[[i]],lib,xvar,tol,log,specsimcut))
}

# PARALLELIZED VERSION
library(snowfall)
addbestlibmatches2=function(recs,lib,xvar="RI",tol=10,log=F,specsimcut=0.8,cores=4) {
  sfInit(parallel=TRUE,cpus=cores,type="SOCK")
  sfExportAll()
  sfLapply(1:length(recs), function(i) addbestlibmatch(recs[[i]],lib,xvar,tol,log,specsimcut))
  sfStop() 
}



# TEST TIMINGS

system.time(addbestlibmatches1(s,lib))
#|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#user  system elapsed 
#83.60    0.06   83.82 

system.time(addbestlibmatches2(s,lib))
#user  system elapsed  - a bit better, but not much
#2.59    0.74   42.37