Implementing a ZARR reader in cytomapper

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: cytomapper
-Version: 1.13.0
+Version: 1.13.1
 Title: Visualization of highly multiplexed imaging data in R
 Description: 
     Highly multiplexed imaging acquires the single-cell expression of

NAMESPACE

@@ -14,6 +14,8 @@ export(measureObjects)
 export(mergeChannels)
 export(plotCells)
 export(plotPixels)
+export(readImagesFromZARR)
+export(readMasksFromZARR)
 export(scaleImages)
 exportClasses(CytoImageList)
 exportMethods("[<-")
@@ -63,6 +65,8 @@ importFrom(EBImage,readImage)
 importFrom(HDF5Array,HDF5Array)
 importFrom(HDF5Array,writeHDF5Array)
 importFrom(RColorBrewer,brewer.pal)
+importFrom(Rarr,read_zarr_array)
+importFrom(Rarr,zarr_overview)
 importFrom(S4Vectors,"mcols<-")
 importFrom(S4Vectors,DataFrame)
 importFrom(S4Vectors,SimpleList)
@@ -97,6 +101,7 @@ importFrom(graphics,rect)
 importFrom(graphics,strheight)
 importFrom(graphics,strwidth)
 importFrom(graphics,text)
+importFrom(jsonlite,fromJSON)
 importFrom(matrixStats,rowRanges)
 importFrom(methods,as)
 importFrom(methods,callNextMethod)
@@ -107,6 +112,7 @@ importFrom(raster,as.raster)
 importFrom(rhdf5,h5delete)
 importFrom(rhdf5,h5ls)
 importFrom(stats,quantile)
+importFrom(stringr,str_split)
 importFrom(svgPanZoom,renderSvgPanZoom)
 importFrom(svgPanZoom,svgPanZoom)
 importFrom(svgPanZoom,svgPanZoomOutput)

NEWS

@@ -183,3 +183,6 @@ Changes in version 1.11.2 (2023-02-02):
 
 Changes in version 1.11.3 (2023-01-23):
 + loadImages: added option to read in single-channel images to multi-channel
+
+Changes in version 1.13.1 (2023-09-18):
++ readZARR function to read in images and masks from ZARR archives

R/ZarrReaders.R

