inst/doc/IPPD.R

In IPPD: Isotopic peak pattern deconvolution for Protein Mass Spectrometry by template matching

### R code from vignette source 'IPPD.Rnw'

###################################################
### code chunk number 1: setup
###################################################

options(prompt = "R> ", continue = " ", warn = -1, width = 90)                                       



###################################################
### code chunk number 2: loadmyo500
###################################################

library(IPPD) 
data(myo500)
x <- myo500[,"mz"]
y <- myo500[,"intensities"]
y <- y[x <= 2500]
x <- x[x <= 2500]


###################################################
### code chunk number 3: plotmyo500first1000
###################################################

layout(matrix(c(1,2), 1, 2))

plot(x[1:1000], y[1:1000], xlab = expression(x[1]~~ldots~~x[1000]), 
     cex.lab = 1.5, cex.axis = 1.25, ylab = expression(y))

plot(x[x >= 804 & x <= 807], y[x >= 804 & x <= 807], 
     xlab = "x: 804 <= x <= 807", 
     cex.lab = 1.5, cex.axis = 1.25, ylab = expression(y), type = "b")

layout(matrix(1))



###################################################
### code chunk number 4: fitGauss
###################################################

fitGauss <- fitModelParameters(mz = x, intensities = y,
                   model = "Gaussian", fitting = "model", formula.sigma = formula(~mz),
                   control = list(window = 6, threshold = 200))


###################################################
### code chunk number 5: fitEMG
###################################################

fitEMG <- fitModelParameters(mz = x, intensities = y,
                   model = "EMG", fitting = "model",
                           formula.alpha = formula(~mz),
                           formula.sigma = formula(~mz),
                           formula.mu = formula(~1),
                           control = list(window = 6, threshold = 200))


###################################################
### code chunk number 6: assessfit
###################################################

show(fitEMG)
mse.EMG <- data.frame(mse = slot(fitEMG,"peakfitresults")[,"rss"] 
                      / slot(fitEMG,"peakfitresults")[,"datapoints"], 
                      peakshape = rep("EMG", nrow( slot(fitEMG,"peakfitresults"))))
mse.Gauss <- data.frame(mse = slot(fitGauss,"peakfitresults")[,"rss"] 
                        / slot(fitGauss,"peakfitresults")[,"datapoints"], 
                        peakshape = rep("Gaussian", nrow( slot(fitGauss,"peakfitresults"))))
mses <- rbind(mse.EMG, mse.Gauss)
with(mses, boxplot(mse ~ peakshape, cex.axis = 1.5, cex.lab = 1.5, ylab = "MSE"))



###################################################
### code chunk number 7: assessfitEMGpeak
###################################################
visualize(fitEMG, type = "peak", cex.lab = 1.5, cex.axis = 1.25)


###################################################
### code chunk number 8: assessfit2
###################################################

visualize(fitEMG, type = "model", modelfit = TRUE, 
          parameters = c("sigma", "alpha"), 
          cex.lab = 1.5, cex.axis = 1.25)



###################################################
### code chunk number 9: getlist
###################################################

EMGlist <- getPeaklist(mz = x, intensities = y, model = "EMG",
model.parameters = fitEMG, 
loss = "L2", trace = FALSE,
control.localnoise = list(factor.place = 2),
control.basis = list(charges = c(1, 2)),
control.postprocessing = list(ppm = 200))

show(EMGlist)



###################################################
### code chunk number 10: showlist
###################################################

threshold(EMGlist, threshold = 3, refit = TRUE, trace = FALSE)



###################################################
### code chunk number 11: plot1
###################################################

visualize(EMGlist, x, y, lower= 963, upper = 973,
           fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
           localnoise = TRUE, quantile = 0.5,
           cutoff.functions = 3)



###################################################
### code chunk number 12: plot2
###################################################


visualize(EMGlist, x, y, lower= 1502, upper = 1510,
           fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
           localnoise = TRUE, quantile = 0.5,
           cutoff.functions = 2)




###################################################
### code chunk number 13: plot3
###################################################

visualize(EMGlist, x, y, lower= 1360, upper = 1364,
          fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
          localnoise = TRUE, quantile = 0.5, 
          cutoff.functions = 2)



###################################################
### code chunk number 14: clear
###################################################

rm(list = ls())



###################################################
### code chunk number 15: mzXML
###################################################

directory <- system.file("data", package = "IPPD")

download.file("http://www.ml.uni-saarland.de/code/IPPD/CytoC_1860-2200_500-600.mzXML",
              destfile = paste(directory, "/samplefile", sep = ""),
              quiet = TRUE)

data <- read.mzXML(paste(directory, "/samplefile", sep = ""))



###################################################
### code chunk number 16: sweepline (eval = FALSE)
###################################################
## 
## processLCMS <- analyzeLCMS(data, 
##                      arglist.getPeaklist = list(control.basis = list(charges = c(1,2,3))),
##                      arglist.threshold = list(threshold = 10),
##                      arglist.sweepline = list(minboxlength = 20))
## 
## boxes <- processLCMS$boxes


###################################################
### code chunk number 17: loadresults
###################################################

directory <- system.file("data", package = "IPPD")
filename <- file.path(directory, "examples_boxes.txt") 
boxes <- as.matrix(read.table(filename, header = TRUE))



###################################################
### code chunk number 18: displaysweepline
###################################################

print(boxes)

rtlist <- lapply(data$scan, function(x)
                 as.numeric(sub("([^0-9]*)([0-9|.]+)([^0-9]*)", "\\2", x$scanAttr)))

rt <- unlist(rtlist)

nscans <- length(rt)

npoints <- length(data$scan[[1]]$mass)

Y <- matrix(unlist(lapply(data$scan, function(x) x$peaks)),
            nrow = nscans, 
            ncol = npoints,
            byrow = TRUE)




contour(rt, data$scan[[1]]$mass, Y, xlab = "t", ylab = "mz", 
        levels = 10^(seq(from = 5, to = 6.75, by = 0.25)),
        drawlabels = FALSE)

for(i in 1:nrow(boxes))
  lines(c(boxes[i,"rt_begin"], boxes[i,"rt_end"]), rep(boxes[i,"loc"], 2), col = "red")