output:

BiocStyle::html_document

## <script type="text/javascript">
## document.addEventListener("DOMContentLoaded", function() {
##   document.querySelector("h1").style.marginTop = "0";
## });
## </script>
## <script type="text/javascript">
## document.addEventListener("DOMContentLoaded", function() {
##   var links = document.links;  
##   for (var i = 0, linksLength = links.length; i < linksLength; i++)
##     if (links[i].hostname != window.location.hostname)
##       links[i].target = '_blank';
## });
## </script>

Visualisation of proteomics data using R and Bioconductor

Laurent Gatto, Lisa Breckels and Sebastian Gibb

Introduction

References

• CRAN Task View: Graphic Displays & Dynamic Graphics & Graphic Devices & Visualization: http://cran.r-project.org/web/views/Graphics.html
• CRAN Task View: Web Technologies and Services: http://cran.r-project.org/web/views/WebTechnologies.html
• ggplot2 book (syntax is slightly outdated) (code), web page and on-line docs
• lattice book and web page
• R Graphics book
• R Cookbook and R Graphics Cookbook

Relevant packages

There are currently 65 Proteomics and 44 MassSpectrometry packages in Bioconductor version 3.0. Other non-Bioconductor packages are described in the RforProteomics vignette.

Ascombe's quartet

tab <- matrix(NA, 5, 4)
colnames(tab) <- 1:4
rownames(tab) <- c("var(x)", "mean(x)",
                   "var(y)", "mean(y)",
                   "cor(x,y)")

for (i in 1:4)
    tab[, i] <- c(var(anscombe[, i]),
                  mean(anscombe[, i]),
                  var(anscombe[, i+4]),
                  mean(anscombe[, i+4]),
                  cor(anscombe[, i], anscombe[, i+4]))

ff <- y ~ x

par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))
for (i in 1:4) {
    ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
    plot(ff, data = anscombe, pch = 19, xlim = c(3, 19), ylim = c(3, 13))
    abline(lm(ff, data = anscombe))
}

The MA plot example

The following code chunk connects to the PXD000001 data set on the ProteomeXchange repository and fetches the mzTab file. After missing values filtering, we extract relevant data (log2 fold-changes and log10 mean expression intensities) into data.frames.

library("rpx")
px1 <- PXDataset("PXD000001")
mztab <- pxget(px1, "PXD000001_mztab.txt")

## Downloading 1 file

library("MSnbase")
qnt <- readMzTabData(mztab, what = "PEP")

## Detected a metadata section
## Detected a peptide section

sampleNames(qnt) <- reporterNames(TMT6)
qnt <- filterNA(qnt)
## may be combineFeatuers

spikes <- c("P02769", "P00924", "P62894", "P00489")
protclasses <- as.character(fData(qnt)$accession)
protclasses[!protclasses %in% spikes] <- "Background"


madata42 <- data.frame(A = rowMeans(log(exprs(qnt[, c(4, 2)]), 10)),
                       M = log(exprs(qnt)[, 4], 2) - log(exprs(qnt)[, 2], 2),
                       data = rep("4vs2", nrow(qnt)),
                       protein = fData(qnt)$accession,
                       class = protclasses)

madata62 <- data.frame(A = rowMeans(log(exprs(qnt[, c(6, 2)]), 10)),
                       M = log(exprs(qnt)[, 6], 2) - log(exprs(qnt)[, 2], 2),
                       data = rep("6vs2", nrow(qnt)),
                       protein = fData(qnt)$accession,
                       class = protclasses)


madata <- rbind(madata42, madata62)

The traditional plotting system

par(mfrow = c(1, 2))
plot(M ~ A, data = madata42, main = "4vs2",
     xlab = "A", ylab = "M", col = madata62$class)
plot(M ~ A, data = madata62, main = "6vs2",
     xlab = "A", ylab = "M", col = madata62$class)

lattice

library("lattice")
latma <- xyplot(M ~ A | data, data = madata,
                groups = madata$class,
                auto.key = TRUE)
print(latma)

ggplot2

library("ggplot2")
ggma <- ggplot(aes(x = A, y = M, colour = class), data = madata,
               colour = class) +
                   geom_point() +
                       facet_grid(. ~ data)
print(ggma)

