diff --git a/DESCRIPTION b/DESCRIPTION index d91aad48..c444cccf 100644 --- a/DESCRIPTION +++ b/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 diff --git a/NAMESPACE b/NAMESPACE index 2820a59d..4cc09aef 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/AllClasses.R b/R/AllClasses.R index 3ab297df..870b718b 100644 --- a/R/AllClasses.R +++ b/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", diff --git a/R/CTutils.R b/R/CTutils.R index 0dbe1155..a185485f 100644 --- a/R/CTutils.R +++ b/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]]) diff --git a/R/crop.R b/R/crop.R index be61d8e6..9dd706ee 100644 --- a/R/crop.R +++ b/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) }) + diff --git a/R/mask.R b/R/mask.R index 4bed9970..d7505058 100644 --- a/R/mask.R +++ b/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, ...) { diff --git a/R/read.R b/R/read.R index 78874855..d59a5215 100644 --- a/R/read.R +++ b/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) diff --git a/R/trans.R b/R/trans.R index 4dd504a3..3e265e77 100644 --- a/R/trans.R +++ b/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") { diff --git a/R/utils.R b/R/utils.R index db86ae46..30627fb1 100644 --- a/R/utils.R +++ b/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)) diff --git a/inst/NEWS b/inst/NEWS index 171df13b..4d989426 100644 --- a/inst/NEWS +++ b/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 diff --git a/tests/testthat/test-mask.R b/tests/testthat/test-mask.R index cf1e938e..84be06be 100644 --- a/tests/testthat/test-mask.R +++ b/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)) }) diff --git a/tests/testthat/test-sync.R b/tests/testthat/test-sync.R new file mode 100644 index 00000000..bdfc0f16 --- /dev/null +++ b/tests/testthat/test-sync.R @@ -0,0 +1,47 @@ +require(SingleCellExperiment, quietly=TRUE) +zs <- file.path("extdata", "blobs.zarr") +zs <- system.file(zs, package="spatialdataR") +sd <- readSpatialData(zs) + +test_that(".sync_tables_on_crop,label", { + i <- region(table(sd)) + e <- element(sd, i) + + # crop around one label + id <- sample(is <- instances(e), 1) + mx <- as(data(e), "dgCMatrix") + ij <- (nz <- summary(mx))[nz$x == id, ] + dx <- range(ij[, 2]); dy <- range(ij[, 1]) + bb <- list(xmin=dx[1], xmax=dx[2], ymin=dy[1], ymax=dy[2]) + x <- crop(sd, bb) + + # verify table integrity + se <- getTable(x, i) + mx <- data(element(x, i)) + ni <- length(is <- unique(mx[mx != 0])) + + expect_equal(ncol(se), ni) # right instances + expect_identical(region(se), i) # regions untouched +}) + +test_that(".sync_tables_on_crop,shape", { + i <- shapeNames(sd)[1] + e <- element(sd, i) + + # mock shape annotation + mx <- matrix(nc=length(e)) + se <- SingleCellExperiment(list(mx)) + sd <- setTable(sd, i, se) + + # crop around one shape + id <- sample(length(e), 1) + xy <- centroids(e)[id, ] + dx <- xy[[1]]+c(-(. <- 1e-3), .); dy <- xy[[2]]+c(-., .) + bb <- list(xmin=dx[1], xmax=dx[2], ymin=dy[1], ymax=dy[2]) + x <- crop(sd, bb) + + # verify table integrity + se <- getTable(x, i) + expect_equal(ncol(se), 1) # right instances + expect_identical(region(se), i) # regions untouched +}) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index e3d1bc22..1ed50e22 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -11,7 +11,6 @@ test_that("centroids,invalid", { expect_error(centroids(image(x)), "supported") }) test_that("centroids,sdLabel", { - # 2D y <- label(x) z <- centroids(y, "data.frame") expect_is(z, "data.frame") @@ -22,15 +21,6 @@ test_that("centroids,sdLabel", { expect_is(.z, "matrix") z$i <- as.integer(as.character(z$i)) expect_identical(.z, as.matrix(z)) - # 3D - n <- 7; m <- 8 - u <- array(runif(2*n*m), c(1, 5, n, m)) - v <- array(1L, c(5, 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)) }) test_that("centroids,sdPoint", { i <- feature_key(y <- point(x))