Merge pull request #314 from lgatto/plotMsData

LaurentGatto · web-flow · commit f44b035cb9f9 · 2018-02-20T18:22:57.000Z
Add plot,MSnExp,missing, type = "XIC" (issue #313)
diff --git a/NAMESPACE b/NAMESPACE
@@ -17,7 +17,7 @@ importFrom(XML, xmlInternalTreeParse, getDefaultNamespace, xpathApply,
            xmlAttrs)
 
 importFrom(graphics, abline, axis, grid, legend, lines, points, text,
-           barplot, layout, par, matplot)
+           barplot, layout, par, matplot, mtext)
 importFrom(stats, ave, median, medpolish, quantile, reorder, setNames,
            weighted.mean, dist)
 importFrom(utils, combn, head, ls.str, object.size,
@@ -49,6 +49,7 @@ importFrom(affy, MAplot, ma.plot, mva.pairs)
 importFrom(mzID, mzID, flatten)
 importClassesFrom(mzID, mzID, mzIDCollection)
 importFrom(digest, digest)
+importFrom(grDevices, topo.colors)
 
 export(MSnSet,
        abstract,
diff --git a/R/functions-MSnExp.R b/R/functions-MSnExp.R
@@ -355,3 +355,30 @@ removeReporters_MSnExp <- function(object, reporters = NULL,
     return(object)
 }
 
+plotXIC_MSnExp <- function(x, ...) {
+    ## Restrict to MS level 1
+    x <- filterMsLevel(x, 1L)
+    if (!length(x))
+        stop("No MS1 data available")
+    fns <- basename(fileNames(x))
+    if (isMSnbaseVerbose())
+        message("Retrieving data ...", appendLF = FALSE)
+    x <- as(x, "data.frame")
+    x <- split(x, x$file)
+    if (isMSnbaseVerbose())
+        message("OK")
+    ## Check if we are greedy and plot a too large area
+    if (any(unlist(lapply(x, nrow)) > 20000))
+        warning("The MS area to be plotted seems rather large. It is suggested",
+                " to restrict the data first using 'filterRt' and 'filterMz'. ",
+                "See also ?chromatogram and ?Chromatogram for more efficient ",
+                "functions to plot a total ion chromatogram or base peak ",
+                "chromatogram.",
+                immediate = TRUE, call = FALSE)
+    ## Define the layout.
+    dots <- list(...)
+    layout(.vertical_sub_layout(length(x)))
+    tmp <- mapply(x, fns, FUN = function(z, main, ...) {
+        .plotXIC(x = z, main = main, layout = NULL, ...)
+    }, MoreArgs = dots)
+}
diff --git a/R/functions-plotting.R b/R/functions-plotting.R
@@ -0,0 +1,76 @@
+#' Takes the values for a single file.
+#'
+#' @param x `data.frame` with columns `"mz"`, `"rt"` and `"i"`.
+#'
+#' @param main `character(1)` with the title of the plot.
+#'
+#' @param col color for the circles.
+#'
+#' @param colramp color ramp to be used for the points' background.
+#'
+#' @param grid.color color to be used for the grid lines (or `NA` if they should
+#'     not be plotted.
+#'
+#' @param pch The plotting character.
+#'
+#' @param layout `matrix` defining the layout of the plot, or `NULL` if
+#'     `layout` was already called.
+#'
+#' @param ... additional parameters to be passed to the `plot` function.
+#' 
+#' @md
+#'
+#' @author Johannes Rainer
+#' 
+#' @noRd
+.plotXIC <- function(x, main = "", col = "grey", colramp = topo.colors,
+                     grid.color = "lightgrey", pch = 21,
+                     layout = matrix(1:2, ncol = 1), ...) {
+    if (is.matrix(layout))
+        layout(layout)
+    ## Chromatogram.
+    bpi <- unlist(lapply(split(x$i, x$rt), max, na.rm = TRUE))
+    brks <- do.breaks(range(x$i), nint = 256)
+    par(mar = c(0, 4, 2, 1))
+    plot(as.numeric(names(bpi)), bpi, xaxt = "n", col = col, main = main, 
+         bg = level.colors(bpi, at = brks, col.regions = colramp), xlab = "",
+         pch = pch, ylab = "", las = 2, ...)
+    mtext(side = 4, line = 0, "Intensity", cex = par("cex.lab"))
+    grid(col = grid.color)
+    par(mar = c(3.5, 4, 0, 1))
+    plot(x$rt, x$mz, main = "", pch = pch, col = col, xlab = "", ylab = "",
+         yaxt = "n", bg = level.colors(x$i, at = brks, col.regions = colramp),
+         ...)
+    axis(side = 2, las = 2)
+    grid(col = grid.color)
+    mtext(side = 1, line = 2.5, "Retention time", cex = par("cex.lab"))
+    mtext(side = 4, line = 0, "m/z", cex = par("cex.lab"))
+}
+
+#' Create a `matrix` to be used for the `layout` function to allow plotting of
+#' vertically arranged *sub-plots* consisting of `sub_plot` plots.
+#'
+#' @param x `integer(1)` with the number of sub-plots.
+#'
+#' @param sub_plot `integer(1)` with the number of sub-plots per cell/plot.
+#' 
+#' @author Johannes Rainer
+#' 
+#' @md
+#'
+#' @noRd
+#'
+#' @examples
+#'
+#' ## Assum we've got 5 *features* to plot and we want to have a two plots for
+#' ## each feature arranged below each other.
+#'
+#' .vertical_sub_layout(5, sub_plot = 2)
+.vertical_sub_layout <- function(x, sub_plot = 2) {
+    sqrt_x <- sqrt(x)
+    ncol <- ceiling(sqrt_x)
+    nrow <- round(sqrt_x)
+    rws <- split(1:(ncol * nrow * sub_plot), f = rep(1:nrow,
+                                                     each = sub_plot * ncol))
+    do.call(rbind, lapply(rws, matrix, ncol = ncol))
+}
diff --git a/R/methods-MSnExp.R b/R/methods-MSnExp.R
@@ -85,7 +85,13 @@ setMethod("show", "MSnExp",
 
 
 setMethod("plot", c("MSnExp","missing"),
-          function(x, y ,...) plot_MSnExp(x, ...))
+          function(x, y , type = c("spectra", "XIC"), ...) {
+              type <- match.arg(type)
+              if (type == "spectra")
+                  plot_MSnExp(x, ...)
+              if (type == "XIC")
+                  plotXIC_MSnExp(x, ...)
+})
 
 setMethod("plot2d", c("MSnExp"),
           function(object, z, alpha = 1/3, plot = TRUE)
@@ -538,8 +544,12 @@ setMethod("chromatogram", "MSnExp", function(object, rt, mz,
 
 setAs("MSnExp", "data.frame", function(from) {
     do.call(rbind, unname(spectrapply(from, function(z) {
-        ## Directly accessing slots is faster than using methods
-        data.frame(file = z@fromFile, rt = z@rt, mz = z@mz, i = z@intensity)
+        if (length(z@mz))
+            ## Directly accessing slots is faster than using methods
+            data.frame(file = z@fromFile, rt = z@rt, mz = z@mz, i = z@intensity)
+        else
+            data.frame(file = integer(), rt = numeric(), mz = numeric(),
+                       i = numeric())
     })))
 })
 as.data.frame.MSnExp <- function(x, row.names = NULL, optional=FALSE, ...)
diff --git a/man/MSnExp-class.Rd b/man/MSnExp-class.Rd
@@ -186,9 +186,9 @@
       the noise in all profile spectra of \code{object}.  See
       \code{\link{estimateNoise}} documentation for more details and
       examples. }
-    \item{plot}{\code{signature(x = "MSnExp", y = "missing")}: Plots all
-      the spectra of the \code{MSnExp} instance. See
-      \code{\link{plot.MSnExp}} documentation for more details. }
+    \item{plot}{\code{signature(x = "MSnExp", y = "missing")}: Plots
+      the \code{MSnExp} instance. See \code{\link{plot.MSnExp}}
+      documentation for more details. }
     \item{plot2d}{\code{signature(object = "MSnExp", ...)}:
       Plots retention time against precursor MZ for \code{MSnExp}
       instances. See \code{\link{plot2d}} documentation for more
diff --git a/man/plot-methods.Rd b/man/plot-methods.Rd
@@ -9,12 +9,17 @@
 \alias{plot.Spectrum.character}
 \alias{plot}
 
-\title{ Plotting 'Spectrum' object(s) }
+\title{ Plotting 'MSnExp' and 'Spectrum' object(s) }
 
 \description{
 
-  These method plot mass spectra MZ values against the intensities. Full
-  spectra (using the \code{full} parameter) or specific peaks of
+  These methods provide the functionality to plot mass spectrometry data
+  provided as \code{\linkS4class{MSnExp}},
+  \code{\linkS4class{OnDiskMSnExp}} or \code{\linkS4class{Spectrum}}
+  objects. Most functions plot mass spectra M/Z values against
+  intensities.
+
+  Full spectra (using the \code{full} parameter) or specific peaks of
   interest can be plotted using the \code{reporters} parameter. If
   \code{reporters} are specified and \code{full} is set to 'TRUE', a
   sub-figure of the reporter ions is inlaid inside the full spectrum.
@@ -24,8 +29,8 @@
   extract spectra of interest using the \code{[} operator or
   \code{\link{extractPrecSpectra}} methods.
 
-  The methods make use the \code{ggplot2} system. An object of class
-  'ggplot' is returned invisibly.
+  Most methods make use the \code{ggplot2} system in which case an
+  object of class 'ggplot' is returned invisibly.
 
   If a single \code{"\linkS4class{Spectrum2}"} and a \code{"character"}
   representing a valid peptide sequence are passed as argument, the
@@ -46,17 +51,22 @@
   \item{reporters}{ An object of class
     \code{"\linkS4class{ReporterIons}"} that defines the peaks to be
     plotted. If not specified, \code{full} must be set to 'TRUE'.} 
+
   \item{full}{Logical indicating whether full spectrum (respectively
     spectra) of  only reporter ions of interest should be
     plotted. Default is 'FALSE', in which case \code{reporters} must be
     defined. }
+
   \item{centroided.}{Logical indicating if spectrum or spectra are in
     centroided mode, in which case peaks are plotted as histograms,
     rather than curves.}
+
   \item{plot}{Logical specifying whether plot should be printed to
     current device. Default is 'TRUE'.}
+
   \item{w1}{Width of sticks for full centroided spectra. Default is to
     use maximum MZ value divided by 500. }
+
   \item{w2}{Width of histogram bars for centroided reporter ions
     plots. Default is 0.01. }
 
@@ -66,14 +76,51 @@
 
 \section{Methods}{
   \describe{
-    \item{\code{signature(x = "MSnExp", y = "missing", reporters =
-	"ReporterIons", full = "logical", plot = "logical")}}{ Plots
-      all the spectra in the \code{MSnExp} object vertically. One of
-      \code{reporters} must be defined or \code{full} set to 'TRUE'. In
-      case of \code{MSnExp} objects, repoter ions are not inlaid when
-      \code{full} is 'TRUE'. 
+    \item{\code{plot(signature(x = "MSnExp", y = "missing"),
+	type = c("spectra", "XIC"), reporters = "ReporterIons",
+	full = "logical", plot = "logical", ...)}}{
+
+      For \code{type = "spectra"}: Plots all the spectra in the
+      \code{MSnExp} object vertically. One of \code{reporters} must be
+      defined or \code{full} set to 'TRUE'. In case of \code{MSnExp}
+      objects, repoter ions are not inlaid when \code{full} is 'TRUE'.
+
+      For \code{type = "XIC"}: Plots a combined plot of retention time
+      against m/z values and retention time against largest signal per
+      spectrum for each file. Data points are colored by intensity. The
+      lower part of the plot represents the location of the individual
+      signals in the retention time - m/z space, the upper part the base
+      peak chromatogram of the data (i.e. the largest signal for each
+      spectrum). This plot type is restricted to MS level 1 data and is
+      most useful for LC-MS data.
+      Ideally, the \code{MSnExp} (or \code{OnDiskMSnExp})
+      object should be filtered first using the \code{\link{filterRt}}
+      and \code{\link{filterMz}} functions to narrow on an ion of
+      interest. See examples below. This plot uses base R
+      plotting. Additional arguments to the \code{plot} function can be
+      passed with \code{...}.
+      
+      Additional arguments for \code{type = "XIC"} are:
+      \describe{
+	\item{\code{col}}{color for the border of the points. Defaults to
+	  \code{col = "grey"}.}
+
+	\item{\code{colramp}}{color function/ramp to be used for the
+	  intensity-dependent background color of data points. Defaults
+	  to \code{colramp = topo.colors}.}
+
+	\item{\code{grid.color}}{color for the grid lines. Defaults to
+	  \code{grid.color = "lightgrey"}; use \code{grid.color = NA} to
+	  disable grid lines altogether.}
+
+	\item{\code{pch}}{point character. Defaults to \code{pch = 21}}.
+
+	\item{\code{...}}{additional parameters for the low-level
+	  \code{plot} function.}
+	
+      }
     }
-    \item{\code{signature(x = "Spectrum", y = "missing", reporters =
+    \item{\code{plot(signature(x = "Spectrum", y = "missing"), reporters =
 	"ReporterIons", full = "logical", centroided. = "logical", plot =
 	"logical", w1, w2)}}{ Displays the MZs against intensities of
       the \code{Spectrum} object as a line plot.  
@@ -85,7 +132,7 @@
       by default.
     }
 
-    \item{\code{signature(x = "Spectrum2", y = "character", orientation
+    \item{\code{plot(signature(x = "Spectrum2", y = "character"), orientation
 	= "numeric", add = "logical", col = "character", pch, xlab =
 	"character", ylab = "character", xlim = "numeric", ylim =
 	"numeric", tolerance = "numeric", relative = "logical", type =
@@ -111,10 +158,12 @@
   \code{\link{calculateFragments}} to calculate ions produced by
   fragmentation and \code{\link{plot.Spectrum.Spectrum}} to plot and
   compare 2 spectra and their shared peaks.
+
+  \code{\link{Chromatogram}} for plotting of chromatographic data.
 }
 
 \author{
-  Laurent Gatto <lg390@cam.ac.uk> and Sebastian Gibb
+  Laurent Gatto <lg390@cam.ac.uk>, Johannes Rainer and Sebastian Gibb
 }
 
 \examples{
@@ -129,6 +178,18 @@ itraqdata2 <- pickPeaks(itraqdata)
 i <- 14
 s <- as.character(fData(itraqdata2)[i, "PeptideSequence"])
 plot(itraqdata2[[i]], s, main = s)
+
+## Load profile-mode LC-MS files
+library(msdata)
+od <- readMSData(dir(system.file("sciex", package = "msdata"),
+                     full.names = TRUE), mode = "onDisk")
+## Restrict the MS data to signal for serine
+serine <- filterMz(filterRt(od, rt = c(175, 190)), mz = c(106.04, 106.06))
+plot(serine, type = "XIC")
+
+## Same plot but using heat.colors, rectangles and no point border
+plot(serine, type = "XIC", pch = 22, colramp = heat.colors, col = NA)
+
 }
 
 \keyword{methods}
diff --git a/tests/testthat/test_MSnExp.R b/tests/testthat/test_MSnExp.R
@@ -506,3 +506,14 @@ test_that("setAs,MSnExp,data.frame works", {
     expect_equal(unlist(mz(res), use.names = FALSE), df$mz)
     expect_equal(unlist(intensity(res), use.names = FALSE), df$i)
 })
+
+
+test_that("plotXIC_MSnExp works", {
+    im <- microtofq_in_mem_ms1
+    expect_warning(plotXIC_MSnExp(filterMz(im, c(600, 680))))
+    plotXIC_MSnExp(filterMz(im, c(610, 615)), pch = 23)
+    ## filter to get only one
+    plotXIC_MSnExp(filterMz(filterRt(im, c(270, 280)), c(610, 615)), cex = 2)
+
+    expect_error(plotXIC_MSnExp(tmt_erwinia_in_mem_ms2))
+})
diff --git a/tests/testthat/test_functions-plotting.R b/tests/testthat/test_functions-plotting.R
@@ -0,0 +1,42 @@
+test_that(".plotXIC works", {
+    x <- data.frame(rt = rep(abs(rnorm(30)), 100),
+                    mz = rep(abs(rnorm(30)), 100),
+                    i = abs(rnorm(300, mean = 200)))
+    od1 <- filterFile(microtofq_on_disk, 1)
+
+    x <- as(filterMz(filterRt(od1, c(270, 290)), c(610, 615)), "data.frame")
+    .plotXIC(x)
+
+    ## Test passing additional arguments
+    .plotXIC(x, cex = 2)
+    .plotXIC(x, pch = 1)
+    .plotXIC(x, pch = 23)
+})
+
+test_that(".vertical_sub_layout works", {
+    exp <- cbind(c(1, 2, 7, 8),
+                 c(3, 4, 9, 10),
+                 c(5, 6, 11, 12))
+    expect_equal(.vertical_sub_layout(5, 2), exp)
+    ## 8 plots (will be 3 columns), 3 plots each.
+    exp <- cbind(c(1, 2, 3, 10, 11, 12, 19, 20, 21),
+                 c(4, 5, 6, 13, 14, 15, 22, 23, 24),
+                 c(7, 8, 9, 16, 17, 18, 25, 26, 27))
+    expect_equal(.vertical_sub_layout(8, 3), exp)
+
+    ## single plot, 5 each
+    exp <- matrix(1:5, ncol = 1)
+    expect_equal(.vertical_sub_layout(1, 5), exp)
+
+    ## 19 plots 1 each
+    exp <- matrix(1:20, ncol = 5, byrow = TRUE)
+    expect_equal(.vertical_sub_layout(19, 1), exp)
+
+    ## 19 plots 2 each
+    exp <- cbind(c(1, 2, 11, 12, 21, 22, 31, 32),
+                 c(3, 4, 13, 14, 23, 24, 33, 34),
+                 c(5, 6, 15, 16, 25, 26, 35, 36),
+                 c(7, 8, 17, 18, 27, 28, 37, 38),
+                 c(9, 10, 19, 20, 29, 30, 39, 40))
+    expect_equal(.vertical_sub_layout(19, 2), exp)
+})