Customization

library("RColorBrewer")
bcols <- brewer.pal(4, "Set1")
cls <- c("Background" = "#12121230",
         "P02769" = bcols[1],
         "P00924" = bcols[2],
         "P62894" = bcols[3],
         "P00489" = bcols[4])

ggma2 <- ggplot(aes(x = A, y = M, colour = class),
                data = madata) + geom_point(shape = 19) +
                    facet_grid(. ~ data) + scale_colour_manual(values = cls) +
                        guides(colour = guide_legend(override.aes = list(alpha = 1)))
print(ggma2)

The MAplot method for MSnSet instances

MAplot(qnt, cex = .8)

An interactive shiny app for MA plots

This app is based on Mike Love's shinyMA application, adapted for a proteomics data. A screen shot is displayed below. To start the application:

shinyMA()

The application is also available online at https://lgatto.shinyapps.io/shinyMA/.

Visualising mass spectrometry data

Direct access to the raw data

library("lattice")
library("mzR")
mzf <- pxget(px1, 6)

## Downloading 1 file

ms <- openMSfile(mzf)

hd <- header(ms)
ms1 <- which(hd$msLevel == 1)

rtsel <- hd$retentionTime[ms1] / 60 > 30 & hd$retentionTime[ms1] / 60 < 35
library("MSnbase")
(M <- MSmap(ms, ms1[rtsel], 521, 523, .005, hd))

## 1

## Object of class "MSmap"
##  Map [75, 401]
##   [1]  Retention time: 30:1 - 34:58 
##   [2]  M/Z: 521 - 523 (res 0.005)

ff <- colorRampPalette(c("yellow", "steelblue"))
trellis.par.set(regions=list(col=ff(100)))
plot(M, aspect = 1, allTicks = FALSE)

M@map[msMap(M) == 0] <- NA
plot3D(M, FALSE)

library("rgl")
plot3D(M, TRUE)

lout <- matrix(NA, ncol = 10, nrow = 8)
lout[1:2, ] <- 1
for (ii in 3:4)
    lout[ii, ] <- c(2, 2, 2, 2, 2, 2, 3, 3, 3, 3)
lout[5, ] <- rep(4:8, each = 2)
lout[6, ] <- rep(4:8, each = 2)
lout[7, ] <- rep(9:13, each = 2)
lout[8, ] <- rep(9:13, each = 2)

i <- ms1[which(rtsel)][1]
j <- ms1[which(rtsel)][2]
ms2 <- (i+1):(j-1)

layout(lout)

par(mar=c(4,2,1,1))
chromatogram(ms)
abline(v = hd[i, "retentionTime"], col = "red")


par(mar = c(3, 2, 1, 0))
plot(peaks(ms, i), type = "l", xlim = c(400, 1000))
legend("topright", bty = "n",
       legend = paste0(
           "Acquisition ", hd[i, "acquisitionNum"],  "\n",
           "Retention time ", formatRt(hd[i, "retentionTime"])))
abline(h = 0)
abline(v = hd[ms2, "precursorMZ"],
       col = c("#FF000080",
           rep("#12121280", 9)))

par(mar = c(3, 0.5, 1, 1))
plot(peaks(ms, i), type = "l", xlim = c(521, 522.5),
     yaxt = "n")
abline(h = 0)
abline(v = hd[ms2, "precursorMZ"], col = "#FF000080")

##par(mar = omar)
par(mar = c(2, 2, 0, 1))
for (ii in ms2) {
    p <- peaks(ms, ii)
    plot(p, xlab = "", ylab = "", type = "h", cex.axis = .6)
    legend("topright", legend = paste0("Prec M/Z\n",
                           round(hd[ii, "precursorMZ"], 2)),
           bty = "n", cex = .8)
}

M2 <- MSmap(ms, i:j, 100, 1000, 1, hd)

## 1

plot3D(M2)

## 1

MS barcoding

par(mar=c(4,1,1,1))
image(t(matrix(hd$msLevel, 1, nrow(hd))),
      xlab="Retention time",
      xaxt="n", yaxt="n", col=c("black","steelblue"))
k <- round(range(hd$retentionTime) / 60)
nk <- 5
axis(side=1, at=seq(0,1,1/nk), labels=seq(k[1],k[2],k[2]/nk))

