find individual peaks from an average spectrum

avgpeaks.R

library("MALDIquant")

data("fiedler2009subset")

## preprocessing
spectra <- transformIntensity(fiedler2009subset, sqrt)
spectra <- transformIntensity(spectra, savitzkyGolay)
spectra <- removeBaseline(spectra)
peaks <- detectPeaks(spectra)

wf <- determineWarpingFunctions(peaks)

wSpectra <- warpMassSpectra(spectra, wf)
wPeaks <- warpMassPeaks(peaks, wf)

avgSpectra <- mergeMassSpectra(wSpectra)
avgPeaks <- detectPeaks(avgSpectra)
## end of preprocessing

createPeakMatrix <- function(avgPeaks, spectra) {
  ## look for peaks (SNR == 0 to detect all possible maxima)
  peaks <- detectPeaks(spectra, SNR=0)
  peaks <- lapply(peaks, function(p) {
    ## find nearest index
    i <- MALDIquant:::.which.closest(mass(avgPeaks), mass(p))
    ## select only nearest peaks
    return(p[unique(i)])
  })
  return(intensityMatrix(binPeaks(peaks)))
}

createPeakMatrix2 <- function(avgPeaks, spectra) {
  peaks <- lapply(spectra, function(s) {
    ## find nearest index
    i <- MALDIquant:::.which.closest(mass(avgPeaks), mass(s))
    ## select only nearest mass and intensity
    s <- s[unique(i)]
    ## create MassPeaks object out of MassSpectrum object
    return(createMassPeaks(mass=mass(s),
                           intensity=intensity(s),
                           metaData=metaData(s)))
  })
  return(intensityMatrix(binPeaks(peaks)))
}

peakMatrix1 <- createPeakMatrix(avgPeaks, spectra)
peakMatrix2 <- createPeakMatrix2(avgPeaks, spectra)