utilities (query, mask, etc.)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SpatialData
 Title: Representation of Python's SpatialData in R
 Depends: R (>= 4.5)
-Version: 0.99.26
+Version: 0.99.27
 Description: Interface to Python's 'SpatialData', currently including:
     reticulate-based use of 'spatialdata-io' for reading of manufacturer
     data and writing to .zarr, on-disk representation of images/labels as

NAMESPACE

@@ -52,6 +52,7 @@ exportMethods(data)
 exportMethods(data_type)
 exportMethods(dim)
 exportMethods(element)
+exportMethods(feature_key)
 exportMethods(getTable)
 exportMethods(hasTable)
 exportMethods(image)
@@ -93,14 +94,13 @@ importFrom(DelayedArray,ConstantArray)
 importFrom(DelayedArray,DelayedArray)
 importFrom(DelayedArray,cbind)
 importFrom(DelayedArray,rbind)
-importFrom(DelayedArray,realize)
 importFrom(Matrix,rowSums)
+importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,sparseVector)
 importFrom(Matrix,t)
 importFrom(RBGL,sp.between)
 importFrom(Rarr,read_zarr_attributes)
 importFrom(Rarr,zarr_overview)
-importFrom(S4Arrays,as.array.Array)
 importFrom(S4Vectors,"metadata<-")
 importFrom(S4Vectors,coolcat)
 importFrom(S4Vectors,isSequence)
@@ -112,6 +112,8 @@ importFrom(SingleCellExperiment,"int_metadata<-")
 importFrom(SingleCellExperiment,SingleCellExperiment)
 importFrom(SingleCellExperiment,int_colData)
 importFrom(SingleCellExperiment,int_metadata)
+importFrom(SummarizedExperiment,"assay<-")
+importFrom(SummarizedExperiment,"assayNames<-")
 importFrom(SummarizedExperiment,"colData<-")
 importFrom(SummarizedExperiment,assay)
 importFrom(SummarizedExperiment,colData)
@@ -150,12 +152,14 @@ importFrom(methods,setReplaceMethod)
 importFrom(reticulate,import)
 importFrom(sf,"st_geometry<-")
 importFrom(sf,st_as_sf)
+importFrom(sf,st_bbox)
 importFrom(sf,st_coordinates)
+importFrom(sf,st_crop)
 importFrom(sf,st_distance)
 importFrom(sf,st_geometry)
 importFrom(sf,st_geometry_type)
-importFrom(sf,st_point)
-importFrom(sf,st_sfc)
+importFrom(sf,st_intersects)
+importFrom(sf,st_polygon)
 importFrom(utils,.DollarNames)
 importFrom(utils,head)
 importFrom(utils,tail)

R/AllGenerics.R

@@ -57,6 +57,10 @@ setGeneric("rotate", \(x, t, ...) standardGeneric("rotate"))
 setGeneric("transform", \(x, ...) standardGeneric("transform"))
 setGeneric("translation", \(x, t, ...) standardGeneric("translation"))
 
+# sda ----
+
+setGeneric("feature_key", \(x, ...) standardGeneric("feature_key"))
+
 # uts ----
 
 setGeneric("layer", \(x, i, ...) standardGeneric("layer"))

R/Zattrs.R

@@ -1,6 +1,8 @@
 #' @name Zattrs
 #' @title The `Zattrs` class
-#'
+#' 
+#' @aliases feature_key
+#' 
 #' @param x list extracted from a OME-NGFF compliant .zattrs file.
 #' @param name character string for extraction (see ?base::`$`).
 #' 
@@ -16,6 +18,8 @@
 #' CTname(z)
 #' CTtype(z)
 #' CTdata(z, "scale")
+#' 
+#' feature_key(point(x))
 #'
 #' @export
 Zattrs <- \(x=list()) {
@@ -75,3 +79,9 @@ setMethod("$", "Zattrs", \(x, name) x[[name]])
     if (!is.null(cs)) coolcat("channels(%d): %s\n", cs)
 }
 setMethod("show", "Zattrs", .showZattrs)
+
+#' @export
+setMethod("feature_key", "Zattrs", \(x) x$spatialdata_attrs$feature_key)
+
+#' @export
+setMethod("feature_key", "SpatialDataElement", \(x) feature_key(meta(x)))