The MSnbase infrastructure

library("MSnbase")
data(itraqdata)
itraqdata2 <- pickPeaks(itraqdata, verbose = FALSE)
plot(itraqdata[[25]], full=TRUE, reporters = iTRAQ4)

par(oma = c(0, 0, 0, 0))
par(mar = c(4, 4, 1, 1))
plot(itraqdata2[[25]], itraqdata2[[28]], sequences = rep("IMIDLDGTENK", 2))

Preprocessing of MALDI-MS spectra

The following code chunks demonstrate the usage of the mass spectrometry preprocessing and plotting routines in the MALDIquant package. MALDIquant uses the traditional graphics system. Therefore MALDIquant overloads the traditional functions plot, lines and points for its own data types. These data types represents spectrum and peak lists as S4 classes. Please see the MALDIquant vignette and the corresponding website for more details.

After loading some example data a simple plot draws the raw spectrum.

library("MALDIquant")

data("fiedler2009subset", package="MALDIquant")

plot(fiedler2009subset[[14]])

After some preprocessing, namely variance stabilization and smoothing, we use lines to draw our baseline estimate in our processed spectrum.

transformedSpectra <- transformIntensity(fiedler2009subset, method = "sqrt")
smoothedSpectra <- smoothIntensity(transformedSpectra, method = "SavitzkyGolay")

plot(smoothedSpectra[[14]])
lines(estimateBaseline(smoothedSpectra[[14]]), lwd = 2, col = "red")

After removing the background removal we could use plot again to draw our baseline corrected spectrum.

rbSpectra <- removeBaseline(smoothedSpectra)
plot(rbSpectra[[14]])

detectPeaks returns a MassPeaks object that offers the same traditional graphics functions. The next code chunk demonstrates how to mark the detected peaks in a spectrum.

cbSpectra <- calibrateIntensity(rbSpectra, method = "TIC")
peaks <- detectPeaks(cbSpectra, SNR = 5)

plot(cbSpectra[[14]])
points(peaks[[14]], col = "red", pch = 4, lwd = 2)

Additional there is a special function labelPeaks that allows to draw the M/Z values above the corresponding peaks. Next we mark the 5 top peaks in the spectrum.

top5 <- intensity(peaks[[14]]) %in% sort(intensity(peaks[[14]]),
                                         decreasing = TRUE)[1:5]
labelPeaks(peaks[[14]], index = top5, avoidOverlap = TRUE)

Often multiple spectra have to be recalibrated to be comparable. Therefore MALDIquant warps the spectra according to so called reference or landmark peaks. For debugging the determineWarpingFunctions function offers some warping plots. Here we show only the last 4 plots:

par(mfrow = c(2, 2))
warpingFunctions <- determineWarpingFunctions(peaks,
                                              tolerance = 0.001,
                                              plot = TRUE,
                                              plotInteractive = TRUE)

par(mfrow = c(1, 1))

warpedSpectra <- warpMassSpectra(cbSpectra, warpingFunctions)
warpedPeaks <- warpMassPeaks(peaks, warpingFunctions)

In the next code chunk we visualise the need and the effect of the recalibration.

sel <- c(2, 10, 14, 16)
xlim <- c(4180, 4240)
ylim <- c(0, 1.9e-3)
lty <- c(1, 4, 2, 6)

par(mfrow = c(1, 2))
plot(cbSpectra[[1]], xlim = xlim, ylim = ylim, type = "n")

for (i in seq(along = sel)) {
  lines(peaks[[sel[i]]], lty = lty[i], col = i)
  lines(cbSpectra[[sel[i]]], lty = lty[i], col = i)
}

plot(cbSpectra[[1]], xlim = xlim, ylim = ylim, type = "n")

for (i in seq(along = sel)) {
  lines(warpedPeaks[[sel[i]]], lty = lty[i], col = i)
  lines(warpedSpectra[[sel[i]]], lty = lty[i], col = i)
}

par(mfrow = c(1, 1))

The code chunks above generate plots that are very similar to the figure 7 in the corresponding paper “Visualisation of proteomics data using R”. Please find the code to exactly reproduce the figure at: https://github.com/sgibb/MALDIquantExamples/blob/master/R/createFigure1_color.R

