R : 맞춤형 거리 함수 및 다중 기준을 기반으로 레코드를 빠르게 일치시킵니다.

Question

R의 일부 기능을 만들어 맞춤 스펙트럼을 기반으로 한 그러한 스펙트럼 라이브러리에 화학 질량 스펙트럼 (정수 질량 및 강도가있는 두 개의 열이있는 행렬)을 대조했습니다. 유사성 기능 및 화합물의 소위 고정 지수 (즉, 용리 시간)의 일치 (여기에 예를 들어, http://webbook.nist.gov/cgi/cbook.cgi?ID=C630035&Mask=200 참조). 이를 위해 각 레코드의 목록 요소 "RI"를 라이브러리의 요소 "RI"와 비교해야하며 절대 편차가 지정된 허용 오차보다 작 으면 내 레코드에 최상의 스펙트럼 라이브러리 일치를 추가해야합니다. 아래 코드는 제가 작성한 코드입니다. 그러나 문제는 그것이 저의 목적에 너무 느리다는 것입니다 (저는 보통 1000 개의 샘플 스펙트럼과 200,000 개의 라이브러리 스펙트럼을가집니다). 나는 병렬 처리를 시도했으나 많은 도움이되지는 않았다. 아마도 아래 코드를 어떻게 효율적으로 만들 수 있는지에 대한 생각. 더 많은 벡터화를 사용하거나 인라인 C 코드를 사용하거나 다른 R 트릭을 사용합니까? 나는이 점에서 일반적인 조언을 알고 있지만,이 경우에는 쉽게 구현하는 법을 보지 못한다. (그리고 나는 아직 불행히도 C에 익숙하지 않다.) ... 어떤 생각이나 충고? 아 그래, sfLapply을 사용할 때 진행률 표시 줄을 어떻게 추가 할 수 있습니까? 스펙트럼 유사성 함수에서 merge 단계를 피하기 위해, 또는 가장 큰/최대 스펙트럼을 고려할 때 스펙트럼을 고려하는 것과 같은 추가 기준을 사용하여, 아마도 제로를 많이 가진 것처럼, 하나의 큰 (제수가 많은 것처럼) 쿼리 스펙트럼에서 가장 강렬한 피크는 라이브러리 스펙트럼과 동일한 질량을 갖습니다 (또는 라이브러리 스펙트럼의 5 개의 최대 피크를 말합니다). 어쨌든,이 작업을 빠르게하는 방법에 대한 생각은 많이 감사하겠습니다!

EDIT : 아직 남아있는 쿼리 중 하나는 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 

Answer 1

모든 코드를 자세하게 살펴 보지 않고도 모든 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 

을 라이브러리가 하나의 큰 매트릭스 (질량 X 스펙트럼) 인 경우에, 당신은 선택의 스펙트럼에 대한 라이브러리를 부분 집합 및 샘플 및 전체 부분 사이의 거리를 계산할 수 있습니다 이와 같이 한 번에 선택된 라이브러리 :

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))의 행렬

. 이 방법으로 매트릭스를 하위 집합으로 지정하고 (빠르게) 일련의 매트릭스 연산을 수행합니다 (빠름).