R/plotAnnotationDiagnosticMultiplot.R

In peakPantheR: Peak Picking and Annotation of High Resolution Experiments

#' Generate a multiplot of all diagnostic plots
#'
#' Generate a multiplot of all diagnostic plots (as generated by \code{annotationDiagnosticPlots()}) for each compound
#' 
#' @param annotationDiagnosticPlotList (list) List of (one per compound) of list of diagnostic plots as generated by \code{annotationDiagnosticPlots()}
#' 
#' @return A list of multiplots (one per compound)
annotationDiagnosticMultiplot   <- function(annotationDiagnosticPlotList) {
  
  # Init
  nbCpd   <- length(annotationDiagnosticPlotList)
  outList <- vector("list", nbCpd)
  
  # Iterate over compounds
  for (cpd in 1:nbCpd) {
    
    ## Check input
    if (!all(c("EICFit", "rtPeakwidthVert", "rtPeakwidthHorzRunOrder", "mzPeakwidthHorzRunOrder", "areaRunOrder", "rtHistogram", "mzHistogram", "areaHistogram", "title") %in% names(annotationDiagnosticPlotList[[cpd]]))) {
      # cannot generate multiplot, pass to next compound
      message('Required plots missing for compound #',cpd)
      outList[[cpd]]  <- NULL
      next
    }
    
    ## Init plot
    tmp_EIC             <- annotationDiagnosticPlotList[[cpd]]$EICFit + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank())
    tmp_rtPeakwidthVert <- annotationDiagnosticPlotList[[cpd]]$rtPeakwidthVert
    tmp_rtPeakwidthHorz <- annotationDiagnosticPlotList[[cpd]]$rtPeakwidthHorzRunOrder + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank())
    tmp_mzPeakwidthHorz <- annotationDiagnosticPlotList[[cpd]]$mzPeakwidthHorzRunOrder + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank())
    tmp_peakAreaHorz    <- annotationDiagnosticPlotList[[cpd]]$areaRunOrder
    tmp_rtHisto         <- annotationDiagnosticPlotList[[cpd]]$rtHistogram + ggplot2::coord_flip() + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank(), axis.title.y=ggplot2::element_blank(), axis.text.y=ggplot2::element_blank())
    tmp_mzHisto         <- annotationDiagnosticPlotList[[cpd]]$mzHistogram + ggplot2::coord_flip() + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank(), axis.title.y=ggplot2::element_blank(), axis.text.y=ggplot2::element_blank())
    tmp_areaHisto       <- annotationDiagnosticPlotList[[cpd]]$areaHistogram + ggplot2::coord_flip() + ggplot2::theme(axis.title.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank(), axis.title.y=ggplot2::element_blank(), axis.text.y=ggplot2::element_blank())
    
    ## link axis (same lim)
    # EIC + rt peakwidth (is rotated, use y in histogram)
    minXEIC             <- min(ggplot2::layer_scales(tmp_EIC)$x$range$range[1], ggplot2::layer_scales(tmp_rtPeakwidthVert)$y$range$range[1])
    maxXEIC             <- max(ggplot2::layer_scales(tmp_EIC)$x$range$range[2], ggplot2::layer_scales(tmp_rtPeakwidthVert)$y$range$range[2])
    tmp_EIC             <- tmp_EIC + ggplot2::xlim(minXEIC, maxXEIC)
    tmp_rtPeakwidthVert <- tmp_rtPeakwidthVert + ggplot2::ylim(minXEIC, maxXEIC)
    # rt peakwidth runOrder + histo (is rotated, use x in histo)
    minYRT              <- min(ggplot2::layer_scales(tmp_rtPeakwidthHorz)$y$range$range[1], ggplot2::layer_scales(tmp_rtHisto)$x$range$range[1])
    maxYRT              <- max(ggplot2::layer_scales(tmp_rtPeakwidthHorz)$y$range$range[2], ggplot2::layer_scales(tmp_rtHisto)$x$range$range[2])
    tmp_rtPeakwidthHorz <- tmp_rtPeakwidthHorz + ggplot2::ylim(minYRT, maxYRT)
    tmp_rtHisto         <- tmp_rtHisto + ggplot2::xlim(minYRT, maxYRT)
    # mz peakwidth runOrder + histo (is rotated, use x in histo)
    minYMZ              <- min(ggplot2::layer_scales(tmp_mzPeakwidthHorz)$y$range$range[1], ggplot2::layer_scales(tmp_mzHisto)$x$range$range[1])
    maxYMZ              <- max(ggplot2::layer_scales(tmp_mzPeakwidthHorz)$y$range$range[2], ggplot2::layer_scales(tmp_mzHisto)$x$range$range[2])
    tmp_mzPeakwidthHorz <- tmp_mzPeakwidthHorz + ggplot2::ylim(minYMZ, maxYMZ)
    tmp_mzHisto         <- tmp_mzHisto + ggplot2::xlim(minYMZ, maxYMZ)
    # Area runOrder + histo (is rotated, use x in histo)
    minYArea            <- min(ggplot2::layer_scales(tmp_peakAreaHorz)$y$range$range[1], ggplot2::layer_scales(tmp_areaHisto)$x$range$range[1])
    maxYArea            <- max(ggplot2::layer_scales(tmp_peakAreaHorz)$y$range$range[2], ggplot2::layer_scales(tmp_areaHisto)$x$range$range[2])
    tmp_peakAreaHorz    <- tmp_peakAreaHorz + ggplot2::ylim(minYArea, maxYArea)
    tmp_areaHisto       <- tmp_areaHisto + ggplot2::xlim(minYArea, maxYArea)
    
    ## Convert to gtables + match plotted width
    p_EIC               <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_EIC))
    p_rtPeakwidthVert   <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_rtPeakwidthVert))
    p_rtPeakwidthHorz   <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_rtPeakwidthHorz))
    p_mzPeakwidthHorz   <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_mzPeakwidthHorz))
    p_peakAreaHorz      <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_peakAreaHorz))
    p_rtHisto           <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_rtHisto))
    p_mzHisto           <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_mzHisto))
    p_areaHisto         <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(tmp_areaHisto))
    # find the widths of each of the plots, calculate the maximum and then apply it to each of them individually. This effectively applies a uniform layout to each of the plots.
    # EIC + rt peakwidth
    maxWidthEIC     <- grid::unit.pmax(p_EIC$widths[2:3], p_rtPeakwidthVert$widths[2:3])
    p_EIC$widths[2:3]               <- maxWidthEIC
    p_rtPeakwidthVert$widths[2:3]   <- maxWidthEIC
    # rt peakwidth + mz peakwidth + peak area (x axis)
    maxWidthRtMzArea                <- grid::unit.pmax(p_rtPeakwidthHorz$widths[2:3], p_mzPeakwidthHorz$widths[2:3], p_peakAreaHorz$widths[2:3])
    p_rtPeakwidthHorz$widths[2:3]   <- maxWidthRtMzArea
    p_mzPeakwidthHorz$widths[2:3]   <- maxWidthRtMzArea
    p_peakAreaHorz$widths[2:3]      <- maxWidthRtMzArea
    # rt + mz + area peakwidth and histo (y axis)
    maxHeight       <- grid::unit.pmin(p_rtPeakwidthHorz$heights[2:3],p_mzPeakwidthHorz$heights[2:3],p_peakAreaHorz$heights[2:3],  p_rtHisto$widths[2:3],p_mzHisto$widths[2:3],p_areaHisto$widths[2:3])
    p_rtPeakwidthHorz$heights[2:3]  <- maxHeight
    p_mzPeakwidthHorz$heights[2:3]  <- maxHeight
    p_peakAreaHorz$heights[2:3]     <- maxHeight
    p_rtHisto$widths[2:3]           <- maxHeight
    p_mzHisto$widths[2:3]           <- maxHeight
    p_areaHisto$widths[2:3]         <- maxHeight
    
    ## Generate and store multiplot
    outList[[cpd]] <-gridExtra::grid.arrange(p_EIC,
                                             p_rtPeakwidthVert, 
                                             p_rtPeakwidthHorz, p_rtHisto,
                                             p_mzPeakwidthHorz, p_mzHisto,
                                             p_peakAreaHorz, p_areaHisto, 
                                             widths = c(6, 1),
                                             heights = c(100, 45, 30, 30, 30, 14),
                                             layout_matrix = rbind(c(1, 1),
                                                                   c(2, 2),
                                                                   c(3, 4),
                                                                   c(5, 6),
                                                                   c(7, 8),
                                                                   c(7, 9)),
                                             top = annotationDiagnosticPlotList[[cpd]]$title)
  }
  return(outList)
}