Genomic and protein sequences

These visualisations originate from the Pbase Pbase-data and mapping vignettes.

Imaging mass spectrometry

The following code chunk downloads a MALDI imaging dataset from a mouse kidney shared by Adrien Nyakas and Stefan Schurch and generates a plot with the mean spectrum and three slices of interesting M/Z regions.

library("MALDIquant")
library("MALDIquantForeign")

spectra <- importBrukerFlex("http://files.figshare.com/1106682/MouseKidney_IMS_testdata.zip", verbose = FALSE)

spectra <- smoothIntensity(spectra, "SavitzkyGolay",  halfWindowSize = 8)
spectra <- removeBaseline(spectra, method = "TopHat", halfWindowSize = 16)
spectra <- calibrateIntensity(spectra, method = "TIC")
avgSpectrum <- averageMassSpectra(spectra)
avgPeaks <- detectPeaks(avgSpectrum, SNR = 5)

avgPeaks <- avgPeaks[intensity(avgPeaks) > 0.0015]

oldPar <- par(no.readonly = TRUE)
layout(matrix(c(1,1,1,2,3,4), nrow = 2, byrow = TRUE))
plot(avgSpectrum, main = "mean spectrum",
     xlim = c(3000, 6000), ylim = c(0, 0.007))
lines(avgPeaks, col = "red")
labelPeaks(avgPeaks, cex = 1)

par(mar = c(0.5, 0.5, 1.5, 0.5))
for (i in seq(along = avgPeaks)) {
  range <- mass(avgPeaks)[i] + c(-1, 1)
  plotImsSlice(spectra, range = range,
               main = paste(round(range, 2), collapse = " - "))
}
par(oldPar)

An interactive shiny app for Imaging mass spectrometry

There is also an interactive MALDIquant IMS shiny app for demonstration purposes. A screen shot is displayed below. To start the application:

library("shiny")
runGitHub("sgibb/ims-shiny")

Spatial proteomics

library("pRoloc")
library("pRolocdata")

data(tan2009r1)

## these params use class weights
fn <- dir(system.file("extdata", package = "pRoloc"),
          full.names = TRUE, pattern = "params2.rda")
load(fn)

setStockcol(NULL)
setStockcol(paste0(getStockcol(), 90))

w <- table(fData(tan2009r1)[, "pd.markers"])
(w <- 1/w[names(w) != "unknown"])

## 
##  Cytoskeleton            ER         Golgi      Lysosome       Nucleus 
##    0.14285714    0.05000000    0.16666667    0.12500000    0.05000000 
##            PM    Peroxisome    Proteasome  Ribosome 40S  Ribosome 60S 
##    0.06666667    0.25000000    0.09090909    0.07142857    0.04000000 
## mitochondrion 
##    0.07142857

tan2009r1 <- svmClassification(tan2009r1, params2,
                               class.weights = w,
                               fcol = "pd.markers")
ptsze <- exp(fData(tan2009r1)$svm.scores) - 1

lout <- matrix(c(1:4, rep(5, 4)), ncol = 4, nrow = 2)
layout(lout)
cls <- getStockcol()
par(mar = c(4, 4, 1, 1))
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "mitochondrion"), ],
         markers = featureNames(tan2009r1)[which(fData(tan2009r1)$markers.orig == "mitochondrion")],
         mcol = cls[5])
legend("topright", legend = "mitochondrion", bty = "n")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ],
         markers = featureNames(tan2009r1)[which(fData(tan2009r1)$markers.orig == "ER")],
         mcol = cls[2])
legend("topright", legend = "ER", bty = "n")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ],
         markers = featureNames(tan2009r1)[which(fData(tan2009r1)$markers.orig == "Golgi")],
         mcol = cls[3])
legend("topright", legend = "Golgi", bty = "n")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "PM"), ],
         markers = featureNames(tan2009r1)[which(fData(tan2009r1)$markers.orig == "PM")],
         mcol = cls[8])
legend("topright", legend = "PM", bty = "n")
plot2D(tan2009r1, fcol = "svm", cex = ptsze, method = "kpca")
addLegend(tan2009r1, where = "bottomleft", fcol = "svm", bty = "n")

