Merge pull request #4 from HelenaLC/multichannel

PR #2

Author: HelenaLC

.Rbuildignore

@@ -1 +1,3 @@
 ^\.github$
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+*.Rproj

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SpatialData.plot
 Title: SpatialData visualization
 Depends: R (>= 4.4), SpatialData
-Version: 0.99.1
+Version: 0.99.2
 Description: Visualization suit for 'SpatialData' (R). Current functionality 
     includes handling of multiscale 'images', visualizing 'labels', 'points', 
     and 'shapes'. For the latter, POINT, POLYGON, and MULTIPOLYGON geometries 
@@ -41,7 +41,8 @@ Imports:
     rlang,
     sf,
     S4Vectors,
-    SingleCellExperiment
+    SingleCellExperiment,
+    Rarr
 Suggests:
     BiocStyle,
     ggnewscale,

NAMESPACE

@@ -10,9 +10,12 @@ exportMethods(scale)
 exportMethods(translation)
 import(SpatialData)
 importFrom(DelayedArray,realize)
+importFrom(Rarr,zarr_overview)
 importFrom(S4Vectors,metadata)
 importFrom(SingleCellExperiment,int_colData)
 importFrom(SingleCellExperiment,int_metadata)
+importFrom(SpatialData,channels)
+importFrom(SpatialData,getZarrArrayPath)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(ggforce,geom_circle)

R/plotImage.R

@@ -10,6 +10,11 @@
 #' @param k index of the scale of an image; by default (NULL), will auto-select 
 #'   scale in order to minimize memory-usage and blurring for a target size of 
 #'   800 x 800px; use Inf to plot the lowest resolution available.
+#' @param ch image channel(s) to be used for plotting (defaults to 
+#'   the first channel(s) available); use \code{channels()} to see 
+#'   which channels are available for a given \code{ImageArray} 
+#'   
+#' @param c plotting aesthetics; color 
 #'
 #' @return ggplot
 #'
@@ -31,6 +36,86 @@ NULL
 #' @export
 plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme 
 
+# merge/manage image channels
+# if no colors and channels defined, return the first channel
+#' @importFrom grDevices col2rgb
+#' @noRd
+.manage_channels <- \(a, ch, c=NULL){
+    if (length(ch) > (n <- length(.DEFAULT_COLORS)) && is.null(c))
+        stop("Only ", n, " default colors available, but",
+            length(ch), " are needed; please specify 'c'")
+    if (!is.null(c) || (is.null(c) && length(ch) > 1)) {
+        if (is.null(c)) c <- .DEFAULT_COLORS[seq_along(ch)] 
+        c <- col2rgb(c)/255
+        b <- array(0, dim=c(3, dim(a)[-1]))
+        for (i in seq_len(dim(a)[1])) {
+            b[1,,] <- b[1,,,drop=FALSE] + a[i,,,drop=FALSE]*c[1,i]
+            b[2,,] <- b[2,,,drop=FALSE] + a[i,,,drop=FALSE]*c[2,i]
+            b[3,,] <- b[3,,,drop=FALSE] + a[i,,,drop=FALSE]*c[3,i]
+        }
+        a <- pmin(b, 1)
+    } else {
+        a <- a[rep(1, 3), , ]
+    }
+    return(a)
+}
+
+# check if an image is rgb or not
+#' @importFrom SpatialData getZarrArrayPath
+#' @importFrom Rarr zarr_overview
+#' @noRd
+.get_image_dtype <- \(a) {
+    pa <- getZarrArrayPath(a)
+    df <- zarr_overview(pa, as_data_frame=TRUE)
+    if (!is.null(dt <- df$data_type)) return(dt)
+}
+
+# normalize the image data given its data type
+#' @noRd
+.normalize_image_array <- \(a, dt){
+  if (dt %in% names(.DTYPE_MAX_VALUES)) {
+      a <- a/.DTYPE_MAX_VALUES[dt]
+  } else if (max(a) > 1) {
+      for (i in seq_len(dim(a)[1])) 
+          a[i,,] <- a[i,,]/max(a[i,,])
+  }
+  return(a)
+}
+
+# check if an image is RGB or not
+# (NOTE: some RGB channels are named 0, 1, 2)
+#' @importFrom methods is
+#' @noRd
+.is_rgb <- \(x) {
+    if (is(x, "ImageArray") &&
+        !is.null(md <- meta(x)))
+        x <- md$omero$channels$label
+    if (!is.vector(x)) stop("invalid 'x'")
+    is_len <- length(x) == 3
+    is_012 <- setequal(x, seq(0, 2))
+    is_rgb <- setequal(x, c("r", "g", "b"))
+    return(is_len && (is_012 || is_rgb))
+}
+  
+# check if channels are indices or channel names
+#' @importFrom SpatialData channels
+#' @noRd
+.ch_idx <- \(x, ch) {
+    if (is.null(ch))
+        return(1)
+    lbs <- channels(x)
+    if (all(ch %in% lbs)) {
+        return(match(ch, lbs))
+    } else if (!any(ch %in% lbs)) {
+        warning("Couldn't find some channels; picking first one(s)!")
+        return(1)
+    } else {
+        warning("Couldn't find channels; picking first one(s)!")
+        return(1)
+    }
+    return(NULL)
+}
+
 .guess_scale <- \(x, w, h) {
     n <- length(dim(x))
     i <- ifelse(n == 3, -1, TRUE)
@@ -47,14 +132,17 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 #' @importFrom methods as
 #' @importFrom grDevices rgb
 #' @importFrom DelayedArray realize
-.df_i <- \(x, k=NULL) {
+.df_i <- \(x, k=NULL, ch=NULL, c=NULL) {
     a <- .get_plot_data(x, k)
-    a <- if (dim(a)[1] == 1) a[rep(1,3),,] else a
+    ch_i <- .ch_idx(x, ch)
+    if (!.is_rgb(x))
+        a <- a[ch_i, , , drop=FALSE]
+    dt <- .get_image_dtype(a)
     a <- realize(as(a, "DelayedArray"))
-    img <- rgb(
-        maxColorValue=max(a),
-        c(a[1,,]), c(a[2,,]), c(a[3,,]))
-    array(img, dim(a)[-1])
+    a <- .normalize_image_array(a, dt)
+    if (!.is_rgb(x))
+        a <- .manage_channels(a, ch_i, c)
+    apply(a, c(2, 3), \(.) do.call(rgb, as.list(.))) 
 }
 
 .get_wh <- \(x, i, j) {
@@ -75,13 +163,13 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 
 #' @rdname plotImage
 #' @export
-setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL) {
-    if (is.numeric(i)) 
+setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=NULL) {
+    if (is.numeric(i))
         i <- imageNames(x)[i]
     y <- image(x, i)
-    if (is.numeric(j)) 
+    if (is.numeric(j))
         j <- CTname(y)[j]
-    df <- .df_i(y, k)
+    df <- .df_i(y, k, ch, c)
     wh <- .get_wh(x, i, j)
     .gg_i(df, wh$w, wh$h)
 })

R/utils.R

@@ -33,4 +33,15 @@
         plot.title=element_text(hjust=0.5),
         axis.text=element_text(color="grey"),
         axis.ticks=element_line(color="grey"))
