#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.