@@ -0,0 +1,300 @@
+#' @title Read data from ZARR files
+#' @name readImagesFromZARR
+#'
+#' @description
+#' Reads in data from a ZARR file
+#'
+#' @param file TODO
+#' @param type TODO
+#' @param resolution TODO
+#' @param x TODO
+#' @param y TODO
+#' @param c TODO
+#' @param z TODO
+#' @param t TODO
+#' @param fov_names TODO
+#'
+#' @return A \linkS4class{CytoImageList} object
+#' 
+#' @examples
+#' # TODO
+#'
+#' @author Nils Eling (\email{nils.eling@@dqbm.uzh.ch}),
+#'
+#' @export
+#' @importFrom jsonlite fromJSON
+#' @importFrom Rarr read_zarr_array zarr_overview
+#' @importFrom stringr str_split
+readImagesFromZARR <- function(file,
+                     type = c("omengff", "spatialdata"),
+                     resolution = NULL,
+                     x = NULL,
+                     y = NULL,
+                     c = NULL,
+                     z = 1,
+                     t = 1,
+                     fov_names = NULL) {
+
+    # Check if file is spatialData or OME-NGFF
+    #.valid.readZARR.input(file)
+
+    type <- match.arg(type)
+
+    if (type == "spatialdata") {
+
+        cur_meta <- fromJSON(file.path(file, "zmetadata"))
+        cur_img_meta <- str_split(names(cur_meta$metadata), "/", simplify = TRUE)
+        cur_img_meta <- cur_img_meta[cur_img_meta[,1] == "images" & !grepl("^\\.", cur_img_meta[,2]),]
+        cur_fov_names <- unique(cur_img_meta[,2])
+
+        if (is.null(fov_names)) {
+            fov_names <- cur_fov_names[1]
+        }
+
+        if (is.null(resolution)) {
+            if (length(fov_names) > 1) {
+                resolution <- lapply(file.path(file, "images", fov_names), 
+                                     function(cur_name) {
+                                         cur_res <- fromJSON(file.path(cur_name, ".zattrs"))
+                                         cur_res <- cur_res$multiscales$datasets[[1]]$path
+                                         return(sort(cur_res)[length(cur_res)])
+                                     })
+                resolution <- unlist(resolution)
+            } else {
+                resolution <- fromJSON(file.path(file, "images", fov_names, ".zattrs"))
+                resolution <- resolution$multiscales$datasets[[1]]$path
+                resolution <- sort(resolution)[length(resolution)]
+            }
+        }
+
+        cur_names <- file.path(file, "images", fov_names, resolution)
+
+        # Read in metadata
+        cur_meta <- lapply(file.path(file, "images", fov_names),
+                           function(cur_m) {
+                               return(fromJSON(file.path(cur_m, ".zattrs")))
+                           })
+
+        cur_channels <- lapply(cur_meta, function(cur_m){
+            if ("omero" %in% names(cur_m)) {
+                return(cur_m$omero$channels$label)
+            } else if ("channels_metadata" %in% names(cur_m)) {
+                return(cur_m$channels_metadata$channels$label)
+            } 
+        })
+
+        if (length(unique(cur_channels)) > 1) {
+            stop("Channel names need to match across images.")
+        }
+
+        cur_channels <- unlist(unique(cur_channels))
+
+        cur_index <- lapply(cur_meta, function(cur_m){
+            if ("coordinateTransformations" %in% names(cur_m$multiscales)) {
+                return(cur_m$multiscales$coordinateTransformations[[1]]$input$axes[[1]]$name)
+            } else {
+                return(cur_m$multiscales$axes[[1]]$name)
+            } 
+        })
+
+        if (length(unique(cur_index)) > 1) {
+            stop("Dimensions need to match across images.")
+        }
+
+        cur_index <- unlist(unique(cur_index))
+
+    } else if (type == "omengff") {
+
+        cur_meta <- fromJSON(file.path(file, ".zattrs"))
+
+        if (is.null(resolution)) {
+            resolution <- cur_meta$multiscales$datasets[[1]]$path
+            resolution <- sort(resolution)[length(resolution)]
+        }
+
+        fov_names <- sub("\\.[^.]*$", "", basename(file))
+
+        cur_names <- file.path(file, resolution)
+
+        if ("omero" %in% names(cur_meta)) {
+            cur_channels <- cur_meta$omero$channels$label
+        } else if ("channels_metadata" %in% names(cur_meta)) {
+            cur_channels <- cur_meta$channels_metadata$channels$label
+        } 
+
+        if ("coordinateTransformations" %in% names(cur_meta$multiscales)) {
+            cur_index <- cur_meta$multiscales$coordinateTransformations[[1]]$input$axes[[1]]$name
+        } else {
+            cur_index <- cur_meta$multiscales$axes[[1]]$name
+        } 
+    }
+
+    message("Resolution of the image: ", resolution)
+    message("Channels of the image: ", cur_channels)
+
+    cur_out <- lapply(cur_names, function(cur_name){
+
+        cur_ind <- list(t, c, z, y, x)
+        names(cur_ind) <-  c("t", "c", "z", "y", "x")
+        cur_ind <- cur_ind[cur_index]
+
+        cur_perm <- c("x", "y", "c", "z", "t")
+        cur_perm <- match(cur_perm, cur_index)
+        cur_perm <- cur_perm[!is.na(cur_perm)]
+
+        cur_arr <- read_zarr_array(cur_name, index = cur_ind)
+        cur_arr <- aperm(cur_arr, cur_perm)
+        cur_arr <- Image(array(cur_arr, dim = dim(cur_arr)[seq_len(3)]))
+
+        return(cur_arr)
+    })
+
+    names(cur_out) <- fov_names
+
+    cur_CIL <- CytoImageList(cur_out)
+
+    # Set channelNames
+    if (exists("cur_channels")) {
+        channelNames(cur_CIL) <- cur_channels
+    } else {
+        channelNames(cur_CIL) <- as.character(seq_len(dim(cur_CIL[[1]])[3]))
+    }
+
+    return(cur_CIL)
+}
+
+#' @title Read data from ZARR files
+#' @name readMasksFromZARR
+#'
+#' @description
+#' Reads in data from a ZARR file
+#'
+#' @param file TODO
+#' @param type TODO
+#' @param resolution TODO
+#' @param x TODO
+#' @param y TODO
+#' @param c TODO
+#' @param z TODO
+#' @param t TODO
+#' @param fov_names TODO
+#'
+#' @return A \linkS4class{CytoImageList} object
+#' 
+#' @examples
+#' # TODO
+#'
+#' @author Nils Eling (\email{nils.eling@@dqbm.uzh.ch}),
+#'
+#' @export
+#' @importFrom jsonlite fromJSON
+#' @importFrom Rarr read_zarr_array zarr_overview
+#' @importFrom stringr str_split
+readMasksFromZARR <- function(file,
+                     type = c("omengff", "spatialdata"),
+                     resolution = NULL,
+                     x = NULL,
+                     y = NULL,
+                     c = NULL,
+                     z = 1,
+                     t = 1,
+                     fov_names = NULL) {
+
+    # Check if file is spatialData or OME-NGFF
+    #.valid.readZARR.input(file)
+
+    type <- match.arg(type)
+
+    if (type == "spatialdata") {
+
+      cur_meta <- fromJSON(file.path(file, "zmetadata"))
+      cur_points_meta <- str_split(names(cur_meta$metadata), "/", simplify = TRUE)
+      cur_points_meta <- cur_points_meta[cur_points_meta[,1] == "labels" & !grepl("^\\.", cur_points_meta[,2]),]
+
+      cur_fov_names <- unique(cur_points_meta[,2])
+
+      if (is.null(fov_names)) {
+            fov_names <- cur_fov_names[1]
+        }
+
+        if (is.null(resolution)) {
+            if (length(fov_names) > 1) {
+                resolution <- lapply(file.path(file, "labels", fov_names), 
+                                     function(cur_name) {
+                                       cur_res <- fromJSON(file.path(cur_name, ".zattrs"))
+                                       cur_res <- cur_res$multiscales$datasets[[1]]$path
+                                       return(sort(cur_res)[length(cur_res)])
+                                     })
+
+                resolution <- unlist(resolution)
+                } else {
+                resolution <- fromJSON(file.path(file, "labels", fov_names, ".zattrs"))
+                resolution <- resolution$multiscales$datasets[[1]]$path
+                resolution <- sort(resolution)[length(resolution)]
+            }
+        }
+
+        cur_names <- file.path(file, "labels", fov_names, resolution)
+
+        # Read in metadata
+        cur_meta <- lapply(file.path(file, "labels", fov_names),
+                           function(cur_m) {
+                               return(fromJSON(file.path(cur_m, ".zattrs")))
+                           })
+
+        cur_index <- lapply(cur_meta, function(cur_m){
+          if ("coordinateTransformations" %in% names(cur_m$multiscales)) {
+            return(cur_m$multiscales$coordinateTransformations[[1]]$input$axes[[1]]$name)
+          } else {
+            return(cur_m$multiscales$axes[[1]]$name)
+          } 
+        })
+
+        if (length(unique(cur_index)) > 1) {
+          stop("Dimensions need to match across images.")
+        }
+
+        cur_index <- unlist(unique(cur_index))
+
+    } else if (type == "omengff") {
+
+        if (is.null(resolution)) {
+          cur_res <- fromJSON(file.path(file, "labels/Cell/.zattrs"))
+          cur_res <- cur_res$multiscales$datasets[[1]]$path
+          resolution <- cur_res[length(cur_res)]
+          resolution <- unlist(resolution)
+        }
+
+        fov_names <- sub("\\.[^.]*$", "", basename(file))
+
+        cur_names <- file.path(file, "labels/Cell", resolution)
+
+        cur_meta <- fromJSON(file.path(file, "labels/Cell/.zattrs"))
+
+        cur_index <- if ("coordinateTransformations" %in% names(cur_meta$multiscales)) {
+            cur_meta$multiscales$coordinateTransformations[[1]]$input$axes[[1]]$name
+          } else {
+            cur_meta$multiscales$axes[[1]]$name
+          } 
+    }
+
+    message("Resolution of the image: ", resolution)
+
+    cur_out <- lapply(cur_names, function(cur_name){
+        cur_ind <- list(t, c, z, y, x)
+        names(cur_ind) <-  c("t", "c", "z", "y", "x")
+        cur_ind <- cur_ind[cur_index]
+
+        cur_arr <- read_zarr_array(cur_name, index = cur_ind)
+
+        cur_arr <- Image(aperm(cur_arr, rev(seq_along(dim(cur_arr)))))
+
+        return(cur_arr)
+    })
+
+    names(cur_out) <- fov_names
+
+    cur_CIL <- CytoImageList(cur_out)
+
+    return(cur_CIL)
+}