diff --git a/DESCRIPTION b/DESCRIPTION index da1c373c..ca319efa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GiottoClass Title: Giotto Suite Object Definitions and Framework -Version: 0.4.12 +Version: 0.5.0 Authors@R: c( person("Ruben", "Dries", email = "rubendries@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7650-7754")), diff --git a/NAMESPACE b/NAMESPACE index 0271dc38..52e58540 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,6 +8,8 @@ S3method(.DollarNames,processParam) S3method(.DollarNames,spatEnrObj) S3method(.DollarNames,spatLocsObj) S3method(.DollarNames,terraVectData) +S3method(as.data.frame,overlapIntensityDT) +S3method(as.data.frame,overlapPointDT) S3method(as.data.table,SpatVector) S3method(as.data.table,giottoPoints) S3method(as.data.table,giottoPolygon) @@ -28,6 +30,7 @@ export(addNetworkLayout) export(addSpatialCentroidLocations) export(addSpatialCentroidLocationsLayer) export(add_img_array_alpha) +export(aggregateFeatures) export(aggregateStacks) export(aggregateStacksExpression) export(aggregateStacksLocations) @@ -297,6 +300,7 @@ exportClasses(giottoLargeImage) exportClasses(giottoPoints) exportClasses(giottoPolygon) exportClasses(nnNetObj) +exportClasses(overlapPointDT) exportClasses(processParam) exportClasses(spatEnrObj) exportClasses(spatLocsObj) diff --git a/NEWS.md b/NEWS.md index 74dbdae7..7e342708 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,22 @@ +# GiottoClass 0.5.0 + +## changes +- `calculateOverlap()` and `overlapToMatrix()` param harmonization + +## new +- `aggregateFeatures()` wrapper for running `calculateOverlap()` and `overlapToMatrix()` +- `overlapPointDT()` and `overlapIntensityDT()` classes to store overlaps relationships efficiently and help with aggregation pipeline + +## bug fixes +- `overlaps()` will now properly find image overlaps + + # GiottoClass 0.4.12 (2025/12/12) ## enhancements - automatic checking for `"count"` column in feature info + ## new - `misc` slot for storing unstructured data diff --git a/R/aggregate.R b/R/aggregate.R index 063fc038..8c286bf6 100644 --- a/R/aggregate.R +++ b/R/aggregate.R @@ -18,6 +18,178 @@ NULL ## calculate overlap between cellular structures and features #### +#' @name aggregateFeatures +#' @title Aggregate Spatial Features Covered by Polygon Geometries +#' @description +#' Aggregate features (either `feat_info` OR `image_names`) with +#' polygons (`spat_info`). Under the hood, this is performed in two steps: +#' +#' 1. Find the overlapped features via the lower-level generic +#' [calculateOverlap()] +#' 2. Summarize the overlapped features as a matrix via [overlapToMatrix()] +#' @param gobject `giotto` object containing spatial data to aggregate. +#' @param spat_info character. Name of polygon information to use to aggregate +#' features with. +#' @param feat_info character. Name of feature point detections to be +#' aggregated. +#' @param image_names character. Name of image(s) containing intensity values +#' to be aggregated. +#' @param new_spat_unit character (optional). Name of spatial unit to assign +#' the expression info to. Default is the same as `spat_info`. +#' @param new_feat_type character (optional). Name of feature type to assign +#' the expression info to. Default is the same as `feat_info` when used. When +#' `image_names` is used instead, default is "protein". +#' @param name character. (default = "raw") Name to assign the output +#' expresssion values information. +#' @param poly_subset_ids character vector. (optional) Specific poly_IDs to use +#' @param feat_subset_column character. (optional) feature info attribute to +#' subset feature points on when performing overlap calculation. +#' @param feat_subset_values (optional) values matched against +#' in `feat_subset_column` in order to subset feature points when performing +#' overlap calculation. +#' @param feat_count_column character. (optional) feature info column with counts +#' information. Useful in cases when more than one detection is reported per +#' point. +#' @param fun character (default = "sum"). A function usable by +#' [exactextractr::exact_extract()] to aggregate image intensity values. +#' @param return_gobject logical (default = TRUE). Whether to return the +#' `giotto` object or just the aggregated expression values as `exprObj` class. +#' @param verbose logical. Whether to be verbose. +#' @param \dots Additional params to pass to the overlap calculation method. +#' None implemented for point overlaps. For intensity overlaps, passes to +#' [exactextractr::exact_extract()] and additionally the function requested +#' with the `fun` param. +#' @returns `giotto` when `return_gobject=TRUE`, `exprObj` when +#' `return_gobject=FALSE` +#' @details `feat_subset_column`, `feat_subset_values`, and `feat_count_column` +#' are specific to overlaps on feature points info, and should not be provided +#' when overlapping image data. +#' @export +aggregateFeatures <- function(gobject, + spat_info = NULL, + feat_info = NULL, + image_names = NULL, + new_spat_unit = NULL, + new_feat_type = NULL, + name = "raw", + poly_subset_ids = NULL, + feat_subset_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, + fun = "sum", + return_gobject = TRUE, + verbose = TRUE, + ... +) { + checkmate::assert_character(spat_info, len = 1L, null.ok = TRUE) + checkmate::assert_character(feat_info, null.ok = TRUE) + checkmate::assert_character(new_spat_unit, len = 1L, null.ok = TRUE) + checkmate::assert_character(new_feat_type, len = 1L, null.ok = TRUE) + checkmate::assert_character(image_names, null.ok = TRUE) + checkmate::assert_character(fun, len = 1L) + checkmate::assert_character(name, len = 1L) + checkmate::assert_logical(return_gobject) + checkmate::assert_character(feat_count_column, null.ok = TRUE) + + # catch improper feature input usage + fun_tag <- "[aggregateFeatures] " + if (!is.null(feat_info) && !is.null(image_names)) { + stop(fun_tag, "Only one of 'feat_info' or 'image_names' may ", + "be provided.\n", call. = FALSE) + } + + # decide polygons to use + spat_info <- spat_info %null% names(gobject@spatial_info)[[1]] + + # decide target spatial unit + new_spat_unit <- new_spat_unit %null% spat_info + # decide target feature type + if (is.null(new_feat_type)) { + if (!is.null(feat_info)) new_feat_type <- feat_info + if (!is.null(image_names)) new_feat_type <- "protein" + } + + # select overlap type (point geoms vs intensity rasters) + overlap_type <- if (!is.null(feat_info)) { + "point" + } else if (!is.null(image_names)) { + "intensity" + } else { + "point" # fallback assumption + } + + # calculate overlap --------------------------------------------------- # + calculate_overlap_params <- list( # common params + x = gobject, + spat_info = spat_info, + feat_info = feat_info, + name_overlap = new_feat_type, + poly_subset_ids = poly_subset_ids, + return_gobject = TRUE, + verbose = verbose, + ... + ) + + # method specific params + switch(overlap_type, + "point" = { + calculate_overlap_params <- c(calculate_overlap_params, + list( + feat_subset_column = feat_subset_column, + feat_subset_values = feat_subset_values, + feat_count_column = feat_count_column + ) + ) + }, + "intensity" = { + calculate_overlap_params <- c(calculate_overlap_params, + list(image_names = image_names) + ) + } + ) + + gobject <- do.call(calculateOverlap, args = calculate_overlap_params) + + # overlap to matrix --------------------------------------------------- # + + overlap_to_matrix_params <- list( + x = gobject, + name = name, + spat_info = spat_info, + feat_info = new_feat_type, # this difference is on purpose + type = overlap_type, + return_gobject = FALSE, + verbose = verbose + ) + + switch(overlap_type, + "point" = { + overlap_to_matrix_params <- c(overlap_to_matrix_params, + list(feat_count_column = feat_count_column) + ) + }, + "intensity" = { + overlap_to_matrix_params <- c(overlap_to_matrix_params, + list(fun = fun) + ) + } + ) + + ex <- do.call(overlapToMatrix, args = overlap_to_matrix_params) + # this is moved out of the gobject since overlapToMatrix doesn't have a + # way to set the spat unit (if it doesn't match spat_info) + spatUnit(ex) <- new_spat_unit + + gobject <- setGiotto(gobject, ex, verbose = verbose) + return(gobject) +} + + + + + + + ### ** raster way #### #' @title Convert polygon to raster @@ -60,6 +232,9 @@ polygon_to_raster <- function(polygon, field = NULL) { } + + + # calculateOverlap methods #### @@ -76,17 +251,19 @@ polygon_to_raster <- function(polygon, field = NULL) { #' @param poly_subset_ids character vector. (optional) Specific poly_IDs to use #' @param feat_subset_column character. (optional) feature info attribute to #' subset feature points on when performing overlap calculation. -#' @param feat_subset_ids (optional) values matched against +#' @param feat_subset_values (optional) values matched against #' in `feat_subset_column` in order to subset feature points when performing #' overlap calculation. -#' @param count_info_column character. (optional) column with count information. +#' @param feat_subset_ids deprecated. Use `feat_subset_values` instead. +#' @param feat_count_column character. (optional) column with count information. #' Useful in cases when more than one detection is reported per point. If a #' column called "count" is present in the feature points data, it will be #' automatically selected. #' @param verbose be verbose +#' @param count_info_column deprecated. Use `feat_count_column` instead. #' @param \dots additional params to pass to methods. -#' @details `feat_subset_column`, `feat_subset_ids`, and `count_info_column` are -#' specific to overlaps on feature points info, and should not be provided +#' @details `feat_subset_column`, `feat_subset_values`, and `feat_count_column` +#' are specific to overlaps on feature points info, and should not be provided #' when overlapping image data. These three params can also be passed to the #' `giotto` method through the `...` param when working with overlaps on feature #' points info. @@ -117,7 +294,7 @@ polygon_to_raster <- function(polygon, field = NULL) { #' # calculate z1 only #' out_z1 <- calculateOverlap(gpoly, gpoints, #' feat_subset_column = "global_z", -#' feat_subset_ids = c(1) +#' feat_subset_values = c(1) #' ) #' overlaps_z1 <- overlaps(out_z1) #' overlaps_z1$rna @@ -131,7 +308,7 @@ polygon_to_raster <- function(polygon, field = NULL) { #' # calculate z0 overlaps and return as gobject #' out_g <- calculateOverlap(g, #' feat_subset_column = "global_z", -#' feat_subset_ids = 0 +#' feat_subset_values = 0 #' ) #' overlaps(getPolygonInfo(out_g, return_giottoPolygon = TRUE)) #' @@ -153,13 +330,19 @@ setMethod( "calculateOverlap", signature(x = "giotto", y = "missing"), function(x, name_overlap = NULL, - spatial_info = NULL, + spat_info = NULL, feat_info = NULL, image_names = NULL, poly_subset_ids = NULL, return_gobject = TRUE, verbose = TRUE, + spatial_info = deprecated(), ...) { + # deprecations + spat_info <- GiottoUtils::deprecate_param(spatial_info, spat_info, + fun = "calculateOverlap", when = "0.4.7" + ) + # 0. guards # # --------- # @@ -215,8 +398,8 @@ setMethod( if (!is.null(image_names)) name_overlap <- "protein" } - if (is.null(spatial_info)) { - spatial_info <- names(x@spatial_info)[[1]] + if (is.null(spat_info)) { + spat_info <- names(x@spatial_info)[[1]] } @@ -226,7 +409,7 @@ setMethod( # ---[polys to overlap with]--- A <- getPolygonInfo( gobject = x, - polygon_name = spatial_info, + polygon_name = spat_info, return_giottoPolygon = TRUE ) @@ -246,8 +429,8 @@ setMethod( for (img_name in image_names) { if (!img_name %in% potential_large_image_names) { warning( - "image with the name ", img_name, - " was not found and will be skipped \n" + "[calculateOverlap] image with the name ", img_name, + " was not found and will be skipped \n", call. = FALSE ) } } @@ -325,35 +508,53 @@ setMethod( name_overlap = NULL, poly_subset_ids = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, return_gpolygon = TRUE, verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated(), ...) { - if ("count" %in% names(y) && is.null(count_info_column)) { + # deprecations + feat_subset_values <- GiottoUtils::deprecate_param( + feat_subset_ids, feat_subset_values, + fun = "calculateOverlap", when = "0.5.0" + ) + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "calculateOverlap", when = "0.5.0" + ) + + # autodetect count columns + if ("count" %in% names(y) && is.null(feat_count_column)) { vmsg(.v = verbose, "[overlap] Found column \"count\" in feature info. - - Using as `count_info_column` - [!] Set count_info_column = FALSE to disable.") - count_info_column <- "count" + - Using as `feat_count_column` + [!] Set feat_count_column = FALSE to disable.") + feat_count_column <- "count" } - if (isFALSE(count_info_column)) { - count_info_column <- NULL + if (isFALSE(feat_count_column)) { + feat_count_column <- NULL } + # return an overlap info object res <- calculateOverlap( x = x[], y = y[], poly_subset_ids = poly_subset_ids, feat_subset_column = feat_subset_column, - feat_subset_ids = feat_subset_ids, - count_info_column = count_info_column, + feat_subset_values = feat_subset_values, + feat_count_column = feat_count_column, verbose = verbose, ... ) if (isTRUE(return_gpolygon)) { + # update schema metadata in overlap object if (is.null(name_overlap)) name_overlap <- objName(y) + prov(res) <- spatUnit(x) + spatUnit(res) <- spatUnit(x) + featType(res) <- name_overlap # ensure centroids calculated if (is.null(centroids(x))) { @@ -404,8 +605,10 @@ setMethod( verbose = TRUE, ...) { aff <- y@affine + # perform the reverse of the image's affine on the polys + inv_aff_poly <- affine(x, aff, inv = TRUE) res <- calculateOverlap( - x = affine(x, aff, inv = TRUE), + x = inv_aff_poly, y = y@raster_object, name_overlap = name_overlap %null% objName(y), poly_subset_ids = poly_subset_ids, @@ -436,6 +639,7 @@ setMethod( stop("calculateOverlap: name_overlap must be given", call. = FALSE) } + # return overlaps information res <- calculateOverlap( x = x[], y = y, @@ -445,6 +649,11 @@ setMethod( ) if (isTRUE(return_gpolygon)) { + # update schema metadata in overlap object + prov(res) <- spatUnit(x) + spatUnit(res) <- spatUnit(x) + featType(res) <- name_overlap + # ensure centroids calculated if (is.null(centroids(x))) { x <- centroids(x, append_gpolygon = TRUE) @@ -473,6 +682,7 @@ setMethod( function(x, y, poly_subset_ids = NULL, verbose = TRUE, + fun = "sum", ...) { checkmate::assert_true(terra::is.polygons(x)) GiottoUtils::package_check("exactextractr") @@ -485,9 +695,6 @@ setMethod( names(y) <- sprintf("channel_%d", seq_len(nchannel)) } - # NSE vars - coverage_fraction <- NULL - # subset polys if needed if (!is.null(poly_subset_ids)) { x <- x[x$poly_ID %in% poly_subset_ids] @@ -498,54 +705,72 @@ setMethod( vmsg(.v = verbose, "Start image extract") - # perform extraction, producing list of results + eer_cname_fun <- function( # how to construct output colnames + values, # name of value layer + weights, # name of weight layer + fun_name, # value of fun + fun_value, # value associated with fun (quantile/frac/wfrac) + nvalues, # number of value layers + nweights # number of weight layers + ) { + values + } + + # perform extraction extract_res <- exactextractr::exact_extract( x = y, y = sf_x, - include_cols = "poly_ID", + append_cols = "poly_ID", + fun = fun, + colname_fun = eer_cname_fun, ... ) + data.table::setDT(extract_res) + extract_res <- .create_overlap_intensity_dt(extract_res) vmsg(.v = verbose, "End image extract") - # rbind and convert output to data.table - dt_exact <- data.table::as.data.table(do.call("rbind", extract_res)) - - missing_ids <- setdiff(x$poly_ID, dt_exact$poly_ID) - if (length(missing_ids) > 0L) { - missing_dt <- data.table::data.table(poly_ID = missing_ids) - # fill out missing_dt with 0 values so it can be used to rbind - for (col in setdiff(names(dt_exact), names(missing_dt))) { - missing_dt[, (col) := 0] - } - dt_exact <- rbind(dt_exact, missing_dt) - } - - # prepare output - colnames(dt_exact)[2:(length(image_names) + 1)] <- image_names - dt_exact[, coverage_fraction := NULL] - - return(dt_exact) + return(extract_res) } ) # * SpatVector SpatVector #### #' @rdname calculateOverlap +#' @param method character. One of `"vector"` or `"raster"`, +#' (default = `"vector"`). Method for polygon-point feature overlap calculation. +#' Can also set as an option: `"giotto.overlap_point_method"` +#' +#' * `"vector"` uses direct spatial extraction (more accurate to input geometry, +#' will double count features in overlapping polygon regions for all overlapping +#' polygons). +#' * `"raster"` uses rasterization (faster, assigns each feature to only one +#' polygon even in overlapping regions as a byproduct of the rasterization). #' @export setMethod( "calculateOverlap", signature(x = "SpatVector", y = "SpatVector"), function(x, y, poly_subset_ids = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, - verbose = TRUE) { + feat_subset_values = NULL, + feat_count_column = NULL, + method = getOption("giotto.overlap_point_method", "vector"), + verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated()) { + method <- match.arg(method, choices = c("vector", "raster")) + feat_subset_values <- GiottoUtils::deprecate_param( + feat_subset_ids, feat_subset_values, + fun = "calculateOverlap", when = "0.4.7" + ) + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "calculateOverlap", when = "0.4.7" + ) + checkmate::assert_true(terra::is.polygons(x)) checkmate::assert_true(terra::is.points(y)) # TODO allow another poly? - if (!is.null(poly_subset_ids)) { - checkmate::assert_character(poly_subset_ids) - } + checkmate::assert_character(poly_subset_ids, null.ok = TRUE) # subset points and polys if needed # * subset x @@ -553,25 +778,94 @@ setMethod( x <- x[x$poly_ID %in% poly_subset_ids] } + # preserve total all feat ids present + feat_ids <- unique(y$feat_ID) + # * subset points if needed # e.g. to select transcripts within a z-plane - if (!is.null(feat_subset_column) && !is.null(feat_subset_ids)) { - bool_vector <- y[[feat_subset_column]][[1]] %in% feat_subset_ids + if (!is.null(feat_subset_column) && !is.null(feat_subset_values)) { + bool_vector <- y[[feat_subset_column]][[1]] %in% feat_subset_values y <- y[bool_vector] } - .calculate_overlap_raster( - spatvec = x, - pointvec = y, - count_info_column = count_info_column, - verbose = verbose + res <- switch(method, + "raster" = .calculate_overlap_raster( + spatvec = x, + pointvec = y, + keep = feat_count_column, + verbose = verbose + ), + "vector" = .calculate_overlap_vector( + spatvec = x, + pointvec = y, + keep = feat_count_column + ) ) + + .create_overlap_point_dt(x, y, res, feat_ids = feat_ids) } ) +#' @param spatvec SpatVector polygon +#' @param pointvec SpatVector points +#' @param keep other col(s) to keep +#' @keywords internal +#' @noRd +.calculate_overlap_vector <- function(spatvec, pointvec, keep = NULL) { + .terra_extract(x = spatvec, y = pointvec, keep = keep) +} +#' @name .calculate_overlap_raster +#' @title Find feature points overlapped by rasterized polygon. +#' @description Core workflow function that accepts simple `SpatVector` inputs, +#' performs rasterization of the polys and then checks for overlaps. +#' @param spatvec `SpatVector` polygon from a `giottoPolygon` object +#' @param pointvec `SpatVector` points from a `giottoPoints` object +#' @param keep column(s) to keep +#' @param verbose be verbose +#' @concept overlap +#' @returns `SpatVector` of overlapped points info +#' @seealso [calculateOverlapRaster()] +#' @keywords internal +.calculate_overlap_raster <- function(spatvec, + pointvec, + keep = NULL, + verbose = TRUE) { + # DT vars + poly_ID <- poly_i <- ID <- x <- y <- feat_ID <- feat_ID_uniq <- NULL + # spatial vector to raster + if (verbose) GiottoUtils::wrap_msg("1. convert polygon to raster \n") + spatrast_res <- polygon_to_raster(spatvec, field = "poly_ID") + spatrast <- spatrast_res[["raster"]] + ID_vector <- spatrast_res[["ID_vector"]] + + ## overlap between raster and point + if (verbose) GiottoUtils::wrap_msg("2. overlap raster and points \n") + .terra_extract(x = spatrast, y = pointvec, keep = keep) +} + +# x is segmentation, y is point +# keep is additional cols of point metadata to keep +# assume output first two cols from `terra::extract()` are: +# (1) point idx, (2) poly idx +# additional cols (if any) are poly or mask attributes and they will be dropped +.terra_extract <- function(x, y, keep, ...) { + checkmate::assert_character(keep, null.ok = TRUE) + res <- terra::extract(x, y, ...)[, 1L:2L] + # 2nd cols can have NA values. NAs denote points that are not overlapped + res <- res[!is.na(res[[2]]),] # drop NAs (sparsify extracted relations) + + # get any needed attributes for `keep` and append them to relations info + if (!is.null(keep)) { + feat_keep <- do.call( + data.frame, terra::as.list(y[][res[[1]], keep]) + ) # list of vectors + res <- cbind(res, feat_keep) + } + return(res) +} #' @title calculateOverlapRaster @@ -585,10 +879,13 @@ setMethod( #' @param poly_ID_names (optional) list of poly_IDs to use #' @param feat_info character. name of feature information #' @param feat_subset_column feature info column to subset features with -#' @param feat_subset_ids ids within feature info column to use for subsetting -#' @param count_info_column column with count information (optional) +#' @param feat_subset_values value(s) within feature info `feat_subset_column` +#' to use for subsetting +#' @param feat_count_column column with count information (optional) #' @param return_gobject return giotto object (default: TRUE) #' @param verbose be verbose +#' @param feat_subset_ids deprecated. Use `feat_subset_values` instead. +#' @param count_info_column deprecated. Use `feat_count_column` instead. #' @returns giotto object or spatVector with overlapping information #' @details Serial overlapping function. #' @concept overlap @@ -605,10 +902,28 @@ calculateOverlapRaster <- function( poly_ID_names = NULL, feat_info = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, return_gobject = TRUE, - verbose = TRUE) { + verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated()) { + deprecate_warn( + when = "0.4.7", + what = "calculateOverlapRaster()", + with = "aggregateFeatures()", + details = "`calculateOverlap()` is another option if only the overlap + step is desired." + ) + feat_subset_values <- GiottoUtils::deprecate_param( + feat_subset_ids, feat_subset_values, + fun = "calculateOverlapRaster", when = "0.4.7" + ) + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "calculateOverlapRaster", when = "0.4.7" + ) + # set defaults if not provided if (is.null(feat_info)) { feat_info <- names(gobject@feat_info)[[1]] @@ -648,8 +963,9 @@ calculateOverlapRaster <- function( # * subset points if needed # e.g. to select transcripts within a z-plane - if (!is.null(feat_subset_column) & !is.null(feat_subset_ids)) { - bool_vector <- pointvec[[feat_subset_column]][[1]] %in% feat_subset_ids + if (!is.null(feat_subset_column) & !is.null(feat_subset_values)) { + bool_vector <- pointvec[[feat_subset_column]][[1]] %in% + feat_subset_values pointvec <- pointvec[bool_vector] } @@ -657,7 +973,7 @@ calculateOverlapRaster <- function( overlap_points <- .calculate_overlap_raster( spatvec = spatvec, pointvec = pointvec, - count_info_column = count_info_column, + keep = feat_count_column, verbose = verbose ) @@ -672,80 +988,86 @@ calculateOverlapRaster <- function( } } +#' @param overlap_data `data.table` of extracted intensity values per poly_ID +#' @noRd +.create_overlap_intensity_dt <- function(overlap_data) { + odt <- new("overlapIntensityDT", data = overlap_data) + odt@nfeats <- ncol(overlap_data) - 1L + odt +} +#' @param x from data (SpatVector) +#' @param y to data (SpatVector) +#' @param overlap_data relationships (data.frame). Expected to be numeric row +#' indices between x and y +#' @param keep additional col(s) in `y` to keep +#' @noRd +.create_overlap_point_dt <- function(x, y, + overlap_data, keep = NULL, feat_ids) { + poly <- feat_idx <- feat <- feat_id_index <- NULL # NSE vars + # cleanup input overlap_data + checkmate::assert_data_frame(overlap_data) + data.table::setDT(overlap_data) + cnames <- colnames(overlap_data) + data.table::setnames(overlap_data, + old = c(cnames[[2]], cnames[[1]]), + new = c("poly", "feat_idx") + ) + # make relationships table sparse by removing non-overlapped features + # these results are indexed by all features, so no need to filter + # non-overlapped polys + overlap_data <- overlap_data[!is.na(poly)] + + # extract needed info from y + keep <- c("feat_ID", "feat_ID_uniq", keep) + ytab <- terra::as.data.frame(y[overlap_data$feat_idx, keep]) + + # initialize overlap object and needed ids + sids <- x$poly_ID + fids <- unique(ytab$feat_ID) + odt <- new("overlapPointDT", + spat_ids = sids, + feat_ids = feat_ids, + nfeats = as.integer(nrow(y)) + ) -#' @name .calculate_overlap_raster -#' @title Find feature points overlapped by rasterized polygon. -#' @description Core workflow function that accepts simple `SpatVector` inputs, -#' performs rasterization of the polys and then checks for overlaps. -#' @param spatvec `SpatVector` polygon from a `giottoPolygon` object -#' @param pointvec `SpatVector` points from a `giottoPoints` object -#' @param count_info_column column with count information (optional) -#' @param verbose be verbose -#' @concept overlap -#' @returns `SpatVector` of overlapped points info -#' @seealso [calculateOverlapRaster()] -#' @keywords internal -.calculate_overlap_raster <- function(spatvec, - pointvec, - count_info_column = NULL, - verbose = TRUE) { - # DT vars - poly_ID <- poly_i <- ID <- x <- y <- feat_ID <- feat_ID_uniq <- NULL - - # spatial vector to raster - if (verbose) GiottoUtils::wrap_msg("1. convert polygon to raster \n") - spatrast_res <- polygon_to_raster(spatvec, field = "poly_ID") - spatrast <- spatrast_res[["raster"]] - ID_vector <- spatrast_res[["ID_vector"]] - - ## overlap between raster and point - if (verbose) GiottoUtils::wrap_msg("2. overlap raster and points \n") - overlap_test <- terra::extract(x = spatrast, y = pointvec) - - # add poly_ID information - if (verbose) GiottoUtils::wrap_msg("3. add polygon information \n") - overlap_test_dt <- data.table::as.data.table(overlap_test) - overlap_test_dt[, poly_ID := ID_vector[poly_i]] - - # add point information - if (verbose) GiottoUtils::wrap_msg("4. add points information \n") - pointvec_dt <- data.table::as.data.table(pointvec, geom = "XY") - - pointvec_dt_x <- pointvec_dt$x - names(pointvec_dt_x) <- pointvec_dt$geom - pointvec_dt_y <- pointvec_dt$y - names(pointvec_dt_y) <- pointvec_dt$geom - pointvec_dt_feat_ID <- pointvec_dt$feat_ID - names(pointvec_dt_feat_ID) <- pointvec_dt$geom - pointvec_dt_feat_ID_uniq <- pointvec_dt$feat_ID_uniq - names(pointvec_dt_feat_ID_uniq) <- pointvec_dt$geom - - overlap_test_dt[, x := pointvec_dt_x[ID]] - overlap_test_dt[, y := pointvec_dt_y[ID]] - overlap_test_dt[, feat_ID := pointvec_dt_feat_ID[ID]] - overlap_test_dt[, feat_ID_uniq := pointvec_dt_feat_ID_uniq[ID]] - - if (!is.null(count_info_column)) { - pointvec_dt_count <- pointvec_dt[[count_info_column]] - names(pointvec_dt_count) <- pointvec_dt$geom - overlap_test_dt[, c(count_info_column) := pointvec_dt_count[ID]] + # Ensure data is stored as integer or integer-based mapping + ## - if poly/feat_idx contents are NOT integer coercible, establish a map # + if (!overlap_data[, checkmate::test_integerish(head(poly, 100))]) { + overlap_data[, poly := match(poly, sids)] + } + if (!overlap_data[, checkmate::test_integerish(head(feat_idx, 100))]) { + overlap_data[, feat_idx := match(feat_idx, fids)] + } + ## -- if still not integer, coerce to integer --------------------------- # + if (!is.integer(overlap_data$poly[1])) { + overlap_data[, poly := as.integer(poly)] + } + if (!is.integer(overlap_data$feat_idx[1])) { + overlap_data[, feat_idx := as.integer(feat_idx)] } - if (verbose) GiottoUtils::wrap_msg("5. create overlap polygon - information \n") - overlap_test_dt_spatvector <- terra::vect( - x = as.matrix(overlap_test_dt[, c("x", "y"), with = FALSE]), - type = "points", - atts = overlap_test_dt[, c( - "poly_ID", "feat_ID", "feat_ID_uniq", - count_info_column - ), with = FALSE] - ) - names(overlap_test_dt_spatvector) <- c( - "poly_ID", "feat_ID", "feat_ID_uniq", count_info_column + # append y attribute info + overlap_data <- cbind(overlap_data, ytab) + data.table::setnames(overlap_data, + old = c("feat_ID_uniq", "feat_ID"), + new = c("feat", "feat_id_index") ) - return(overlap_test_dt_spatvector) + if (!is.integer(overlap_data$feat[1])) { + overlap_data[, feat := as.integer(feat)] + } + # add feat_ID map + overlap_data[, feat_id_index := match(feat_id_index, odt@feat_ids)] + # remove feat_idx which may not be reliable after feature subsets + overlap_data[, feat_idx := NULL] + # set indices + data.table::setkeyv(overlap_data, "feat") + data.table::setindexv(overlap_data, "poly") + data.table::setcolorder(overlap_data, c("poly", "feat", "feat_id_index")) + # add to object + odt@data <- overlap_data + + odt } @@ -863,8 +1185,8 @@ calculateOverlapPolygonImages <- function(gobject, for (img_name in image_names) { if (!img_name %in% potential_large_image_names) { warning( - "image with the name ", img_name, - " was not found and will be skipped \n" + "[calculateOverlap] image with the name ", img_name, + " was not found and will be skipped \n", call. = FALSE ) } } @@ -878,8 +1200,8 @@ calculateOverlapPolygonImages <- function(gobject, if (!img_name %in% potential_large_image_names) { warning( - "image with the name ", img_name, - " was not found and will be skipped \n" + "[calculateOverlap] image with the name ", img_name, + " was not found and will be skipped \n", call. = FALSE ) } @@ -1233,9 +1555,10 @@ calculateOverlapParallel <- function(gobject, #' @param x object containing overlaps info. Can be giotto object or SpatVector #' points or data.table of overlaps generated from `calculateOverlap` #' @param name name for the overlap count matrix -#' @param count_info_column column with count information. If a +#' @param feat_count_column column with count information. If a #' column called "count" is present in the feature points data, it will be #' automatically selected. +#' @param count_info_column deprecated. Use `feat_count_column` instead. #' @param \dots additional params to pass to methods #' @concept overlap #' @returns giotto object or count matrix @@ -1256,81 +1579,77 @@ NULL # * gobject #### #' @rdname overlapToMatrix -#' @param poly_info character. Polygon information to use +#' @param spat_info character. Polygon information to use #' @param feat_info character. Feature information to use #' @param type character. Type of overlap data (either 'point' or 'intensity') #' @param return_gobject return giotto object (default: TRUE) #' @param verbose be verbose +#' @param poly_info deprecated. Please use spat_info. #' @export setMethod( "overlapToMatrix", signature("giotto"), function(x, name = "raw", - poly_info = NULL, + spat_info = NULL, feat_info = NULL, type = c("point", "intensity"), - count_info_column = NULL, - aggr_function = "sum", + feat_count_column = NULL, + fun = "sum", return_gobject = TRUE, verbose = TRUE, + aggr_function = deprecated(), + poly_info = deprecated(), + count_info_column = deprecated(), ...) { + # deprecations + spat_info <- GiottoUtils::deprecate_param(poly_info, spat_info, + fun = "overlapToMatrix", when = "0.4.7" + ) + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "overlapToMatrix", when = "0.4.7" + ) + fun <- GiottoUtils::deprecate_param(aggr_function, fun, + fun = "overlapToMatrix", when = "0.4.7" + ) + type <- match.arg(type, choices = c("point", "intensity")) checkmate::assert_character(name, len = 1L) - if (!is.null(count_info_column) && !isFALSE(count_info_column)) { - checkmate::assert_character(count_info_column, len = 1L) + if (!is.null(feat_count_column) && !isFALSE(feat_count_column)) { + checkmate::assert_character(feat_count_column, len = 1L) } checkmate::assert_logical(return_gobject) - poly_info <- set_default_spat_unit( + spat_info <- set_default_spat_unit( gobject = x, - spat_unit = poly_info + spat_unit = spat_info ) feat_info <- set_default_feat_type( gobject = x, - spat_unit = poly_info, + spat_unit = spat_info, feat_type = feat_info ) # get data gpoly <- getPolygonInfo( gobject = x, - polygon_name = poly_info, + polygon_name = spat_info, return_giottoPolygon = TRUE, verbose = verbose ) o2m_args <- list( x = gpoly, - col_names = spatIDs(x, spat_unit = poly_info), - row_names = featIDs(x, feat_type = feat_info), feat_info = feat_info, - count_info_column = count_info_column, - aggr_function = aggr_function, - # output = 'Matrix', # Do not specify here. methods must return - # something that operates similarly to a [matrix] - # object by default. + feat_count_column = feat_count_column, + output = "exprobj", + sort = TRUE, type = type, - verbose = verbose, ... ) # pass to giottoPolygon method - overlapmatrix <- do.call(overlapToMatrix, args = o2m_args) - - # order matrix row/col - mat_r_names <- rownames(overlapmatrix) - mat_c_names <- colnames(overlapmatrix) - overlapmatrix <- overlapmatrix[ - match(mixedsort(mat_r_names), mat_r_names), - match(mixedsort(mat_c_names), mat_c_names) - ] - - overlapExprObj <- create_expr_obj( - name = name, - exprMat = overlapmatrix, - spat_unit = poly_info, - feat_type = feat_info, - provenance = poly_info - ) + overlapExprObj <- do.call(overlapToMatrix, args = o2m_args) + objName(overlapExprObj) <- name if (isTRUE(return_gobject)) { centroidsDT <- centroids(gpoly) %>% @@ -1345,8 +1664,8 @@ setMethod( spatlocs <- createSpatLocsObj( coordinates = centroidsDT_loc, name = name, - spat_unit = poly_info, - provenance = poly_info, + spat_unit = spat_info, + provenance = spat_info, verbose = FALSE ) @@ -1378,15 +1697,21 @@ setMethod( # * giottoPolygon #### #' @rdname overlapToMatrix -#' @param output data format/class to return the results as +#' @param output data format/class to return the results as. Default is "Matrix" #' @export setMethod( "overlapToMatrix", signature("giottoPolygon"), function(x, feat_info = "rna", type = c("point", "intensity"), - count_info_column = NULL, - output = c("Matrix", "data.table"), + feat_count_column = NULL, + count_info_column = deprecated(), ...) { + # deprecations + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "overlapToMatrix", when = "0.4.7" + ) + type <- match.arg(type, choices = c("point", "intensity")) overlaps_data <- switch(type, @@ -1404,14 +1729,13 @@ setMethod( argslist <- list( x = overlaps_data, - count_info_column = count_info_column, - output = output, + feat_count_column = feat_count_column, ... ) # remove args not accepted by specific method if (type == "intensity") { - argslist$count_info_column <- NULL + argslist$feat_count_column <- NULL argslist$col_names <- NULL argslist$row_names <- NULL argslist$verbose <- NULL @@ -1435,10 +1759,17 @@ setMethod( "overlapToMatrix", signature("SpatVector"), function(x, col_names = NULL, row_names = NULL, - count_info_column = NULL, + feat_count_column = NULL, output = c("Matrix", "data.table"), verbose = TRUE, + count_info_column = deprecated(), ...) { + # deprecations + feat_count_column <- GiottoUtils::deprecate_param( + count_info_column, feat_count_column, + fun = "overlapToMatrix", when = "0.4.7" + ) + output <- match.arg( toupper(output), choices = c("MATRIX", "DATA.TABLE") @@ -1454,33 +1785,32 @@ setMethod( # 2. Perform aggregation to counts DT - # autodetect counts col - if ("count" %in% names(dtoverlap) && is.null(count_info_column)) { + if ("count" %in% names(dtoverlap) && is.null(feat_count_column)) { vmsg(.v = verbose, "[overlap] Found column \"count\" in feature info. - - Using as `count_info_column` - [!] Set count_info_column = FALSE to disable.") - count_info_column <- "count" + - Using as `feat_count_column` + [!] Set feat_count_column = FALSE to disable.") + feat_count_column <- "count" } - if (isFALSE(count_info_column)) { - count_info_column <- NULL + if (isFALSE(feat_count_column)) { + feat_count_column <- NULL } - if (!is.null(count_info_column)) { # if there is a counts col + if (!is.null(feat_count_column)) { # if there is a counts col selected - if (!count_info_column %in% colnames(dtoverlap)) { - .gstop("count_info_column ", count_info_column, + if (!feat_count_column %in% colnames(dtoverlap)) { + .gstop("feat_count_column ", feat_count_column, " does not exist", .n = 2L ) } # aggregate counts of features - dtoverlap[, c(count_info_column) := as.numeric( - get(count_info_column) + dtoverlap[, c(feat_count_column) := as.numeric( + get(feat_count_column) )] - aggr_dtoverlap <- dtoverlap[, base::sum(get(count_info_column)), + aggr_dtoverlap <- dtoverlap[, base::sum(get(feat_count_column)), by = c("poly_ID", "feat_ID") ] data.table::setnames(aggr_dtoverlap, "V1", "N") @@ -1551,12 +1881,20 @@ setMethod( # * data.frame #### # images #' @rdname overlapToMatrix -#' @param aggr_function function to aggregate image information (default = sum) +#' @param fun character. Function to aggregate image information +#' (default = "sum") +#' @param aggr_function deprecated. Use `fun` instead. #' @export setMethod( "overlapToMatrix", signature("data.table"), function(x, - aggr_function = "sum", - output = c("Matrix", "data.table")) { + fun = "sum", + output = c("Matrix", "data.table"), + aggr_function = deprecated()) { + + fun <- GiottoUtils::deprecate_param(aggr_function, fun, + fun = "overlapToMatrix", when = "0.4.7" + ) + output <- match.arg( toupper(output), choices = c("MATRIX", "DATA.TABLE") @@ -1571,7 +1909,7 @@ setMethod( variable.name = "feat_ID" ) - aggr_fun <- get(aggr_function) + aggr_fun <- get(fun) aggr_comb <- melt_image_info[, aggr_fun(value), by = .(poly_ID, feat_ID) ] @@ -1593,7 +1931,58 @@ setMethod( } ) +# * overlapPointDT #### +#' @rdname overlapToMatrix +#' @param sort logical (default = TRUE). Whether to perform a mixed sort on +#' output matrix row and col names. +#' @export +setMethod("overlapToMatrix", signature("overlapPointDT"), + function(x, + name = "raw", + sort = TRUE, + feat_count_column = NULL, + output = c("Matrix", "exprObj"), + ...) { + output <- match.arg(tolower(output), choices = c("matrix", "exprobj")) + m <- as.matrix(x, feat_count_column = feat_count_column, ...) + if (isTRUE(sort)) m <- .mixedsort_rowcols(m) + switch(output, + "matrix" = m, + "exprobj" = createExprObj( + expression_data = m, + name = name, + spat_unit = spatUnit(x), + feat_type = featType(x), + provenance = prov(x) + ) + ) +}) + +# * overlapIntensityDT #### + +#' @rdname overlapToMatrix +#' @export +setMethod("overlapToMatrix", signature("overlapIntensityDT"), + function(x, + name = "raw", + sort = TRUE, + output = c("Matrix", "exprObj"), + ...) { + output <- match.arg(tolower(output), choices = c("matrix", "exprobj")) + m <- as.matrix(x, ...) + if (isTRUE(sort)) m <- .mixedsort_rowcols(m) + switch(output, + "matrix" = m, + "exprobj" = createExprObj( + expression_data = m, + name = name, + spat_unit = spatUnit(x), + feat_type = featType(x), + provenance = prov(x) + ) + ) +}) @@ -1719,7 +2108,16 @@ overlapToMatrixMultiPoly <- function(gobject, } - +.mixedsort_rowcols <- function(x) { + # order matrix row/col + mat_r_names <- rownames(x) + mat_c_names <- colnames(x) + x[ + match(mixedsort(mat_r_names), mat_r_names), + match(mixedsort(mat_c_names), mat_c_names), + drop = FALSE + ] +} #' @title overlapImagesToMatrix @@ -1748,6 +2146,15 @@ overlapImagesToMatrix <- function(gobject, image_names = NULL, spat_locs_name = "raw", return_gobject = TRUE) { + + deprecate_warn( + when = "0.4.7", + what = "overlapImagesToMatrix()", + with = "aggregateFeatures()", + details = "`overlapToMatrix()` is another option if only the matrix + construction from overlaps information step is desired." + ) + # data.table vars value <- poly_ID <- feat_ID <- x <- y <- NULL @@ -2304,24 +2711,24 @@ aggregateStacksPolygonOverlaps <- function(gobject, for (i in seq_len(length(spat_units))) { spat_unit <- spat_units[i] - vecDT <- gobject@spatial_info[[spat_unit]]@overlaps[[feat_type]] + ovlp <- getPolygonInfo(gobject, + polygon_name = spat_unit, + polygon_overlap = feat_type + ) + # vecDT <- gobject@spatial_info[[spat_unit]]@overlaps[[feat_type]] - if (!is.null(vecDT)) { - vecDT <- .spatvector_to_dt(vecDT) - vecDT[, "stack" := i] - polygon_list[[spat_unit]] <- vecDT + if (!is.null(ovlp)) { + ovlp@data[, "stack" := i] + polygon_list[[spat_unit]] <- ovlp } } if (length(polygon_list) == 0) { wrap_msg("No feature overlaps found for stack aggregation \n") } else { - polygon_DT <- data.table::rbindlist(polygon_list) - polygon <- .dt_to_spatvector_points( - dt = polygon_DT, - include_values = TRUE - ) - gobject@spatial_info[[new_spat_unit]]@overlaps[[feat_type]] <- polygon + comb_ovlp <- do.call(rbind, polygon_list) + gobject@spatial_info[[new_spat_unit]]@overlaps[[feat_type]] <- + comb_ovlp } return(gobject) diff --git a/R/classes.R b/R/classes.R index b8c13fe9..af14a9e9 100644 --- a/R/classes.R +++ b/R/classes.R @@ -489,6 +489,48 @@ updateGiottoObject <- function(gobject) { # -------------------------------------------------------------------------# + # subobject updates + if (!is.null(attr(gobject, "feat_info"))) { + info_list <- get_feature_info_list(gobject) + # update S4 object if needed + info_list <- lapply(info_list, function(info) { + try_val <- try(validObject(info), silent = TRUE) + if (inherits(try_val, "try-error")) { + info <- updateGiottoPointsObject(info) + } + return(info) + }) + ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### + gobject <- setFeatureInfo( + gobject = gobject, + x = info_list, + verbose = FALSE, + initialize = FALSE + ) + ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### + } + if (!is.null(attr(gobject, "spatial_info"))) { + info_list <- get_polygon_info_list(gobject) + # update S4 object if needed + info_list <- lapply(info_list, function(info) { + try_val <- try(validObject(info), silent = TRUE) + if (inherits(try_val, "try-error") || + .gversion(gobject) <= "0.4.7") { + info <- updateGiottoPolygonObject(info) + } + return(info) + }) + ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### + gobject <- setPolygonInfo( + gobject = gobject, + x = info_list, + verbose = FALSE, + centroids_to_spatlocs = FALSE, + initialize = FALSE + ) + ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### + } + # finally, set updated version number .gversion(gobject) <- packageVersion("GiottoClass") @@ -1407,7 +1449,7 @@ giottoPolygon <- setClass( #' @export updateGiottoPolygonObject <- function(gpoly) { if (!inherits(gpoly, "giottoPolygon")) { - stop("This function is only for giottoPoints") + stop("This function is only for giottoPolygon") } # 3.2.X adds cacheing of IDs @@ -1417,6 +1459,10 @@ updateGiottoPolygonObject <- function(gpoly) { ) } + # 0.4.7 changes overlaps representation + # intersection `SpatVector` -> `overlapInfo`-inheriting classes + gpoly <- .update_overlaps(gpoly) + gpoly } @@ -1536,10 +1582,185 @@ setClass( +## overlapInfo #### +setClass("overlapInfo", + contains = c("spatFeatData", "VIRTUAL"), + slots = list(data = "ANY") +) +setClass("overlapPoint", contains = c("overlapInfo", "VIRTUAL")) +setClass("overlapIntensity", contains = c("overlapInfo", "VIRTUAL")) + +#' @name overlapPointDT-class +#' @title Polygon and Point Relationships +#' @description +#' Utility class for storing overlaps relationships between polygons and points +#' in a sparse `data.table` format. Retrieve the unique ID index of overlapped +#' points `[i, ]`. Get indices of which polys are overlapping specific feature +#' species using `[, j]`. +#' +#' Subsetting with `ids = FALSE` and `[i, j]` indexing is also supported. +#' +#' Supports `as.matrix` for conversion to `dgCMatrix`. Contained poly and +#' feature names simplify rownames/colnames and empty row/col creation. +#' +#' @slot data data.table. Table containing 3 integer cols: +#' +#' * `poly` - polygon index. Maps to `spat_ids` slot. +#' * `feat` - feat_ID_uniq (unique integer identifier) of a point detection +#' * `feat_id_index` - index of feature name mapping in `@feat_ids` slot. +#' @slot spat_unit character. Spatial unit (usually name of polygons information) +#' @slot feat_type character. Feature type (usually name of points information) +#' @slot provenance character. provenance information +#' @slot spat_ids character. Polygon names +#' @slot feat_ids character. Feature names +#' @slot nfeats integer (optional metadata). How many feature points were +#' used in overlap operation. Gives an idea of sparsity, but has no effect on +#' processing. +#' +#' @param x object +#' @param i numeric, character, logical. Index of or name of poly in overlapping +#' polygons +#' @param j numeric, character, logical. Index of or name of feature being +#' overlapped. +#' @param use_names logical (default = `FALSE`). Whether to return as integer +#' indices or with character ids. +#' @param ids logical (default = `TRUE`). Whether to return the requested +#' integer indices (`TRUE`) or the subset overlap object (`FALSE`). +#' @param drop not used. +#' @param \dots additional params to pass (none implemented) +#' @returns integer or character if only `i` or `j` provided, depending on +#' `use_names`. A subset `overlapPointDT` if both `i` and `j` are used. +#' @examples +#' g <- GiottoData::loadGiottoMini("vizgen") +#' poly <- g[["spatial_info", "z0"]][[1]] +#' ovlp <- overlaps(poly, "rna") +#' ovlp +#' +#' as.matrix(ovlp) +#' +#' dim(ovlp) +#' nrow(ovlp) # number of relationships +#' +#' # get feature unique IDs overlapped by nth poly +#' ovlp[1] # check one (no overlaps returns integer(0)) +#' ovlp[1:5] # check multiple +#' ovlp[1:5, use_names = TRUE] # returns feature names, but no longer unique +#' +#' # get integer index of poly(s) overlapping particular feature species +#' ovlp[, 1] +#' ovlp[, "Mlc1"] # this is the same +#' +#' # get a subset of overlap object +#' ovlp[1:10, ids = FALSE] # subset to first 10 polys +#' ovlp[, 1:10, ids = FALSE] # subset to first 10 feature species +#' ovlp[1:10, 1:10] # subset to first 10 polys and first 10 features species +#' @exportClass overlapPointDT +setClass("overlapPointDT", + contains = "overlapPoint", + slots = list( + spat_ids = "character", # spat_ids are unique, no need to record npoly + feat_ids = "character", + nfeats = "integer" + ) +) +setClass("overlapIntensityDT", + contains = "overlapIntensity", + slots = list( + nfeats = "integer", # skip spat_ids/feat_ids. This data is not sparse + fun = "character" + ) +) +# update old overlaps information to new `overlapInfo` +.update_overlaps <- function(x, ...) { + + if (inherits(x, "giottoPolygon")) { + res <- .update_overlaps(x@overlaps, + poly_ids = x$poly_ID, + spat_unit = spatUnit(x), + ... + ) + x@overlaps <- res + return(x) + } + if (inherits(x, "list")) { + list_names <- names(x) + list_res <- lapply(list_names, function(ovlp_name) { + if (ovlp_name == "intensity") { + # recurse over list of intensity overlaps (can't set feat_type) + intensity_res <- .update_overlaps(x[["intensity"]], ...) + return(intensity_res) + } else { + updated_res <- .update_overlaps(x[[ovlp_name]], + feat_type = ovlp_name, + ... + ) + return(updated_res) + } + }) + names(list_res) <- list_names + return(list_res) + } + if (inherits(x, "SpatVector")) { + return(.update_overlaps_points(x, ...)) + } + if (inherits(x, "data.table")) { + return(.update_overlaps_intensity(x, ...)) + } + x # allow passthrough if not matching either signature +} + +#' @param x (SpatVector) the old overlaps representation to convert +#' @param poly_ids the `spatIDs()` of the giottoPolygon. This is since the +#' ordering of the polygon IDs within the overlaps data usually does not match. +#' @param spat_unit,feat_type spat_unit / feat_type +#' @noRd +.update_overlaps_points <- function(x, poly_ids, + spat_unit = NA_character_, feat_type = NA_character_, ...) { + checkmate::assert_character(poly_ids) + data <- terra::as.data.frame(x) + data.table::setDT(data) + odt <- new("overlapPointDT", + spat_unit = spat_unit, + feat_type = feat_type, + provenance = spat_unit, + nfeats = as.integer(nrow(data)) + ) + odt@feat_ids <- unique(data$feat_ID) + odt@spat_ids <- unique(poly_ids) + data[, feat := as.integer(feat_ID_uniq)] + data[, feat_ID_uniq := NULL] + data.table::setnames(data, + old = c("poly_ID", "feat_ID"), + new = c("poly", "feat_id_index") + ) + data <- data[!is.na(poly) & !is.na(feat),] # drop NAs + # Ensure data is stored as integer-based mapping + data[, poly := match(poly, odt@spat_ids)] + data[, feat_id_index := match(feat_id_index, odt@feat_ids)] + data.table::setkeyv(data, "feat") + data.table::setindexv(data, "poly") + data.table::setcolorder(data, c("poly", "feat", "feat_id_index")) + # add to object + odt@data <- data + odt +} +.update_overlaps_intensity <- function(x, + spat_unit = NA_character_, feat_type = NA_character_, ...) { + odt <- new("overlapIntensityDT", + spat_unit = spat_unit, + feat_type = feat_type, + provenance = spat_unit, + nfeats = as.integer(ncol(x) - 1) + ) + fids <- setdiff(names(x), "poly_ID") + x <- x[, lapply(.SD, sum), by = "poly_ID", .SDcols = fids] + odt@data <- x + odt +} ## featureNetwork class #### diff --git a/R/combine_metadata.R b/R/combine_metadata.R index bd40b029..40e5f49a 100644 --- a/R/combine_metadata.R +++ b/R/combine_metadata.R @@ -475,15 +475,21 @@ combineFeatureOverlapData <- function(gobject, # overlap poly and feat info poly_list <- list() for (poly in poly_info) { - feat_overlap_info_spatvec <- getPolygonInfo( + feat_overlap <- getPolygonInfo( gobject = gobject, polygon_name = poly, polygon_overlap = feat ) - feat_overlap_info <- .spatvector_to_dt( - feat_overlap_info_spatvec + + pts <- getFeatureInfo(gobject, + feat_type = feat_type, + return_giottoPoints = FALSE ) + # extract overlapped points + pts <- pts[pts$feat_ID_uniq %in% feat_overlap@data$feat,] + feat_overlap_info <- .spatvector_to_dt(pts) + if (!is.null(sel_feats[[feat]])) { selected_features <- sel_feats[[feat]] feat_overlap_info <- feat_overlap_info[ diff --git a/R/methods-coerce.R b/R/methods-coerce.R index 512688d9..83636ccd 100644 --- a/R/methods-coerce.R +++ b/R/methods-coerce.R @@ -31,6 +31,9 @@ NULL #' @description Coerce to matrix #' @param x object to coerce #' @param id_rownames logical. Retain the spatial IDs as the rownames +#' @param feat_count_column character. If provided, column in overlaps info +#' that contains count information to take into account when generating matrix. +#' This is important when point detections represent more than one count. #' @param \dots additional params to pass (none implemented) #' @examples #' sl <- GiottoData::loadSubObjectMini("spatLocsObj") @@ -147,8 +150,23 @@ as.data.table.giottoPoints <- function(x, ...) { as.data.table(x[], geomtype = "points", ...) } +#' @rdname as.data.table +#' @method as.data.frame overlapPointDT +#' @export +as.data.frame.overlapPointDT <- function(x, ...) { + data.frame( + poly_ID = x@spat_ids[x@data$poly], + feat_ID = x@feat_ids[x@data$feat_id_index], + feat_ID_uniq = x@data$feat # these map from feat_ID_uniq + ) +} - +#' @rdname as.data.table +#' @method as.data.frame overlapIntensityDT +#' @export +as.data.frame.overlapIntensityDT <- function(x, ...) { + x[] +} # to matrix #### @@ -169,6 +187,44 @@ setMethod("as.matrix", signature("spatLocsObj"), function(x, id_rownames = TRUE, return(m) }) +# not using `setAs()` since we need more params +#' @rdname as.matrix +#' @export +setMethod("as.matrix", signature("overlapPointDT"), + function(x, feat_count_column = NULL, ...) { + if (!is.null(feat_count_column)) { + feat_count <- x@data[[feat_count_column]] + } else { + feat_count <- 1 + } + Matrix::sparseMatrix( + i = x@data$feat_id_index, + j = x@data$poly, + x = feat_count, + dimnames = list( + x@feat_ids, + x@spat_ids + ), + dims = c( + length(x@feat_ids), + length(x@spat_ids) + ), + ... + ) + }) + +#' @rdname as.matrix +#' @export +setMethod("as.matrix", signature("overlapIntensityDT"), function(x, ...) { + m <- Matrix::Matrix(as.matrix( + x[][,-1]), + dimnames = list( + x[]$poly_ID, + colnames(x[][, -1]) + ) + ) + t(m) +}) #' @rdname as.matrix #' @param attr Either NULL or `character` providing the name of an edge diff --git a/R/methods-crop.R b/R/methods-crop.R index a0467508..6691fbc3 100644 --- a/R/methods-crop.R +++ b/R/methods-crop.R @@ -233,13 +233,7 @@ setMethod( # iterate through all overlaps, removing cell_ids that were removed in the # crop. - for (feat in names(x@overlaps)) { - cell_id_bool <- terra::as.list( - x@overlaps[[feat]] - )$poly_ID %in% x@unique_ID_cache - x@overlaps[[feat]] <- x@overlaps[[feat]][cell_id_bool] - } - + x <- .subset_overlaps_poly(x, x@unique_ID_cache) x } ) diff --git a/R/methods-dims.R b/R/methods-dims.R index 3c5d36af..ee31ee3c 100644 --- a/R/methods-dims.R +++ b/R/methods-dims.R @@ -69,7 +69,13 @@ setMethod("nrow", signature("enrData"), function(x) nrow(x@enrichDT)) #' @export setMethod("nrow", signature("dimObj"), function(x) nrow(x@coordinates)) +#' @rdname dims-generic +#' @export +setMethod("nrow", signature("overlapPointDT"), function(x) nrow(x@data)) +#' @rdname dims-generic +#' @export +setMethod("nrow", signature("overlapIntensityDT"), function(x) nrow(x@data)) # ncol #### @@ -103,8 +109,13 @@ setMethod("ncol", signature("enrData"), function(x) ncol(x@enrichDT)) #' @export setMethod("ncol", signature("dimObj"), function(x) ncol(x@coordinates)) +#' @rdname dims-generic +#' @export +setMethod("ncol", signature("overlapPointDT"), function(x) ncol(x@data)) - +#' @rdname dims-generic +#' @export +setMethod("ncol", signature("overlapIntensityDT"), function(x) ncol(x@data)) ## dim() generic #### @@ -145,3 +156,11 @@ setMethod("dim", signature("giottoPolygon"), function(x) dim(x[])) #' @rdname dims-generic #' @export setMethod("dim", signature("giottoPoints"), function(x) dim(x[])) + +#' @rdname dims-generic +#' @export +setMethod("dim", signature("overlapPointDT"), function(x) dim(x@data)) + +#' @rdname dims-generic +#' @export +setMethod("dim", signature("overlapIntensityDT"), function(x) dim(x@data)) diff --git a/R/methods-extract.R b/R/methods-extract.R index cbd45936..3417a503 100644 --- a/R/methods-extract.R +++ b/R/methods-extract.R @@ -1059,17 +1059,151 @@ setMethod( if (is.null(x@overlaps)) { return(x) } # if no overlaps, skip following + x <- .subset_overlaps_poly(x, i = i) + x + } +) - for (feat in names(x@overlaps)) { - cell_id_bool <- terra::as.list( - x@overlaps[[feat]] - )$poly_ID %in% x@unique_ID_cache - x@overlaps[[feat]] <- x@overlaps[[feat]][cell_id_bool] +#' @description +#' Perform a feature subset to keep only those features referenced by `i` +#' across all `feat_type` requested. +#' ** If subsetting for paired sets of `i` and `feat_type`, then this function +#' should be looped. +#' @param x giottoPolygon to use +#' @param i features (within `feat_type`) to keep +#' (character, numeric, or logical) +#' @param feat_type feature types to subset features within +#' @examples +#' # this is pseudocode # +#' overlaps(gpoly) +#' gpoly <- .subset_overlaps_feat(gpoly, c("Selplg", "Mlc1"), "rna") +#' overlaps(gpoly) +#' @noRd +.subset_overlaps_feat <- function(x, i, feat_type) { + if (length(x@overlaps) == 0L) return(x) # return early if no overlaps + checkmate::assert_character(feat_type) + is_pnt_feat <- feat_type %in% names(x@overlaps) + is_int_feat <- feat_type %in% names(x@overlaps$intensity) + # return early if no relevant feat_types + if (sum(is_pnt_feat, is_int_feat) == 0) return(x) + + .do_subset <- function(o, ftypes) { + for (feat in ftypes) { + if (!inherits(o[[feat]], "overlapInfo")) { + warning("[subset overlaps feats] unrecognized overlap type") + } + o[[feat]] <- o[[feat]][, i, ids = FALSE] } + o + } - x + if (any(is_pnt_feat)) { # point overlaps + x@overlaps <- .do_subset(x@overlaps, + ftypes = feat_type[is_pnt_feat] + ) } -) + if (any(is_int_feat)) { # intensity overlaps + x@overlaps$intensity <- .do_subset(x@overlaps$intensity, + ftypes = feat_type[is_int_feat] + ) + } + + x +} + +#' @param x giottoPolygon to use +#' @param i polys to keep (character, numeric, or logical) +#' @examples +#' # this is pseudocode # +#' overlaps(gpoly) +#' gpoly <- .subset_overlaps_poly(gpoly, 1:10) +#' overlaps(gpoly) +#' @noRd +.subset_overlaps_poly <- function(x, i) { + if (length(x@overlaps) == 0L) return(x) # return early if no overlaps + for (feat in names(x@overlaps)) { + if (feat == "intensity") { + for (feat_intensity in names(x@overlaps$intensity)) { + x@overlaps$intensity[[feat_intensity]] <- + x@overlaps$intensity[[feat_intensity]][i] + } + next + } + x@overlaps[[feat]] <- x@overlaps[[feat]][i, ids = FALSE] + } + x +} + +# i should be poly ids to keep +.subset_overlap_point_dt_i <- function(x, i) { + if (is.numeric(i) || is.logical(i)) { + i <- x@spat_ids[i] + } + + poly <- NULL # NSE vars + idx <- match(i, x@spat_ids) # poly indices to keep + idx <- idx[!is.na(idx)] # drop unmatched NAs + x@spat_ids <- x@spat_ids[x@spat_ids %in% i] # replace spatial ids + + x@data <- x@data[poly %in% idx] + x@data[, poly := match(poly, idx)] + data.table::setkeyv(x@data, "feat") + data.table::setindex(x@data, "poly") + x +} + +.subset_overlap_point_dt_j <- function(x, j) { + # ---- convert j to numerical index ---- # + if (is.logical(j)) { + if (length(j) != length(x@feat_ids)) { + # recycle logical if needed + j <- rep(j, length.out = length(x@feat_ids)) + } + j <- which(j) + } + if (is.character(j)) { + j <- match(j, x@feat_ids) + } + + x@feat_ids <- x@feat_ids[j] # replace feature ids + + # subset on feat_id_index matches + x@data <- x@data[feat_id_index %in% j] + x@data[, feat_id_index := match(feat_id_index, j)] + data.table::setkeyv(x@data, "feat") + data.table::setindex(x@data, "poly") + x +} + +.subset_overlap_intensity_dt_i <- function(x, i) { + if (is.character(i)) { + x@data <- x@data[match(i, poly_ID)] + return(x) + } + + poly_ID <- NULL # NSE vars + nr <- nrow(x@data) + if (is.logical(i) && length(i) != nr) { + i <- rep(i, length.out = nr) # manual logical recycling for DT + } + x@data <- x@data[i,] + x +} + +.subset_overlap_intensity_dt_j <- function(x, j) { + # convert j to char col reference + if (is.numeric(j)) { + j <- j + 1L + } else if (is.logical(j)) { + j <- c(TRUE, j) + } + if (!is.character(j)) { + j <- colnames(x@data)[j] + } + j <- unique(c("poly_ID", j)) + x@data <- x@data[, .SD, .SDcols = j] + x +} # this behavior is different from normal spatvectors # SpatVector defaults to col subsetting when character is provided to i @@ -1142,6 +1276,147 @@ setMethod( } ) + + + +# * overlapPointDT #### + +#' @rdname overlapPointDT-class +#' @export +setMethod("[", + signature( + x = "overlapPointDT", + i = "gIndex", + j = "missing", + drop = "missing" + ), function(x, i, j, ..., use_names = FALSE, ids = TRUE, drop) { + if (!isTRUE(ids)) { + res <- .subset_overlap_point_dt_i(x, i) + return(res) + } + .select_overlap_point_dt_i(x, i, use_names = use_names) + } +) + +#' @rdname overlapPointDT-class +#' @export +setMethod("[", + signature( + x = "overlapPointDT", + i = "missing", + j = "gIndex", + drop = "missing" + ), function(x, i, j, ..., use_names = FALSE, ids = TRUE, drop) { + if (!isTRUE(ids)) { + res <- .subset_overlap_point_dt_j(x, j) + return(res) + } + .select_overlap_point_dt_j(x, j, use_names = use_names) + } +) + +#' @rdname overlapPointDT-class +#' @export +setMethod("[", + signature( + x = "overlapPointDT", + i = "gIndex", + j = "gIndex", + drop = "missing" + ), + function(x, i, j, ..., use_names = FALSE, drop) { + x <- .subset_overlap_point_dt_j(x, j) + .subset_overlap_point_dt_i(x, i) + } +) + +.select_overlap_point_dt_i <- function(x, i, use_names = FALSE) { + # coerce inputs to numeric poly indexing + if (is.character(i)) { + i <- unique(i) + i <- match(i, x@spat_ids) + } else if (is.logical(i)) { + npoly <- length(x@spat_ids) + if (length(i) != npoly) { + i <- rep(i, length.out = npoly) + } + i <- which(i) + } + # numeric indexing only below + if (isTRUE(use_names)) { + return(x@feat_ids[x@data[poly %in% i, feat_id_index]]) + } + x@data[poly %in% i, feat] +} + +.select_overlap_point_dt_j <- function(x, j, use_names = FALSE) { + # coerce inputs to numeric feat_id_index + if (is.character(j)) { + j <- unique(j) + j <- match(j, x@feat_ids) + } else if (is.logical(j)) { + nfeat <- length(x@feat_ids) + if (length(j) != nfeat) { + j <- rep(j, length.out = nfeat) + } + j <- which(j) + } + # numeric indexing only below + pids <- x@data[feat_id_index %in% j, unique(poly)] + if (isTRUE(use_names)) { + return(x@spat_ids[pids]) + } + pids +} + +# * overlapIntensityDT #### + +setMethod("[", + signature( + x = "overlapIntensityDT", + i = "missing", + j = "missing", + drop = "missing" + ), + function(x, i, j, ..., drop) { + x@data + } +) + +setMethod("[", + signature( + x = "overlapIntensityDT", + i = "gIndex", + j = "missing", + drop = "missing" + ), + function(x, i, j, ..., drop) { + .subset_overlap_intensity_dt_i(x, i) + } +) + +setMethod("[", + signature(x = "overlapIntensityDT", + i = "missing", + j = "gIndex", + drop = "missing"), + function(x, i, j, ..., drop) { + .subset_overlap_intensity_dt_j(x, j) +}) + +setMethod("[", + signature(x = "overlapIntensityDT", + i = "gIndex", + j = "gIndex", + drop = "missing"), +function(x, i, j, ..., drop) { + .subset_overlap_intensity_dt_j(x, j) + .subset_overlap_intensity_dt_i(x, i) +}) + + + + # * giottoLargeImage #### #' @rdname subset_bracket #' @export diff --git a/R/methods-initialize.R b/R/methods-initialize.R index 48082055..4b7386bb 100644 --- a/R/methods-initialize.R +++ b/R/methods-initialize.R @@ -409,58 +409,6 @@ setMethod("initialize", signature("giottoAffineImage"), function(.Object, ...) { avail_se <- list_spatial_enrichments(.Object) - ## Perform any subobject updates ## - ## ----------------------------- ## - - # Feature Info # - if (!is.null(avail_fi)) { - info_list <- get_feature_info_list(.Object) - # update S4 object if needed - info_list <- lapply(info_list, function(info) { - try_val <- try(validObject(info), silent = TRUE) - if (inherits(try_val, "try-error")) { - info <- updateGiottoPointsObject(info) - } - return(info) - }) - ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### - .Object <- setFeatureInfo( - gobject = .Object, - x = info_list, - verbose = FALSE, - initialize = FALSE - ) - ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### - } - - - # Spatial Info # - if (!is.null(avail_si)) { - info_list <- get_polygon_info_list(.Object) - - # update S4 object if needed - info_list <- lapply(info_list, function(info) { - try_val <- try(validObject(info), silent = TRUE) - if (inherits(try_val, "try-error")) { - info <- updateGiottoPolygonObject(info) - } - return(info) - }) - ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### - .Object <- setPolygonInfo( - gobject = .Object, - x = info_list, - verbose = FALSE, - centroids_to_spatlocs = FALSE, - initialize = FALSE - ) - ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### - } - - - - - ## Set active/default spat_unit and feat_type ## ## ------------------------------------------ ## diff --git a/R/methods-overlaps.R b/R/methods-overlaps.R index 33510f25..cbaba19f 100644 --- a/R/methods-overlaps.R +++ b/R/methods-overlaps.R @@ -27,7 +27,10 @@ setMethod( return(x@overlaps) } else { # return named entry - return(x@overlaps[[name]]) + o <- x@overlaps[[name]] + # if still null, look in intensity overlaps + if (is.null(o)) o <- x@overlaps$intensity[[name]] } + o } ) diff --git a/R/methods-rbind.R b/R/methods-rbind.R index d31a47d7..51f8d1b4 100644 --- a/R/methods-rbind.R +++ b/R/methods-rbind.R @@ -90,6 +90,36 @@ setMethod( } ) +#' @rdname rbind-generic +#' @export +setMethod("rbind2", signature("overlapPointDT", "overlapPointDT"), + function(x, y, ...) { + comb_spat <- unique(c(x@spat_ids, y@spat_ids)) + comb_feat <- unique(c(x@feat_ids, y@feat_ids)) + + x_spat_map <- match(x@spat_ids, comb_spat) + y_spat_map <- match(y@spat_ids, comb_spat) + x_feat_map <- match(x@feat_ids, comb_feat) + y_feat_map <- match(y@feat_ids, comb_feat) + + # replace id dictionaries for x (output object) + x@spat_ids <- comb_spat + x@feat_ids <- comb_feat + + # remap indices + x@data[, poly := x_spat_map[poly]] + y@data[, poly := y_spat_map[poly]] + x@data[, feat_id_index := x_feat_map[feat_id_index]] + y@data[, feat_id_index := y_feat_map[feat_id_index]] + + x@data <- rbind(x@data, y@data) + x@nfeats <- x@nfeats + y@nfeats + + data.table::setkeyv(x@data, "feat") + data.table::setindex(x@data, "poly") + x + } +) if (!isGeneric("rbind")) setGeneric("rbind", signature = "...") @@ -113,7 +143,14 @@ setMethod("rbind", "spatLocsObj", function(..., deparse.level = 1) { } }) - +setMethod("rbind", "overlapPointDT", function(..., deparse.level = 1) { + if (nargs() <= 2L) { + rbind2(...) + } else { + xs <- list(...) + rbind2(xs[[1L]], do.call(Recall, xs[-1L])) + } +}) # internals #### @@ -170,10 +207,9 @@ rbind2_giotto_polygon_homo <- function(x, y) { ) } - if (is.null(slot(x, "overlaps"))) { - slot(x, "overlaps") <- slot(y, "overlaps") - } else { - slot(x, "overlaps") <- rbind(slot(x, "overlaps"), slot(y, "overlaps")) + if (!is.null(slot(x, "overlaps")) || !is.null(slot(y, "overlaps"))) { + .rbind2_giotto_polygon_overlap_message() + slot(x, "overlaps") <- NULL } slot(x, "unique_ID_cache") <- unique(c(spatIDs(x), spatIDs(y))) @@ -201,7 +237,8 @@ rbind2_giotto_polygon_hetero <- function(x, y, new_name, add_list_ID = TRUE) { return(gpoly) } - null_xsv <- null_xsvc <- null_xovlp <- FALSE + # init flags + null_xsv <- null_xsvc <- FALSE # Add list_ID if (!is.null(slot(x, "spatVector"))) { @@ -230,19 +267,6 @@ rbind2_giotto_polygon_hetero <- function(x, y, new_name, add_list_ID = TRUE) { } } - if (!is.null(slot(x, "overlaps"))) { - if (!"list_ID" %in% names(slot(x, "overlaps"))) { - slot(x, "overlaps")$list_ID <- slot(x, "name") - } - } else { - null_xovlp <- TRUE - } - if (!is.null(y@overlaps)) { - if (!"list_ID" %in% names(slot(y, "overlaps"))) { - slot(y, "overlaps")$list_ID <- slot(y, "name") - } - } - # Perform rbinds if (isTRUE(null_xsv)) { new_sv <- slot(y, "spatVector") @@ -259,18 +283,21 @@ rbind2_giotto_polygon_hetero <- function(x, y, new_name, add_list_ID = TRUE) { ) } - if (isTRUE(null_xovlp)) { - new_ovlp <- slot(y, "overlaps") - } else { - new_ovlp <- rbind(slot(x, "overlaps"), slot(y, "overlaps")) + if (!is.null(slot(x, "overlaps")) || !is.null(slot(y, "overlaps"))) { + .rbind2_giotto_polygon_overlap_message() } new_poly <- create_giotto_polygon_object( name = new_name, spatVector = new_sv, spatVectorCentroids = new_svc, - overlaps = new_ovlp, + overlaps = NULL, unique_IDs = unique(c(spatIDs(x), spatIDs(y))) ) new_poly } + +.rbind2_giotto_polygon_overlap_message <- function() { + vmsg("[rbind giottoPolygon] Overlap information removed. + Please recalculate with calculateOverlap() if needed.") +} diff --git a/R/methods-show.R b/R/methods-show.R index 9e276c37..0c94420b 100644 --- a/R/methods-show.R +++ b/R/methods-show.R @@ -667,8 +667,23 @@ setMethod( - - +## overlapPointDT #### +setMethod("show", signature("overlapPointDT"), function(object) { + cat(sprintf("<%s>\n", class(object))) + .show_spat_and_feat(object) + .show_prov(object) + cat(sprintf("* polygons : %d\n", length(object@spat_ids))) + cat(sprintf("* features : %d\n", object@nfeats)) + cat(sprintf("* relations: %d\n", nrow(object@data))) +}) +## overlapIntensityDT #### +setMethod("show", signature("overlapIntensityDT"), function(object) { + cat(sprintf("<%s>\n", class(object))) + .show_spat_and_feat(object) + .show_prov(object) + print(object@data) + cat(color_blue("* use `[]` to drop to `data.table`")) +}) diff --git a/R/save_load.R b/R/save_load.R index 8ee35101..05e217e1 100644 --- a/R/save_load.R +++ b/R/save_load.R @@ -177,42 +177,6 @@ saveGiotto <- function( filename = filename, overwrite = TRUE ) } - - # overlap information - if (!is.null(gobject@spatial_info[[spatinfo]]@overlaps)) { - for (feature in names( - gobject@spatial_info[[spatinfo]]@overlaps - )) { - if (feature == "intensity") next - # intensities are stored as data.table - # They are already saveable with the rest of the gobject. - # Skip. - - # write names of spatvector - spatvecnames <- names( - gobject@spatial_info[[spatinfo]]@overlaps[[feature]] - ) - filename_names <- paste0( - spatinfo_dir, "/", feature, "_", - spatinfo, "_spatInfo_spatVectorOverlaps_names.txt" - ) - write.table( - x = spatvecnames, file = filename_names, - col.names = FALSE, row.names = FALSE - ) - - # write spatvector - filename <- paste0( - spatinfo_dir, "/", feature, "_", - spatinfo, - "_spatInfo_spatVectorOverlaps.shp" - ) - terra::writeVector( - gobject@spatial_info[[spatinfo]]@overlaps[[feature]], - filename = filename, overwrite = TRUE - ) - } - } } } @@ -655,6 +619,11 @@ loadGiotto <- function(path_to_folder, # load and append to gobject the polygons overlaps information .load_giotto_spatial_info_overlaps <- function(gobject, manifest, verbose = NULL) { + # objects from GiottoClass v0.5 and onwards do not need this for overlaps + if ("versions" %in% names(attributes(gobject))) { + if (.gversion(gobject) >= "0.5.0") return(gobject) + } + ## 3.3. overlaps vmsg(.v = verbose, "3.3 read Giotto spatial overlap information \n") diff --git a/R/slot_show.R b/R/slot_show.R index b104c871..d41b8abd 100644 --- a/R/slot_show.R +++ b/R/slot_show.R @@ -1104,12 +1104,10 @@ showGiottoImageNames <- function(gobject) { -#' @title Print abbreviated matrix -#' @name .abbrev_mat -#' @description print abbreviated matrix exprObj. Works for Matrix pkg denseMatrix, -#' matrix, data.frame and classes that inherit them. -#' @keywords internal -#' @returns abbreviated matrix exprObj + +# print abbreviated matrix exprObj. Works for Matrix pkg denseMatrix, +# matrix, data.frame and classes that inherit them. +# returns abbreviated matrix exprObj .abbrev_mat <- function(exprObj, nrows, ncols, print_prov = TRUE, header = TRUE ) { diff --git a/R/subset.R b/R/subset.R index 621a3439..41882b51 100644 --- a/R/subset.R +++ b/R/subset.R @@ -617,24 +617,18 @@ res_list[[spat_info]] <- spat_subset } else { # even if the spatial info is not one selected directly through - # poly_info, - # still subset subset any existing feature overlaps matching the - # feat_type - # for the feat_ids + # `poly_info`, still subset any existing feature overlaps matching + # the `feat_type` for the `feat_ids` if (!is.null(si@overlaps) && !is.null(feat_ids)) { if (isTRUE(feat_type) == ":all:") { feat_type <- names(si@overlaps) } - for (feat in names(si@overlaps)) { - if (isTRUE(feat %in% feat_type)) { - feat_id_bool <- terra::as.list( - si@overlaps[[feat]] - )$feat_ID %in% feat_ids - si@overlaps[[feat]] <- si@overlaps[[feat]][feat_id_bool] - } - } + si <- .subset_overlaps_feat(si, + i = feat_ids, + feat_type = feat_type + ) } res_list[[spat_info]] <- si @@ -1715,17 +1709,11 @@ subsetGiottoLocsSubcellular <- function( return(x) } - gpoly_overlap_names <- names(x@overlaps) - for (overlap_feat in gpoly_overlap_names) { - if (isTRUE(overlap_feat %in% names(cropped_feats))) { - feat_id_bool <- terra::as.list( - x@overlaps[[overlap_feat]] - )$feat_ID %in% - cropped_feats[[overlap_feat]] - x@overlaps[[overlap_feat]] <- x@overlaps[[ - overlap_feat - ]][feat_id_bool] - } + for (ftype in names(cropped_feats)) { + x <- .subset_overlaps_feat(x, + i = cropped_feats[[ftype]], + feat_type = ftype + ) } x }) @@ -1905,33 +1893,26 @@ subsetGiottoLocsSubcellular <- function( } } - # cell ID and feat ID subsets - if (!is.null(gpolygon@overlaps)) { - if (isTRUE(feat_type) == ":all:") feat_type <- names(gpolygon@overlaps) - - for (feat in names(gpolygon@overlaps)) { - # TODO check this for intensity image overlaps - if (!is.null(cell_ids)) { - cell_id_bool <- terra::as.list( - gpolygon@overlaps[[feat]] - )$poly_ID %in% cell_ids - gpolygon@overlaps[[feat]] <- gpolygon@overlaps[[ - feat - ]][cell_id_bool] - } + # overlap subsets + if (is.null(overlaps(gpolygon))) { + return(gpolygon) # return early if none + } - if (!is.null(feat_ids) && - isTRUE(feat %in% feat_type)) { - feat_id_bool <- terra::as.list( - gpolygon@overlaps[[feat]] - )$feat_ID %in% feat_ids - gpolygon@overlaps[[feat]] <- gpolygon@overlaps[[ - feat - ]][feat_id_bool] - } - } + if (feat_type == ":all:") { + feat_type <- names(gpolygon@overlaps) } - return(gpolygon) + + if (!is.null(cell_ids)) { + gpolygon <- .subset_overlaps_poly(gpolygon, i = cell_ids) + } + if (!is.null(feat_ids)) { + gpolygon <- .subset_overlaps_feat(gpolygon, + i = feat_ids, + feat_type = feat_type + ) + } + + gpolygon } diff --git a/R/zzz.R b/R/zzz.R index ff2c7c1a..588d8d60 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -15,4 +15,5 @@ init_option("giotto.update_param", TRUE) init_option("giotto.no_python_warn", FALSE) init_option("giotto.init_check_severity", "stop") + init_option("giotto.overlap_point_method", "vector") } diff --git a/man/aggregateFeatures.Rd b/man/aggregateFeatures.Rd new file mode 100644 index 00000000..51525bb9 --- /dev/null +++ b/man/aggregateFeatures.Rd @@ -0,0 +1,90 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/aggregate.R +\name{aggregateFeatures} +\alias{aggregateFeatures} +\title{Aggregate Spatial Features Covered by Polygon Geometries} +\usage{ +aggregateFeatures( + gobject, + spat_info = NULL, + feat_info = NULL, + image_names = NULL, + new_spat_unit = NULL, + new_feat_type = NULL, + name = "raw", + poly_subset_ids = NULL, + feat_subset_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, + fun = "sum", + return_gobject = TRUE, + verbose = TRUE, + ... +) +} +\arguments{ +\item{gobject}{\code{giotto} object containing spatial data to aggregate.} + +\item{spat_info}{character. Name of polygon information to use to aggregate +features with.} + +\item{feat_info}{character. Name of feature point detections to be +aggregated.} + +\item{image_names}{character. Name of image(s) containing intensity values +to be aggregated.} + +\item{new_spat_unit}{character (optional). Name of spatial unit to assign +the expression info to. Default is the same as \code{spat_info}.} + +\item{new_feat_type}{character (optional). Name of feature type to assign +the expression info to. Default is the same as \code{feat_info} when used. When +\code{image_names} is used instead, default is "protein".} + +\item{name}{character. (default = "raw") Name to assign the output +expresssion values information.} + +\item{poly_subset_ids}{character vector. (optional) Specific poly_IDs to use} + +\item{feat_subset_column}{character. (optional) feature info attribute to +subset feature points on when performing overlap calculation.} + +\item{feat_subset_values}{(optional) values matched against +in \code{feat_subset_column} in order to subset feature points when performing +overlap calculation.} + +\item{feat_count_column}{character. (optional) feature info column with counts +information. Useful in cases when more than one detection is reported per +point.} + +\item{fun}{character (default = "sum"). A function usable by +\code{\link[exactextractr:exact_extract]{exactextractr::exact_extract()}} to aggregate image intensity values.} + +\item{return_gobject}{logical (default = TRUE). Whether to return the +\code{giotto} object or just the aggregated expression values as \code{exprObj} class.} + +\item{verbose}{logical. Whether to be verbose.} + +\item{\dots}{Additional params to pass to the overlap calculation method. +None implemented for point overlaps. For intensity overlaps, passes to +\code{\link[exactextractr:exact_extract]{exactextractr::exact_extract()}} and additionally the function requested +with the \code{fun} param.} +} +\value{ +\code{giotto} when \code{return_gobject=TRUE}, \code{exprObj} when +\code{return_gobject=FALSE} +} +\description{ +Aggregate features (either \code{feat_info} OR \code{image_names}) with +polygons (\code{spat_info}). Under the hood, this is performed in two steps: +\enumerate{ +\item Find the overlapped features via the lower-level generic +\code{\link[=calculateOverlap]{calculateOverlap()}} +\item Summarize the overlapped features as a matrix via \code{\link[=overlapToMatrix]{overlapToMatrix()}} +} +} +\details{ +\code{feat_subset_column}, \code{feat_subset_values}, and \code{feat_count_column} +are specific to overlaps on feature points info, and should not be provided +when overlapping image data. +} diff --git a/man/as.data.table.Rd b/man/as.data.table.Rd index abeb9767..3b9e0169 100644 --- a/man/as.data.table.Rd +++ b/man/as.data.table.Rd @@ -5,6 +5,8 @@ \alias{as.data.table.SpatVector} \alias{as.data.table.giottoPolygon} \alias{as.data.table.giottoPoints} +\alias{as.data.frame.overlapPointDT} +\alias{as.data.frame.overlapIntensityDT} \title{Coerce to data.table} \usage{ \method{as.data.table}{SpatVector}( @@ -19,6 +21,10 @@ \method{as.data.table}{giottoPolygon}(x, ...) \method{as.data.table}{giottoPoints}(x, ...) + +\method{as.data.frame}{overlapPointDT}(x, ...) + +\method{as.data.frame}{overlapIntensityDT}(x, ...) } \arguments{ \item{x}{The object to coerce} diff --git a/man/as.matrix.Rd b/man/as.matrix.Rd index 96a33da3..ef1f22b6 100644 --- a/man/as.matrix.Rd +++ b/man/as.matrix.Rd @@ -3,11 +3,17 @@ \name{as.matrix} \alias{as.matrix} \alias{as.matrix,spatLocsObj-method} +\alias{as.matrix,overlapPointDT-method} +\alias{as.matrix,overlapIntensityDT-method} \alias{as.matrix,nnNetObj-method} \title{Coerce to matrix} \usage{ \S4method{as.matrix}{spatLocsObj}(x, id_rownames = TRUE, ...) +\S4method{as.matrix}{overlapPointDT}(x, feat_count_column = NULL, ...) + +\S4method{as.matrix}{overlapIntensityDT}(x, ...) + \S4method{as.matrix}{nnNetObj}(x, attr = NULL, ...) } \arguments{ @@ -17,6 +23,10 @@ \item{\dots}{additional params to pass (none implemented)} +\item{feat_count_column}{character. If provided, column in overlaps info +that contains count information to take into account when generating matrix. +This is important when point detections represent more than one count.} + \item{attr}{Either NULL or \code{character} providing the name of an edge attribute to include in the ajacency matrix. The edge attribute to use must be either \code{logical} or \code{numeric}.} diff --git a/man/calculateOverlap.Rd b/man/calculateOverlap.Rd index 1facc0c6..6753e46d 100644 --- a/man/calculateOverlap.Rd +++ b/man/calculateOverlap.Rd @@ -14,12 +14,13 @@ \S4method{calculateOverlap}{giotto,missing}( x, name_overlap = NULL, - spatial_info = NULL, + spat_info = NULL, feat_info = NULL, image_names = NULL, poly_subset_ids = NULL, return_gobject = TRUE, verbose = TRUE, + spatial_info = deprecated(), ... ) @@ -29,10 +30,12 @@ name_overlap = NULL, poly_subset_ids = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, return_gpolygon = TRUE, verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated(), ... ) @@ -66,16 +69,26 @@ ... ) -\S4method{calculateOverlap}{SpatVector,SpatRaster}(x, y, poly_subset_ids = NULL, verbose = TRUE, ...) +\S4method{calculateOverlap}{SpatVector,SpatRaster}( + x, + y, + poly_subset_ids = NULL, + verbose = TRUE, + fun = "sum", + ... +) \S4method{calculateOverlap}{SpatVector,SpatVector}( x, y, poly_subset_ids = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, - verbose = TRUE + feat_subset_values = NULL, + feat_count_column = NULL, + method = getOption("giotto.overlap_point_method", "vector"), + verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated() ) } \arguments{ @@ -85,8 +98,6 @@ polygons. Can also be a \code{giotto} object} \item{name_overlap}{name for the overlap results (default to feat_info parameter)} -\item{spatial_info}{character. Name polygon information} - \item{feat_info}{character. Name of vector feature information to overlap} \item{image_names}{character vector. Name(s) of the image feature information @@ -98,6 +109,8 @@ to overlap} \item{verbose}{be verbose} +\item{spatial_info}{character. Name polygon information} + \item{\dots}{additional params to pass to methods.} \item{y}{Object with features to overlap: \code{giottoPoints}, \code{giottoLargeImage}, @@ -106,11 +119,11 @@ to overlap} \item{feat_subset_column}{character. (optional) feature info attribute to subset feature points on when performing overlap calculation.} -\item{feat_subset_ids}{(optional) values matched against +\item{feat_subset_values}{(optional) values matched against in \code{feat_subset_column} in order to subset feature points when performing overlap calculation.} -\item{count_info_column}{character. (optional) column with count information. +\item{feat_count_column}{character. (optional) column with count information. Useful in cases when more than one detection is reported per point. If a column called "count" is present in the feature points data, it will be automatically selected.} @@ -118,6 +131,21 @@ automatically selected.} \item{return_gpolygon}{default = TRUE. Whether to return the entire giottoPolygon provided to \code{x}, but with the overlaps information appended or as a bare terra \code{SpatVector}} + +\item{feat_subset_ids}{deprecated. Use \code{feat_subset_values} instead.} + +\item{count_info_column}{deprecated. Use \code{feat_count_column} instead.} + +\item{method}{character. One of \code{"vector"} or \code{"raster"}, +(default = \code{"vector"}). Method for polygon-point feature overlap calculation. +Can also set as an option: \code{"giotto.overlap_point_method"} +\itemize{ +\item \code{"vector"} uses direct spatial extraction (more accurate to input geometry, +will double count features in overlapping polygon regions for all overlapping +polygons). +\item \code{"raster"} uses rasterization (faster, assigns each feature to only one +polygon even in overlapping regions as a byproduct of the rasterization). +}} } \value{ Usually an object of the same class as \code{x}, with the overlaps @@ -131,8 +159,8 @@ polygon annotations. This provides a summary of the spatial data overlapped by the polygon which can be further processed to become an expression matrix. } \details{ -\code{feat_subset_column}, \code{feat_subset_ids}, and \code{count_info_column} are -specific to overlaps on feature points info, and should not be provided +\code{feat_subset_column}, \code{feat_subset_values}, and \code{feat_count_column} +are specific to overlaps on feature points info, and should not be provided when overlapping image data. These three params can also be passed to the \code{giotto} method through the \code{...} param when working with overlaps on feature points info. @@ -160,7 +188,7 @@ overlaps_all$rna # calculate z1 only out_z1 <- calculateOverlap(gpoly, gpoints, feat_subset_column = "global_z", - feat_subset_ids = c(1) + feat_subset_values = c(1) ) overlaps_z1 <- overlaps(out_z1) overlaps_z1$rna @@ -174,7 +202,7 @@ overlaps_img$intensity # calculate z0 overlaps and return as gobject out_g <- calculateOverlap(g, feat_subset_column = "global_z", - feat_subset_ids = 0 + feat_subset_values = 0 ) overlaps(getPolygonInfo(out_g, return_giottoPolygon = TRUE)) diff --git a/man/calculateOverlapRaster.Rd b/man/calculateOverlapRaster.Rd index 7db26aed..476f27ea 100644 --- a/man/calculateOverlapRaster.Rd +++ b/man/calculateOverlapRaster.Rd @@ -11,10 +11,12 @@ calculateOverlapRaster( poly_ID_names = NULL, feat_info = NULL, feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, + feat_subset_values = NULL, + feat_count_column = NULL, return_gobject = TRUE, - verbose = TRUE + verbose = TRUE, + feat_subset_ids = deprecated(), + count_info_column = deprecated() ) } \arguments{ @@ -31,13 +33,18 @@ results (default to feat_info parameter)} \item{feat_subset_column}{feature info column to subset features with} -\item{feat_subset_ids}{ids within feature info column to use for subsetting} +\item{feat_subset_values}{value(s) within feature info \code{feat_subset_column} +to use for subsetting} -\item{count_info_column}{column with count information (optional)} +\item{feat_count_column}{column with count information (optional)} \item{return_gobject}{return giotto object (default: TRUE)} \item{verbose}{be verbose} + +\item{feat_subset_ids}{deprecated. Use \code{feat_subset_values} instead.} + +\item{count_info_column}{deprecated. Use \code{feat_count_column} instead.} } \value{ giotto object or spatVector with overlapping information diff --git a/man/dims-generic.Rd b/man/dims-generic.Rd index ae95bd5e..9d4817a2 100644 --- a/man/dims-generic.Rd +++ b/man/dims-generic.Rd @@ -11,11 +11,15 @@ \alias{nrow,spatialNetworkObj-method} \alias{nrow,enrData-method} \alias{nrow,dimObj-method} +\alias{nrow,overlapPointDT-method} +\alias{nrow,overlapIntensityDT-method} \alias{ncol,giotto-method} \alias{ncol,exprData-method} \alias{ncol,metaData-method} \alias{ncol,enrData-method} \alias{ncol,dimObj-method} +\alias{ncol,overlapPointDT-method} +\alias{ncol,overlapIntensityDT-method} \alias{dim,giotto-method} \alias{dim,spatLocsObj-method} \alias{dim,exprData-method} @@ -24,6 +28,8 @@ \alias{dim,giottoLargeImage-method} \alias{dim,giottoPolygon-method} \alias{dim,giottoPoints-method} +\alias{dim,overlapPointDT-method} +\alias{dim,overlapIntensityDT-method} \title{Dimensions of giotto objects} \usage{ \S4method{nrow}{giotto}(x) @@ -44,6 +50,10 @@ \S4method{nrow}{dimObj}(x) +\S4method{nrow}{overlapPointDT}(x) + +\S4method{nrow}{overlapIntensityDT}(x) + \S4method{ncol}{giotto}(x) \S4method{ncol}{exprData}(x) @@ -54,6 +64,10 @@ \S4method{ncol}{dimObj}(x) +\S4method{ncol}{overlapPointDT}(x) + +\S4method{ncol}{overlapIntensityDT}(x) + \S4method{dim}{giotto}(x) \S4method{dim}{spatLocsObj}(x) @@ -69,6 +83,10 @@ \S4method{dim}{giottoPolygon}(x) \S4method{dim}{giottoPoints}(x) + +\S4method{dim}{overlapPointDT}(x) + +\S4method{dim}{overlapIntensityDT}(x) } \arguments{ \item{x}{object to check dimensions of} diff --git a/man/dot-abbrev_mat.Rd b/man/dot-abbrev_mat.Rd deleted file mode 100644 index 57f23510..00000000 --- a/man/dot-abbrev_mat.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/slot_show.R -\name{.abbrev_mat} -\alias{.abbrev_mat} -\title{Print abbreviated matrix} -\usage{ -.abbrev_mat(exprObj, nrows, ncols, print_prov = TRUE, header = TRUE) -} -\value{ -abbreviated matrix exprObj -} -\description{ -print abbreviated matrix exprObj. Works for Matrix pkg denseMatrix, -matrix, data.frame and classes that inherit them. -} -\keyword{internal} diff --git a/man/dot-calculate_overlap_raster.Rd b/man/dot-calculate_overlap_raster.Rd index d4dc38b0..3caec6fb 100644 --- a/man/dot-calculate_overlap_raster.Rd +++ b/man/dot-calculate_overlap_raster.Rd @@ -4,19 +4,14 @@ \alias{.calculate_overlap_raster} \title{Find feature points overlapped by rasterized polygon.} \usage{ -.calculate_overlap_raster( - spatvec, - pointvec, - count_info_column = NULL, - verbose = TRUE -) +.calculate_overlap_raster(spatvec, pointvec, keep = NULL, verbose = TRUE) } \arguments{ \item{spatvec}{\code{SpatVector} polygon from a \code{giottoPolygon} object} \item{pointvec}{\code{SpatVector} points from a \code{giottoPoints} object} -\item{count_info_column}{column with count information (optional)} +\item{keep}{column(s) to keep} \item{verbose}{be verbose} } diff --git a/man/overlapPointDT-class.Rd b/man/overlapPointDT-class.Rd new file mode 100644 index 00000000..57e8ceee --- /dev/null +++ b/man/overlapPointDT-class.Rd @@ -0,0 +1,100 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/classes.R, R/methods-extract.R +\docType{class} +\name{overlapPointDT-class} +\alias{overlapPointDT-class} +\alias{[,overlapPointDT,gIndex,missing,missing-method} +\alias{[,overlapPointDT,missing,gIndex,missing-method} +\alias{[,overlapPointDT,gIndex,gIndex,missing-method} +\title{Polygon and Point Relationships} +\usage{ +\S4method{[}{overlapPointDT,gIndex,missing,missing}(x, i, j, ..., use_names = FALSE, ids = TRUE, drop) + +\S4method{[}{overlapPointDT,missing,gIndex,missing}(x, i, j, ..., use_names = FALSE, ids = TRUE, drop) + +\S4method{[}{overlapPointDT,gIndex,gIndex,missing}(x, i, j, ..., use_names = FALSE, drop) +} +\arguments{ +\item{x}{object} + +\item{i}{numeric, character, logical. Index of or name of poly in overlapping +polygons} + +\item{j}{numeric, character, logical. Index of or name of feature being +overlapped.} + +\item{\dots}{additional params to pass (none implemented)} + +\item{use_names}{logical (default = \code{FALSE}). Whether to return as integer +indices or with character ids.} + +\item{ids}{logical (default = \code{TRUE}). Whether to return the requested +integer indices (\code{TRUE}) or the subset overlap object (\code{FALSE}).} + +\item{drop}{not used.} +} +\value{ +integer or character if only \code{i} or \code{j} provided, depending on +\code{use_names}. A subset \code{overlapPointDT} if both \code{i} and \code{j} are used. +} +\description{ +Utility class for storing overlaps relationships between polygons and points +in a sparse \code{data.table} format. Retrieve the unique ID index of overlapped +points \verb{[i, ]}. Get indices of which polys are overlapping specific feature +species using \verb{[, j]}. + +Subsetting with \code{ids = FALSE} and \verb{[i, j]} indexing is also supported. + +Supports \code{as.matrix} for conversion to \code{dgCMatrix}. Contained poly and +feature names simplify rownames/colnames and empty row/col creation. +} +\section{Slots}{ + +\describe{ +\item{\code{data}}{data.table. Table containing 3 integer cols: +\itemize{ +\item \code{poly} - polygon index. Maps to \code{spat_ids} slot. +\item \code{feat} - feat_ID_uniq (unique integer identifier) of a point detection +\item \code{feat_id_index} - index of feature name mapping in \verb{@feat_ids} slot. +}} + +\item{\code{spat_unit}}{character. Spatial unit (usually name of polygons information)} + +\item{\code{feat_type}}{character. Feature type (usually name of points information)} + +\item{\code{provenance}}{character. provenance information} + +\item{\code{spat_ids}}{character. Polygon names} + +\item{\code{feat_ids}}{character. Feature names} + +\item{\code{nfeats}}{integer (optional metadata). How many feature points were +used in overlap operation. Gives an idea of sparsity, but has no effect on +processing.} +}} + +\examples{ +g <- GiottoData::loadGiottoMini("vizgen") +poly <- g[["spatial_info", "z0"]][[1]] +ovlp <- overlaps(poly, "rna") +ovlp + +as.matrix(ovlp) + +dim(ovlp) +nrow(ovlp) # number of relationships + +# get feature unique IDs overlapped by nth poly +ovlp[1] # check one (no overlaps returns integer(0)) +ovlp[1:5] # check multiple +ovlp[1:5, use_names = TRUE] # returns feature names, but no longer unique + +# get integer index of poly(s) overlapping particular feature species +ovlp[, 1] +ovlp[, "Mlc1"] # this is the same + +# get a subset of overlap object +ovlp[1:10, ids = FALSE] # subset to first 10 polys +ovlp[, 1:10, ids = FALSE] # subset to first 10 feature species +ovlp[1:10, 1:10] # subset to first 10 polys and first 10 features species +} diff --git a/man/overlapToMatrix.Rd b/man/overlapToMatrix.Rd index c5a31209..5d6af34f 100644 --- a/man/overlapToMatrix.Rd +++ b/man/overlapToMatrix.Rd @@ -6,18 +6,23 @@ \alias{overlapToMatrix,giottoPolygon-method} \alias{overlapToMatrix,SpatVector-method} \alias{overlapToMatrix,data.table-method} +\alias{overlapToMatrix,overlapPointDT-method} +\alias{overlapToMatrix,overlapIntensityDT-method} \title{overlapToMatrix} \usage{ \S4method{overlapToMatrix}{giotto}( x, name = "raw", - poly_info = NULL, + spat_info = NULL, feat_info = NULL, type = c("point", "intensity"), - count_info_column = NULL, - aggr_function = "sum", + feat_count_column = NULL, + fun = "sum", return_gobject = TRUE, verbose = TRUE, + aggr_function = deprecated(), + poly_info = deprecated(), + count_info_column = deprecated(), ... ) @@ -25,8 +30,8 @@ x, feat_info = "rna", type = c("point", "intensity"), - count_info_column = NULL, - output = c("Matrix", "data.table"), + feat_count_column = NULL, + count_info_column = deprecated(), ... ) @@ -34,13 +39,36 @@ x, col_names = NULL, row_names = NULL, - count_info_column = NULL, + feat_count_column = NULL, output = c("Matrix", "data.table"), verbose = TRUE, + count_info_column = deprecated(), + ... +) + +\S4method{overlapToMatrix}{data.table}( + x, + fun = "sum", + output = c("Matrix", "data.table"), + aggr_function = deprecated() +) + +\S4method{overlapToMatrix}{overlapPointDT}( + x, + name = "raw", + sort = TRUE, + feat_count_column = NULL, + output = c("Matrix", "exprObj"), ... ) -\S4method{overlapToMatrix}{data.table}(x, aggr_function = "sum", output = c("Matrix", "data.table")) +\S4method{overlapToMatrix}{overlapIntensityDT}( + x, + name = "raw", + sort = TRUE, + output = c("Matrix", "exprObj"), + ... +) } \arguments{ \item{x}{object containing overlaps info. Can be giotto object or SpatVector @@ -48,30 +76,40 @@ points or data.table of overlaps generated from \code{calculateOverlap}} \item{name}{name for the overlap count matrix} -\item{poly_info}{character. Polygon information to use} +\item{spat_info}{character. Polygon information to use} \item{feat_info}{character. Feature information to use} \item{type}{character. Type of overlap data (either 'point' or 'intensity')} -\item{count_info_column}{column with count information. If a +\item{feat_count_column}{column with count information. If a column called "count" is present in the feature points data, it will be automatically selected.} -\item{aggr_function}{function to aggregate image information (default = sum)} +\item{fun}{character. Function to aggregate image information +(default = "sum")} \item{return_gobject}{return giotto object (default: TRUE)} \item{verbose}{be verbose} -\item{\dots}{additional params to pass to methods} +\item{aggr_function}{deprecated. Use \code{fun} instead.} -\item{output}{data format/class to return the results as} +\item{poly_info}{deprecated. Please use spat_info.} + +\item{count_info_column}{deprecated. Use \code{feat_count_column} instead.} + +\item{\dots}{additional params to pass to methods} \item{col_names, row_names}{character vector. (optional) Set of row and col names that are expected to exist. This fixes the dimensions of the matrix since the overlaps information does not directly report rows and cols where no values were detected.} + +\item{output}{data format/class to return the results as. Default is "Matrix"} + +\item{sort}{logical (default = TRUE). Whether to perform a mixed sort on +output matrix row and col names.} } \value{ giotto object or count matrix diff --git a/man/rbind-generic.Rd b/man/rbind-generic.Rd index 8b89e3bf..aa8cf2cc 100644 --- a/man/rbind-generic.Rd +++ b/man/rbind-generic.Rd @@ -6,6 +6,7 @@ \alias{rbind2,featMetaObj,featMetaObj-method} \alias{rbind2,spatLocsObj,spatLocsObj-method} \alias{rbind2,giottoPolygon,giottoPolygon-method} +\alias{rbind2,overlapPointDT,overlapPointDT-method} \title{Combine objects by rows (Giotto-related)} \usage{ \S4method{rbind2}{cellMetaObj,cellMetaObj}(x, y, ...) @@ -15,6 +16,8 @@ \S4method{rbind2}{spatLocsObj,spatLocsObj}(x, y, ...) \S4method{rbind2}{giottoPolygon,giottoPolygon}(x, y, add_list_ID = TRUE, ...) + +\S4method{rbind2}{overlapPointDT,overlapPointDT}(x, y, ...) } \arguments{ \item{x}{item 1 to rbind} diff --git a/tests/testthat/test-aggregate.R b/tests/testthat/test-aggregate.R new file mode 100644 index 00000000..27de3c9d --- /dev/null +++ b/tests/testthat/test-aggregate.R @@ -0,0 +1,245 @@ +random_pts_names <- function(n, species = 20) { + GiottoUtils::local_seed(1234) + sampleset <- c(letters, LETTERS, seq(from = 0, to = 9)) + nameset <- vapply(seq_len(species), FUN = function(i) { + paste(sample(sampleset, size = 8, replace = TRUE), collapse = "") + + }, FUN.VALUE = character(1L)) + sample(nameset, replace = TRUE, size = n, prob = runif(species)) +} + +random_points_gen <- function(n = 500, extent = ext(gpoly), count = TRUE) { + GiottoUtils::local_seed(1234) + evect <- as.numeric(ext(extent)[]) + cts <- abs(round(rnorm(n, 0, sd = 0.8))) + 1 + d <- data.table::data.table( + id = random_pts_names(n), + x = runif(n, min = evect[[1]], max = evect[[2]]), + y = runif(n, min = evect[[3]], max = evect[[4]]) + ) + if (count) { + d[, count := cts] + } + d +} + +g <- test_data$viz +gpoly <- g[["spatial_info", "aggregate"]][[1]] +gpoly@overlaps = NULL +gpoly@spatVectorCentroids <- NULL +gpts <- createGiottoPoints(random_points_gen(80000, count = FALSE), + verbose = FALSE) +gpts_cts <- createGiottoPoints(random_points_gen(80000, count = TRUE), + verbose = FALSE) +imglist <- g[["images",]] +img <- imglist[[1]] + +# calculateOverlap #### +## --- poly vs pts level tests + +# these tests can change if the source test dataset changes + +test_that("calculateOverlap works for points", { + # [NO COUNTS] + # 1. test raster method + res_rast <- calculateOverlap(gpoly, gpts, verbose = FALSE, method = "raster") + expect_identical(names(res_rast@overlaps), "rna") + ovlp_rast <- overlaps(res_rast, "rna") + checkmate::expect_class(ovlp_rast, "overlapInfo") + expect_equal(nrow(ovlp_rast@data), 12383) + expect_identical(as.numeric(ovlp_rast@data[100,]), c(385, 685, 12)) + + # 2. test vector method (default) + res_vect <- calculateOverlap(gpoly, gpts, verbose = FALSE, method = "vector") + + ovlp_vect <- overlaps(res_vect, "rna") + expect_equal(nrow(ovlp_vect@data), 12311) + expect_identical(as.numeric(ovlp_vect@data[100,]), c(12, 671, 3)) + + # 3. check expected col names in output + expect_identical( + names(ovlp_vect@data), + c("poly", "feat", "feat_id_index") + ) + + # [WITH COUNTS] + # counts info now automatically included if a `count` column is available + # in points input + + # 1. test raster method + res_rast <- calculateOverlap(gpoly, gpts_cts, verbose = FALSE, method = "raster") + expect_identical(names(res_rast@overlaps), "rna") + ovlp_rast <- overlaps(res_rast, "rna") + checkmate::expect_class(ovlp_rast, "overlapInfo") + expect_equal(nrow(ovlp_rast@data), 12383) + expect_identical(as.numeric(ovlp_rast@data[100,]), c(385, 685, 12, 2)) + + # 2. test vector method (default) + res_vect <- calculateOverlap(gpoly, gpts_cts, verbose = FALSE, method = "vector") + + ovlp_vect <- overlaps(res_vect, "rna") + expect_equal(nrow(ovlp_vect@data), 12311) + expect_identical(as.numeric(ovlp_vect@data[100,]), c(12, 671, 3, 2)) + + # 3. check expected col names in output + expect_identical( + names(ovlp_vect@data), + c("poly", "feat", "feat_id_index", "count") + ) +}) + +test_that("calculateOverlap works for basic images", { + res <- calculateOverlap(gpoly, img, name = "image", verbose = FALSE, ) + o <- overlaps(res, "image") + vals <- c(514207.812, 576698.000, 593017.562) + expect_equal(head(o@data[[2]], 3), vals) +}) + +test_that("calculateOverlap works for affine images", { + res_base <- calculateOverlap(gpoly, img) + aimg <- spin(imglist[[1]], 45, x0 = 0, y0 = 0) + spin_gpoly <- spin(gpoly, 45, x0 = 0, y0 = 0) + res_aff <- calculateOverlap(spin_gpoly, aimg, + name = "affine_image", verbose = FALSE + ) + ovlp_base <- overlaps(res_base, "dapi_z0") + ovlp_aff <- overlaps(res_aff, "affine_image") + # original vs affine calculates same value + expect_equal(ovlp_base@data$mini_dataset_dapi_z0, + ovlp_aff@data$mini_dataset_dapi_z0) +}) + + +# overlapToMatrix #### +## --- poly level tests + +# test data + + + +test_that("overlapToMatrix works for point overlaps", { + res_vect <- calculateOverlap(gpoly, gpts, + verbose = FALSE, + method = "vector" + ) + ovlp_vect <- overlaps(res_vect, "rna") + expect_identical(names(ovlp_vect@data), + c("poly", "feat", "feat_id_index")) + + # with a counts column summation + res_vect_cts <- calculateOverlap(gpoly, gpts_cts, + feat_count_column = "count", + verbose = FALSE, + method = "vector" + # count col should be autodetected + ) + ovlp_vect_cts <- overlaps(res_vect_cts, "rna") + expect_identical(names(ovlp_vect_cts@data), + c("poly", "feat", "feat_id_index", "count")) + + # matrix without counts + mat <- as.matrix(ovlp_vect) + mat_val <- mat["Bp5vKRUi", "104929374280991324709110264935409912418"] + expect_equal(mat_val, 4) + # matrix with counts + mat <- as.matrix(ovlp_vect_cts, feat_count_column = "count") + mat_val <- mat["Bp5vKRUi", "104929374280991324709110264935409912418"] + expect_equal(mat_val, 6) + # matrix with counts, calculated differently + ovlp_val_table <- ovlp_vect_cts@data[ + poly == match("104929374280991324709110264935409912418", + ovlp_vect_cts@spat_ids) & + feat_id_index == match("Bp5vKRUi", ovlp_vect_cts@feat_ids) + ] + ovlp_val <- sum(ovlp_val_table$count) + expect_equal(ovlp_val, 6) +}) + +test_that("overlapToMatrix works for intensity overlaps", { + res_rast <- calculateOverlap(gpoly, img, verbose = FALSE) + ovlp_rast <- overlaps(res_rast, "dapi_z0") + m <- as.matrix(ovlp_rast) + expect_equal(m[1,10], 536529.94) +}) + + +# overlap point object tests #### + +res_vect <- calculateOverlap(gpoly, gpts, + verbose = FALSE, + method = "vector" +) +ovlp <- overlaps(res_vect, "rna") + +test_that("overlap `[]` subset works", { + # expect the feat_ID_uniq overlapped by poly 9 + fuid9 <- ovlp[9] + expected_fuid9 <- c( + 1949, 4234, 6934, 9360, 11891, 12037, 13671, 13766, 14633, 14732, + 22670, 25184, 27910, 30742, 37026, 44881, 45922, 50053, 53519, 67618, + 68997, 69003, 69911, 70336, 72325, 74205, 76817 + ) + expect_equal(fuid9, expected_fuid9) + # expect the feat_ID_uniq overlapped by poly 40 + fuid40 <- ovlp[40] + expected_fuid40 <- c( + 2184, 3710, 3870, 5862, 12664, 19589, 24959, 25368, 26082, 29026, + 36516, 45774, 49727, 50094, 50514, 53664, 53834, 54290, 54365, 54534, + 55520, 59178, 66336, 66812, 70350, 73361, 75951 + ) + expect_equal(fuid40, expected_fuid40) + + # expect poly overlapping specific features + feat_idx <- match("Bp5vKRUi", ovlp@feat_ids) + poly_idx_by_fidx <- head(ovlp[, feat_idx]) + poly_idx_by_fname <- head(ovlp[, "Bp5vKRUi"]) + expect_identical(poly_idx_by_fidx, poly_idx_by_fname) + expected_idx <- c(459, 215, 11, 30, 294, 301) + expect_equal(poly_idx_by_fidx, expected_idx) +}) + +test_that("overlap `as.data.frame` works", { + res <- as.data.frame(ovlp) + expect_equal(ncol(res), 3) + expect_equal(nrow(res), 12311) + checkmate::expect_character(res$poly_ID) + checkmate::expect_character(res$feat_ID) + checkmate::expect_integer(res$feat_ID_uniq) +}) + +# aggregateFeatures #### +## --- gobject-level checks + +test_that("aggregateFeatures works and generates exprObj", { + g <- aggregateFeatures(g, + spat_info = "z0", + feat_info = "rna", + new_spat_unit = "test_spat", + new_feat_type = "test_feat", + name = "test_mat", + verbose = FALSE, + return_gobject = TRUE + ) + e <- g@expression$test_spat$test_feat$test_mat + expect_s4_class(e, "exprObj") + m <- e[] + expect_s4_class(m, "dgCMatrix") + expect_identical(head(colnames(m)), c( + "1132915601442719251817312578799507532", + "1900302715660571356090444774019116326", + "2467888014556719520437642348850497467", + "2582675971475731682260721390861103474", + "3151102251621248215463444891373423271", + "3188909098022617910378369321966273882" + )) + expect_identical(head(rownames(m)), c( + "Abcc9", "Ackr1", "Ackr3", "Adcyap1r1", "Adgra1", "Adgra2" + )) + expect_identical(objName(e), "test_mat") + expect_identical(spatUnit(e), "test_spat") + expect_identical(featType(e), "test_feat") + expect_equal(dim(e), c(337, 498)) +}) + + + diff --git a/tests/testthat/test-create_mini_vizgen.R b/tests/testthat/test-create_mini_vizgen.R index da2de8dc..1a9e7f26 100644 --- a/tests/testthat/test-create_mini_vizgen.R +++ b/tests/testthat/test-create_mini_vizgen.R @@ -141,7 +141,7 @@ test_that("giottoLargeImages are created", { vizsubc <- addGiottoImage( gobject = vizsubc, - largeImages = imagelist + images = imagelist ) test_that("images were added", { @@ -159,20 +159,20 @@ test_that("images were added", { # we can set a global option or specify this for each command # options('giotto.spat_unit' = 'z1') # now you don't need to think about setting spat_unit each time -vizsubc <- calculateOverlapRaster( +vizsubc <- calculateOverlap( vizsubc, - spatial_info = "z0", + spat_info = "z0", feat_info = "rna", feat_subset_column = "global_z", - feat_subset_ids = 0 + feat_subset_values = 0 ) -vizsubc <- calculateOverlapRaster( +vizsubc <- calculateOverlap( vizsubc, - spatial_info = "z1", + spat_info = "z1", feat_info = "rna", feat_subset_column = "global_z", - feat_subset_ids = 1 + feat_subset_values = 1 ) @@ -194,22 +194,20 @@ nfeats <- length(feats) test_that("overlaps are calculated", { - expect_class(overlaps(z0_gpoly)$rna, "SpatVector") - expect_class(overlaps(z0_gpoly)$rna, "SpatVector") - expect_identical(terra::geomtype(overlaps(z0_gpoly)$rna), "points") - expect_identical(terra::geomtype(overlaps(z0_gpoly)$rna), "points") + expect_class(overlaps(z0_gpoly)$rna, "overlapPointDT") + expect_class(overlaps(z0_gpoly)$rna, "overlapPointDT") }) vizsubc <- overlapToMatrix( vizsubc, - poly_info = "z0", + spat_info = "z0", feat_info = "rna", name = "raw" ) vizsubc <- overlapToMatrix( vizsubc, - poly_info = "z1", + spat_info = "z1", feat_info = "rna", name = "raw" )