refactoring

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: spatialdataR
 Title: Representation of Python's spatialdata in R
 Depends: R (>= 4.6)
-Version: 0.99.42
+Version: 0.99.43
 Description: R interface to Python/scverse's 'spatialdata' framework for 
     unified spatial omics data handling. Adheres to OME-NGFF standards, 
     providing lazy, on-disk representations for multiscale images and 

NAMESPACE

@@ -150,6 +150,7 @@ importFrom(SummarizedExperiment,"assay<-")
 importFrom(SummarizedExperiment,"assayNames<-")
 importFrom(SummarizedExperiment,"colData<-")
 importFrom(SummarizedExperiment,assay)
+importFrom(SummarizedExperiment,assayNames)
 importFrom(SummarizedExperiment,colData)
 importFrom(ZarrArray,ZarrArray)
 importFrom(ZarrArray,path)
@@ -194,6 +195,7 @@ importFrom(methods,"slot<-")
 importFrom(methods,as)
 importFrom(methods,callNextMethod)
 importFrom(methods,is)
+importFrom(methods,new)
 importFrom(methods,setClass)
 importFrom(methods,setClassUnion)
 importFrom(methods,setMethod)

R/AllClasses.R

@@ -1,5 +1,5 @@
 #' @importFrom S4Vectors SimpleList
