0

R でいくつかの関数を作成して、化学質量スペクトル (整数の質量と強度を持つ 2 つの列を持つマトリックス) を、カスタム スペクトル類似度関数と化合物のいわゆる保持インデックスのマッチングに基づいて、そのようなスペクトルのライブラリに一致させました (つまり、溶出時間) (例については、http://webbook.nist.gov/cgi/cbook.cgi? ID=C630035&Mask=200 を参照)）。そのためには、各レコードのリスト要素「RI」をライブラリ内の要素と比較する必要があり、絶対偏差が所定の許容範囲よりも小さい場合、スペクトル ライブラリの最適な一致をレコードに追加する必要があります。以下は、これを行うために私が書いたコードですが、問題は、私の目的には遅すぎることです (通常、約 1000 のサンプル スペクトルと 200 000 のライブラリ スペクトルがあります)。私はそれを並列化しようとしましたが、それもあまり役に立たないようです。以下のコードをより効率的にする方法について、おそらく何か考えはありますか? たとえば、より多くのベクトル化を使用したり、インライン C コードを使用したり、その他の R トリックを使用したりしますか? この点に関する一般的なアドバイスは知っていますが、この場合に簡単に実装する方法がわかりません (残念ながら、私はまだ C に精通していません)...何か考えやアドバイスはありますか? そうそう、sfLapply? おそらく、スペクトル類似度関数のステップを回避するために、最初にスペクトルを1つの大きな（ゼロがたくさんあるため、スパース）マトリックスに入れるmergeか、最大/最大の場合にのみスペクトルを考慮するなどの追加の基準を使用するのに役立ちますかクエリ スペクトルの強いピークは、ライブラリ スペクトルと同じ質量を持っていますか (またはライブラリ スペクトルの 5 つの最大ピークのセットに含まれていますか)? とにかく、このタスクをスピードアップする方法についての考えは大歓迎です!

編集: 私がまだ持っている残りのクエリの 1 つは、関数 addbestlibmatches1 でサンプル レコード recs の完全なコピーを作成することを回避する方法ですが、ライブラリの一致がある場所のレコードのみを変更する方法です。また、保持インデックスの一致があるライブラリ レコードの選択を渡すことは、おそらく効率的ではありません (関数 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
4

1 に答える 1

2

すべてのコードを詳細に調べなくても、まだすべて C++ に移行しなくても、SpecSim 関数には改善の余地があると思います。マトリックスをdata.framesに強制するマージを使用しています。それは常にパフォーマンスに悪影響を及ぼします。コード時間のほとんどは、おそらく merge() とサブセット化に費やされています。data.frame サブセット化は遅く、行列またはベクトルは高速です。

SpecSim2 <- function (ms1,ms2,log=F) {
  i <- unique(c(ms1[,1], ms2[,1]))
  y <- x <- numeric(length(i))
  x[match(ms1[,1], i)] <- ms1[, 2]
  y[match(ms2[,1], i)] <- ms2[, 2]
  x <- 100 * x / max(x)
  y <- 100 * y / max(y)
  if (log){
    x <- log2(x + 1)
    y <- log2(y + 1)
  }
  similarity.score <- x %*% y / (sqrt(sum(x^2)) * sqrt(sum(y^2)))
  if (is.na(similarity.score)){
    return(-1)
  } else {
    return(similarity.score)
  }
}

それぞれ、書き直しとオリジナルのタイミングを次に示します。

> system.time(addbestlibmatches1(s,lib))
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
   user  system elapsed 
   4.16    0.00    4.17

> system.time(addbestlibmatches1(s,lib))
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
   user  system elapsed 
  34.25    0.02   34.34

必要な速度が得られない場合がありますが、8 倍以上の改善が得られます...

addbestlibmatches() の処理方法については、マトリックスの問題として再考する必要があると思います。リストを使用してライブラリを保持するのではなく、ライブラリとサンプルの両方の RI 値にベクトルを使用します。次に、次のように初期画面を実行できます。

selected <- outer(sRI, libRI, FUN = '-') < 10

ライブラリが 1 つの大きなマトリックス (質量 x スペクトル) である場合、選択したスペクトルのライブラリをサブセット化し、サンプルと 1 回のパスで選択されたライブラリの全体との間の距離を次のように計算できます。

SpecSimMat <- function(x, lib, log = F){
  stopifnot(require(matrixStats))
  x <- 100 * x / max(x)
  lib <- sweep(lib, 2, colMaxs(lib))
  x %*% lib / (sqrt(sum(x^2)) * sqrt(colSums(lib^2)))

}

ここで、x はサンプル、lib は選択したスペクトル (質量 x スペクトル) の行列です。このようにして、行列をサブセット化し (高速)、一連の行列演算を実行します (高速)。

于 2014-07-07T17:50:25.153 に答える