-)
+)
+
+# default colors (from ImageJ/Fiji)
+.DEFAULT_COLORS <- c("red", "green", "blue", "gray", "cyan", "magenta", "yellow")
+
+# image data type factors (max values)
+# TODO: add more cases from other data types
+# https://doc.embedded-wizard.de/uint-type
+.DTYPE_MAX_VALUES <- c("uint8" = 255,
+                       "uint16" = 65535,
+                       "uint32" = 4294967295,
+                       "uint64" = 2^64 - 1)

inst/NEWS

@@ -1,3 +1,8 @@
+changes in version 0.99.2
+
+- in 'plotImage', added support to visualize channels of choice
+- updated vignette to include corresponding examples
+
 changes in version 0.99.1
 
 - various fixes related to moving 'instance/region_key' to 'int_colData'

man/plotImage.Rd

@@ -0,0 +1,30 @@
+require(SpatialData, quietly=TRUE)
+x <- file.path("extdata", "blobs.zarr")
+x <- system.file(x, package="SpatialData")
+x <- readSpatialData(x, tables=FALSE)
+
+test_that(".is_rgb()", {
+    # valid integer vector
+    expect_false(.is_rgb(c(0, 1, 1)))
+    expect_true(.is_rgb(. <- seq(0, 2)))
+    expect_true(.is_rgb(rev(.)))
+    # valid character vector
+    expect_false(.is_rgb(c("r", "g", "g")))
+    expect_true(.is_rgb(. <- c("r", "g", "b")))
+    expect_true(.is_rgb(rev(.)))
+    # only works for 'ImageArray'
+    expect_true(.is_rgb(image(x, 1)))
+    expect_error(.is_rgb(label(x, 1)))
+})
+
+test_that(".ch_idx()", {
+  # get indices of channels
+  expect_equal(.ch_idx(image(x,1), ch=c(2,0,1)), c(3,1,2))
+  # return first if no matching channel
+  expect_warning(expect_equal(.ch_idx(image(x,1), ch=99), 1)) 
+})
+
+# TODO: any tests for image array normalization ? 
+test_that(".normalize_image_array", {
+    skip()
+})

tests/testthat/test-plotImage.R

@@ -238,6 +238,53 @@ wrap_plots(nrow=1, lapply(seq(3), \(.)
     plot_layout(guides="collect")
 ```
 
+## CyCIF (MCMICRO)
+
+Small lung adenocarcinoma, 250 MB; 1 image, 2 labels, 2 tables.
+
+```{r mcmicro-read}
+dir.create(td <- tempfile())
+pa <- unzip_spd_demo(
+    zipname="mcmicro_io.zip", 
+    dest=td, source="biocOSN")
+(x <- readSpatialData(pa, anndataR=FALSE))
+```
+
+Getting channel names for the image:
+
+```{r mcmicro-channels}
+channels(image(x))
+```
+
+Plotting with multiple image channels:
+
+```{r mcmicro-plot}
+plotSpatialData() + plotImage(x,
+    ch=c("DNA_6", "CD45", "CD57"), 
+    c=c("blue", "cyan", "yellow"))
+```
+
+## IMC (Steinbock)
+
+4 different cancers (SCCHN, BCC, NSCLC, CRC), 820 MB; 14 images, 14 labels, 1 table.
+
+```{r steinbock-read}
+dir.create(td <- tempfile())
+pa <- unzip_spd_demo(
+    zipname="steinbock_io.zip", 
+    dest=td, source="biocOSN")
+x <- readSpatialData(pa, anndataR=FALSE)
+```
+
+Plotting with multiple image channels.
+
+```{r steinbock-plot}
+plotSpatialData() + plotImage(x, 
+    i="Patient3_003_image",
+    ch=c(6, 22, 39),
+    c=c("blue", "cyan", "yellow"))
+```
+
 # Masking
 
 Back to blobs...