R/coord_utils.R

@@ -1,6 +1,6 @@
 #' @name coord-utils
 #' @title Coordinate transformations
-#' @aliases axes CTname CTtype CTdata CTpath CTgraph addCT rmvCT
+#' @aliases axes CTlist CTname CTtype CTdata CTpath CTgraph addCT rmvCT
 #' 
 #' @param x \code{SpatialData}, an element, or \code{Zattrs}.
 #' @param i for \code{CTpath}, source node label; else, string or 

R/mask.R

@@ -1,104 +1,166 @@
 #' @name mask
 #' @title Masking
 #'
-#' @description ...
+#' @description 
+#' Masking operations serve to aggregate data across layers, e.g., 
+#' counting points in shapes, averaging image channels by labels, etc.
+#' For added flexibility, these may be carried out directly between elements,
+#' or using an input \code{SpatialData} object and specifying element names.
 #' 
 #' @param x \code{\link{SpatialData}} object.
 #' @param i,j character string; names of elements to mask,
 #'   specifically, \code{i} will be masked by \code{j},
 #'   adding a \code{table} for \code{j} in \code{x}.
+#' @param how character string; statistic to use for masking.
+#' @param name function use to generate the new \code{table}'s name.
 #' @param ... optional arguments passed to and from other methods.
 #'
-#' @return \code{\link{SingleCellExperiment}}
+#' @return Input \code{SpatialData} object \code{x} with an additional table.
 #'
 #' @examples
 #' library(SingleCellExperiment)
 #' x <- file.path("extdata", "blobs.zarr")
 #' x <- system.file(x, package="SpatialData")
 #' x <- readSpatialData(x, tables=FALSE)
 #'
-#' # count points in circles
-#' x <- mask(x, "blobs_points", "blobs_circles")
-#' x <- mask(x, "blobs_image", "blobs_labels")
-#' tables(x)
+#' # count points in shapes
+#' y <- mask(x, "blobs_points", "blobs_circles")
+#' tail(tables(y), 1)
+#' 
+#' # average image channels by labels
+#' y <- mask(x, "blobs_image", "blobs_labels")
+#' tail(tables(y), 1)
 #'
+#' library(SpatialData.data)
+#' x <- get_demo_SDdata("merfish")
+#' x <- readSpatialData(x)
+#' 
+#' # sum table counts by shapes
+#' y <- mask(x, "cells", "anatomical")
+#' tail(tables(y), 1)
+#' 
 #' @export
 NULL
 
-# TODO: table from point + shape, image + label etc. etc. etc.
+.check_ij <- \(x, .) stopifnot(length(.) == 1, is.character(.), . %in% unlist(colnames(x)))
 
 #' @rdname mask
+#' @importFrom methods as
+#' @importFrom SummarizedExperiment assay assay<-
 #' @importFrom SingleCellExperiment int_colData int_colData<- int_metadata<-
 #' @export