See the pRoloc-tutorial vignette (pdf) from the pRoloc package for details about spatial proteomics data analysis and visualisation.

Session information

print(sessionInfo(), locale = FALSE)

## R version 3.1.2 (2014-10-31)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggplot2_1.0.0         lattice_0.20-29       pRolocdata_1.4.1     
##  [4] pRoloc_1.6.0          MLInterfaces_1.46.0   cluster_1.15.3       
##  [7] annotate_1.44.0       XML_3.98-1.1          AnnotationDbi_1.28.1 
## [10] GenomeInfoDb_1.2.3    IRanges_2.0.0         S4Vectors_0.4.0      
## [13] MALDIquantForeign_0.9 MALDIquant_1.11       RColorBrewer_1.0-5   
## [16] xtable_1.7-4          rpx_1.2.0             knitr_1.8            
## [19] BiocInstaller_1.16.1  RforProteomics_1.4.1  MSnbase_1.14.1       
## [22] BiocParallel_1.0.0    mzR_2.0.0             Rcpp_0.11.3          
## [25] Biobase_2.26.0        BiocGenerics_0.12.1   BiocStyle_1.4.1      
## 
## loaded via a namespace (and not attached):
##  [1] BBmisc_1.8                   BatchJobs_1.5               
##  [3] BradleyTerry2_1.0-5          Category_2.32.0             
##  [5] DBI_0.3.1                    FNN_1.1                     
##  [7] GSEABase_1.28.0              MASS_7.3-35                 
##  [9] Matrix_1.1-4                 R.methodsS3_1.6.1           
## [11] R.oo_1.18.0                  R.utils_1.34.0              
## [13] R6_2.0.1                     RBGL_1.42.0                 
## [15] RCurl_1.95-4.3               RJSONIO_1.3-0               
## [17] RSQLite_1.0.0                RUnit_0.4.27                
## [19] affy_1.44.0                  affyio_1.34.0               
## [21] base64enc_0.1-2              biocViews_1.34.1            
## [23] brew_1.0-6                   brglm_0.5-9                 
## [25] car_2.0-21                   caret_6.0-37                
## [27] checkmate_1.5.0              class_7.3-11                
## [29] codetools_0.2-9              colorspace_1.2-4            
## [31] digest_0.6.4                 doParallel_1.0.8            
## [33] downloader_0.3               e1071_1.6-4                 
## [35] evaluate_0.5.5               fail_1.2                    
## [37] foreach_1.4.2                formatR_1.0                 
## [39] gdata_2.13.3                 genefilter_1.48.1           
## [41] graph_1.44.0                 grid_3.1.2                  
## [43] gridSVG_1.4-0                gtable_0.1.2                
## [45] gtools_3.4.1                 highr_0.4                   
## [47] htmltools_0.2.6              httpuv_1.3.2                
## [49] impute_1.40.0                interactiveDisplay_1.4.0    
## [51] interactiveDisplayBase_1.4.0 iterators_1.0.7             
## [53] kernlab_0.9-19               labeling_0.3                
## [55] limma_3.22.1                 lme4_1.1-7                  
## [57] lpSolve_5.6.10               mclust_4.4                  
## [59] mime_0.2                     minqa_1.2.4                 
## [61] munsell_0.4.2                mvtnorm_1.0-1               
## [63] mzID_1.4.1                   nlme_3.1-118                
## [65] nloptr_1.0.4                 nnet_7.3-8                  
## [67] pcaMethods_1.56.0            pls_2.4-3                   
## [69] plyr_1.8.1                   preprocessCore_1.28.0       
## [71] proto_0.3-10                 proxy_0.4-13                
## [73] randomForest_4.6-10          rda_1.0.2-2                 
## [75] readBrukerFlexData_1.8.1     readMzXmlData_2.8           
## [77] reshape2_1.4                 rpart_4.1-8                 
## [79] sampling_2.6                 scales_0.2.4                
## [81] sendmailR_1.2-1              sfsmisc_1.0-26              
## [83] shiny_0.10.2.1               splines_3.1.2               
## [85] stringr_0.6.2                survival_2.37-7             
## [87] tools_3.1.2                  vsn_3.34.0                  
## [89] zlibbioc_1.12.0