add PDM

R/Peak.R

@@ -54,3 +54,78 @@ integration <- function(x,yf){
   integral <- 0.5*sum((x[2:n] - x[1:(n-1)]) * (yf[2:n] + yf[1:(n-1)]))
   return(integral)
 }
+
+PDM <- function(pic, min_snr, level, pval, iter){
+  library(matrix)
+  library(IRanges)
+  vec <- pic[,2]
+  rts <- pic[,1]
+  mzs <- pic[,3]
+  peaks <- peak_detection(vec, min_snr, level)
+
+  vec.smooth <- WhittakerSmooth(vec, 2)
+  helf.height <- vec.smooth[peaks$peakIndex]*0.5
+  split.point <- sapply(1:(length(peaks$peakIndex)-1), function(s){
+    peaks$peakIndex[s]+which.min(vec.smooth[peaks$peakIndex[s]:peaks$peakIndex[s+1]])
+  })
+  split.point <- unique(c(1, split.point, length(vec)))
+
+  peak.ranges <- lapply(1:length(peaks$peakIndex),function(s){
+    split.point[s]+which(vec.smooth[split.point[s]:split.point[s+1]]>helf.height[s])-1
+  })
+
+  peak.mzs <- lapply(peak.ranges, function(peak.range){
+    mzs[peak.range]
+  })
+
+  pvals <- sapply(1:(length(peak.mzs)-1),function(s){
+    t.test(peak.mzs[[s]], peak.mzs[[s+1]])$p.value
+  })
+
+  TP <- rep(TRUE, length(peaks$peakIndex))
+  for (i in 1:length(pvals)){
+    if (pvals[i] > pval){
+      if (peaks$signals[i] > peaks$signals[i+1]){
+        TP[i+1] <- FALSE
+      } else {
+        TP[i] <- FALSE
+      }
+    }
+  }
+  peaks$peakIndex <- peaks$peakIndex[TP]
+  peaks$peakScale <- peaks$peakScale[TP]
+  peaks$snr <- peaks$snr[TP]
+  peaks$signals <- peaks$signals[TP]
+
+  res <- .PICfit(peaks, pic, iter)
+  plot.resolve(pic, res)
+  return(res)
+}
+
+WhittakerSmooth <- function(y,lambda){
+  M <- length(y)
+  E <- sparseMatrix(i=1:M,j=1:M,x=1)
+  D <- Matrix::diff(E)
+  C <- chol(E+lambda*Matrix::t(D)%*%D)
+  z <- solve(C,solve(t(C),y))
+  return(as.numeric(z))
+}
+
+plot.resolve <- function(pic, res){
+  library(plotly)
+  rts <- pic[,1]
+  raw.vec <- pic[,2]
+  fit.pics <- do.call(rbind, res$fitpics)
+  sum.vec <- colSums(fit.pics)
+
+  p <- plot_ly(x = rts, y = raw.vec, type = 'scatter', mode = 'lines', name = 'raw') %>% 
+    add_trace(x = rts, y = sum.vec, name = 'fitted') %>%
+    layout(xaxis = list(title = 'Retention Time (s)'),
+           yaxis = list (title = 'Intensity'))
+
+  for (i in 1:nrow(fit.pics)) {
+    name <- paste('peak ', i)
+    p <- add_trace(p, x = rts, y = fit.pics[i,],  line = list(dash = 'dash'), name = name)
+  }
+  p
+}