Reimplement plotTIS()/plotRIC() on the Set classes; add single-sample methods

R/align-GCIMSSample.R

@@ -102,9 +102,8 @@ methods::setMethod(
 #' @return The modified [GCIMSSample]
 #' @export
 alignDt <- function(object, rip_ref_ms) {
-  tis <- getTIS(object)
   dt <- dtime(object)
-  rip_pos_ms <- dt[which.max(tis)]
+  rip_pos_ms <- dt[rip_position(object)]
   Kcorr <- rip_ref_ms/rip_pos_ms
 
   int_mat <- intensity(object)

R/getTIS_getRIC-GCIMSDataset.R

@@ -35,12 +35,13 @@ setMethod("getRIC", "GCIMSDataset", function(object) {
 #' @param object A [GCIMSDataset] object
 #' @inheritParams dt_rt_range_normalization
 #' @param sample A number or a string with the sample index or name. If `NULL`, all samples are returned
+#' @param color_by The name of a `pData(object)` column (or `"SampleID"`) used to color the spectra
 #' @return The plot of the TIS
 #' @export
 setMethod(
   "plotTIS",
   "GCIMSDataset",
-  function(object, dt_range = NULL, sample = NULL) {
+  function(object, dt_range = NULL, sample = NULL, color_by = "SampleID") {
     tis <- getTIS(object)
     dt <- dtime(object)
     sample_names <- sampleNames(object)
@@ -49,22 +50,20 @@ setMethod(
     }
     sample_idx <- sample_name_or_number_to_both(sample, sample_names)
     idx <- dt_rt_range_normalization(dt = dt, dt_range = dt_range)
-    tis_long <- reshape2::melt(tis[sample_idx$idx, idx$dt_logical, drop = FALSE], value.name = "TIS")
-    gplt <- ggplot2::ggplot() +
-      ggplot2::geom_line(
-        data = tis_long,
-        mapping = ggplot2::aes(
-          x = .data$drift_time_ms,
-          y = .data$TIS,
-          color = .data$SampleID
-        )
-      ) +
-      ggplot2::labs(
-        x = "Drift time (ms)",
-        y = "TIS Intensity (a.u.)",
-        color = "SampleID"
-      )
-    gplt
+    tis_sub <- tis[sample_idx$idx, idx$dt_logical, drop = FALSE]
+    dt_sub <- dt[idx$dt_logical]
+
+    pd <- pData(object)
+    spectra <- stats::setNames(
+      purrr::map(rownames(tis_sub), function(sample_id) {
+        GCIMSSpectrum(drift_time = dt_sub, intensity = unname(tis_sub[sample_id, ]), description = sample_id)
+      }),
+      rownames(tis_sub)
+    )
+    spectrum_set <- GCIMSSpectrumSet(spectra = spectra, pData = pd[pd$SampleID %in% rownames(tis_sub), , drop = FALSE])
+
+    plot(spectrum_set, color_by = color_by) +
+      ggplot2::labs(y = "TIS Intensity (a.u.)")
   }
 )
 
@@ -73,12 +72,13 @@ setMethod(
 #' @param object A [GCIMSDataset] object
 #' @inheritParams dt_rt_range_normalization
 #' @param sample A number or a string with the sample index or name. If `NULL`, all samples are returned
+#' @param color_by The name of a `pData(object)` column (or `"SampleID"`) used to color the chromatograms
 #' @return A plot
 #' @export
 setMethod(
   "plotRIC",
   "GCIMSDataset",
-  function(object, rt_range = NULL, sample = NULL) {
+  function(object, rt_range = NULL, sample = NULL, color_by = "SampleID") {
     ric <- getRIC(object)
     rt <- rtime(object)
     sample_names <- sampleNames(object)
@@ -87,23 +87,20 @@ setMethod(
     }
     sample_idx <- sample_name_or_number_to_both(sample, sample_names)
     idx <- dt_rt_range_normalization(rt = rt, rt_range = rt_range)
-    ric_long <- reshape2::melt(ric[sample_idx$idx, idx$rt_logical, drop = FALSE], value.name = "RIC")
-
-    gplt <- ggplot2::ggplot() +
-      ggplot2::geom_line(
-        data = ric_long,
-        mapping = ggplot2::aes(
-          x = .data$retention_time_s,
-          y = .data$RIC,
-          color = .data$SampleID
-        )
-      ) +
-      ggplot2::labs(
-        x = "Retention time (s)",
-        y = "RIC Intensity (a.u.)",
-        color = "SampleID"
-      )
-    gplt
+    ric_sub <- ric[sample_idx$idx, idx$rt_logical, drop = FALSE]
+    rt_sub <- rt[idx$rt_logical]
+
+    pd <- pData(object)
+    chromatograms <- stats::setNames(
+      purrr::map(rownames(ric_sub), function(sample_id) {
+        GCIMSChromatogram(retention_time = rt_sub, intensity = unname(ric_sub[sample_id, ]), description = sample_id)
+      }),
+      rownames(ric_sub)
+    )
+    chromatogram_set <- GCIMSChromatogramSet(chromatograms = chromatograms, pData = pd[pd$SampleID %in% rownames(ric_sub), , drop = FALSE])
+
+    plot(chromatogram_set, color_by = color_by) +
+      ggplot2::labs(y = "RIC Intensity (a.u.)")
   }
 )
 

R/getTIS_getRIC-GCIMSSample.R

@@ -28,12 +28,63 @@ setMethod("getTIS", "GCIMSSample", function(object) {
 #' ric <- getRIC(s)
 setMethod("getRIC", "GCIMSSample", function(object) {
   intmat <- intensity(object)
-  tis <- rowSums(intmat)
-  ric_pos <- which.max(tis)
+  ric_pos <- rip_position(object)
   ric <- intmat[ric_pos, ]
   ric <- max(ric) - ric
   ric <- ric/sum(ric)
   ric
 })
 
+#' Drift time index of the Reactant Ion Peak (RIP)
+#'
+#' The RIP is the drift time with the highest total ion signal, summed
+#' across all retention times (the `which.max()` of [getTIS()]).
+#'
+#' @param object A [GCIMSSample] object
+#' @return An integer, the drift time index of the RIP
+#' @noRd
+rip_position <- function(object) {
+  which.max(rowSums(intensity(object)))
+}
+
+#' Plot the total ion spectrum of a sample
+#' @param object A [GCIMSSample] object
+#' @inheritParams dt_rt_range_normalization
+#' @return A ggplot2 plot object
+#' @export
+setMethod("plotTIS", "GCIMSSample", function(object, dt_range = NULL) {
+  dt <- dtime(object)
+  rt <- rtime(object)
+  idx <- dt_rt_range_normalization(dt, rt, dt_range = dt_range)
+  tis <- getTIS(object)
+  spec <- GCIMSSpectrum(
+    drift_time = dt[idx[["dt_logical"]]],
+    intensity = tis[idx[["dt_logical"]]],
+    retention_time_idx = unique(c(idx[["rt_idx_min"]], idx[["rt_idx_max"]])),
+    retention_time_s = unique(c(idx[["rt_s_min"]], idx[["rt_s_max"]])),
+    description = object@description
+  )
+  plot(spec) + ggplot2::labs(y = "TIS Intensity (a.u.)")
+})
+
+#' Plot the reverse ion chromatogram of a sample
+#' @param object A [GCIMSSample] object
+#' @inheritParams dt_rt_range_normalization
+#' @return A ggplot2 plot object
+#' @export
+setMethod("plotRIC", "GCIMSSample", function(object, rt_range = NULL) {
+  dt <- dtime(object)
+  rt <- rtime(object)
+  ric <- getRIC(object)
+  ric_pos <- rip_position(object)
+  idx <- dt_rt_range_normalization(dt, rt, dt_idx = ric_pos, rt_range = rt_range)
+  chrom <- GCIMSChromatogram(
+    retention_time = rt[idx[["rt_logical"]]],
+    intensity = ric[idx[["rt_logical"]]],
+    drift_time_idx = unique(c(idx[["dt_idx_min"]], idx[["dt_idx_max"]])),
+    drift_time_ms = unique(c(idx[["dt_ms_min"]], idx[["dt_ms_max"]])),
+    description = object@description
+  )
+  plot(chrom) + ggplot2::labs(y = "RIC Intensity (a.u.)")
+})
 

tests/testthat/test-getTIS_getRIC-GCIMSDataset.R

@@ -10,7 +10,8 @@ make_dataset <- function() {
   }
   s1 <- make_s(2, 25)
   s2 <- make_s(2.2, 25)
-  GCIMSDataset$new_from_list(samples = list(s1 = s1, s2 = s2), on_ram = TRUE, scratch_dir = NULL)
+  pd <- data.frame(SampleID = c("s1", "s2"), Group = c("A", "B"))
+  GCIMSDataset$new_from_list(samples = list(s1 = s1, s2 = s2), pData = pd, on_ram = TRUE, scratch_dir = NULL)
 }
 
 test_that("getTIS()/getRIC() extract per-sample TIS and RIC matrices, keyed by SampleID", {
@@ -39,17 +40,17 @@ test_that("plotTIS() returns a ggplot with one line per sample over the full dri
   expect_s3_class(p, "ggplot")
   expect_equal(p$labels$x, "Drift time (ms)")
   expect_equal(p$labels$y, "TIS Intensity (a.u.)")
-  expect_setequal(as.character(unique(p$layers[[1]]$data$SampleID)), c("s1", "s2"))
-  expect_equal(range(p$layers[[1]]$data$drift_time_ms), c(0, 4))
+  expect_setequal(as.character(unique(p$data$SampleID)), c("s1", "s2"))
+  expect_equal(range(p$data$drift_time_ms), c(0, 4))
 })
 
 test_that("plotTIS() restricts to dt_range and a single sample when given", {
   ds <- make_dataset()
 
   p <- plotTIS(ds, dt_range = c(1, 3), sample = "s1")
 
-  expect_equal(as.character(unique(p$layers[[1]]$data$SampleID)), "s1")
-  expect_equal(range(p$layers[[1]]$data$drift_time_ms), c(1, 3))
+  expect_equal(as.character(unique(p$data$SampleID)), "s1")
+  expect_equal(range(p$data$drift_time_ms), c(1, 3))
 })
 
 test_that("plotRIC() returns a ggplot with one line per sample over the full retention time range by default", {
@@ -60,14 +61,24 @@ test_that("plotRIC() returns a ggplot with one line per sample over the full ret
   expect_s3_class(p, "ggplot")
   expect_equal(p$labels$x, "Retention time (s)")
   expect_equal(p$labels$y, "RIC Intensity (a.u.)")
-  expect_setequal(as.character(unique(p$layers[[1]]$data$SampleID)), c("s1", "s2"))
+  expect_setequal(as.character(unique(p$data$SampleID)), c("s1", "s2"))
 })
 
 test_that("plotRIC() restricts to rt_range and a single sample when given", {
   ds <- make_dataset()
 
   p <- plotRIC(ds, rt_range = c(10, 30), sample = 1)
 
-  expect_equal(as.character(unique(p$layers[[1]]$data$SampleID)), "s1")
-  expect_equal(range(p$layers[[1]]$data$retention_time_s), c(10, 30))
+  expect_equal(as.character(unique(p$data$SampleID)), "s1")
+  expect_equal(range(p$data$retention_time_s), c(10, 30))
+})
+
+test_that("plotTIS()/plotRIC() can color by a pData column instead of SampleID", {
+  ds <- make_dataset()
+
+  p_tis <- plotTIS(ds, color_by = "Group")
+  p_ric <- plotRIC(ds, color_by = "Group")
+
+  expect_equal(sort(unique(p_tis$data$Group)), c("A", "B"))
+  expect_equal(sort(unique(p_ric$data$Group)), c("A", "B"))
 })

tests/testthat/test-getTIS_getRIC-GCIMSSample.R

@@ -0,0 +1,55 @@
+make_rip_sample <- function() {
+  dt <- seq(1, 10, by = 0.5) # 19 pts, dt[10] == 5.5
+  rt <- seq(1, 60, by = 1)
+  data <- matrix(1, nrow = length(dt), ncol = length(rt))
+  data[10, ] <- 1000 # unambiguous RIP row at dt = 5.5
+  GCIMSSample(drift_time = dt, retention_time = rt, data = data, description = "s1")
+}
+
+test_that("plotTIS(GCIMSSample) plots the full spectrum by default, with an accurate subtitle", {
+  s <- make_rip_sample()
+
+  p <- plotTIS(s)
+
+  expect_s3_class(p, "ggplot")
+  expect_equal(p$labels$x, "Drift time (ms)")
+  expect_equal(p$labels$y, "TIS Intensity (a.u.)")
+  expect_equal(p$labels$title, "s1")
+  expect_match(p$labels$subtitle, "Retention time 1 - 60 s")
+  expect_equal(unname(p$data$y), unname(getTIS(s)))
+})
+
+test_that("plotTIS(GCIMSSample) restricts the displayed drift time range without changing values", {
+  s <- make_rip_sample()
+
+  p_full <- plotTIS(s)
+  p_sub <- plotTIS(s, dt_range = c(2, 8))
+
+  expect_equal(range(p_sub$data$x), c(2, 8))
+  # Values within the window are identical to the unrestricted plot's values:
+  common <- merge(p_full$data, p_sub$data, by = "x")
+  expect_equal(common$y.x, common$y.y)
+})
+
+test_that("plotRIC(GCIMSSample) has an accurate single-value drift time subtitle (the RIP position)", {
+  s <- make_rip_sample()
+
+  p <- plotRIC(s)
+
+  expect_s3_class(p, "ggplot")
+  expect_equal(p$labels$x, "Retention time (s)")
+  expect_equal(p$labels$y, "RIC Intensity (a.u.)")
+  expect_equal(p$labels$title, "s1")
+  expect_equal(p$labels$subtitle, "Drift time 5.5 ms")
+})
+
+test_that("plotRIC(GCIMSSample) restricts the displayed retention time range without changing values", {
+  s <- make_rip_sample()
+
+  p_full <- plotRIC(s)
+  p_sub <- plotRIC(s, rt_range = c(10, 30))
+
+  expect_equal(range(p_sub$data$x), c(10, 30))
+  common <- merge(p_full$data, p_sub$data, by = "x")
+  expect_equal(common$y.x, common$y.y)
+})