-#' @importFrom methods setClass setClassUnion setOldClass
+#' @importFrom methods new setClass setClassUnion setOldClass
 
 .sdLayerList <- setClass(
     Class="sdLayerList", 

R/CTutils.R

@@ -74,7 +74,7 @@ setMethod("CTlist", "SpatialDataAttrs", \(x, ...) {
 #' @rdname CTutils
 #' @export
 setMethod("CTdata", "SpatialDataAttrs", \(x, i=1, ...) {
-    i <- .resolve_id(i, CTname(x))
+    i <- .val_id(i, CTname(x))
     t <- CTtype(x)[i]
     if (t != "sequence") 
         return(CTlist(x)[[i]][[t]])

R/crop.R

@@ -107,7 +107,7 @@ NULL
         y=c(y$ymin, y$ymin, y$ymax, y$ymax, y$ymin),
         id=seq_len(5))
     # get transformation for space 'j'
-    j <- .resolve_id(j, CTname(x))
+    j <- .val_id(j, CTname(x))
     ct <- CTlist(x)[[j]]
     # helper to adapt transformation data to spatial (XY) dims
     axs <- axes(x)
@@ -253,35 +253,8 @@ setMethod("crop", "SpatialData", \(x, y, j=1, ...) {
     z <- .lapplyElement(z, \(.) if (length(.) > 0) .)
     z <- do.call("SpatialData", z)
     tables(z) <- tables(x)
-    # filter tables for remaining region(s)/instance(s)
-    rs <- unlist(colnames(z))
-    ts <- lapply(tables(z), \(t) {
-        # filter for remaining element(s)
-        t <- t[, regions(t) %in% rs]
-        region(t) <- intersect(region(t), rs)
-        # table's regions-instances
-        df <- data.frame(
-            r=regions(t), 
-            i=instances(t),
-            keep=seq_len(ncol(t)))
-        # for each annotated element
-        rs <- intersect(region(t), unlist(colnames(z)))
-        is <- lapply(rs, \(r) {
-            # subset look-up
-            df <- df[df$r == r, ]
-            e <- element(z, r)
-            if (is(e, "SpatialDataShape")) {
-                # element's regions-instances
-                ik <- instance_key(t)
-                i <- if (ik %in% names(e)) e[[ik]] else seq_along(e)
-                fd <- data.frame(r, i)
-                # return table indices in element
-                right_join(df, fd, names(fd))$keep
-            } else df$keep
-        })
-        # subset table instances
-        t <- t[, unlist(is)]
-    })
-    tables(z) <- ts
+    # filter table instances
+    z <- .sync_tables_on_crop(z)
     return(z)
 })
+

R/mask.R

@@ -56,7 +56,7 @@ setMethod("mask", c("SpatialData", "ANY", "ANY"), \(x, i, j, k,
     if (!length(ct)) stop(
         "can't mask; found no common ",
         "coordinates between 'i' and 'j'")
-    k <- if (missing(k)) 1 else .resolve_id(k, ct)
+    k <- if (missing(k)) 1 else .val_id(k, ct)
     .i <- transform(.i, ct[k])
     .j <- transform(.j, ct[k])
     t <- tryCatch(error=\(.) NULL, getTable(x, i))
@@ -79,8 +79,8 @@ setGeneric("mask_i_by_j", \(i, j, ...) standardGeneric("mask_i_by_j"))
 #' @noRd
 #' @importFrom methods as
 #' @importFrom Matrix sparseVector
-#' @importFrom SummarizedExperiment assayNames<-
 #' @importFrom SingleCellExperiment SingleCellExperiment
+#' @importFrom SummarizedExperiment assayNames assayNames<-
 setMethod("mask_i_by_j", 
     c("SpatialDataImage", "SpatialDataLabel"), 
     \(i, j, how=NULL, ...) {

R/read.R

@@ -142,7 +142,7 @@ readSpatialData <- function(x,
         opt <- args[[l]]
         if (!isTRUE(opt)) {
             # validate each requested element
-            j <- j[vapply(opt, .resolve_id, integer(1), ok=nms, nm=l)]
+            j <- j[vapply(opt, .val_id, integer(1), ok=nms, nm=l)]
         }
         f <- paste0("read", toupper(substr(l, 1, 1)), substr(l, 2, nchar(l)-1))
         lapply(j, f)

R/trans.R

@@ -52,7 +52,7 @@ NULL
 #' @rdname trans
 #' @importFrom BiocGenerics transform
 setMethod("transform", "SpatialDataElement", \(x, i=1, ...) {
-    i <- .resolve_id(i, CTname(x))
+    i <- .val_id(i, CTname(x))
     f <- CTtype(x)[i]
     t <- CTdata(x, i)
     if (f == "sequence") {

R/utils.R

@@ -1,37 +1,8 @@
-# get/make DuckDB connection
-#' @importFrom DBI dbIsValid
-#' @importFrom duckspatial ddbs_create_conn
-.conn <- \() {
-    nm <- ".SpatialData_DuckDB_conn"
-    if (!exists(nm, envir=.GlobalEnv) ||
-            !dbIsValid(.GlobalEnv[[nm]])) {
-        .GlobalEnv[[nm]] <- ddbs_create_conn()
-    }
-    .GlobalEnv[[nm]]
-}
-
-# internal helper to resolve name/index to integer index
-.resolve_id <- \(i, ok, nm=deparse1(substitute(i))) {
-    nm <- sprintf("'%s'", nm)
-    if (is.character(i)) {
-        i <- match.arg(i, ok)
-        return(match(i, ok))
-    }
-    if (is.numeric(i) && i == round(i) && length(i) == 1) {
-        if (i < 1 || i > length(ok)) {
-            stop(sprintf("invalid %s index: %d (max: %d)", nm, i, length(ok)))
-        }
-        return(as.integer(i))
-    }
-    stop(sprintf("invalid %s; expected character or integer index", nm))
-}
-
 # internal helper for null-coalescing
 `%||%` <- \(a, b) if (is.null(a)) b else a
 
 # internal helpers for object-wide iteration 
 # across spatial elements (excluding tables)
-
 .ls <- .LAYERS[.LAYERS != "tables"]
 
 .lapplyLayer <- \(x, FUN, ...) {
@@ -47,6 +18,37 @@
     return(x)
 }
 
+# get/make DuckDB connection
+#' @importFrom DBI dbIsValid
+#' @importFrom duckspatial ddbs_create_conn
+.conn <- \() {
+    nm <- ".SpatialData_DuckDB_conn"
+    if (!exists(nm, envir=.GlobalEnv) ||
+            !dbIsValid(.GlobalEnv[[nm]])) {
+        .GlobalEnv[[nm]] <- ddbs_create_conn()
+    }
+    .GlobalEnv[[nm]]
+}
+
+# tables ----
+
+.sync_shapes_on_drop <- \(x, i) {
+    # skip when there aren't any shapes
+    if (!length(shapes(x))) return(x)
+    t <- table(x, i)
+    for (j in region(t)) {
+        # skip non-shape elements
+        if (layer(x, j) != "shapes") next
+        # get element 'y' annotated by table 't'
+        y <- element(x, j)
+        # match instances between them
+        y <- y[match(instances(t), instances(y), nomatch=0)] 
+        # return matching shape instances
+        shape(x, j) <- y
+    }
+    return(x)
+}
+
 #' @importFrom methods slot<-
 .sync_tables_sdattrs <- \(x, old, new) {
     if (!length(ts <- tables(x))) return(x)
@@ -73,23 +75,6 @@
     return(x)
 }
 
-.sync_shapes_on_drop <- \(x, i) {
-    # skip when there aren't any shapes
-    if (!length(shapes(x))) return(x)
-    t <- table(x, i)
-    for (j in region(t)) {
-        # skip non-shape elements
-        if (layer(x, j) != "shapes") next
-        # get element 'y' annotated by table 't'
-        y <- element(x, j)
-        # match instances between them
-        y <- y[match(instances(t), instances(y), nomatch=0)] 
-        # return matching shape instances
-        shape(x, j) <- y
-    }
-    return(x)
-}
-
 #' @importFrom methods slot<-
 .sync_tables_on_drop <- \(x) {
     if (!length(ts <- tables(x))) return(x)
@@ -120,22 +105,68 @@
     return(x)
 }
 
+#' @importFrom dplyr right_join
+.sync_tables_on_crop <- \(x) {
+    # filter tables for remaining region(s)/instance(s)
+    rs <- unlist(colnames(x))
+    ts <- lapply(tables(x), \(t) {
+        # filter for remaining element(s)
+        t <- t[, regions(t) %in% rs]
+        region(t) <- intersect(region(t), rs)
+        # table's regions-instances
+        df <- data.frame(
+            r=regions(t), 
+            i=instances(t),
+            keep=seq_len(ncol(t)))
+        # for each annotated element
+        rs <- intersect(region(t), unlist(colnames(x)))
+        is <- lapply(rs, \(r) {
+            # subset look-up
+            e <- element(x, r)
+            df <- df[df$r == r, ]
+            # keep all for labels
+            lb <- is(e, "SpatialDataLabel")
+            if (lb) return(df$keep)
+            # element's regions-instances
+            ik <- instance_key(t)
+            i <- if (ik %in% names(e)) e[[ik]] else seq_along(e)
+            fd <- data.frame(r, i)
+            # return table indices in element
+            right_join(df, fd, names(fd))$keep
+        })
+        # subset table instances
+        t <- t[, unlist(is)]
+    })
+    tables(x) <- ts
+    return(x)
+}
+
 # internal helper to resolve spatial (XY) axis indices
 .get_xy_axes <- \(x) {
     nm <- axes(x, "name")
     ix <- match("x", nm)
     iy <- match("y", nm)
-    # fallback: OME-NGFF usually places spatial dimensions at the end (YX)
-    if (is.na(ix) || is.na(iy)) {
-        n <- length(nm)
-        ix <- n
-        iy <- n-1
-    }
     return(list(x=ix, y=iy))
 }
 
 # validation ----
 
+# internal helper to verify & resolve name/index to index
+.val_id <- \(i, ok, nm=deparse1(substitute(i))) {
+    nm <- sprintf("'%s'", nm)
+    if (is.character(i)) {
+        i <- match.arg(i, ok)
+        return(match(i, ok))
+    }
+    if (is.numeric(i) && i == round(i) && length(i) == 1) {
+        if (i < 1 || i > length(ok)) {
+            stop(sprintf("invalid %s index: %d (max: %d)", nm, i, length(ok)))
+        }
+        return(as.integer(i))
+    }
+    stop(sprintf("invalid %s; expected character or integer index", nm))
+}
+
 # validate OME version
 .val_ome_ver <- \(v) {
     ok <- length(v) == 1 && is.character(v) && (v <- gsub("-.*", "", v)) %in% sprintf("0.%d", seq_len(6))

inst/NEWS

@@ -1,7 +1,15 @@
+changes in version 0.99.43
+
+- revised validity to allow 4D data (time dimension)
+- revised 'mask()' to support 4D, returning one sheet per timepoint
+- revised 'SpatialDataAttrs()' to generate attributes for any dimensions
+- refactored functions for spatial axiis, multiscales & tables handling
+- added unit tests for 'show()', >2D support, table synchronization
+
 changes in version 0.99.42
 
+- revised subsetting to support >2D arrays (t & z dims. are protected)
 - revised 'crop()' to adjust for array multiscales with scale factor != 1
-- revised subsetting to support >2/3D arrays (t & z dims. are protected)
 
 changes in version 0.99.41
 

tests/testthat/test-mask.R

@@ -52,6 +52,36 @@ test_that("mask,sdImage,sdLabel", {
     .i@data <- lapply(.i@data, \(.) .[,,-1])
     .x <- x; image(.x, i) <- .i
     expect_error(mask(.x, i, j))
+
+    # 3D
+    z <- 5; n <- 7; m <- 8
+    u <- array(runif(z*n*m), c(1, z, n, m))
+    v <- array(1L, c(z, n, m))
+    i <- SpatialDataImage(u, SpatialDataAttrs(dim=3, nch=1))
+    l <- SpatialDataLabel(v, SpatialDataAttrs(type="label", dim=3))
+    sd <- SpatialData(images=list(a=i), labels=list(b=l))
+    sd <- expect_silent(mask(sd, "a", "b", how="mean"))
+    expect_identical(as.numeric(assay(table(sd))), mean(u))
+
+    # 4D
+    t <- 4; z <- 5; n <- 7; m <- 8
+    u <- array(runif(t*z*n*m), c(t, 1, z, n, m))
+    v <- array(sample(9L, t*z*n*m, TRUE), c(t, z, n, m))
+    i <- SpatialDataImage(u, SpatialDataAttrs(dim=4, nch=1))
+    l <- SpatialDataLabel(v, SpatialDataAttrs(type="label", dim=4))
+    sd <- SpatialData(images=list(a=i), labels=list(b=l))
+    sd <- expect_silent(mask(sd, "a", "b", how="mean"))
+    se <- table(sd)
+    # should get one sheet per timepoint
+    expect_length(assays(se), t)
+    expect_equal(dim(se), c(1,9))
+    # verify that aggregation went right
+    se <- table(mask(sd, "a", "b", how="sum"))
+    for (t in seq_along(assays(se)))
+        for (i in seq_len(ncol(se))) {
+            n <- (v[t,,,] == i) %*% drop(u)[t,,,]
+            expect_equal(as.numeric(n), assay(se, t)[i])
+        }
 })
 
 test_that("mask w/ transform", {
@@ -179,4 +209,7 @@ test_that("mask,sdShape,sdShape", {
         expect_equal(dim(sf), c(7,2))
         expect_identical(assay(sf)[,"1"], fun(mx))
     }
+
+    # default to 'sum' with a message
+    expect_message(mask(y, i, j))
 })