#' This function creates Figure 1 in Gibb and Strimmer 2012 in a colorized #' version. #' #' @title Figure 1 #' #' @references #' S. Gibb and K. Strimmer. 2012. MALDIquant: a versatile R package for the #' analysis of mass spectrometry data. Bioinformatics 28: 2270-2271 #' #' @seealso \code{\link[MALDIquantExamples]{createFigure1}} #' #' @export #' @examples #' \dontrun{ #' library("MALDIquantExamples") #' pdfWidth <- 7 #' pdfHeight <- pdfWidth*0.35 #' #' pdf(file="figure1.pdf", height=pdfHeight, width=pdfWidth) #' createFigure1Color() #' dev.off() #' } createFigure1Color <- function() { ## load example spectra data("fiedler2009subset", package="MALDIquant") ## some preprocessing ## sqrt transform (for variance stabilization) tSpectra <- transformIntensity(fiedler2009subset, method="sqrt") ## simple 5 point moving average for smoothing spectra tSpectra <- smoothIntensity(tSpectra, method="MovingAverage", halfWindowSize=2) ## remove baseline rbSpectra <- removeBaseline(tSpectra) ## calibrate intensity values by "total ion current" cbSpectra <- calibrateIntensity(rbSpectra, method="TIC") ## run peak detection peaks <- detectPeaks(cbSpectra, SNR=5) ### warping reference <- referencePeaks(peaks) warpingFunctions <- determineWarpingFunctions(peaks, reference=reference, tolerance=0.001) ## warp spectra warpedSpectra <- warpMassSpectra(cbSpectra, warpingFunctions) ## warp peaks warpedPeaks <- warpMassPeaks(peaks, warpingFunctions) ## merge technical replicates mergedSpectra <- averageMassSpectra(warpedSpectra, rep(1:8, each=2)) binnedPeaks <- binPeaks(warpedPeaks) mergedPeaks <- mergeMassPeaks(binnedPeaks, rep(1:8, each=2)) ## helper function to mark plots as LETTERS[1:4] labelPlot <- function(char, cex=1.5) { usr <- par("usr") text(x=usr[2]-(cex*strwidth(char)), y=usr[4]-(cex*strheight(char)), labels=char, cex=cex) } par(mfrow=c(2, 3)) par(cex=0.4) par(yaxt="n") par(mar=c(2.5, 1, 1, 1)) # bottom, left, top, right xlim <- c(1e3, 1e4) ## select 1 spectra for plot A/B AB <- 14 ## first row ## plot A plot(fiedler2009subset[[AB]], lwd=0.25, sub="", main="", ylab="", xlab="", xlim=xlim) lines(estimateBaseline(fiedler2009subset[[AB]]), lwd=0.75, col=2) labelPlot("A") ## plot B plot(cbSpectra[[AB]], lwd=0.25, sub="", main="", ylab="", xlab="", xlim=xlim) points(peaks[[AB]], pch=4, lwd=0.25, col=4) labelPlot("B") ## plot C par(yaxt="s") par(mar=c(2.5, 2, 1, 1)) # bottom, left, top, right determineWarpingFunctions(peaks[[10]], reference=reference, tolerance=0.001, plot=TRUE, plotInteractive=TRUE, ylim=c(-2, 4), lwd=0.5, xlab="", ylab="", main="", sub="") labelPlot("C") ## second row par(mar=c(4, 1, 1, 1)) # bottom, left, top, right par(yaxt="n") ## select 4 spectra for plot D/E DE <- c(2, 10, 14, 16) ## limits for plot D/E xlimDE <- c(4180, 4240) ylimDE <- c(0, 1.9e-3) ## line types lty <- c(1, 4, 2, 6) ## plot D plot(cbSpectra[[1]], xlim=xlimDE, ylim=ylimDE, type="n", main="", xlab="", ylab="") labelPlot("D") for (i in seq(along=DE)) { lines(peaks[[DE[i]]], lty=lty[i], lwd=0.5, col=i) lines(cbSpectra[[DE[i]]], lty=lty[i], lwd=0.5, col=i) } ## plot E plot(cbSpectra[[1]], xlim=xlimDE, ylim=ylimDE, type="n", main="", xlab="", ylab="") labelPlot("E") for (i in seq(along=DE)) { lines(warpedPeaks[[DE[i]]], lty=lty[i], lwd=0.5, col=i) lines(warpedSpectra[[DE[i]]], lty=lty[i], lwd=0.5, col=i) } ## plot F F <- 7 plot(mergedSpectra[[F]], lwd=0.25, sub="", main="", ylab="", xlab="", xlim=xlim, ylim=c(0, (max(intensity(mergedSpectra[[F]]))*1.05))) labelPlot("F") points(mergedPeaks[[F]], lwd=0.25, pch=4, col=4) ## label highest peaks top <- intensity(mergedPeaks[[F]]) %in% sort(intensity(mergedPeaks[[F]]), decreasing=TRUE)[1:10] labelPeaks(mergedPeaks[[F]], index=top, underline=TRUE, cex=0.8, lwd=0.25) }