From 9df43381656649533744e18f5efaeff696c53d4f Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 6 Feb 2024 01:46:26 -0500 Subject: [PATCH 01/24] new: `createNetwork()` generic - make generic function `createNetwork()` - new method for all matrices to create sNN and kNN networks. Flexibly allows inclusion of weight, distance cols and name attributes. Allows return as data.table or igraph --- R/NN_network.R | 220 ++++++++++++++++++++++++++++++++++++----- R/classes.R | 7 ++ R/generics.R | 2 + R/package_imports.R | 1 + R/spatial_structures.R | 6 +- 5 files changed, 207 insertions(+), 29 deletions(-) diff --git a/R/NN_network.R b/R/NN_network.R index fe63521e..45fa1e14 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -1,3 +1,170 @@ + + + + + +setMethod( + "createNetwork", + signature("allMatrix"), + function( + x, + type = c("sNN", "kNN"), + k = 30, + minimum_shared = 5, + top_shared = 3, + node_ids = NULL, + include_weight = TRUE, + include_distance = TRUE, + as.igraph = TRUE, + ... + ) + { + + # check params + k <- as.integer(k) + minimum_shared <- as.integer(minimum_shared) + top_shared <- as.integer(top_shared) + type <- match.arg(type, choices = c("sNN", "kNN")) + + if (k >= nrow(x)) { + k <- (nrow(x) - 1L) + vmsg(.v = verbose, "k is higher than total number of cells. + Adjusted to (total number of cells - 1)") + } + + # get common params + alist <- list( + x = x, + k = k, + node_ids = node_ids, + include_weight = include_weight, + include_distance = include_distance, + ... + ) + + # generate network data.table + network_dt <- switch( + type, + "kNN" = do.call(.net_dt_knn, args = alist), + "sNN" = do.call( + .net_dt_snn, + args = c( + alist, list(top_shared = top_shared, minimum_shared = minimum_shared) + ) + ) + ) + + + # return early if igraph not required + if (!as.igraph) return(network_dt) + + + ## convert to igraph object + + # cols to include in igraph object + igraph_cols <- if (is.null(node_ids)) { + c("from", "to") + } else { + c("from_cell_ID", "to_cell_ID") + } + + all_index <- network_dt[, unique(unlist(.SD)), .SDcols = igraph_cols] + if (include_weight) igraph_cols <- c(igraph_cols, "weight") + if (include_distance) igraph_cols <- c(igraph_cols, "distance") + + if (type == "sNN") igraph_cols <- c(igraph_cols, "shared", "rank") + + out <- igraph::graph_from_data_frame( + network_dt[, igraph_cols, with = FALSE], + directed = TRUE, + vertices = all_index + ) + + return(out) + } +) + + +# x input is a matrix +.net_dt_knn <- function( + x, k, node_ids, include_weight = TRUE, include_distance = TRUE, ... +) { + # NSE vars + from <- to <- NULL + + nn_network <- dbscan::kNN(x = x, k = k, sort = TRUE, ...) + + nn_network_dt <- data.table::data.table( + from = rep(1:nrow(nn_network$id), k), + to = as.vector(nn_network$id) + ) + + if (!is.null(node_ids)) { + # if node_ids are provided, include as extra cols + names(node_ids) <- seq_along(node_ids) + nn_network_dt[, "from_cell_ID" := node_ids[from]] + nn_network_dt[, "to_cell_ID" := node_ids[to]] + } + + # optional info + if (include_weight) { + nn_network_dt[, "weight" := 1 / (1 + as.vector(nn_network$dist))] + } + if (include_distance) { + nn_network_dt[, "distance" := as.vector(nn_network$dist)] + } + + return(nn_network_dt) +} + +# x input is a matrix +.net_dt_snn <- function( + x, k, node_ids, include_weight = TRUE, include_distance = TRUE, + top_shared = 3, minimum_shared = 5, ... +) { + # NSE vars + from <- to <- shared <- NULL + + nn_network <- dbscan::kNN(x = x, k = k, sort = TRUE, ...) + snn_network <- dbscan::sNN(x = nn_network, k = k, kt = NULL, ...) + + snn_network_dt <- data.table::data.table( + from = rep(1:nrow(snn_network$id), k), + to = as.vector(snn_network$id), + shared = as.vector(snn_network$shared) + ) + snn_network_dt <- snn_network_dt[stats::complete.cases(snn_network_dt)] + + if (!is.null(node_ids)) { + # if node_ids are provided, include as extra cols + names(node_ids) <- seq_along(node_ids) + snn_network_dt[, "from_cell_ID" := node_ids[from]] + snn_network_dt[, "to_cell_ID" := node_ids[to]] + } + + # optional info + if (include_weight) { + snn_network_dt[, "weight" := 1 / (1 + as.vector(snn_network$dist))] + } + if (include_distance) { + snn_network_dt[, "distance" := as.vector(snn_network$dist)] + } + + + # rank snn + data.table::setorder(snn_network_dt, from, -shared) + snn_network_dt[, rank := 1:.N, by = from] + + # filter snn + snn_network_dt <- snn_network_dt[rank <= top_shared | shared >= minimum_shared] + + return(snn_network_dt) +} + + + + + #' @title createNearestNetwork #' @name createNearestNetwork #' @description create a nearest neighbour (NN) network @@ -49,28 +216,25 @@ #' } #' #' @export -createNearestNetwork <- function(gobject, - spat_unit = NULL, - feat_type = NULL, - type = c("sNN", "kNN"), - dim_reduction_to_use = "pca", - dim_reduction_name = NULL, - dimensions_to_use = 1:10, - feats_to_use = NULL, - genes_to_use = NULL, - expression_values = c("normalized", "scaled", "custom"), - name = NULL, - return_gobject = TRUE, - k = 30, - minimum_shared = 5, - top_shared = 3, - verbose = T, - ...) { - ## deprecated arguments - if (!is.null(genes_to_use)) { - feats_to_use <- genes_to_use - warning("genes_to_use is deprecated, use feats_to_use in the future \n") - } +createNearestNetwork <- function( + gobject, + spat_unit = NULL, + feat_type = NULL, + type = c("sNN", "kNN"), + dim_reduction_to_use = "pca", + dim_reduction_name = NULL, + dimensions_to_use = 1:10, + feats_to_use = NULL, + expression_values = c("normalized", "scaled", "custom"), + name = NULL, + return_gobject = TRUE, + k = 30, + minimum_shared = 5, + top_shared = 3, + verbose = TRUE, + ... +) +{ # Set feat_type and spat_unit spat_unit <- set_default_spat_unit( @@ -105,7 +269,11 @@ createNearestNetwork <- function(gobject, ) if (!dim_reduction_name %in% dim_red_names) { - stop("\n dimension reduction: ", dim_reduction_to_use, " or dimension reduction name: ", dim_reduction_name, " is not available \n") + stop(sprintf( + "\n dimension reduction: %s or dimension reduction name: %s is not available \n", + dim_reduction_to_use, + dim_reduction_name + )) } # check = gobject@dimension_reduction[['cells']][[spat_unit]][[dim_reduction_to_use]][[dim_reduction_name]] @@ -160,7 +328,8 @@ createNearestNetwork <- function(gobject, ## run nearest-neighbour algorithm ## if (k >= nrow(matrix_to_use)) { k <- (nrow(matrix_to_use) - 1) - if (verbose == TRUE) cat("\n k is higher than total number of cells, adjusted to (total number of cells - 1) \n") + vmsg(.v = verbose, "k is higher than total number of cells. + Adjusted to (total number of cells - 1)") } nn_network <- dbscan::kNN(x = matrix_to_use, k = k, sort = TRUE, ...) @@ -206,6 +375,7 @@ createNearestNetwork <- function(gobject, if (type == "kNN") { nn_network_igraph <- igraph::graph_from_data_frame(nn_network_dt[, .(from_cell_ID, to_cell_ID, weight, distance)], directed = TRUE, vertices = all_index) } else if (type == "sNN") { + # TODO never returned? missing_indices <- all_index[!all_index %in% unique(snn_network_dt$from)] nn_network_igraph <- igraph::graph_from_data_frame(snn_network_dt[, .(from_cell_ID, to_cell_ID, weight, distance, shared, rank)], directed = TRUE, vertices = all_index) } @@ -242,8 +412,6 @@ createNearestNetwork <- function(gobject, nn_network = nnObj ) - # gobject@nn_network[[spat_unit]][[type]][[name]][['igraph']] = nn_network_igraph - ## update parameters used ## gobject <- update_giotto_params(gobject, description = "_nn_network") diff --git a/R/classes.R b/R/classes.R index afb04f55..a4a1775d 100644 --- a/R/classes.R +++ b/R/classes.R @@ -30,6 +30,13 @@ setClassUnion("nullOrDatatable", c("NULL", "data.table")) #' @noRd setClassUnion("gIndex", c("numeric", "logical", "character")) +#' @title Matrix classes +#' @description combined class definition for matrix representations for use in +#' dispatching. +#' @keywords internal +#' @noRd +setClassUnion("allMatrix", c("matrix", "Matrix", "DelayedArray")) + diff --git a/R/generics.R b/R/generics.R index 8c807466..0364b79d 100644 --- a/R/generics.R +++ b/R/generics.R @@ -52,6 +52,8 @@ setGeneric("copy", function(x) standardGeneric("copy"), useAsDefault = data.tabl setGeneric("calculateOverlap", function(x, y, ...) standardGeneric("calculateOverlap")) setGeneric("overlapToMatrix", function(x, ...) standardGeneric("overlapToMatrix")) +# networks +setGeneric("createNetwork", function(x, ...) standardGeneric("createNetwork")) # Giotto subnesting #### # All methods and documentations found in methods-nesting.R diff --git a/R/package_imports.R b/R/package_imports.R index 99653da7..8320719a 100644 --- a/R/package_imports.R +++ b/R/package_imports.R @@ -27,6 +27,7 @@ #' @importMethodsFrom terra nrow ncol #' @importClassesFrom terra SpatExtent #' @importClassesFrom terra SpatVector +#' @importClassesFrom DelayedArray DelayedArray #' @import GiottoUtils #' @import data.table #' @import utils diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 0f6c1f4a..2aedda96 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -290,8 +290,9 @@ spat_net_to_igraph <- function(spatialNetworkObj, attr = NULL) { ## create delaunay network delaunay_triangle <- geometry::delaunayn( - p = spatial_locations[, c(sdimx, sdimy), with = F], - options = options, ... + p = spatial_locations[, c(sdimx, sdimy), with = FALSE], + options = options, + ... ) ## save delaunay network object @@ -384,7 +385,6 @@ spat_net_to_igraph <- function(spatialNetworkObj, attr = NULL) { delaunay_edges_dedup2 <- igraph::get.data.frame(igraph_obj2) delaunay_edges_dedup <- data.table::as.data.table(delaunay_edges_dedup2) ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### - ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### xbegin_name <- paste0(sdimx, "_begin") ybegin_name <- paste0(sdimy, "_begin") From 3016f20354c4e188b35c1d793eadc52b0de596cc Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Sun, 11 Feb 2024 17:55:19 -0500 Subject: [PATCH 02/24] feat: WIP createNetwork() hub function Still needs further work aligning all the params and making obvious things to edit in order to get similar behavior to default Giotto operations --- DESCRIPTION | 2 +- NAMESPACE | 3 + R/NN_network.R | 342 ++++++++++++++++++++++++++++++------ R/package_imports.R | 1 + R/spatial_structures.R | 1 - man/createNearestNetwork.Rd | 5 +- man/createNetwork.Rd | 75 ++++++++ 7 files changed, 367 insertions(+), 62 deletions(-) create mode 100644 man/createNetwork.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 93eae3de..3b832b44 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -79,10 +79,10 @@ Suggests: Remotes: drieslab/GiottoUtils Config/testthat/edition: 3 Collate: - 'NN_network.R' 'package_imports.R' 'classes.R' 'generics.R' + 'NN_network.R' 'aggregate.R' 'auxilliary.R' 'combine_metadata.R' diff --git a/NAMESPACE b/NAMESPACE index 4f3c1943..6ba0e54f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -329,6 +329,7 @@ exportMethods(colnames) exportMethods(copy) exportMethods(createGiottoPoints) exportMethods(createGiottoPolygon) +exportMethods(createNetwork) exportMethods(crop) exportMethods(dim) exportMethods(dimnames) @@ -358,6 +359,8 @@ exportMethods(wrap) import(GiottoUtils) import(data.table) import(utils) +importClassesFrom(DelayedArray,DelayedArray) +importClassesFrom(Matrix,Matrix) importClassesFrom(terra,SpatExtent) importClassesFrom(terra,SpatVector) importFrom(grDevices,dev.size) diff --git a/R/NN_network.R b/R/NN_network.R index 45fa1e14..fc33c3d1 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -1,20 +1,64 @@ +#' @include generics.R +NULL + + + + + +#' @name createNetwork +#' @title Create a network +#' @description Create networks from node values. +#' @param x data to treat as nodes. (`matrices` or `data.frame` inputs accepted) +#' `matrix` type inputs are assumed to be for creation of nearest neighbors. +#' `data.frame` input is assumed to be spatial type info, made for creation of +#' kNN or delaunay triangulation networks. +#' @param type specific type of network to create, given the input +#' @param method method used to create the type of network requested +#' @param max_distance maximum distance allowed to be connected +#' @param node_ids character. Node ID values to assign. If none are provided, +#' only integer indices will be used as node IDs. +#' @param include_weight include edge weight attribute in output +#' @param include_distance include edge distance attribute in output +#' @param weight_fun function to calculate weights based on distance if +#' `include_distance = TRUE`. Default is \eqn{weight = 1 / (1 + distance}) +#' @param as.igraph logical. Whether to return as `igraph`. Otherwise returns +#' as `data.table` +#' @param ... additional params to pass. See details section +#' @returns Either `igraph` if `as.igraph = TRUE` and `data.table` otherwise. +#' @details +#' This is a hub function for many different methods of finding nearest +#' neighbors. The `...` param can be used to pass additional params relevant +#' to specific methods. Additionally, some params used in only specific +#' functions are mentioned below. +#' - **k** number of neighbors to find for sNN and kNN type networks. Default is 30 +#' +NULL - - - - +#' @rdname createNetwork +#' @param minimum_shared minimum shared neighbors +#' @param top_shared keep at ... +#' @examples +#' e <- GiottoData::loadSubObjectMini("exprObj")[] # dgCMatrix +#' expr_igraph <- createNetwork(t(e), node_ids = colnames(e)) +#' +#' pca <- GiottoData::loadSubObjectMini("dimObj")[] +#' pca_igraph <- createNetwork(pca, node_ids = rownames(pca)) +#' @export setMethod( "createNetwork", signature("allMatrix"), function( x, - type = c("sNN", "kNN"), - k = 30, + type = c("sNN", "kNN", "delaunay"), + method = c("dbscan", "geometry", "RTriangle", "deldir"), # TODO + max_distance = NULL, # TODO + minimum_k = 0L, # TODO minimum_shared = 5, top_shared = 3, node_ids = NULL, - include_weight = TRUE, include_distance = TRUE, + include_weight = TRUE, + weight_fun = function(d) 1 / (1 + d), as.igraph = TRUE, ... ) @@ -24,7 +68,16 @@ setMethod( k <- as.integer(k) minimum_shared <- as.integer(minimum_shared) top_shared <- as.integer(top_shared) - type <- match.arg(type, choices = c("sNN", "kNN")) + type <- match.arg(type, choices = c("sNN", "kNN", "delaunay")) + method <- switch( + type, + "sNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), + "kNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), + "delaunay" = match.arg( + method, choices = c("geometry", "RTriangle", "deldir"), + several.ok = TRUE + ) + ) if (k >= nrow(x)) { k <- (nrow(x) - 1L) @@ -35,8 +88,6 @@ setMethod( # get common params alist <- list( x = x, - k = k, - node_ids = node_ids, include_weight = include_weight, include_distance = include_distance, ... @@ -44,38 +95,49 @@ setMethod( # generate network data.table network_dt <- switch( - type, - "kNN" = do.call(.net_dt_knn, args = alist), - "sNN" = do.call( + sprintf("%s:%s", type, method), + "kNN:dbscan" = do.call(.net_dt_knn, args = c(alist, k = k)), + "sNN:dbscan" = do.call( .net_dt_snn, args = c( - alist, list(top_shared = top_shared, minimum_shared = minimum_shared) + alist, + list(k = k, top_shared = top_shared, minimum_shared = minimum_shared) ) - ) + ), + "delaunay:deldir" = do.call(.net_dt_del_deldir, + args = alist)$delaunay_network_DT, + "delaunay:RTriangle" = do.call(.net_dt_del_rtriangle, + args = alist)$delaunay_network_DT, + "delaunay:geometry" = do.call(.net_dt_del_geometry, + args = alist)$delaunay_network_DT ) + # replace indices with node_ids if desired + if (!is.null(node_ids)) { + # if node_ids are provided, include as extra cols + names(node_ids) <- seq_along(node_ids) + network_dt[, "from" := node_ids[from]] + network_dt[, "to" := node_ids[to]] + } - # return early if igraph not required - if (!as.igraph) return(network_dt) - + ## outputs ## - ## convert to igraph object + # cols to include in output + keep_cols <- c("from", "to") - # cols to include in igraph object - igraph_cols <- if (is.null(node_ids)) { - c("from", "to") - } else { - c("from_cell_ID", "to_cell_ID") - } + all_index <- network_dt[, unique(unlist(.SD)), .SDcols = keep_cols] + if (include_weight) keep_cols <- c(keep_cols, "weight") + if (include_distance) keep_cols <- c(keep_cols, "distance") - all_index <- network_dt[, unique(unlist(.SD)), .SDcols = igraph_cols] - if (include_weight) igraph_cols <- c(igraph_cols, "weight") - if (include_distance) igraph_cols <- c(igraph_cols, "distance") + if (type == "sNN") keep_cols <- c(keep_cols, "shared", "rank") - if (type == "sNN") igraph_cols <- c(igraph_cols, "shared", "rank") + # return early if igraph not required + network_dt_final <- network_dt[, keep_cols, with = FALSE] + if (!as.igraph) return(network_dt_final) + ## convert to igraph object out <- igraph::graph_from_data_frame( - network_dt[, igraph_cols, with = FALSE], + network_dt_final, directed = TRUE, vertices = all_index ) @@ -87,7 +149,7 @@ setMethod( # x input is a matrix .net_dt_knn <- function( - x, k, node_ids, include_weight = TRUE, include_distance = TRUE, ... + x, k = 30, include_weight = TRUE, include_distance = TRUE, ... ) { # NSE vars from <- to <- NULL @@ -99,28 +161,22 @@ setMethod( to = as.vector(nn_network$id) ) - if (!is.null(node_ids)) { - # if node_ids are provided, include as extra cols - names(node_ids) <- seq_along(node_ids) - nn_network_dt[, "from_cell_ID" := node_ids[from]] - nn_network_dt[, "to_cell_ID" := node_ids[to]] - } - # optional info - if (include_weight) { - nn_network_dt[, "weight" := 1 / (1 + as.vector(nn_network$dist))] - } - if (include_distance) { + if (include_distance || include_weight) { nn_network_dt[, "distance" := as.vector(nn_network$dist)] } + if (include_weight) { + nn_network_dt[, "weight" := 1 / (1 + distance)] + } return(nn_network_dt) } # x input is a matrix .net_dt_snn <- function( - x, k, node_ids, include_weight = TRUE, include_distance = TRUE, - top_shared = 3, minimum_shared = 5, ... + x, k = 30, include_weight = TRUE, include_distance = TRUE, + top_shared = 3, minimum_shared = 5, weight_fun = function(d) 1 / (1 + d), + ... ) { # NSE vars from <- to <- shared <- NULL @@ -135,20 +191,13 @@ setMethod( ) snn_network_dt <- snn_network_dt[stats::complete.cases(snn_network_dt)] - if (!is.null(node_ids)) { - # if node_ids are provided, include as extra cols - names(node_ids) <- seq_along(node_ids) - snn_network_dt[, "from_cell_ID" := node_ids[from]] - snn_network_dt[, "to_cell_ID" := node_ids[to]] - } - # optional info - if (include_weight) { - snn_network_dt[, "weight" := 1 / (1 + as.vector(snn_network$dist))] - } - if (include_distance) { + if (include_distance || include_weight) { snn_network_dt[, "distance" := as.vector(snn_network$dist)] } + if (include_weight) { + snn_network_dt[, "weight" := weight_fun(distance)] + } # rank snn @@ -161,6 +210,188 @@ setMethod( return(snn_network_dt) } +.net_dt_del_geometry <- function( + x, include_weight = TRUE, include_distance = TRUE, options = "Pp", + weight_fun = function(d) 1 / (1 + d), ... +) { + package_check("geometry", repository = "CRAN:geometry") + + # data.table variables + from <- to <- NULL + + delaunay_simplex_mat <- geometry::delaunayn( + p = x, options = options, ... + ) + + geometry_obj <- list("delaunay_simplex_mat" = delaunay_simplex_mat) + edge_combs <- utils::combn(x = ncol(delaunay_simplex_mat), m = 2L) + delaunay_edges <- data.table::as.data.table(apply( + edge_combs, MARGIN = 1L, function(comb) delaunay_simplex_mat[, comb] + )) + + ### making sure of no duplication ### + delaunay_edges_dedup <- unique(delaunay_edges) + igraph_obj <- igraph::graph_from_edgelist(as.matrix(delaunay_edges_dedup)) + adj_obj <- igraph::as_adjacency_matrix(igraph_obj) + igraph_obj2 <- igraph::graph.adjacency(adj_obj) + delaunay_edges_dedup2 <- igraph::get.data.frame(igraph_obj2) + delaunay_network_dt <- data.table::as.data.table(delaunay_edges_dedup2) + delaunay_network_dt[, from := as.integer(from)] + delaunay_network_dt[, to := as.integer(to)] + data.table::setorder(delaunay_network_dt, from, to) + + # optional cols + if (include_distance || include_weight) { + delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), + .SDcols = c("from", "to")] + } + if (include_weight) { + delaunay_network_dt[, "weight" := weight_fun(distance)] + } + + out_object <- list( + "geometry_obj" = geometry_obj, + "delaunay_network_DT" = delaunay_network_dt + ) + return(out_object) +} + +.net_dt_del_rtriangle <- function( + x, include_weight = TRUE, include_distance = TRUE, + Y = TRUE, j = TRUE, S = 0, + weight_fun = function(d) 1 / (1 + d), + ... +) { + # NSE vars + from <- to <- NULL + + package_check("RTriangle", repository = "CRAN:RTriangle") + + rtriangle_obj <- RTriangle::triangulate( + RTriangle::pslg(x), + Y = Y, j = j, S = S, + ... + ) + + delaunay_network_dt <- data.table::data.table( + from = rtriangle_obj$E[, 1], + to = rtriangle_obj$E[, 2] + ) + + data.table::setorder(delaunay_network_dt, from, to) + + # optional cols + if (include_distance || include_weight) { + delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), + .SDcols = c("from", "to")] + } + if (include_weight) { + delaunay_network_dt[, "weight" := weight_fun(distance)] + } + + out_object <- list( + "RTriangle_obj" = rtriangle_obj, + "delaunay_network_DT" = delaunay_network_dt + ) + return(out_object) +} + +.net_dt_del_deldir <- function( + x, include_weight = TRUE, include_distance = TRUE, + weight_fun = function(d) 1 / (1 + d), + ... +) { + # NSE variables + from <- to <- NULL + + if (ncol(x) > 2L) { + .gstop("\'deldir\' delaunay method only applies to 2D data. + use method \'geometry\' or \'RTriangle\' instead") + } + + deldir_obj <- deldir::deldir(x = x, ...) + + delaunay_network_dt <- data.table::data.table( + from = deldir_obj$delsgs$ind1, + to = deldir_obj$delsgs$ind2 + ) + + data.table::setorder(delaunay_network_dt, from, to) + + # optional cols + if (include_distance || include_weight) { + delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), + .SDcols = c("from", "to")] + } + if (include_weight) { + delaunay_network_dt[, "weight" := weight_fun(distance)] + } + + out_object <- list( + "deldir_obj" = deldir_obj, + "delaunay_network_DT" = delaunay_network_dt + ) + return(out_object) +} + + + + +# Nodes row order is assumed to be the same as the network indices +#' @name edge_coords_array +#' @param x matrix of nodes info with coords +#' @param y network data.table +#' @param x_node_ids if y is indexed by character in from and to cols, then the +#' node IDs that apply to the coords in x must be supplied as a character vector +#' @keywords internal +.edge_coords_array <- function(x, y, x_node_ids = NULL) { + + checkmate::assert_matrix(x) + checkmate::assert_data_table(y) + + # if indexed by character + if (y[, is.character(from) && is.character(to)]) { + # try to match against the cell_ID col in nodes info + if (is.null(x_node_ids)) { + .gstop("y is indexed by node ID. + Node IDs for x must be provided as a vector to 'x_node_ids'") + } + # convert to int indexing (should match x by row) + y <- data.table::copy(y) + y[, from := match(from, x_node_ids)] + y[, to := match(to, x_node_ids)] + } + + edge_coords_array <- array( + dim = c(nrow(y), ncol(x), 2), + dimnames = list( + c(), + paste0("dim_", seq(ncol(x))), + c("start", "end") + ) + ) + + edge_coords_array[, ,1] <- x[y$from, ] + edge_coords_array[, ,2] <- x[y$to, ] + edge_coords_array <- aperm(edge_coords_array, perm = c(3,2,1)) + class(edge_coords_array) <- c("edge_coords_array", class(edge_coords_array)) + + return(edge_coords_array) +} + +# x should be an edge_coords_array +.calc_edge_dist <- function(x, method = "euclidean", ...) { + + checkmate::assert_class(x, "edge_coords_array") + + sapply( + seq(dim(x)[3L]), + function(pair_i) stats::dist(x[, , pair_i], method = method, ...) + ) +} + + + @@ -176,7 +407,6 @@ setMethod( #' @param name arbitrary name for NN network. Defaults to #' \[type\].\[dim_reduction_to_use\] #' @param feats_to_use if dim_reduction_to_use = NULL, which genes to use -#' @param genes_to_use deprecated, use feats_to_use #' @param expression_values expression values to use #' @param return_gobject boolean: return giotto object (default = TRUE) #' @param k number of k neighbors to use diff --git a/R/package_imports.R b/R/package_imports.R index 8320719a..2214ca23 100644 --- a/R/package_imports.R +++ b/R/package_imports.R @@ -28,6 +28,7 @@ #' @importClassesFrom terra SpatExtent #' @importClassesFrom terra SpatVector #' @importClassesFrom DelayedArray DelayedArray +#' @importClassesFrom Matrix Matrix #' @import GiottoUtils #' @import data.table #' @import utils diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 2aedda96..4aa41e39 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -97,7 +97,6 @@ convert_to_reduced_spatial_network <- function(full_spatial_network_DT) { stop("parameter networkDT can not be NULL \n") } - # number of spatial dimensions TODO: chech with Huipeng! # d2_or_d3 = match.arg(d2_or_d3, choices = c(2,3)) if (d2_or_d3 == 3) { diff --git a/man/createNearestNetwork.Rd b/man/createNearestNetwork.Rd index 5c74b183..b091ef1f 100644 --- a/man/createNearestNetwork.Rd +++ b/man/createNearestNetwork.Rd @@ -13,14 +13,13 @@ createNearestNetwork( dim_reduction_name = NULL, dimensions_to_use = 1:10, feats_to_use = NULL, - genes_to_use = NULL, expression_values = c("normalized", "scaled", "custom"), name = NULL, return_gobject = TRUE, k = 30, minimum_shared = 5, top_shared = 3, - verbose = T, + verbose = TRUE, ... ) } @@ -41,8 +40,6 @@ createNearestNetwork( \item{feats_to_use}{if dim_reduction_to_use = NULL, which genes to use} -\item{genes_to_use}{deprecated, use feats_to_use} - \item{expression_values}{expression values to use} \item{name}{arbitrary name for NN network. Defaults to diff --git a/man/createNetwork.Rd b/man/createNetwork.Rd new file mode 100644 index 00000000..7aef4896 --- /dev/null +++ b/man/createNetwork.Rd @@ -0,0 +1,75 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/NN_network.R +\name{createNetwork} +\alias{createNetwork} +\alias{createNetwork,allMatrix-method} +\title{Create a network} +\usage{ +\S4method{createNetwork}{allMatrix}( + x, + type = c("sNN", "kNN", "delaunay"), + method = c("dbscan", "geometry", "RTriangle", "deldir"), + max_distance = NULL, + minimum_k = 0L, + minimum_shared = 5, + top_shared = 3, + node_ids = NULL, + include_weight = TRUE, + include_distance = TRUE, + as.igraph = TRUE, + ... +) +} +\arguments{ +\item{x}{data to treat as nodes. (\code{matrices} or \code{data.frame} inputs accepted) +\code{matrix} type inputs are assumed to be for creation of nearest neighbors. +\code{data.frame} input is assumed to be spatial type info, made for creation of +kNN or delaunay triangulation networks.} + +\item{type}{specific type of network to create, given the input} + +\item{method}{method used to create the type of network requested} + +\item{max_distance}{maximum distance allowed to be connected} + +\item{minimum_shared}{minimum shared neighbors} + +\item{top_shared}{keep at ...} + +\item{node_ids}{character. Node ID values to assign. If none are provided, +only integer indices will be used as node IDs.} + +\item{include_weight}{include edge weight attribute in output} + +\item{include_distance}{include edge distance attribute in output} + +\item{as.igraph}{logical. Whether to return as \code{igraph}. Otherwise returns +as \code{data.table}} + +\item{...}{additional params to pass. See details section} + +\item{weight_fun}{function to calculate weights based on distance if +\code{include_distance = TRUE}. Default is \eqn{weight = 1 / (1 + distance})} +} +\value{ +Either \code{igraph} if \code{as.igraph = TRUE} and \code{data.table} otherwise. +} +\description{ +Create networks from node values. +} +\details{ +This is a hub function for many different methods of finding nearest +neighbors. The \code{...} param can be used to pass additional params relevant +to specific methods. Additionally, some params used in only specific +functions are mentioned below. +\itemize{ +\item \strong{k} number of neighbors to find for sNN and kNN type networks. Default is +} +} +\examples{ +e <- GiottoData::loadSubObjectMini("exprObj")[] # dgCMatrix +expr_igraph <- createNetwork(t(e), node_ids = colnames(e)) + +pca <- GiottoData::loadSubObjectMini("dimObj")[] +pca_igraph <- createNetwork(pca, node_ids = rownames(pca)) +} From 1092bcf38019feff66c1ccf6c7cd53ec4f395167 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 01:00:41 -0500 Subject: [PATCH 03/24] chore: remove stray browser() --- R/subset.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/subset.R b/R/subset.R index f517516a..eb89892f 100644 --- a/R/subset.R +++ b/R/subset.R @@ -942,7 +942,7 @@ - # browser() + # # mirai # mirai_res = get_mirai_list(gobject) # if (length(mirai_res > 0)) { From 807839dd3945a42526d9d4386282959975f7a736 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 01:01:57 -0500 Subject: [PATCH 04/24] chore: use setDT instead of as.data.table --- R/spatial_structures.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 4aa41e39..8461d1a9 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -1071,7 +1071,7 @@ create_KNNnetwork_dbscan <- function(spatial_locations, distance = as.vector(knn_spatial$dist) ) nw_sptial.norm <- igraph::graph_from_data_frame(knn_sptial.norm, directed = FALSE) - network_DT <- data.table::as.data.table(knn_sptial.norm) + network_DT <- data.table::setDT(knn_sptial.norm) # spatial_network_DT[, `:=`(from, cell_ID_vec[from])] From 4ae1cf32b42caa3971d999dabb8885daafb1bef6 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 01:03:40 -0500 Subject: [PATCH 05/24] chore: improve docs and comments for `convert_to_full_spatial_network()` --- R/spatial_structures.R | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 8461d1a9..a71360cd 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -3,7 +3,13 @@ #' @title convert_to_full_spatial_network #' @name convert_to_full_spatial_network -#' @description convert to a full spatial network +#' @description Convert to a full spatial network, ie ensuring that all edges +#' that may currently only be represented as \eqn{a -> b} also have the reverse +#' \eqn{b -> a}. The entries are then made unique, after which all interactions +#' are ranked by distance, where rank increases from smaller to larger distances. +#' This rank is appended to the `data.table` as a `rank_int` column. Another +#' `rnk_src_trgt` column is added with the IDs of \eqn{a} and \eqn{b} pasted +#' together #' @param reduced_spatial_network_DT reduced spatial network in data.table format #' @keywords internal #' @export @@ -11,20 +17,22 @@ convert_to_full_spatial_network <- function(reduced_spatial_network_DT) { # data.table variables distance <- rank_int <- NULL - # find location coordinates - coordinates <- grep("sdim", colnames(reduced_spatial_network_DT), value = T) + # find location coordinates cols + coordinates <- grep("sdim", colnames(reduced_spatial_network_DT), value = TRUE) - begin_coordinates <- grep("begin", coordinates, value = T) + # convert names from sdimx_being and sdimy_begin to source_x and source_y + begin_coordinates <- grep("begin", coordinates, value = TRUE) new_begin_coordinates <- gsub(x = begin_coordinates, pattern = "_begin", replacement = "") new_begin_coordinates <- gsub(x = new_begin_coordinates, pattern = "sdim", replacement = "source_") - end_coordinates <- grep("end", coordinates, value = T) + # convert names from sdimx_end and sdimy_end to target_x and target_y + end_coordinates <- grep("end", coordinates, value = TRUE) new_end_coordinates <- gsub(x = end_coordinates, pattern = "_end", replacement = "") new_end_coordinates <- gsub(x = new_end_coordinates, pattern = "sdim", replacement = "target_") # create normal source --> target part1 <- data.table::copy(reduced_spatial_network_DT) - part1 <- part1[, c("from", "to", begin_coordinates, end_coordinates, "distance", "weight"), with = F] + part1 <- part1[, c("from", "to", begin_coordinates, end_coordinates, "distance", "weight"), with = FALSE] colnames(part1) <- c("source", "target", new_begin_coordinates, new_end_coordinates, "distance", "weight") # revert order target (now source) --> source (now target) @@ -35,11 +43,12 @@ convert_to_full_spatial_network <- function(reduced_spatial_network_DT) { full_spatial_network_DT <- rbind(part1, part2) full_spatial_network_DT <- unique(full_spatial_network_DT) - # create ranking of interactions + # create ranking of interactions by distance per source + # the lower the ranking, the shorter the distance data.table::setorder(full_spatial_network_DT, source, distance) full_spatial_network_DT[, rank_int := 1:.N, by = "source"] - # create unified column + # create unified column for source and target as rnk_src_trgt full_spatial_network_DT <- dt_sort_combine_two_columns(full_spatial_network_DT, "source", "target", "rnk_src_trgt") return(full_spatial_network_DT) From c4c25f757d4ee05395ac3bca0970c4ce6c91e04e Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 01:15:05 -0500 Subject: [PATCH 06/24] chore: improve docs for `convert_to_full_spatial_network()` --- R/spatial_structures.R | 14 +++++++------- man/convert_to_full_spatial_network.Rd | 10 ++++++++-- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/R/spatial_structures.R b/R/spatial_structures.R index a71360cd..81f7e504 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -4,13 +4,13 @@ #' @title convert_to_full_spatial_network #' @name convert_to_full_spatial_network #' @description Convert to a full spatial network, ie ensuring that all edges -#' that may currently only be represented as \eqn{a -> b} also have the reverse -#' \eqn{b -> a}. The entries are then made unique, after which all interactions -#' are ranked by distance, where rank increases from smaller to larger distances. -#' This rank is appended to the `data.table` as a `rank_int` column. Another -#' `rnk_src_trgt` column is added with the IDs of \eqn{a} and \eqn{b} pasted -#' together -#' @param reduced_spatial_network_DT reduced spatial network in data.table format +#' that may currently only be represented as \eqn{a} -> \eqn{b} also have the +#' reverse \eqn{b} -> \eqn{a}. The entries are then made unique, after which all +#' interactions are ranked by distance, where rank increases from smaller to +#' larger distances. This rank is appended to the `data.table` as a `rank_int` +#' column. Another `rnk_src_trgt` column is added with the IDs of \eqn{a} and +#' \eqn{b} pasted together +#' @param reduced_spatial_network_DT reduced spatial network in `data.table` format #' @keywords internal #' @export convert_to_full_spatial_network <- function(reduced_spatial_network_DT) { diff --git a/man/convert_to_full_spatial_network.Rd b/man/convert_to_full_spatial_network.Rd index 426b8408..59308aaf 100644 --- a/man/convert_to_full_spatial_network.Rd +++ b/man/convert_to_full_spatial_network.Rd @@ -7,9 +7,15 @@ convert_to_full_spatial_network(reduced_spatial_network_DT) } \arguments{ -\item{reduced_spatial_network_DT}{reduced spatial network in data.table format} +\item{reduced_spatial_network_DT}{reduced spatial network in \code{data.table} format} } \description{ -convert to a full spatial network +Convert to a full spatial network, ie ensuring that all edges +that may currently only be represented as \eqn{a} -> \eqn{b} also have the +reverse \eqn{b} -> \eqn{a}. The entries are then made unique, after which all +interactions are ranked by distance, where rank increases from smaller to +larger distances. This rank is appended to the \code{data.table} as a \code{rank_int} +column. Another \code{rnk_src_trgt} column is added with the IDs of \eqn{a} and +\eqn{b} pasted together } \keyword{internal} From 6f635cc37839207bca648d9bc642ab09adb276e3 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 01:43:20 -0500 Subject: [PATCH 07/24] feat: add compatibility for spatial network DTs that do not have coords info in `convert_to_reduced_spatial_network()` --- R/spatial_structures.R | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 81f7e504..6b22817c 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -68,6 +68,22 @@ convert_to_reduced_spatial_network <- function(full_spatial_network_DT) { reduced_spatial_network_DT <- full_spatial_network_DT[!duplicated(rnk_src_trgt)] reduced_spatial_network_DT[, c("rank_int", "rnk_src_trgt") := NULL] # don't make sense in a reduced network + # TODO moving forward, coords info start/end may not be included 24.02.12 + has_coords <- any(grepl( + "source_|target_", + x = colnames(reduced_spatial_network_DT) + )) + + if (!has_coords) { + data.table::setnames( + reduced_spatial_network_DT, + old = c("source", "target"), + new = c("from", "to") + ) + return(reduced_spatial_network_DT) + } + + # return col names to sdimx/sdimy naming scheme # convert to names for a reduced network source_coordinates <- grep("source_", colnames(reduced_spatial_network_DT), value = T) new_source_coordinates <- gsub(x = source_coordinates, pattern = "source_", replacement = "sdim") From 1ea0e514526de3fed1934cb15e94282e2ce3591e Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 13:42:51 -0500 Subject: [PATCH 08/24] chore: document `convert_to_reduced_spatial_network()` --- R/spatial_structures.R | 4 +++- man/convert_to_reduced_spatial_network.Rd | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 6b22817c..715449ae 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -56,7 +56,9 @@ convert_to_full_spatial_network <- function(reduced_spatial_network_DT) { #' @title convert_to_reduced_spatial_network #' @name convert_to_reduced_spatial_network -#' @description convert to a reduced spatial network +#' @description Convert to a reduced spatial network. Specifically, removes +#' the duplicated connections so that only \eqn{a} -> \eqn{b} interactions +#' remain. #' @param full_spatial_network_DT full spatial network in data.table format #' @keywords internal #' @export diff --git a/man/convert_to_reduced_spatial_network.Rd b/man/convert_to_reduced_spatial_network.Rd index 2c74387b..80b73693 100644 --- a/man/convert_to_reduced_spatial_network.Rd +++ b/man/convert_to_reduced_spatial_network.Rd @@ -10,6 +10,8 @@ convert_to_reduced_spatial_network(full_spatial_network_DT) \item{full_spatial_network_DT}{full spatial network in data.table format} } \description{ -convert to a reduced spatial network +Convert to a reduced spatial network. Specifically, removes +the duplicated connections so that only \eqn{a} -> \eqn{b} interactions +remain. } \keyword{internal} From 3b0e6d2460d2d4248645b9b9930d217fce9b12d6 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 13:45:42 -0500 Subject: [PATCH 09/24] feat: edge distance calculation `edge_distances()` Add function to calculate euclidean distances for any m x n matrix when provided a data.table of `from` and `to` interactions --- NAMESPACE | 1 + R/NN_network.R | 35 +++++++++++++++++++++++++++++++++-- man/dot-edge_coords_array.Rd | 24 ++++++++++++++++++++++++ man/edge_distances.Rd | 29 +++++++++++++++++++++++++++++ 4 files changed, 87 insertions(+), 2 deletions(-) create mode 100644 man/dot-edge_coords_array.Rd create mode 100644 man/edge_distances.Rd diff --git a/NAMESPACE b/NAMESPACE index 6ba0e54f..d5786a86 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -100,6 +100,7 @@ export(create_spat_locs_obj) export(create_spat_net_obj) export(cropGiottoLargeImage) export(distGiottoImage) +export(edge_distances) export(estimateImageBg) export(evaluate_input) export(fDataDT) diff --git a/R/NN_network.R b/R/NN_network.R index fc33c3d1..5ed678bd 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -336,11 +336,41 @@ setMethod( +# distances calculation #### + +#' @name edge_distances +#' @title Calculate network edge euclidean distances +#' @param x matrix of nodes info with coords. Rows should be samples, Cols +#' should be variables +#' @param y network data.table with `from` and `to` cols. Usually integer +#' indices matching the rows of x. +#' @param x_node_ids if y is indexed by character in from and to cols, then the +#' node IDs that apply to the coords in x must be supplied as a character vector +#' @examples +#' m <- matrix(c(0,0,0, 1,1,1, 3,2,4), byrow = TRUE, nrow = 3) +#' edges <- data.table( +#' from = c(1, 1), +#' to = c(2, 3) +#' ) +#' edge_distances(m, edges) +#' @export +edge_distances <- function(x, y, x_node_ids = NULL) { + .calc_edge_dist(.edge_coords_array(x, y)) +} + + # Nodes row order is assumed to be the same as the network indices -#' @name edge_coords_array +#' @title Numerical array of edge start and end +#' @name .edge_coords_array +#' @description +#' Generate a \eqn{2} x \eqn{j} x \eqn{k} numerical array of edge start and end +#' coordinates. Rows correspond to start and end. Cols are for each variable +#' ie x, y, (z) or whatever other variable is used to measure sample location +#' in graph space. The third dim is for each sample. This layout makes it easy +#' to iterate across matrix slices of this array with `[stats::dist()]`. #' @param x matrix of nodes info with coords -#' @param y network data.table +#' @param y network data.table with `from` and `to` cols #' @param x_node_ids if y is indexed by character in from and to cols, then the #' node IDs that apply to the coords in x must be supplied as a character vector #' @keywords internal @@ -394,6 +424,7 @@ setMethod( +# original implementations #### #' @title createNearestNetwork diff --git a/man/dot-edge_coords_array.Rd b/man/dot-edge_coords_array.Rd new file mode 100644 index 00000000..6ff753a7 --- /dev/null +++ b/man/dot-edge_coords_array.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/NN_network.R +\name{.edge_coords_array} +\alias{.edge_coords_array} +\title{Numerical array of edge start and end} +\usage{ +.edge_coords_array(x, y, x_node_ids = NULL) +} +\arguments{ +\item{x}{matrix of nodes info with coords} + +\item{y}{network data.table with \code{from} and \code{to} cols} + +\item{x_node_ids}{if y is indexed by character in from and to cols, then the +node IDs that apply to the coords in x must be supplied as a character vector} +} +\description{ +Generate a \eqn{2} x \eqn{j} x \eqn{k} numerical array of edge start and end +coordinates. Rows correspond to start and end. Cols are for each variable +ie x, y, (z) or whatever other variable is used to measure sample location +in graph space. The third dim is for each sample. This layout makes it easy +to iterate across matrix slices of this array with \verb{[stats::dist()]}. +} +\keyword{internal} diff --git a/man/edge_distances.Rd b/man/edge_distances.Rd new file mode 100644 index 00000000..a458ddaf --- /dev/null +++ b/man/edge_distances.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/NN_network.R +\name{edge_distances} +\alias{edge_distances} +\title{Calculate network edge euclidean distances} +\usage{ +edge_distances(x, y, x_node_ids = NULL) +} +\arguments{ +\item{x}{matrix of nodes info with coords. Rows should be samples, Cols +should be variables} + +\item{y}{network data.table with \code{from} and \code{to} cols. Usually integer +indices matching the rows of x.} + +\item{x_node_ids}{if y is indexed by character in from and to cols, then the +node IDs that apply to the coords in x must be supplied as a character vector} +} +\description{ +Calculate network edge euclidean distances +} +\examples{ +m <- matrix(c(0,0,0, 1,1,1, 3,2,4), byrow = TRUE, nrow = 3) +edges <- data.table( + from = c(1, 1), + to = c(2, 3) +) +edge_distances(m, edges) +} From 734b63c310d5b9580388a9be492e6712bf28beda Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 13:49:45 -0500 Subject: [PATCH 10/24] feat: `createNetwork()` Refactors Giotto's spatial and NN networks creation into a hub function that allows starting directly from matrix data. Added tests that currently only ensure parity with original implementation. --- R/NN_network.R | 278 ++++++++++++++++------- R/spatial_structures.R | 2 +- man/createNetwork.Rd | 145 +++++++++--- tests/testthat/test-networks.R | 197 ++++++++++++++++ tests/testthat/test-spatial_structures.R | 31 --- 5 files changed, 510 insertions(+), 143 deletions(-) create mode 100644 tests/testthat/test-networks.R delete mode 100644 tests/testthat/test-spatial_structures.R diff --git a/R/NN_network.R b/R/NN_network.R index 5ed678bd..bd6996bb 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -3,46 +3,130 @@ NULL - +# createNetwork #### #' @name createNetwork #' @title Create a network -#' @description Create networks from node values. -#' @param x data to treat as nodes. (`matrices` or `data.frame` inputs accepted) -#' `matrix` type inputs are assumed to be for creation of nearest neighbors. -#' `data.frame` input is assumed to be spatial type info, made for creation of -#' kNN or delaunay triangulation networks. -#' @param type specific type of network to create, given the input -#' @param method method used to create the type of network requested -#' @param max_distance maximum distance allowed to be connected -#' @param node_ids character. Node ID values to assign. If none are provided, -#' only integer indices will be used as node IDs. -#' @param include_weight include edge weight attribute in output -#' @param include_distance include edge distance attribute in output -#' @param weight_fun function to calculate weights based on distance if -#' `include_distance = TRUE`. Default is \eqn{weight = 1 / (1 + distance}) +#' @description Create networks from node values. This is a hub function for +#' many different methods of finding nearest neighbors. See details for additional +#' params important for generating specific types of networks. +#' @param x matrix. Data to treat as nodes. Examples include expression +#' information, PCA matrix, spatial locations xy(z) coordinates. +#' @param type type of network to create. Currently: "sNN", "kNN", or "delaunay". +#' @param method method used to create the type of network requested. +#' One of "dbscan" for sNN and kNN or "geometry", "RTriangle", or "deldir" +#' for delaunay. +#' @param node_ids character. Node ID values to assign. If NULL, integer indices +#' will be used as node IDs. +#' @param include_weight logical. include edge weight attribute in output +#' @param include_distance logical. include edge distance attribute in output #' @param as.igraph logical. Whether to return as `igraph`. Otherwise returns #' as `data.table` #' @param ... additional params to pass. See details section #' @returns Either `igraph` if `as.igraph = TRUE` and `data.table` otherwise. #' @details -#' This is a hub function for many different methods of finding nearest -#' neighbors. The `...` param can be used to pass additional params relevant -#' to specific methods. Additionally, some params used in only specific -#' functions are mentioned below. -#' - **k** number of neighbors to find for sNN and kNN type networks. Default is 30 +#' Additional params are described below. Items in parenthesis refer to which +#' network types and/or methods the params are specific to. +#' - \[**`k`**\] numeric. (*sNN, kNN*) number of neighbors to find. Default is 30 +#' - \[**`minimum_shared`**\] numeric. (*sNN*) minimum shared neighbors allowed +#' per edge +#' - \[**`top_shared`**\] numeric. (*sNN*) keep at least this many edges per node, +#' where kept edges are top ranked in terms of number of shared neighbors. +#' - \[**`filter`**\] logical. (*kNN*) whether to filter for only unique +#' edges and apply `minimum_k` and `maximum_distance` filters. Should be set +#' `TRUE` when generating a spatial kNN network. Default is `FALSE.` +#' - \[**`minimum_k`**\] (*delaunay, kNN*) minimum nearest neighbours if +#' `maximum_distance != NULL` +#' - \[**`maximum_distance`**\] (*delaunay, kNN*) edge maximum euclidean +#' distance allowed +#' - \[**`Y`**\] (*RTriangle*) If TRUE prohibits the insertion of Steiner points +#' on the mesh boundary. Default is TRUE +#' - \[**`j`**\] (*RTriangle*) If TRUE jettisons vertices that are not part of +#' the final triangulation from the output. Default is TRUE +#' - \[**`S`**\] (*RTriangle*) Specifies the maximum number of added Steiner +#' points. Default is 0 +#' - \[**`options`**\] (*geometry*) default is "Pp". See [geometry::delaunayn] +#' - \[**`weight_fun`**\] function to calculate weights based on distance if +#' `include_weight = TRUE`. Default is \eqn{weight = 1 / (1 + distance)} for +#' `"kNN"` and `"sNN"` types and \eqn{weight = 1 / distance} for `delaunay` +#' type +#' networks +#' @examples +#' \dontrun{ +#' pca <- GiottoData::loadSubObjectMini("dimObj")[] +#' sl <- GiottoData::loadSubObjectMini("spatLocsObj")[] +#' +#' # Delaunay via geometry::delaunayn() +#' del_geom <- createNetwork( +#' x = as.matrix(sl[, .(sdimx, sdimy)]), +#' type = "delaunay", +#' method = "geometry", +#' include_weight = TRUE, +#' weight_fun = function(d) 1 / d, +#' as.igraph = FALSE, +#' node_ids = sl$cell_ID +#' ) +#' +#' # Delaunay via RTriangle::triangulate() +#' del_rt <- createNetwork( +#' x = as.matrix(sl[, .(sdimx, sdimy)]), +#' type = "delaunay", +#' method = "RTriangle", +#' include_weight = TRUE, +#' weight_fun = function(d) 1 / d, +#' as.igraph = FALSE, +#' node_ids = sl$cell_ID +#' ) +#' +#' # Delaunay via deldir::deldir() +#' del_dd <- createNetwork( +#' x = as.matrix(sl[, .(sdimx, sdimy)]), +#' type = "delaunay", +#' method = "deldir", +#' include_weight = TRUE, +#' weight_fun = function(d) 1 / d, +#' as.igraph = FALSE, +#' node_ids = sl$cell_ID +#' ) +#' +#' # kNN spatial network +#' kNN_spat <- createNetwork( +#' x = as.matrix(sl[, .(sdimx, sdimy)]), +#' type = "kNN", +#' method = "dbscan", +#' include_weight = TRUE, +#' weight_fun = function(d) 1 / d, # not the default +#' as.igraph = FALSE, +#' node_ids = sl$cell_ID, +#' k = 4L, +#' maximum_distance = NULL, +#' minimum_k = 0L +#' ) +#' +#' # kNN NN network +#' kNN <- createNetwork( +#' pca[, 1:10], +#' type = "kNN", +#' method = "dbscan", +#' node_ids = rownames(pca), +#' as.igraph = TRUE +#' ) +#' +#' # sNN NN network +#' sNN <- createNetwork( +#' pca[, 1:10], +#' type = "sNN", +#' method = "dbscan", +#' node_ids = rownames(pca), +#' as.igraph = TRUE +#' ) #' +#' # using defaults for sNN with index IDs to create igraph +#' sNN_idx <- createNetwork(pca[, 1:10]) +#' } NULL #' @rdname createNetwork -#' @param minimum_shared minimum shared neighbors -#' @param top_shared keep at ... -#' @examples -#' e <- GiottoData::loadSubObjectMini("exprObj")[] # dgCMatrix -#' expr_igraph <- createNetwork(t(e), node_ids = colnames(e)) -#' -#' pca <- GiottoData::loadSubObjectMini("dimObj")[] -#' pca_igraph <- createNetwork(pca, node_ids = rownames(pca)) #' @export setMethod( "createNetwork", @@ -50,24 +134,16 @@ setMethod( function( x, type = c("sNN", "kNN", "delaunay"), - method = c("dbscan", "geometry", "RTriangle", "deldir"), # TODO - max_distance = NULL, # TODO - minimum_k = 0L, # TODO - minimum_shared = 5, - top_shared = 3, + method = c("dbscan", "geometry", "RTriangle", "deldir"), node_ids = NULL, include_distance = TRUE, include_weight = TRUE, - weight_fun = function(d) 1 / (1 + d), as.igraph = TRUE, ... ) { # check params - k <- as.integer(k) - minimum_shared <- as.integer(minimum_shared) - top_shared <- as.integer(top_shared) type <- match.arg(type, choices = c("sNN", "kNN", "delaunay")) method <- switch( type, @@ -79,31 +155,18 @@ setMethod( ) ) - if (k >= nrow(x)) { - k <- (nrow(x) - 1L) - vmsg(.v = verbose, "k is higher than total number of cells. - Adjusted to (total number of cells - 1)") - } - # get common params alist <- list( x = x, include_weight = include_weight, - include_distance = include_distance, ... ) # generate network data.table network_dt <- switch( sprintf("%s:%s", type, method), - "kNN:dbscan" = do.call(.net_dt_knn, args = c(alist, k = k)), - "sNN:dbscan" = do.call( - .net_dt_snn, - args = c( - alist, - list(k = k, top_shared = top_shared, minimum_shared = minimum_shared) - ) - ), + "kNN:dbscan" = do.call(.net_dt_knn, args = alist), + "sNN:dbscan" = do.call(.net_dt_snn, args = alist), "delaunay:deldir" = do.call(.net_dt_del_deldir, args = alist)$delaunay_network_DT, "delaunay:RTriangle" = do.call(.net_dt_del_rtriangle, @@ -149,11 +212,23 @@ setMethod( # x input is a matrix .net_dt_knn <- function( - x, k = 30, include_weight = TRUE, include_distance = TRUE, ... + x, k = 30L, include_weight = TRUE, include_distance = TRUE, filter = FALSE, + maximum_distance = NULL, minimum_k = 0L, weight_fun = function(d) 1 / (1 + d), + ... ) { # NSE vars from <- to <- NULL + k <- as.integer(k) + + if (k >= nrow(x)) { + k <- (nrow(x) - 1L) + vmsg(.v = verbose, "k is higher than total number of cells. + Adjusted to (total number of cells - 1)") + } + # distances must be calculated when a limit is set + if (!is.null(maximum_distance)) include_distance <- TRUE + nn_network <- dbscan::kNN(x = x, k = k, sort = TRUE, ...) nn_network_dt <- data.table::data.table( @@ -163,10 +238,29 @@ setMethod( # optional info if (include_distance || include_weight) { - nn_network_dt[, "distance" := as.vector(nn_network$dist)] + if (!is.null(maximum_distance)) { + # maximum_distance flag treated as a flag to use this function for + # spatial network purposes. + # + # Use the input matrix coords instead of those exported from dbscan + # needed for filtering + nn_network_dt[, "distance" := edge_distances(x, .SD), + .SDcols = c("from", "to")] + } else { + nn_network_dt[, "distance" := as.vector(nn_network$dist)] + } } if (include_weight) { - nn_network_dt[, "weight" := 1 / (1 + distance)] + nn_network_dt[, "weight" := weight_fun(distance)] + } + + # filtering by distance and min k is done when maximum_distance is not NULL + if (filter) { + nn_network_dt <- .filter_network( + networkDT = nn_network_dt, + maximum_distance = maximum_distance, + minimum_k = minimum_k + ) } return(nn_network_dt) @@ -174,13 +268,23 @@ setMethod( # x input is a matrix .net_dt_snn <- function( - x, k = 30, include_weight = TRUE, include_distance = TRUE, - top_shared = 3, minimum_shared = 5, weight_fun = function(d) 1 / (1 + d), + x, k = 30L, include_weight = TRUE, include_distance = TRUE, + top_shared = 3L, minimum_shared = 5L, weight_fun = function(d) 1 / (1 + d), ... ) { # NSE vars from <- to <- shared <- NULL + k <- as.integer(k) + top_shared <- as.integer(top_shared) + minimum_shared <- as.integer(minimum_shared) + + if (k >= nrow(x)) { + k <- (nrow(x) - 1L) + vmsg(.v = verbose, "k is higher than total number of cells. + Adjusted to (total number of cells - 1)") + } + nn_network <- dbscan::kNN(x = x, k = k, sort = TRUE, ...) snn_network <- dbscan::sNN(x = nn_network, k = k, kt = NULL, ...) @@ -200,19 +304,23 @@ setMethod( } - # rank snn + # rank snn. LOWER ranking means MORE shared per source data.table::setorder(snn_network_dt, from, -shared) snn_network_dt[, rank := 1:.N, by = from] # filter snn + # keep at at least `top_shared` - 1 interactions where the ones selected should + # have more connections than the cutoff. Also keep any interactions with more + # shared than `minimum_shared` snn_network_dt <- snn_network_dt[rank <= top_shared | shared >= minimum_shared] return(snn_network_dt) } .net_dt_del_geometry <- function( - x, include_weight = TRUE, include_distance = TRUE, options = "Pp", - weight_fun = function(d) 1 / (1 + d), ... + x, include_weight = TRUE, options = "Pp", maximum_distance = "auto", + minimum_k = 0L, weight_fun = function(d) 1 / d, + ... ) { package_check("geometry", repository = "CRAN:geometry") @@ -240,15 +348,21 @@ setMethod( delaunay_network_dt[, to := as.integer(to)] data.table::setorder(delaunay_network_dt, from, to) + # needed for filtering + delaunay_network_dt[, "distance" := edge_distances(x, .SD), + .SDcols = c("from", "to")] + # optional cols - if (include_distance || include_weight) { - delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), - .SDcols = c("from", "to")] - } if (include_weight) { delaunay_network_dt[, "weight" := weight_fun(distance)] } + delaunay_network_dt <- .filter_network( + networkDT = delaunay_network_dt, + maximum_distance = maximum_distance, + minimum_k = minimum_k + ) + out_object <- list( "geometry_obj" = geometry_obj, "delaunay_network_DT" = delaunay_network_dt @@ -257,9 +371,8 @@ setMethod( } .net_dt_del_rtriangle <- function( - x, include_weight = TRUE, include_distance = TRUE, - Y = TRUE, j = TRUE, S = 0, - weight_fun = function(d) 1 / (1 + d), + x, include_weight = TRUE, maximum_distance = "auto", minimum_k = 0L, + Y = TRUE, j = TRUE, S = 0, weight_fun = function(d) 1 / d, ... ) { # NSE vars @@ -280,15 +393,21 @@ setMethod( data.table::setorder(delaunay_network_dt, from, to) + # needed for filtering + delaunay_network_dt[, "distance" := edge_distances(x, .SD), + .SDcols = c("from", "to")] + # optional cols - if (include_distance || include_weight) { - delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), - .SDcols = c("from", "to")] - } if (include_weight) { delaunay_network_dt[, "weight" := weight_fun(distance)] } + delaunay_network_dt <- .filter_network( + networkDT = delaunay_network_dt, + maximum_distance = maximum_distance, + minimum_k = minimum_k + ) + out_object <- list( "RTriangle_obj" = rtriangle_obj, "delaunay_network_DT" = delaunay_network_dt @@ -297,8 +416,8 @@ setMethod( } .net_dt_del_deldir <- function( - x, include_weight = TRUE, include_distance = TRUE, - weight_fun = function(d) 1 / (1 + d), + x, include_weight = TRUE, maximum_distance = "auto", minimum_k = 0L, + weight_fun = function(d) 1 / d, ... ) { # NSE variables @@ -318,15 +437,20 @@ setMethod( data.table::setorder(delaunay_network_dt, from, to) + delaunay_network_dt[, "distance" := edge_distances(x, .SD), + .SDcols = c("from", "to")] + # optional cols - if (include_distance || include_weight) { - delaunay_network_dt[, "distance" := .calc_edge_dist(.edge_coords_array(x, .SD)), - .SDcols = c("from", "to")] - } if (include_weight) { delaunay_network_dt[, "weight" := weight_fun(distance)] } + delaunay_network_dt <- .filter_network( + networkDT = delaunay_network_dt, + maximum_distance = maximum_distance, + minimum_k = minimum_k + ) + out_object <- list( "deldir_obj" = deldir_obj, "delaunay_network_DT" = delaunay_network_dt diff --git a/R/spatial_structures.R b/R/spatial_structures.R index 715449ae..a829e126 100644 --- a/R/spatial_structures.R +++ b/R/spatial_structures.R @@ -205,7 +205,7 @@ get_distance <- function(networkDT, temp_fullnetwork <- convert_to_full_spatial_network(networkDT) ## filter based on distance or minimum number of neighbors - if (maximum_distance == "auto") { + if (isTRUE(maximum_distance == "auto")) { temp_fullnetwork <- temp_fullnetwork[distance <= grDevices::boxplot.stats(temp_fullnetwork$distance)$stats[5] | rank_int <= minimum_k] } else if (!is.null(maximum_distance)) { temp_fullnetwork <- temp_fullnetwork[distance <= maximum_distance | rank_int <= minimum_k] diff --git a/man/createNetwork.Rd b/man/createNetwork.Rd index 7aef4896..cdf5251c 100644 --- a/man/createNetwork.Rd +++ b/man/createNetwork.Rd @@ -9,67 +9,144 @@ x, type = c("sNN", "kNN", "delaunay"), method = c("dbscan", "geometry", "RTriangle", "deldir"), - max_distance = NULL, - minimum_k = 0L, - minimum_shared = 5, - top_shared = 3, node_ids = NULL, - include_weight = TRUE, include_distance = TRUE, + include_weight = TRUE, as.igraph = TRUE, ... ) } \arguments{ -\item{x}{data to treat as nodes. (\code{matrices} or \code{data.frame} inputs accepted) -\code{matrix} type inputs are assumed to be for creation of nearest neighbors. -\code{data.frame} input is assumed to be spatial type info, made for creation of -kNN or delaunay triangulation networks.} - -\item{type}{specific type of network to create, given the input} - -\item{method}{method used to create the type of network requested} +\item{x}{matrix. Data to treat as nodes. Examples include expression +information, PCA matrix, spatial locations xy(z) coordinates.} -\item{max_distance}{maximum distance allowed to be connected} +\item{type}{type of network to create. Currently: "sNN", "kNN", or "delaunay".} -\item{minimum_shared}{minimum shared neighbors} +\item{method}{method used to create the type of network requested. +One of "dbscan" for sNN and kNN or "geometry", "RTriangle", or "deldir" +for delaunay.} -\item{top_shared}{keep at ...} +\item{node_ids}{character. Node ID values to assign. If NULL, integer indices +will be used as node IDs.} -\item{node_ids}{character. Node ID values to assign. If none are provided, -only integer indices will be used as node IDs.} +\item{include_distance}{logical. include edge distance attribute in output} -\item{include_weight}{include edge weight attribute in output} - -\item{include_distance}{include edge distance attribute in output} +\item{include_weight}{logical. include edge weight attribute in output} \item{as.igraph}{logical. Whether to return as \code{igraph}. Otherwise returns as \code{data.table}} \item{...}{additional params to pass. See details section} - -\item{weight_fun}{function to calculate weights based on distance if -\code{include_distance = TRUE}. Default is \eqn{weight = 1 / (1 + distance})} } \value{ Either \code{igraph} if \code{as.igraph = TRUE} and \code{data.table} otherwise. } \description{ -Create networks from node values. +Create networks from node values. This is a hub function for +many different methods of finding nearest neighbors. See details for additional +params important for generating specific types of networks. } \details{ -This is a hub function for many different methods of finding nearest -neighbors. The \code{...} param can be used to pass additional params relevant -to specific methods. Additionally, some params used in only specific -functions are mentioned below. +Additional params are described below. Items in parenthesis refer to which +network types and/or methods the params are specific to. \itemize{ -\item \strong{k} number of neighbors to find for sNN and kNN type networks. Default is +\item [\strong{\code{k}}] numeric. (\emph{sNN, kNN}) number of neighbors to find. Default is 30 +\item [\strong{\code{minimum_shared}}] numeric. (\emph{sNN}) minimum shared neighbors allowed +per edge +\item [\strong{\code{top_shared}}] numeric. (\emph{sNN}) keep at least this many edges per node, +where kept edges are top ranked in terms of number of shared neighbors. +\item [\strong{\code{filter}}] logical. (\emph{kNN}) whether to filter for only unique +edges and apply \code{minimum_k} and \code{maximum_distance} filters. Should be set +\code{TRUE} when generating a spatial kNN network. Default is \code{FALSE.} +\item [\strong{\code{minimum_k}}] (\emph{delaunay, kNN}) minimum nearest neighbours if +\code{maximum_distance != NULL} +\item [\strong{\code{maximum_distance}}] (\emph{delaunay, kNN}) edge maximum euclidean +distance allowed +\item [\strong{\code{Y}}] (\emph{RTriangle}) If TRUE prohibits the insertion of Steiner points +on the mesh boundary. Default is TRUE +\item [\strong{\code{j}}] (\emph{RTriangle}) If TRUE jettisons vertices that are not part of +the final triangulation from the output. Default is TRUE +\item [\strong{\code{S}}] (\emph{RTriangle}) Specifies the maximum number of added Steiner +points. Default is 0 +\item [\strong{\code{options}}] (\emph{geometry}) default is "Pp". See \link[geometry:delaunayn]{geometry::delaunayn} +\item [\strong{\code{weight_fun}}] function to calculate weights based on distance if +\code{include_weight = TRUE}. Default is \eqn{weight = 1 / (1 + distance)} for +\code{"kNN"} and \code{"sNN"} types and \eqn{weight = 1 / distance} for \code{delaunay} +type +networks } } \examples{ -e <- GiottoData::loadSubObjectMini("exprObj")[] # dgCMatrix -expr_igraph <- createNetwork(t(e), node_ids = colnames(e)) - +\dontrun{ pca <- GiottoData::loadSubObjectMini("dimObj")[] -pca_igraph <- createNetwork(pca, node_ids = rownames(pca)) +sl <- GiottoData::loadSubObjectMini("spatLocsObj")[] + +# Delaunay via geometry::delaunayn() +del_geom <- createNetwork( + x = as.matrix(sl[, .(sdimx, sdimy)]), + type = "delaunay", + method = "geometry", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID +) + +# Delaunay via RTriangle::triangulate() +del_rt <- createNetwork( + x = as.matrix(sl[, .(sdimx, sdimy)]), + type = "delaunay", + method = "RTriangle", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID +) + +# Delaunay via deldir::deldir() +del_dd <- createNetwork( + x = as.matrix(sl[, .(sdimx, sdimy)]), + type = "delaunay", + method = "deldir", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID +) + +# kNN spatial network +kNN_spat <- createNetwork( + x = as.matrix(sl[, .(sdimx, sdimy)]), + type = "kNN", + method = "dbscan", + include_weight = TRUE, + weight_fun = function(d) 1 / d, # not the default + as.igraph = FALSE, + node_ids = sl$cell_ID, + k = 4L, + maximum_distance = NULL, + minimum_k = 0L +) + +# kNN NN network +kNN <- createNetwork( + pca[, 1:10], + type = "kNN", + method = "dbscan", + node_ids = rownames(pca), + as.igraph = TRUE +) + +# sNN NN network +sNN <- createNetwork( + pca[, 1:10], + type = "sNN", + method = "dbscan", + node_ids = rownames(pca), + as.igraph = TRUE +) + +# using defaults for sNN with index IDs to create igraph +sNN_idx <- createNetwork(pca[, 1:10]) +} } diff --git a/tests/testthat/test-networks.R b/tests/testthat/test-networks.R new file mode 100644 index 00000000..290caad1 --- /dev/null +++ b/tests/testthat/test-networks.R @@ -0,0 +1,197 @@ +# Ignore internal usage of deprecated accessors +lifecycle_opt <- getOption("lifecycle_verbosity") +options("lifecycle_verbosity" = "quiet") + +# ignore conda +options("giotto.has_conda" = FALSE) + + +# load data to test +g <- GiottoData::loadGiottoMini("viz") +sn <- GiottoData::loadSubObjectMini("spatialNetworkObj") + + +test_that("full spatial network can be created and reverted", { + # reduced to full + n_edges <- nrow(sn) + full <- convert_to_full_spatial_network(sn[]) + expect_equal(n_edges * 2, nrow(full)) + + # revert from full to reduced + reduced <- convert_to_reduced_spatial_network(full) + expect_equal(n_edges, nrow(reduced)) +}) + +test_that("spatial weight matrix can be created", { + rlang::local_options(lifecycle_verbosity = "quiet") + test <- createSpatialWeightMatrix(g, spat_unit = "aggregate", return_gobject = TRUE) + mat <- getSpatialNetwork(test, spat_unit = "aggregate", name = "kNN_network")@misc$weight_matrix$spat_weights + + expect_true(inherits(mat, c("matrix", "Matrix"))) +}) + + +# network generation #### + +e <- GiottoData::loadSubObjectMini("exprObj") +pca <- GiottoData::loadSubObjectMini("dimObj") +sl <- GiottoData::loadSubObjectMini("spatLocsObj") + +# make gobject +options("giotto.use_conda" = FALSE) +# ignore warnings & messages from no py usage + obj setup +g <- suppressMessages(suppressWarnings({ + giotto() %>% + setGiotto(e) %>% + setGiotto(pca) %>% + setGiotto(sl) +})) +options("giotto.use_conda" = TRUE) + +# createNetwork is current separately implemented from Giotto's +# gobject-based network generation functions. These tests currently check +# for parity with the gobject-based functions which have been battle-tested. +# +# TODO update tests after switch to the createNetwork methods for giotto +# functions + +test_that("delaunay [geometry] produces expected results", { + # gobject + g2 <- createSpatialNetwork(g, delaunay_method = "delaunayn_geometry") + del_geom <- getSpatialNetwork(g2)[] + + del_geom2 <- createNetwork( + x = as.matrix(sl[][, .(sdimx, sdimy)]), + type = "delaunay", + method = "geometry", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID + ) + + expect_identical(del_geom$from, del_geom2$from) + expect_identical(del_geom$to, del_geom2$to) + expect_identical(del_geom$weight, del_geom2$weight) + expect_identical(del_geom$distance, del_geom2$distance) +}) + +test_that("delaunay [RTriangle] produces expected results", { + # gobject + g2 <- createSpatialNetwork(g, delaunay_method = "RTriangle") + del_rt <- getSpatialNetwork(g2)[] + + del_rt2 <- createNetwork( + x = as.matrix(sl[][, .(sdimx, sdimy)]), + type = "delaunay", + method = "RTriangle", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID + ) + + expect_identical(del_rt$from, del_rt2$from) + expect_identical(del_rt$to, del_rt2$to) + expect_identical(del_rt$weight, del_rt2$weight) + expect_identical(del_rt$distance, del_rt2$distance) +}) + +test_that("delaunay [deldir] produces expected results", { + # gobject + g2 <- createSpatialNetwork(g, delaunay_method = "deldir") + del_dd <- getSpatialNetwork(g2)[] + + del_dd2 <- createNetwork( + x = as.matrix(sl[][, .(sdimx, sdimy)]), + type = "delaunay", + method = "deldir", + include_weight = TRUE, + weight_fun = function(d) 1 / d, + as.igraph = FALSE, + node_ids = sl$cell_ID + ) + + expect_identical(del_dd$from, del_dd2$from) + expect_identical(del_dd$to, del_dd2$to) + + # deldir original implementation generates distances from x and y values + # that are reported by the deldir object. The new createNetwork() + # implementation uses the spatiallocations like the other delaunay methods do. + # A very slight difference less than 1e-4 is expected + cutoff <- 1e-4 + expect_true(all((del_dd$weight - del_dd2$weight) < cutoff)) + expect_true(all((del_dd$distance - del_dd2$distance) < cutoff)) +}) + +test_that("knn [dbscan] produces expected results", { + # gobject + g2 <- createSpatialNetwork( + g, + method = "kNN", + knn_method = "dbscan", + k = 8L, + maximum_distance_knn = NULL, # (900 by itself) + minimum_k = 0L # no change by itself, 1488 with max_dist_knn + ) + kNN_spat <- getSpatialNetwork(g2)[] + + kNN_spat2 <- createNetwork( + x = as.matrix(sl[][, .(sdimx, sdimy)]), + type = "kNN", + method = "dbscan", + include_weight = TRUE, + weight_fun = function(d) 1 / d, # not the default + as.igraph = FALSE, + node_ids = sl$cell_ID, + k = 8L, + filter = TRUE, + maximum_distance = NULL, + minimum_k = 0L + ) + + data.table::setorder(kNN_spat, from, to) + data.table::setorder(kNN_spat2, from, to) + + expect_identical(kNN_spat$from, kNN_spat2$from) + expect_identical(kNN_spat$to, kNN_spat2$to) + # expect_identical(kNN_spat$weight, kNN_spat2$weight) # expected to be different + expect_identical(kNN_spat$distance, kNN_spat2$distance) +}) + +test_that("kNN [dbscan] produces expected results", { + # gobject + g2 <- createNearestNetwork(g, type = "kNN", dimensions_to_use = 1:10) + kNN_g <- suppressMessages(getNearestNetwork(g2)[]) + + # low level + kNN_cn <- createNetwork( + pca[][, 1:10], + type = "kNN", + method = "dbscan", + node_ids = rownames(pca), + as.igraph = TRUE, + ) + + expect_true(igraph::identical_graphs(kNN_g, kNN_cn)) +}) + +test_that("sNN [dbscan] produces expected results", { + # gobject + g2 <- createNearestNetwork(g, type = "sNN", dimensions_to_use = 1:10) + sNN_g <- suppressMessages(getNearestNetwork(g2)[]) + + # low level + sNN_cn <- createNetwork( + pca[][, 1:10], + type = "sNN", + method = "dbscan", + node_ids = rownames(pca), + as.igraph = TRUE, + ) + + expect_true(igraph::identical_graphs(sNN_g, sNN_cn)) +}) + + + diff --git a/tests/testthat/test-spatial_structures.R b/tests/testthat/test-spatial_structures.R deleted file mode 100644 index 7a934f80..00000000 --- a/tests/testthat/test-spatial_structures.R +++ /dev/null @@ -1,31 +0,0 @@ -# Ignore internal usage of deprecated accessors -lifecycle_opt <- getOption("lifecycle_verbosity") -options("lifecycle_verbosity" = "quiet") - -# ignore conda -options("giotto.has_conda" = FALSE) - - -# load data to test -g <- GiottoData::loadGiottoMini("viz") -sn <- GiottoData::loadSubObjectMini("spatialNetworkObj") - - -test_that("full spatial network can be created and reverted", { - # reduced to full - n_edges <- nrow(sn) - full <- convert_to_full_spatial_network(sn[]) - expect_equal(n_edges * 2, nrow(full)) - - # revert from full to reduced - reduced <- convert_to_reduced_spatial_network(full) - expect_equal(n_edges, nrow(reduced)) -}) - -test_that("spatial weight matrix can be created", { - rlang::local_options(lifecycle_verbosity = "quiet") - test <- createSpatialWeightMatrix(g, spat_unit = "aggregate", return_gobject = TRUE) - mat <- getSpatialNetwork(test, spat_unit = "aggregate", name = "kNN_network")@misc$weight_matrix$spat_weights - - expect_true(inherits(mat, c("matrix", "Matrix"))) -}) From 3fc568f31d8ac6a982897152d6a3b9422951cadc Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 16:35:20 -0500 Subject: [PATCH 11/24] chore: run document --- NAMESPACE | 3 --- man/aggregateStacksExpression.Rd | 4 ++-- man/aggregateStacksLocations.Rd | 4 ++-- man/aggregateStacksPolygonOverlaps.Rd | 4 ++-- man/aggregateStacksPolygons.Rd | 4 ++-- man/seuratToGiottoV4.Rd | 5 ++++- man/seuratToGiottoV5.Rd | 5 ++++- 7 files changed, 16 insertions(+), 13 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 0c8bbbc4..0d996a9a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -367,9 +367,6 @@ importClassesFrom(Matrix,Matrix) importClassesFrom(SeuratObject,SpatialImage) importClassesFrom(terra,SpatExtent) importClassesFrom(terra,SpatVector) -importFrom(GiottoUtils,getDistinctColors) -importFrom(GiottoUtils,getMonochromeColors) -importFrom(GiottoUtils,getRainbowColors) importFrom(grDevices,dev.size) importFrom(graphics,legend) importFrom(graphics,par) diff --git a/man/aggregateStacksExpression.Rd b/man/aggregateStacksExpression.Rd index 94091fd0..e32de879 100644 --- a/man/aggregateStacksExpression.Rd +++ b/man/aggregateStacksExpression.Rd @@ -37,9 +37,9 @@ aggregate expression matrices from different z-stacks } \seealso{ Other aggregate stacks: -\code{\link{aggregateStacks}()}, \code{\link{aggregateStacksLocations}()}, \code{\link{aggregateStacksPolygonOverlaps}()}, -\code{\link{aggregateStacksPolygons}()} +\code{\link{aggregateStacksPolygons}()}, +\code{\link{aggregateStacks}()} } \concept{aggregate stacks} diff --git a/man/aggregateStacksLocations.Rd b/man/aggregateStacksLocations.Rd index ef667912..6bd4fbe9 100644 --- a/man/aggregateStacksLocations.Rd +++ b/man/aggregateStacksLocations.Rd @@ -31,9 +31,9 @@ aggregate expression matrices from different z-stacks } \seealso{ Other aggregate stacks: -\code{\link{aggregateStacks}()}, \code{\link{aggregateStacksExpression}()}, \code{\link{aggregateStacksPolygonOverlaps}()}, -\code{\link{aggregateStacksPolygons}()} +\code{\link{aggregateStacksPolygons}()}, +\code{\link{aggregateStacks}()} } \concept{aggregate stacks} diff --git a/man/aggregateStacksPolygonOverlaps.Rd b/man/aggregateStacksPolygonOverlaps.Rd index 6e282ef4..032dceac 100644 --- a/man/aggregateStacksPolygonOverlaps.Rd +++ b/man/aggregateStacksPolygonOverlaps.Rd @@ -28,9 +28,9 @@ aggregate polygons overlap information from different z-stacks } \seealso{ Other aggregate stacks: -\code{\link{aggregateStacks}()}, \code{\link{aggregateStacksExpression}()}, \code{\link{aggregateStacksLocations}()}, -\code{\link{aggregateStacksPolygons}()} +\code{\link{aggregateStacksPolygons}()}, +\code{\link{aggregateStacks}()} } \concept{aggregate stacks} diff --git a/man/aggregateStacksPolygons.Rd b/man/aggregateStacksPolygons.Rd index 2cac98eb..0077887a 100644 --- a/man/aggregateStacksPolygons.Rd +++ b/man/aggregateStacksPolygons.Rd @@ -31,9 +31,9 @@ aggregate polygons from different z-stacks } \seealso{ Other aggregate stacks: -\code{\link{aggregateStacks}()}, \code{\link{aggregateStacksExpression}()}, \code{\link{aggregateStacksLocations}()}, -\code{\link{aggregateStacksPolygonOverlaps}()} +\code{\link{aggregateStacksPolygonOverlaps}()}, +\code{\link{aggregateStacks}()} } \concept{aggregate stacks} diff --git a/man/seuratToGiottoV4.Rd b/man/seuratToGiottoV4.Rd index f5ae5e81..8600a295 100644 --- a/man/seuratToGiottoV4.Rd +++ b/man/seuratToGiottoV4.Rd @@ -8,7 +8,10 @@ seuratToGiottoV4( sobject, spatial_assay = "Spatial", dim_reduction = c("pca", "umap"), - subcellular_assay = "Vizgen" + subcellular_assay = "Vizgen", + sp_network = NULL, + nn_network = NULL, + verbose = TRUE ) } \arguments{ diff --git a/man/seuratToGiottoV5.Rd b/man/seuratToGiottoV5.Rd index d11670a2..7f0a416e 100644 --- a/man/seuratToGiottoV5.Rd +++ b/man/seuratToGiottoV5.Rd @@ -8,7 +8,10 @@ seuratToGiottoV5( sobject, spatial_assay = "Spatial", dim_reduction = c("pca", "umap"), - subcellular_assay = "Vizgen" + subcellular_assay = "Vizgen", + sp_network = NULL, + nn_network = NULL, + verbose = TRUE ) } \arguments{ From 7481276d584d61eb097d8303199b484f46c85c68 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 18:48:01 -0500 Subject: [PATCH 12/24] chore: change `createNetwork()` to non-generic function Avoids needing to import matrix class definitions. --- NAMESPACE | 4 +- R/NN_network.R | 118 +++++++++++++++++++++---------------------- R/classes.R | 7 --- R/generics.R | 2 - R/package_imports.R | 2 - man/createNetwork.Rd | 3 +- 6 files changed, 59 insertions(+), 77 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 0d996a9a..018df106 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -77,6 +77,7 @@ export(createGiottoPolygonsFromMask) export(createMetafeats) export(createNearestNetObj) export(createNearestNetwork) +export(createNetwork) export(createSpatEnrObj) export(createSpatLocsObj) export(createSpatNetObj) @@ -332,7 +333,6 @@ exportMethods(colnames) exportMethods(copy) exportMethods(createGiottoPoints) exportMethods(createGiottoPolygon) -exportMethods(createNetwork) exportMethods(crop) exportMethods(dim) exportMethods(dimnames) @@ -362,8 +362,6 @@ exportMethods(wrap) import(GiottoUtils) import(data.table) import(utils) -importClassesFrom(DelayedArray,DelayedArray) -importClassesFrom(Matrix,Matrix) importClassesFrom(SeuratObject,SpatialImage) importClassesFrom(terra,SpatExtent) importClassesFrom(terra,SpatVector) diff --git a/R/NN_network.R b/R/NN_network.R index bd6996bb..d0fc76b8 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -128,10 +128,7 @@ NULL #' @rdname createNetwork #' @export -setMethod( - "createNetwork", - signature("allMatrix"), - function( +createNetwork <- function( x, type = c("sNN", "kNN", "delaunay"), method = c("dbscan", "geometry", "RTriangle", "deldir"), @@ -140,74 +137,73 @@ setMethod( include_weight = TRUE, as.igraph = TRUE, ... - ) - { - - # check params - type <- match.arg(type, choices = c("sNN", "kNN", "delaunay")) - method <- switch( - type, - "sNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), - "kNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), - "delaunay" = match.arg( - method, choices = c("geometry", "RTriangle", "deldir"), - several.ok = TRUE - ) - ) +) +{ - # get common params - alist <- list( - x = x, - include_weight = include_weight, - ... + # check params + type <- match.arg(type, choices = c("sNN", "kNN", "delaunay")) + method <- switch( + type, + "sNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), + "kNN" = match.arg(method, choices = c("dbscan"), several.ok = TRUE), + "delaunay" = match.arg( + method, choices = c("geometry", "RTriangle", "deldir"), + several.ok = TRUE ) + ) - # generate network data.table - network_dt <- switch( - sprintf("%s:%s", type, method), - "kNN:dbscan" = do.call(.net_dt_knn, args = alist), - "sNN:dbscan" = do.call(.net_dt_snn, args = alist), - "delaunay:deldir" = do.call(.net_dt_del_deldir, - args = alist)$delaunay_network_DT, - "delaunay:RTriangle" = do.call(.net_dt_del_rtriangle, - args = alist)$delaunay_network_DT, - "delaunay:geometry" = do.call(.net_dt_del_geometry, - args = alist)$delaunay_network_DT - ) + # get common params + alist <- list( + x = x, + include_weight = include_weight, + ... + ) - # replace indices with node_ids if desired - if (!is.null(node_ids)) { - # if node_ids are provided, include as extra cols - names(node_ids) <- seq_along(node_ids) - network_dt[, "from" := node_ids[from]] - network_dt[, "to" := node_ids[to]] - } + # generate network data.table + network_dt <- switch( + sprintf("%s:%s", type, method), + "kNN:dbscan" = do.call(.net_dt_knn, args = alist), + "sNN:dbscan" = do.call(.net_dt_snn, args = alist), + "delaunay:deldir" = do.call(.net_dt_del_deldir, + args = alist)$delaunay_network_DT, + "delaunay:RTriangle" = do.call(.net_dt_del_rtriangle, + args = alist)$delaunay_network_DT, + "delaunay:geometry" = do.call(.net_dt_del_geometry, + args = alist)$delaunay_network_DT + ) - ## outputs ## + # replace indices with node_ids if desired + if (!is.null(node_ids)) { + # if node_ids are provided, include as extra cols + names(node_ids) <- seq_along(node_ids) + network_dt[, "from" := node_ids[from]] + network_dt[, "to" := node_ids[to]] + } - # cols to include in output - keep_cols <- c("from", "to") + ## outputs ## - all_index <- network_dt[, unique(unlist(.SD)), .SDcols = keep_cols] - if (include_weight) keep_cols <- c(keep_cols, "weight") - if (include_distance) keep_cols <- c(keep_cols, "distance") + # cols to include in output + keep_cols <- c("from", "to") - if (type == "sNN") keep_cols <- c(keep_cols, "shared", "rank") + all_index <- network_dt[, unique(unlist(.SD)), .SDcols = keep_cols] + if (include_weight) keep_cols <- c(keep_cols, "weight") + if (include_distance) keep_cols <- c(keep_cols, "distance") - # return early if igraph not required - network_dt_final <- network_dt[, keep_cols, with = FALSE] - if (!as.igraph) return(network_dt_final) + if (type == "sNN") keep_cols <- c(keep_cols, "shared", "rank") - ## convert to igraph object - out <- igraph::graph_from_data_frame( - network_dt_final, - directed = TRUE, - vertices = all_index - ) + # return early if igraph not required + network_dt_final <- network_dt[, keep_cols, with = FALSE] + if (!as.igraph) return(network_dt_final) - return(out) - } -) + ## convert to igraph object + out <- igraph::graph_from_data_frame( + network_dt_final, + directed = TRUE, + vertices = all_index + ) + + return(out) +} # x input is a matrix diff --git a/R/classes.R b/R/classes.R index a4a1775d..afb04f55 100644 --- a/R/classes.R +++ b/R/classes.R @@ -30,13 +30,6 @@ setClassUnion("nullOrDatatable", c("NULL", "data.table")) #' @noRd setClassUnion("gIndex", c("numeric", "logical", "character")) -#' @title Matrix classes -#' @description combined class definition for matrix representations for use in -#' dispatching. -#' @keywords internal -#' @noRd -setClassUnion("allMatrix", c("matrix", "Matrix", "DelayedArray")) - diff --git a/R/generics.R b/R/generics.R index 0364b79d..8c807466 100644 --- a/R/generics.R +++ b/R/generics.R @@ -52,8 +52,6 @@ setGeneric("copy", function(x) standardGeneric("copy"), useAsDefault = data.tabl setGeneric("calculateOverlap", function(x, y, ...) standardGeneric("calculateOverlap")) setGeneric("overlapToMatrix", function(x, ...) standardGeneric("overlapToMatrix")) -# networks -setGeneric("createNetwork", function(x, ...) standardGeneric("createNetwork")) # Giotto subnesting #### # All methods and documentations found in methods-nesting.R diff --git a/R/package_imports.R b/R/package_imports.R index 2214ca23..99653da7 100644 --- a/R/package_imports.R +++ b/R/package_imports.R @@ -27,8 +27,6 @@ #' @importMethodsFrom terra nrow ncol #' @importClassesFrom terra SpatExtent #' @importClassesFrom terra SpatVector -#' @importClassesFrom DelayedArray DelayedArray -#' @importClassesFrom Matrix Matrix #' @import GiottoUtils #' @import data.table #' @import utils diff --git a/man/createNetwork.Rd b/man/createNetwork.Rd index cdf5251c..63e6a187 100644 --- a/man/createNetwork.Rd +++ b/man/createNetwork.Rd @@ -2,10 +2,9 @@ % Please edit documentation in R/NN_network.R \name{createNetwork} \alias{createNetwork} -\alias{createNetwork,allMatrix-method} \title{Create a network} \usage{ -\S4method{createNetwork}{allMatrix}( +createNetwork( x, type = c("sNN", "kNN", "delaunay"), method = c("dbscan", "geometry", "RTriangle", "deldir"), From 38293d23450b6f3cf7303186687ec82b61f140dc Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 22:59:55 -0500 Subject: [PATCH 13/24] chore: fix network test script --- .github/workflows/covr.yml | 2 +- .github/workflows/dev_check.yml | 2 +- .github/workflows/main_check.yml | 2 +- tests/testthat/test-networks.R | 12 ++++++++++++ 4 files changed, 15 insertions(+), 3 deletions(-) diff --git a/.github/workflows/covr.yml b/.github/workflows/covr.yml index a5ae4864..c1baec4c 100644 --- a/.github/workflows/covr.yml +++ b/.github/workflows/covr.yml @@ -35,7 +35,7 @@ jobs: _R_CHECK_RD_XREFS: false with: dependencies: '"hard"' # do not use suggested dependencies - extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::remotes, any::covr, any::sp, any::stars, any::raster, any::sf + extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::remotes, any::covr, any::sp, any::stars, any::raster, any::sf, any::RTriangle, any::geometry needs: coverage - name: Set up dependencies (GiottoData) diff --git a/.github/workflows/dev_check.yml b/.github/workflows/dev_check.yml index 98e00617..24d4d08d 100644 --- a/.github/workflows/dev_check.yml +++ b/.github/workflows/dev_check.yml @@ -46,7 +46,7 @@ jobs: _R_CHECK_RD_XREFS: false with: dependencies: '"hard"' # do not use suggested dependencies - extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::knitr, any::rmarkdown, any::sp, any::stars, any::raster, any::sf, any::scattermore, any::exactextractr, github::drieslab/GiottoData + extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::knitr, any::rmarkdown, any::sp, any::stars, any::raster, any::sf, any::scattermore, any::exactextractr, any::RTriangle, any::geometry, github::drieslab/GiottoData - name: Run R CMD check uses: r-lib/actions/check-r-package@v2 diff --git a/.github/workflows/main_check.yml b/.github/workflows/main_check.yml index 6545b3fb..fdcc921e 100644 --- a/.github/workflows/main_check.yml +++ b/.github/workflows/main_check.yml @@ -62,7 +62,7 @@ jobs: _R_CHECK_RD_XREFS: false with: dependencies: '"hard"' # do not use suggested dependencies - extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::remotes, any::knitr, any::rmarkdown, any::sp, any::stars, any::raster, any::sf, any::scattermore, any::exactextractr, github::drieslab/GiottoData + extra-packages: any::rcmdcheck, any::testthat, any::rlang, any::R.utils, any::remotes, any::knitr, any::rmarkdown, any::sp, any::stars, any::raster, any::sf, any::scattermore, any::exactextractr, any::RTriangle, any::geometry, github::drieslab/GiottoData - name: Test python env build run: | diff --git a/tests/testthat/test-networks.R b/tests/testthat/test-networks.R index 290caad1..d2bd2108 100644 --- a/tests/testthat/test-networks.R +++ b/tests/testthat/test-networks.R @@ -56,6 +56,8 @@ options("giotto.use_conda" = TRUE) # functions test_that("delaunay [geometry] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createSpatialNetwork(g, delaunay_method = "delaunayn_geometry") del_geom <- getSpatialNetwork(g2)[] @@ -77,6 +79,8 @@ test_that("delaunay [geometry] produces expected results", { }) test_that("delaunay [RTriangle] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createSpatialNetwork(g, delaunay_method = "RTriangle") del_rt <- getSpatialNetwork(g2)[] @@ -98,6 +102,8 @@ test_that("delaunay [RTriangle] produces expected results", { }) test_that("delaunay [deldir] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createSpatialNetwork(g, delaunay_method = "deldir") del_dd <- getSpatialNetwork(g2)[] @@ -125,6 +131,8 @@ test_that("delaunay [deldir] produces expected results", { }) test_that("knn [dbscan] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createSpatialNetwork( g, @@ -160,6 +168,8 @@ test_that("knn [dbscan] produces expected results", { }) test_that("kNN [dbscan] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createNearestNetwork(g, type = "kNN", dimensions_to_use = 1:10) kNN_g <- suppressMessages(getNearestNetwork(g2)[]) @@ -177,6 +187,8 @@ test_that("kNN [dbscan] produces expected results", { }) test_that("sNN [dbscan] produces expected results", { + rlang::local_options(lifecycle_verbosity = "quiet") + # gobject g2 <- createNearestNetwork(g, type = "sNN", dimensions_to_use = 1:10) sNN_g <- suppressMessages(getNearestNetwork(g2)[]) From 7bc7afce632235b95dc0295b965ea5af27f48cfc Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Tue, 13 Feb 2024 23:01:16 -0500 Subject: [PATCH 14/24] chore: change deprecated igraph function names --- R/data_evaluation.R | 4 +- R/interoperability.R | 66 +++++++++++++++--------------- tests/testthat/test-createObject.R | 10 ++--- 3 files changed, 40 insertions(+), 40 deletions(-) diff --git a/R/data_evaluation.R b/R/data_evaluation.R index 448f6286..2be1fa6b 100644 --- a/R/data_evaluation.R +++ b/R/data_evaluation.R @@ -511,8 +511,8 @@ evaluate_input <- function(type, x, ...) { nn_network[] <- .evaluate_nearest_networks(nn_network = nn_network[]) return(nn_network) } else if (inherits(nn_network, "igraph")) { - v_attr <- igraph::list.vertex.attributes(nn_network) - e_attr <- igraph::list.edge.attributes(nn_network) + v_attr <- igraph::vertex_attr_names(nn_network) + e_attr <- igraph::edge_attr_names(nn_network) if (!"name" %in% v_attr) { .gstop( diff --git a/R/interoperability.R b/R/interoperability.R index 995c50a5..6dc47f1d 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -1497,7 +1497,7 @@ giottoToSeuratV5 <- function(gobject, for (i in seq(sobj@assays)) { sobj@assays[[i]]@meta.data <- meta_genes } - + # dim reduction # note: Seurat requires assay name specification for each dim reduc @@ -1619,10 +1619,10 @@ giottoToSeuratV5 <- function(gobject, } } } - + all_x <- NULL all_y <- NULL - + if(length(gobject@largeImages) > 0){ for (i in seq_along(gobject@largeImages)){ # spatVec <- terra::as.points(gobject@largeImages[[i]]@raster_object) @@ -1634,25 +1634,25 @@ giottoToSeuratV5 <- function(gobject, img <- terra::as.array(gobject@largeImages[[i]]@raster_object) img[, , 1:3] <- img[, , 1:3] / 255 coord <- data.frame(imagerow = imagerow, imagecol = imagecol) - + scalefactors <- Seurat::scalefactors(spot = gobject@largeImages[[i]]@scale_factor, fiducial = gobject@largeImages[[i]]@resolution, hires = max(img), lowres = min(img)) - - + + newV1 <- new(Class = "VisiumV1", image = img, scale.factors = scalefactors, coordinates = coord) - + sobj@images[[gobject@largeImages[[i]]@name]] <- newV1 - - + + } - + } - + return(sobj) @@ -1784,25 +1784,25 @@ seuratToGiottoV4 <- function(sobject, } } } - + #Networks - + #spatial_network - + if (!is.null(sp_network)) { for (i in seq(sp_network)) { if (verbose) message("Copying spatial networks") stempgraph <- Seurat::as.Graph(sobject@graphs[[sp_network[i]]]) m <- as.matrix(stempgraph) sobjIgraph <- igraph::graph.adjacency(m, weighted = TRUE) - + if (verbose) message("Calculating distances") for (j in seq_along(sobject@graphs[[sp_network[i]]]@x)){ dist_matrix <- sobject@graphs[[sp_network[i]]]@x[[j]] sobjIgraph <- igraph::set_edge_attr(graph = sobjIgraph, name = "distance", value = dist_matrix, index = j) } - e_attr <- igraph::list.edge.attributes(sobjIgraph) + e_attr <- igraph::edge_attr_names(sobjIgraph) if ("weight" %in% e_attr) { igDT <- data.table::setDT(igraph::as_data_frame(sobjIgraph)) igDT[, weight := 1 / (1 + distance)] @@ -1810,7 +1810,7 @@ seuratToGiottoV4 <- function(sobject, sobjIgraph <- igraph::graph_from_data_frame(igDT) } DT <- data.table() - + edges <- igraph::ends(sobjIgraph, es = E(sobjIgraph), names = TRUE) DT$from <- edges[, 1] DT$to <- edges[, 2] @@ -1833,28 +1833,28 @@ seuratToGiottoV4 <- function(sobject, stempgraph <- Seurat::as.Graph(sobject@graphs[[nn_network[i]]]) m <- as.matrix(stempgraph) sobjIgraph <- igraph::graph.adjacency(m, weighted = TRUE) - + if (verbose) message("Calculating distances") for (j in seq_along(sobject@graphs[[nn_network[i]]]@x)){ dist_matrix <- sobject@graphs[[nn_network[i]]]@x[[j]] - + sobjIgraph <- igraph::set_edge_attr(graph = sobjIgraph, name = "distance", value = dist_matrix, index = j) } - e_attr <- igraph::list.edge.attributes(sobjIgraph) + e_attr <- igraph::edge_attr_names(sobjIgraph) if ("weight" %in% e_attr) { igDT <- data.table::setDT(igraph::as_data_frame(sobjIgraph)) igDT[, weight := 1 / (1 + distance)] data.table::setcolorder(igDT, neworder = c("from", "to", "weight", "distance")) sobjIgraph <- igraph::graph_from_data_frame(igDT) } - + if (verbose) message("Copying nearest neighbour networks") nnNetObj <- GiottoClass::createNearestNetObj(name = names(sobject@graphs)[i], network = sobjIgraph) return(nnNetObj) }) - + for (i in seq_along(nnNetObj_list)) { gobject <- GiottoClass::set_NearestNetwork(gobject = gobject, nn_network = nnNetObj_list[[i]]) @@ -2063,7 +2063,7 @@ seuratToGiottoV5 <- function(sobject, r <- as.matrix(sobject@images[[i]]@image[, , 1]) g <- as.matrix(sobject@images[[i]]@image[, , 2]) b <- as.matrix(sobject@images[[i]]@image[, , 3]) - + r <- round(r * 255) g <- round(g * 255) b <- round(b * 255) @@ -2095,25 +2095,25 @@ seuratToGiottoV5 <- function(sobject, gobject <- set_expression_values(gobject = gobject, values = exprObj, set_defaults = FALSE) # gobject@expression$cell$rna$normalized = normexp } - + #Networks - + #spatial_network - + if (!is.null(sp_network)) { for (i in seq(sp_network)) { if (verbose) message("Copying spatial networks") stempgraph <- Seurat::as.Graph(sobject@graphs[[sp_network[i]]]) m <- as.matrix(stempgraph) sobjIgraph <- igraph::graph.adjacency(m, weighted = TRUE) - + if (verbose) message("Calculating distances") for (j in seq_along(sobject@graphs[[sp_network[i]]]@x)){ dist_matrix <- sobject@graphs[[sp_network[i]]]@x[[j]] sobjIgraph <- igraph::set_edge_attr(graph = sobjIgraph, name = "distance", value = dist_matrix, index = j) } - e_attr <- igraph::list.edge.attributes(sobjIgraph) + e_attr <- igraph::edge_attr_names(sobjIgraph) if ("weight" %in% e_attr) { igDT <- data.table::setDT(igraph::as_data_frame(sobjIgraph)) igDT[, weight := 1 / (1 + distance)] @@ -2121,7 +2121,7 @@ seuratToGiottoV5 <- function(sobject, sobjIgraph <- igraph::graph_from_data_frame(igDT) } DT <- data.table() - + edges <- igraph::ends(sobjIgraph, es = E(sobjIgraph), names = TRUE) DT$from <- edges[, 1] DT$to <- edges[, 2] @@ -2145,15 +2145,15 @@ seuratToGiottoV5 <- function(sobject, stempgraph <- Seurat::as.Graph(sobject@graphs[[nn_network[i]]]) m <- as.matrix(stempgraph) sobjIgraph <- igraph::graph.adjacency(m, weighted = TRUE) - + if (verbose) message("Calculating distances") for (j in seq_along(sobject@graphs[[nn_network[i]]]@x)){ dist_matrix <- sobject@graphs[[nn_network[i]]]@x[[j]] - + sobjIgraph <- igraph::set_edge_attr(graph = sobjIgraph, name = "distance", value = dist_matrix, index = j) } - e_attr <- igraph::list.edge.attributes(sobjIgraph) + e_attr <- igraph::edge_attr_names(sobjIgraph) if ("weight" %in% e_attr) { igDT <- data.table::setDT(igraph::as_data_frame(sobjIgraph)) igDT[, weight := 1 / (1 + distance)] @@ -2166,7 +2166,7 @@ seuratToGiottoV5 <- function(sobject, network = sobjIgraph) return(nnNetObj) }) - + for (i in seq_along(nnNetObj_list)) { gobject <- GiottoClass::set_NearestNetwork(gobject = gobject, nn_network = nnNetObj_list[[i]]) diff --git a/tests/testthat/test-createObject.R b/tests/testthat/test-createObject.R index 3b60e98d..f6a59433 100644 --- a/tests/testthat/test-createObject.R +++ b/tests/testthat/test-createObject.R @@ -142,8 +142,8 @@ test_that("Eval with missing info throws error", { test_that("Eval with minimum info works", { ig_new <- expect_no_error(.evaluate_nearest_networks(nnDT_min)) expect_s3_class(ig_new, "igraph") - expect_true(all(c("distance", "weight") %in% igraph::list.edge.attributes(ig_new))) - expect_true("name" %in% igraph::list.vertex.attributes(ig_new)) + expect_true(all(c("distance", "weight") %in% igraph::edge_attr_names(ig_new))) + expect_true("name" %in% igraph::vertex_attr_names(ig_new)) }) @@ -172,9 +172,9 @@ test_that("Eval of minimal igraph adds weight attr", { ig_min <- igraph::delete_edge_attr(ig, "weight") ig_min <- .evaluate_nearest_networks(ig_min) expect_s3_class(ig_min, "igraph") - expect_true(all(c("distance", "weight") %in% igraph::list.edge.attributes(ig_min))) - expect_true("name" %in% igraph::list.vertex.attributes(ig_min)) - expect_true("weight" %in% igraph::list.edge.attributes(ig_min)) + expect_true(all(c("distance", "weight") %in% igraph::edge_attr_names(ig_min))) + expect_true("name" %in% igraph::vertex_attr_names(ig_min)) + expect_true("weight" %in% igraph::edge_attr_names(ig_min)) }) From ccaa1ac76bb833476789dd537fd0e2729f773e6c Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 01:14:16 -0500 Subject: [PATCH 15/24] chore: fixes for rcmdcheck --- R/NN_network.R | 19 ++++++++++++------- R/interoperability.R | 4 ++-- man/createNetwork.Rd | 3 +++ 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/R/NN_network.R b/R/NN_network.R index d0fc76b8..220819c6 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -22,6 +22,7 @@ NULL #' @param include_distance logical. include edge distance attribute in output #' @param as.igraph logical. Whether to return as `igraph`. Otherwise returns #' as `data.table` +#' @param verbose be verbose. Default = NULL (uses giotto verbosity option) #' @param ... additional params to pass. See details section #' @returns Either `igraph` if `as.igraph = TRUE` and `data.table` otherwise. #' @details @@ -136,6 +137,7 @@ createNetwork <- function( include_distance = TRUE, include_weight = TRUE, as.igraph = TRUE, + verbose = NULL, ... ) { @@ -210,10 +212,10 @@ createNetwork <- function( .net_dt_knn <- function( x, k = 30L, include_weight = TRUE, include_distance = TRUE, filter = FALSE, maximum_distance = NULL, minimum_k = 0L, weight_fun = function(d) 1 / (1 + d), - ... + verbose, ... ) { # NSE vars - from <- to <- NULL + from <- to <- distance <- NULL k <- as.integer(k) @@ -266,10 +268,10 @@ createNetwork <- function( .net_dt_snn <- function( x, k = 30L, include_weight = TRUE, include_distance = TRUE, top_shared = 3L, minimum_shared = 5L, weight_fun = function(d) 1 / (1 + d), - ... + verbose, ... ) { # NSE vars - from <- to <- shared <- NULL + from <- to <- shared <- distance <- NULL k <- as.integer(k) top_shared <- as.integer(top_shared) @@ -321,7 +323,7 @@ createNetwork <- function( package_check("geometry", repository = "CRAN:geometry") # data.table variables - from <- to <- NULL + from <- to <- distance <- NULL delaunay_simplex_mat <- geometry::delaunayn( p = x, options = options, ... @@ -372,7 +374,7 @@ createNetwork <- function( ... ) { # NSE vars - from <- to <- NULL + from <- to <- distance <- NULL package_check("RTriangle", repository = "CRAN:RTriangle") @@ -417,7 +419,7 @@ createNetwork <- function( ... ) { # NSE variables - from <- to <- NULL + from <- to <- distance <- NULL if (ncol(x) > 2L) { .gstop("\'deldir\' delaunay method only applies to 2D data. @@ -496,6 +498,9 @@ edge_distances <- function(x, y, x_node_ids = NULL) { #' @keywords internal .edge_coords_array <- function(x, y, x_node_ids = NULL) { + # NSE vars + from <- to <- NULL + checkmate::assert_matrix(x) checkmate::assert_data_table(y) diff --git a/R/interoperability.R b/R/interoperability.R index 6dc47f1d..004931d4 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -1811,7 +1811,7 @@ seuratToGiottoV4 <- function(sobject, } DT <- data.table() - edges <- igraph::ends(sobjIgraph, es = E(sobjIgraph), names = TRUE) + edges <- igraph::ends(sobjIgraph, es = igraph::E(sobjIgraph), names = TRUE) DT$from <- edges[, 1] DT$to <- edges[, 2] ed_attr <- igraph::edge.attributes(sobjIgraph) @@ -2122,7 +2122,7 @@ seuratToGiottoV5 <- function(sobject, } DT <- data.table() - edges <- igraph::ends(sobjIgraph, es = E(sobjIgraph), names = TRUE) + edges <- igraph::ends(sobjIgraph, es = igraph::E(sobjIgraph), names = TRUE) DT$from <- edges[, 1] DT$to <- edges[, 2] ed_attr <- igraph::edge.attributes(sobjIgraph) diff --git a/man/createNetwork.Rd b/man/createNetwork.Rd index 63e6a187..f26c2549 100644 --- a/man/createNetwork.Rd +++ b/man/createNetwork.Rd @@ -12,6 +12,7 @@ createNetwork( include_distance = TRUE, include_weight = TRUE, as.igraph = TRUE, + verbose = NULL, ... ) } @@ -35,6 +36,8 @@ will be used as node IDs.} \item{as.igraph}{logical. Whether to return as \code{igraph}. Otherwise returns as \code{data.table}} +\item{verbose}{be verbose. Default = NULL (uses giotto verbosity option)} + \item{...}{additional params to pass. See details section} } \value{ From fdcd664c695eed87fc081951c3921e1b8b0769ad Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 01:16:55 -0500 Subject: [PATCH 16/24] fix: alternative method of creating VisiumV1 Try to avoid import of SeuratObject --- DESCRIPTION | 1 - NAMESPACE | 3 --- R/interoperability.R | 20 ++++++++++++++++---- R/visiumImage.R | 19 ------------------- man/VisiumV1-class.Rd | 26 -------------------------- 5 files changed, 16 insertions(+), 53 deletions(-) delete mode 100644 R/visiumImage.R delete mode 100644 man/VisiumV1-class.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 6aa34dcf..3b832b44 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -137,7 +137,6 @@ Collate: 'stitch_coordinates.R' 'subset.R' 'suite_reexports.R' - 'visiumImage.R' 'zzz.R' URL: https://drieslab.github.io/GiottoClass/ VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 018df106..eb0b3781 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,7 +10,6 @@ S3method(as.data.table,giottoPoints) S3method(as.data.table,giottoPolygon) S3method(t,spatLocsObj) S3method(t,spatialNetworkObj) -export(VisiumV1) export(addCellMetadata) export(addFeatMetadata) export(addGiottoImage) @@ -291,7 +290,6 @@ export(updateGiottoPointsObject) export(updateGiottoPolygonObject) export(update_giotto_params) export(writeGiottoLargeImage) -exportClasses(VisiumV1) exportClasses(cellMetaObj) exportClasses(dimObj) exportClasses(exprObj) @@ -362,7 +360,6 @@ exportMethods(wrap) import(GiottoUtils) import(data.table) import(utils) -importClassesFrom(SeuratObject,SpatialImage) importClassesFrom(terra,SpatExtent) importClassesFrom(terra,SpatVector) importFrom(grDevices,dev.size) diff --git a/R/interoperability.R b/R/interoperability.R index 004931d4..4c74f0ec 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -1640,11 +1640,23 @@ giottoToSeuratV5 <- function(gobject, hires = max(img), lowres = min(img)) + senv <- rlang::ns_env("SeuratObject") + + newV1 <- do.call( + new, + args = list( + Class = "VisiumV1", + image = img, + scale.factors = scalefactors, + coordinates = coord + ), + envir = senv + ) - newV1 <- new(Class = "VisiumV1", - image = img, - scale.factors = scalefactors, - coordinates = coord) + # newV1 <- new(Class = "VisiumV1", + # image = img, + # scale.factors = scalefactors, + # coordinates = coord) sobj@images[[gobject@largeImages[[i]]@name]] <- newV1 diff --git a/R/visiumImage.R b/R/visiumImage.R deleted file mode 100644 index 5ded2cd7..00000000 --- a/R/visiumImage.R +++ /dev/null @@ -1,19 +0,0 @@ -#' Create VisiumV1 Class -#' VisiumV1 -#' @include interoperability.R -#' @description Creates a VisiumV1 object -#' @slot image array. -#' @slot scale.factors scalefactors. -#' @slot coordinates data.frame. -#' @slot spot.radius numeric. -#' @importClassesFrom SeuratObject SpatialImage -#' @return visium object -#' @export - - VisiumV1 <- setClass(Class = 'VisiumV1', - contains = 'SpatialImage', - slots = list('image' = 'array', - 'scale.factors' = 'scalefactors', - 'coordinates' = 'data.frame', - 'spot.radius' = 'numeric' )) - \ No newline at end of file diff --git a/man/VisiumV1-class.Rd b/man/VisiumV1-class.Rd deleted file mode 100644 index d0c10278..00000000 --- a/man/VisiumV1-class.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/visiumImage.R -\docType{class} -\name{VisiumV1-class} -\alias{VisiumV1-class} -\alias{VisiumV1} -\title{Create VisiumV1 Class -VisiumV1} -\value{ -visium object -} -\description{ -Creates a VisiumV1 object -} -\section{Slots}{ - -\describe{ -\item{\code{image}}{array.} - -\item{\code{scale.factors}}{scalefactors.} - -\item{\code{coordinates}}{data.frame.} - -\item{\code{spot.radius}}{numeric.} -}} - From e31d3b1a77efb5ba3e98fd01a8c4068112a70195 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 01:27:23 -0500 Subject: [PATCH 17/24] chore: fix `edge_distances()` example --- R/NN_network.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/NN_network.R b/R/NN_network.R index 220819c6..9e13e611 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -470,7 +470,7 @@ createNetwork <- function( #' node IDs that apply to the coords in x must be supplied as a character vector #' @examples #' m <- matrix(c(0,0,0, 1,1,1, 3,2,4), byrow = TRUE, nrow = 3) -#' edges <- data.table( +#' edges <- data.table::data.table( #' from = c(1, 1), #' to = c(2, 3) #' ) From 0395a95f0f2fb9c539b28307d2a384d7c08a160d Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 02:28:57 -0500 Subject: [PATCH 18/24] fix: VisiumV1 definition --- R/interoperability.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/interoperability.R b/R/interoperability.R index 4c74f0ec..8660537e 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -1640,7 +1640,7 @@ giottoToSeuratV5 <- function(gobject, hires = max(img), lowres = min(img)) - senv <- rlang::ns_env("SeuratObject") + senv <- rlang::ns_env("Seurat") newV1 <- do.call( new, From 5c2dd6c49e8c7c7a026197f9220a6c85f6575db6 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 14:40:44 -0500 Subject: [PATCH 19/24] chore: add verbose default value --- R/NN_network.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/NN_network.R b/R/NN_network.R index 9e13e611..c983ed81 100644 --- a/R/NN_network.R +++ b/R/NN_network.R @@ -212,7 +212,7 @@ createNetwork <- function( .net_dt_knn <- function( x, k = 30L, include_weight = TRUE, include_distance = TRUE, filter = FALSE, maximum_distance = NULL, minimum_k = 0L, weight_fun = function(d) 1 / (1 + d), - verbose, ... + verbose = NULL, ... ) { # NSE vars from <- to <- distance <- NULL @@ -268,7 +268,7 @@ createNetwork <- function( .net_dt_snn <- function( x, k = 30L, include_weight = TRUE, include_distance = TRUE, top_shared = 3L, minimum_shared = 5L, weight_fun = function(d) 1 / (1 + d), - verbose, ... + verbose = NULL, ... ) { # NSE vars from <- to <- shared <- distance <- NULL From 74b6acb8a638c2c96c7ff073f964be3abce3d103 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 14:53:08 -0500 Subject: [PATCH 20/24] chore: remove extra code in favor of `loadNamespace()` --- R/interoperability.R | 23 +++++------------------ 1 file changed, 5 insertions(+), 18 deletions(-) diff --git a/R/interoperability.R b/R/interoperability.R index 8660537e..7f269a83 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -1381,7 +1381,7 @@ giottoToSeuratV5 <- function(gobject, # verify if optional package is installed package_check(pkg_name = "Seurat", repository = "CRAN") - requireNamespace("Seurat") + loadNamespace("Seurat") # check whether any raw data exist -- required for Seurat avail_expr <- list_expression(gobject = gobject, spat_unit = spat_unit) @@ -1640,23 +1640,10 @@ giottoToSeuratV5 <- function(gobject, hires = max(img), lowres = min(img)) - senv <- rlang::ns_env("Seurat") - - newV1 <- do.call( - new, - args = list( - Class = "VisiumV1", - image = img, - scale.factors = scalefactors, - coordinates = coord - ), - envir = senv - ) - - # newV1 <- new(Class = "VisiumV1", - # image = img, - # scale.factors = scalefactors, - # coordinates = coord) + newV1 <- new(Class = "VisiumV1", + image = img, + scale.factors = scalefactors, + coordinates = coord) sobj@images[[gobject@largeImages[[i]]@name]] <- newV1 From 8414c8ac23b1e812f847692cc9345023dd8b146b Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 14:57:50 -0500 Subject: [PATCH 21/24] chore: document example change for `edge_distances()` --- man/edge_distances.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/edge_distances.Rd b/man/edge_distances.Rd index a458ddaf..7d4ea5c5 100644 --- a/man/edge_distances.Rd +++ b/man/edge_distances.Rd @@ -21,7 +21,7 @@ Calculate network edge euclidean distances } \examples{ m <- matrix(c(0,0,0, 1,1,1, 3,2,4), byrow = TRUE, nrow = 3) -edges <- data.table( +edges <- data.table::data.table( from = c(1, 1), to = c(2, 3) ) From f7b3dc5ddc03d2c02e4f3c446be1b32f397c45d9 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 22:21:19 -0500 Subject: [PATCH 22/24] fix: Bugfix and improve docs for `calculateOverlap` --- R/aggregate.R | 85 +++++++++++++++++++++++++++++----------- man/calculateOverlap.Rd | 87 +++++++++++++++++++++++++++++++---------- 2 files changed, 128 insertions(+), 44 deletions(-) diff --git a/R/aggregate.R b/R/aggregate.R index 2e434c8f..15c5ca9c 100644 --- a/R/aggregate.R +++ b/R/aggregate.R @@ -58,16 +58,72 @@ polygon_to_raster <- function(polygon, field = NULL) { #' @title Calculate features overlapped by polygons #' @name calculateOverlap +#' @description +#' Calculate subcellular points/feature info or image values overlapped by polygon +#' annotations. This provides a summary of the spatial data overlapped by the +#' polygon which can be further processed to become an expression matrix. #' @param x Object with spatial annotations: `giottoPolygon`, or `SpatVector` #' polygons. Can also be a `giotto` object #' @param y Object with features to overlap: `giottoPoints`, `giottoLargeImage`, #' `SpatVector` points or `SpatRaster` -#' @param poly_subset_ids (optional) character vector of poly_IDs to use -#' @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 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 in `feat_subset_column` +#' in order to subset feature points when performing overlap calculation. +#' @param count_info_column character. (optional) column with count information. +#' Useful in cases when more than one detection is reported per point. #' @param verbose be verbose -#' @param \dots additional params to pass to methods +#' @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 +#' 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. +#' @returns Usually an object of the same class as `x`, with the overlaps +#' information appended. `return_*` logical params usually allow return of +#' a lower-level representation of the results instead. Only the +#' `SpatVector,SpatRaster` method is different in that it returns a `data.table`. +#' @examples +#' g <- GiottoData::loadGiottoMini("vizgen") +#' gpoly <- getPolygonInfo(g, polygon_name = "aggregate", +#' return_giottoPolygon = TRUE) +#' gpoints <- getFeatureInfo(g, return_giottoPoints = TRUE) +#' gimg <- getGiottoImage(g, image_type = "largeImage") +#' +#' gpoly@overlaps <- NULL +#' overlaps(gpoly) # Should now be NULL +#' +#' # detections from 2 z-layers are provided +#' table(gpoints$global_z) +#' +#' # calculate all transcripts overlapped +#' out_all <- calculateOverlap(gpoly, gpoints) +#' overlaps_all <- overlaps(out_all) +#' overlaps_all$rna +#' +#' # calculate z1 only +#' out_z1 <- calculateOverlap(gpoly, gpoints, +#' feat_subset_column = "global_z", +#' feat_subset_ids = c(1)) +#' overlaps_z1 <- overlaps(out_z1) +#' overlaps_z1$rna +#' +#' # overlap image to get sum intensities per cell +#' out_img <- calculateOverlap(gpoly, gimg) +#' overlaps_img <- overlaps(out_img) +#' overlaps_img$intensity +#' +#' # giotto method +#' # calculate z0 overlaps and return as gobject +#' out_g <- calculateOverlap(g, feat_subset_column = "global_z", +#' feat_subset_ids = 0) +#' overlaps(getPolygonInfo(out_g, return_giottoPolygon = TRUE)) +#' +#' # note that z0 and z1 nrows match that from the table of global z values. +#' # With points overlaps, all points are returned, but non-overlapped points +#' # only have an `NA` value for the `poly_ID` column. Overlapped points will +#' have the `poly_ID` of their overlapping polygon. NULL # * gobject #### @@ -86,16 +142,12 @@ setMethod( feat_info = NULL, image_names = NULL, poly_subset_ids = NULL, - feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, return_gobject = TRUE, verbose = TRUE, ...) { # 0. guards # # --------- # - if (!is.null(count_info_column)) checkmate::assert_character(count_info_column) if (!is.null(feat_info)) checkmate::assert_character(feat_info) if (!is.null(image_names)) checkmate::assert_character(image_names) if (!is.null(feat_info) && !is.null(image_names)) { @@ -204,12 +256,9 @@ setMethod( overlap_args_list <- list( x = A, - y = B, + y = B, # points or c(images) name_overlap = name_overlap, poly_subset_ids = poly_subset_ids, - feat_subset_column = feat_subset_column, - feat_subset_ids = feat_subset_ids, - count_info_column = count_info_column, verbose = verbose, return_gpolygon = isTRUE(return_gobject), ... @@ -244,16 +293,6 @@ setMethod( #' @param return_gpolygon default = TRUE. Whether to return the entire #' giottoPolygon provided to `x`, but with the overlaps information appended or #' as a bare terra `SpatVector` -#' @examples -#' \dontrun{ -#' x <- GiottoData::loadSubObjectMini("giottoPolygon") -#' y <- GiottoData::loadSubObjectMini("giottoPoints") -#' x@overlaps <- NULL -#' overlaps(x) # Should now be NULL -#' -#' a <- calculateOverlap(x, y) -#' overlaps(a) -#' } #' @export setMethod( "calculateOverlap", signature(x = "giottoPolygon", y = "giottoPoints"), diff --git a/man/calculateOverlap.Rd b/man/calculateOverlap.Rd index e28f4657..3d335cae 100644 --- a/man/calculateOverlap.Rd +++ b/man/calculateOverlap.Rd @@ -17,9 +17,6 @@ feat_info = NULL, image_names = NULL, poly_subset_ids = NULL, - feat_subset_column = NULL, - feat_subset_ids = NULL, - count_info_column = NULL, return_gobject = TRUE, verbose = TRUE, ... @@ -83,38 +80,86 @@ polygons. Can also be a \code{giotto} object} \item{image_names}{character vector. Name(s) of the image feature information to overlap} -\item{poly_subset_ids}{(optional) character vector of poly_IDs to use} - -\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{count_info_column}{column with count information (optional)} +\item{poly_subset_ids}{character vector. (optional) Specific poly_IDs to use} \item{return_gobject}{return giotto object (default: TRUE)} \item{verbose}{be verbose} -\item{\dots}{additional params to pass to methods} +\item{\dots}{additional params to pass to methods.} \item{y}{Object with features to overlap: \code{giottoPoints}, \code{giottoLargeImage}, \code{SpatVector} points or \code{SpatRaster}} +\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 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. +Useful in cases when more than one detection is reported per point.} + \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}} } +\value{ +Usually an object of the same class as \code{x}, with the overlaps +information appended. \verb{return_*} logical params usually allow return of +a lower-level representation of the results instead. Only the +\verb{SpatVector,SpatRaster} method is different in that it returns a \code{data.table}. +} \description{ -Calculate features overlapped by polygons +Calculate subcellular points/feature info or image values overlapped by polygon +annotations. This provides a summary of the spatial data overlapped by the +polygon which can be further processed to become an expression matrix. } -\examples{ -\dontrun{ -x <- GiottoData::loadSubObjectMini("giottoPolygon") -y <- GiottoData::loadSubObjectMini("giottoPoints") -x@overlaps <- NULL -overlaps(x) # Should now be NULL - -a <- calculateOverlap(x, y) -overlaps(a) +\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 +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. } +\examples{ +g <- GiottoData::loadGiottoMini("vizgen") +gpoly <- getPolygonInfo(g, polygon_name = "aggregate", + return_giottoPolygon = TRUE) +gpoints <- getFeatureInfo(g, return_giottoPoints = TRUE) +gimg <- getGiottoImage(g, image_type = "largeImage") + +gpoly@overlaps <- NULL +overlaps(gpoly) # Should now be NULL + +# detections from 2 z-layers are provided +table(gpoints$global_z) + +# calculate all transcripts overlapped +out_all <- calculateOverlap(gpoly, gpoints) +overlaps_all <- overlaps(out_all) +overlaps_all$rna + +# calculate z1 only +out_z1 <- calculateOverlap(gpoly, gpoints, + feat_subset_column = "global_z", + feat_subset_ids = c(1)) +overlaps_z1 <- overlaps(out_z1) +overlaps_z1$rna + +# overlap image to get sum intensities per cell +out_img <- calculateOverlap(gpoly, gimg) +overlaps_img <- overlaps(out_img) +overlaps_img$intensity + +# giotto method +# calculate z0 overlaps and return as gobject +out_g <- calculateOverlap(g, feat_subset_column = "global_z", + feat_subset_ids = 0) +overlaps(getPolygonInfo(out_g, return_giottoPolygon = TRUE)) + +# note that z0 and z1 nrows match that from the table of global z values. +# With points overlaps, all points are returned, but non-overlapped points +# only have an `NA` value for the `poly_ID` column. Overlapped points will +have the `poly_ID` of their overlapping polygon. } From 571aa35fdc57e85510dc49f1c5dd4dcc91e7d645 Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 22:27:48 -0500 Subject: [PATCH 23/24] chore: fix formatting in example --- R/aggregate.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/aggregate.R b/R/aggregate.R index 15c5ca9c..65007ba9 100644 --- a/R/aggregate.R +++ b/R/aggregate.R @@ -123,7 +123,7 @@ polygon_to_raster <- function(polygon, field = NULL) { #' # note that z0 and z1 nrows match that from the table of global z values. #' # With points overlaps, all points are returned, but non-overlapped points #' # only have an `NA` value for the `poly_ID` column. Overlapped points will -#' have the `poly_ID` of their overlapping polygon. +#' # have the `poly_ID` of their overlapping polygon. NULL # * gobject #### From 4e8ea2640a14a64d892ea3f4ec55d962df67e83a Mon Sep 17 00:00:00 2001 From: jiajic <72078254+jiajic@users.noreply.github.com> Date: Wed, 14 Feb 2024 22:35:53 -0500 Subject: [PATCH 24/24] chore: document --- man/calculateOverlap.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/calculateOverlap.Rd b/man/calculateOverlap.Rd index 3d335cae..6ad5f3b0 100644 --- a/man/calculateOverlap.Rd +++ b/man/calculateOverlap.Rd @@ -161,5 +161,5 @@ overlaps(getPolygonInfo(out_g, return_giottoPolygon = TRUE)) # note that z0 and z1 nrows match that from the table of global z values. # With points overlaps, all points are returned, but non-overlapped points # only have an `NA` value for the `poly_ID` column. Overlapped points will -have the `poly_ID` of their overlapping polygon. +# have the `poly_ID` of their overlapping polygon. }