-setMethod("mask", "SpatialData", \(x, i, j, ...) {
-    stopifnot(length(i) == 1, is.character(i), i %in% unlist(colnames(x)))
-    stopifnot(length(j) == 1, is.character(j), j %in% unlist(colnames(x)))
-    # get element types
-    ls <- vapply(
-        list(i, j), \(e) rownames(x)[vapply(colnames(x), 
-            \(es) e %in% es, logical(1))], character(1))
-    a <- element(x, ls[[1]], i)
-    b <- element(x, ls[[2]], j)
-    t <- .mask(a, b, ...)
-    md <- list(region=j, 
-        region_key="region", 
-        instance_key="instance")
-    int_metadata(t)$spatialdata_attrs <- md
-    cd <- data.frame(region=j, instance=colnames(t))
-    int_colData(t) <- cbind(int_colData(t), cd)
-    nm <- paste0(i, "_masked_by_", j)
-    `table<-`(x, nm, value=t)
+setMethod("mask", c("SpatialData", "ANY", "ANY"), \(x, i, j, 
+    how=NULL, name=\(i, j) sprintf("%s_by_%s", i, j), ...) {
+    .check_ij(x, i); .check_ij(x, j)
+    #if (!is.null(how)) how <- match.arg(how, c("sum", "mean"))
+    ok <- is.character(name) && length(name) == 1 && !name %in% tableNames(x)
+    nm <- if (is.function(name)) name(i, j) else if (ok) name else stop(
+        "Invalid 'name'; should be a function or a ",
+        "character string not yet in 'tableNames(x)'")
+    f <- \(i) names(which(rapply(colnames(x), \(.) i %in% ., "character")))
+    .i <- element(x, f(i), i)
+    .j <- element(x, f(j), j)
+    t <- tryCatch(error=\(.) NULL, getTable(x, i))
+    se <- .mask(.i, .j, how=how, table=t, ...)
+    md <- list(region=j, region_key="region", instance_key="instance")
+    int_metadata(se)$spatialdata_attrs <- md
+    assay(se) <- as(assay(se), "dgCMatrix")
+    cd <- int_colData(se)
+    cd$region <- j
+    cd$instance <- colnames(se)
+    int_colData(se) <- cd
+    `table<-`(x, nm, value=se)
 })
 
-setGeneric(".mask", \(a, b, ...) standardGeneric(".mask"))
+setGeneric(".mask", \(i, j, ...) standardGeneric(".mask"))
 
+#' @noRd
 #' @importFrom methods as
-#' @importFrom Matrix rowSums sparseVector t
+#' @importFrom Matrix sparseVector
+#' @importFrom SummarizedExperiment assayNames<-
 #' @importFrom SingleCellExperiment SingleCellExperiment
-#' @importFrom sf st_as_sf st_geometry_type st_sfc st_point st_distance
-setMethod(".mask", c("PointFrame", "ShapeFrame"), \(a, b) {
-    n <- nrow(b <- st_as_sf(data(b)))
-    fk <- meta(a)$spatialdata_attrs$feature_key
-    switch(paste(st_geometry_type(b)[1]),
-        POINT={
-            # realize one feature at a time
-            is <- split(seq_len(length(a)), a[[fk]])
-            ns <- lapply(is, \(.) {
-                # make points 'sf'-compliant
-                xy <- as.data.frame(a[., c("x", "y")])
-                ps <- st_sfc(lapply(asplit(xy, 1), st_point))
-                # for each circle, count points within radius
-                z <- rowSums(st_distance(b, ps) < b$radius)
-                # sparsify counts
-                sv <- sparseVector(z[i <- z > 0], which(i), n)
-                sm <- as(sv, "sparseMatrix")
-            })
-            # collect intro matrix w/ dim. features x circles
-            ns <- t(as(do.call(cbind, ns), "dgCMatrix"))
-            rownames(ns) <- names(is)
-            colnames(ns) <- seq(ncol(ns))
-        })
-    SingleCellExperiment(list(counts=ns))
+setMethod(".mask", c("ImageArray", "LabelArray"), \(i, j, how=NULL, ...) {
+    if (is.null(how)) { how <- "mean"; message("Missing 'how'; defaulting to 'mean'") }
+    stopifnot(dim(i)[-1] == dim(j))
+    .j <- as(data(j), "sparseVector")
+    .j <- as.vector(.j[ok <- .j > 0])
+    mx <- apply(data(i), 1, \(.i) {
+        .i <- as(.i, "sparseVector")
+        .i <- as.vector(.i[ok])
+        tapply(.i, .j, how)
+    })
+    colnames(mx) <- channels(i)
+    se <- SingleCellExperiment(list(t(mx)))
+    assayNames(se) <- how
+    return(se)
 })
 
+#' @noRd
 #' @importFrom methods as
-#' @importFrom DelayedArray realize
-#' @importFrom S4Arrays as.array.Array
+#' @importFrom Matrix t rowSums sparseVector sparseMatrix
 #' @importFrom SingleCellExperiment SingleCellExperiment
-setMethod(".mask", c("ImageArray", "LabelArray"), \(a, b, fun=mean) {
-    # TODO: somehow rewrite w/o realizing everything 
-    # at once (maybe w/ 'DelayedArray::blockApply'?)
-    .a2v <- \(.) as.vector(as.array.Array(.))
-    stopifnot(dim(a)[-1] == dim(b))
-    w <- .a2v(data(b)); w[w == 0] <- NA 
-    n <- length(i <- unique(w[!is.na(w)]))
-    ns <- vapply(seq_len(dim(a)[1]), \(.) {
-        v <- .a2v(data(a, 1)[., , ])
-        tapply(v, w, sum, na.rm=TRUE)
-    }, numeric(n))
-    ns <- t(as(ns, "dgCMatrix"))
-    dimnames(ns) <- list(seq(dim(a)[1]), i)
+#' @importFrom sf st_as_sf st_geometry_type st_distance
+setMethod(".mask", c("PointFrame", "ShapeFrame"), \(i, j, how=NULL, ...) {
+    if (!is.null(how)) warning("Can only count when masking points; ignoring 'how'")
+    n <- nrow(j <- st_as_sf(data(j)))
+    fun <- switch(as.character(st_geometry_type(j[1, ])),
+        POINT=\(i, j) rowSums(st_distance(j, i) <= j$radius),
+        \(i, j) vapply(st_intersects(j, i), length, integer(1)))
+    # realize one feature at i time
+    is <- split(seq_len(length(i)), i[[feature_key(i)]])
+    ns <- lapply(is, \(.) {
+        # make points 'sf'-compliant
+        i <- as.data.frame(i[., c("x", "y")])
+        i <- st_as_sf(i, coords=c("x", "y"))
+        # for each shape, count intersecting points
+        z <- fun(i, j)
+        # sparsify counts
+        sv <- sparseVector(z[i <- z > 0], which(i), n)
+        sm <- as(sv, "sparseMatrix")
+    })    
+    # collect into matrix w/ dim. features x shapes
+    ns <- t(do.call(cbind, ns))
+    rownames(ns) <- names(is)
+    colnames(ns) <- seq(ncol(ns))
     SingleCellExperiment(list(counts=ns))
 })
-setMethod(".mask", c("ANY", "ANY"), \(a, b) 
-    stop("'mask'ing between these element types not supported."))
+
+#' @noRd
+#' @importFrom methods as
+#' @importFrom Matrix sparseMatrix
+#' @importFrom SummarizedExperiment assay
+#' @importFrom SingleCellExperiment SingleCellExperiment
+setMethod(".mask", c("ShapeFrame", "ShapeFrame"), \(i, j, how=NULL, table=NULL, value=NULL, assay=1, ...) {
+    # validity
+    if (is.null(table)) stop("Missing 'table'; can't mask shapes without")
+    ok <- is.null(value) || (is.character(value) && all(value %in% rownames(table)))
+    if (!ok) stop("Invalid 'value'; should be in 'rownames(table(x, i))'")
+    if (is.null(how)) { how <- "sum"; message("Missing 'how'; defaulting to 'sum'") }
+    if (is.character(how)) how <- match.arg(how, c("sum", "mean", "detected", "prop.detected"))
+    # grouping
+    js <- st_intersects(st_as_sf(data(j)), st_as_sf(data(i)))
+    is <- factor(integer(nrow(i)), seq(0, nrow(j)))
+    is[unlist(js)] <- rep(seq_along(js), lengths(js))
+    ns <- tabulate(is, ni <- nlevels(is))
+    # aggregation
+    mx <- assay(table, assay)
+    if (grepl("detected$", how)) mx <- mx > 0
+    my <- sparseMatrix(
+        x=rep(1, length(is)), 
+        i=seq_along(is), j=is, 
+        dims=c(ncol(table), ni))
+    mx <- mx %*% my
+    if (grepl("mean|prop", how)) mx <- t(t(mx)/ns)
+    # wrangling
+    mx <- as(mx, "dgCMatrix")
+    colnames(mx) <- levels(is)
+    mx <- list(mx); names(mx) <- how
+    se <- SingleCellExperiment(mx)
+    nm <- paste0("n_", meta(table)$region)
+    se[[nm]] <- ns
+    return(se)
+})
+
+#' @noRd
+setMethod(".mask", c("ANY", "ANY"), \(i, j, ...) 
+    stop("'mask'ing between these element types not yet supported"))