From 7b2b283cea7474a6f803ff3e5d2be53b48ed716e Mon Sep 17 00:00:00 2001 From: Crystal Shin <54536735+sjcshin@users.noreply.github.com> Date: Thu, 30 May 2024 15:54:50 -0400 Subject: [PATCH] SpatialData conversions (both directions) + error fixes --- .RData | Bin 0 -> 49 bytes NAMESPACE | 2 + R/interoperability.R | 916 +++++++++++++++++++++---------- inst/python/g2ad.py | 2 +- inst/python/g2sd.py | 47 ++ inst/python/sd2g.py | 372 +++++++++++++ man/anndataToGiotto.Rd | 20 +- man/giottoToAnnData.Rd | 5 - man/giottoToSeuratV5.Rd | 14 +- man/giottoToSpatialData.Rd | 43 ++ man/giottoToSpatialExperiment.Rd | 6 + man/seuratToGiottoV5.Rd | 6 - man/spatialExperimentToGiotto.Rd | 8 +- man/spatialdataToGiotto.Rd | 58 ++ 14 files changed, 1160 insertions(+), 339 deletions(-) create mode 100644 .RData create mode 100644 inst/python/g2sd.py create mode 100644 inst/python/sd2g.py create mode 100644 man/giottoToSpatialData.Rd create mode 100644 man/spatialdataToGiotto.Rd diff --git a/.RData b/.RData new file mode 100644 index 0000000000000000000000000000000000000000..72af5a548858b470167ed95252bc62dad57d7e6e GIT binary patch literal 49 zcmb2|=3oE=X6~X+gJ)e2k`fXU(h?HW(h|~GjU*$So$r+BHuTgpoIQH#iULq20AR2W A{Qv*} literal 0 HcmV?d00001 diff --git a/NAMESPACE b/NAMESPACE index eca152c0..ef56593c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -149,6 +149,7 @@ export(giottoToAnnData) export(giottoToSeurat) export(giottoToSeuratV4) export(giottoToSeuratV5) +export(giottoToSpatialData) export(giottoToSpatialExperiment) export(hexVertices) export(installGiottoEnvironment) @@ -272,6 +273,7 @@ export(spatQueryGiottoPolygons) export(spatValues) export(spat_net_to_igraph) export(spatialExperimentToGiotto) +export(spatialdataToGiotto) export(standardise_flex) export(stitchFieldCoordinates) export(stitchGiottoLargeImage) diff --git a/R/interoperability.R b/R/interoperability.R index c89fbe79..74466136 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -21,9 +21,9 @@ #' See SAW pipeline for additional information about the gef file. #' @returns giotto object #' @export -gefToGiotto <- function( - gef_file, bin_size = "bin100", verbose = FALSE, - h5_file = NULL) { + +gefToGiotto <- function(gef_file, bin_size = "bin100", verbose = FALSE, + h5_file = NULL) { # data.table vars genes <- gene_idx <- x <- y <- sdimx <- sdimy <- cell_ID <- bin_ID <- count <- i.bin_ID <- NULL @@ -106,6 +106,7 @@ gefToGiotto <- function( #' @title Check Scanpy Installation #' @name check_py_for_scanpy +#' @import reticulate #' @returns character #' @description checks current python environment for scanpy 1.9.0 #' @keywords internal @@ -129,8 +130,8 @@ check_py_for_scanpy <- function() { \n ", errWidth = TRUE)) } else if (module_test == FALSE && genv_in_use) { - wrap_msg("Python module scanpy is required for conversion. - Installing scanpy now in the Giotto Miniconda Environment.") + cat("Python module scanpy is required for conversion. + Installing scanpy now in the Giotto Miniconda Environment.\n") conda_path <- reticulate::miniconda_path() py_ver <- reticulate::py_config()$version_string @@ -148,8 +149,8 @@ check_py_for_scanpy <- function() { python_version = py_ver ) } else { - wrap_msg("Required Python module scanpy has been previously installed. - Proceeding with conversion.") + cat("Required Python module scanpy has been previously installed. + Proceeding with conversion.\n") } } @@ -164,19 +165,19 @@ check_py_for_scanpy <- function() { #' neighbors(). If multiple spatial networks are in the anndata object, a list #' of key_added terms may be provided. If converting an anndata object from #' giottoToAnnData, a .txt file may be provided, which was generated in that -#' function, i.e. \{spat_unit\}_\{feat_type\}_nn_network_keys_added.txt. Cannot +#' function, i.e. {spat_unit}_{feat_type}_nn_network_keys_added.txt. Cannot #' be "spatial". This becomes the name of the nearest network in the gobject. #' @param spatial_n_key_added equivalent of "key_added" argument from #' squidpy.gr.spatial_neighbors. If multiple spatial networks are in the #' anndata object, a list of key_added terms may be provided. If converting an #' anndata object from giottoToAnnData, a .txt file may be provided, which was #' generated in that function, -#' i.e. \{spat_unit\}_\{feat_type\}_spatial_network_keys_added.txt +#' i.e. {spat_unit}_{feat_type}_spatial_network_keys_added.txt #' Cannot be the same as n_key_added. #' @param delaunay_spat_net binary parameter for spatial network. If TRUE, the #' spatial network is a delaunay network. -#' @param spat_unit desired spatial unit to use for conversion, default NULL -#' @param feat_type desired feature type to use for conversion, default NULL +#' @param spat_unit desired spatial unit for conversion, default NULL +#' @param feat_type desired feature type for conversion, default NULL #' @param h5_file name to create and on-disk HDF5 file #' @param python_path path to python executable within a conda/miniconda #' environment @@ -188,17 +189,6 @@ check_py_for_scanpy <- function() { #' See \code{\link{changeGiottoInstructions}} to modify instructions after #' creation. #' @returns Giotto object -#' @examples -#' g <- GiottoData::loadGiottoMini("visium") -#' # create saved object to test -#' temp_directory <- tempdir() -#' giottoToAnnData(g, save_directory = paste0(temp_directory, "/")) -#' -#' g2 <- anndataToGiotto(anndata_path = file.path( -#' temp_directory, -#' "/cell_rna_converted_gobject.h5ad" -#' )) -#' force(g2) #' @export anndataToGiotto <- function( anndata_path = NULL, @@ -209,8 +199,7 @@ anndataToGiotto <- function( feat_type = NULL, h5_file = NULL, python_path = NULL, - env_name = "giotto_env" -) { + env_name = "giotto_env") { # Preliminary file checks and guard clauses if (is.null(anndata_path)) { stop("Please provide a path to an AnnData .h5ad file for conversion.\n") @@ -237,7 +226,7 @@ anndataToGiotto <- function( # should trigger a stop() downstream if not installed # Import ad2g, a python module for parsing anndata - ad2g_path <- system.file("python", "ad2g.py", package = "Giotto") + ad2g_path <- system.file("python", "ad2g.py", package = "GiottoClass") reticulate::source_python(ad2g_path) adata <- read_anndata_from_path(anndata_path) @@ -606,17 +595,14 @@ anndataToGiotto <- function( #' The save_directory will be created if it does not already exist. #' The default save_directory is the working directory. #' @returns vector containing .h5ad file path(s) -#' @examples -#' g <- GiottoData::loadGiottoMini("visium") -#' -#' giottoToAnnData(g, save_directory = paste0(tempdir(), "/")) #' @export -giottoToAnnData <- function(gobject = NULL, - spat_unit = NULL, - feat_type = NULL, - python_path = NULL, - env_name = "giotto_env", - save_directory = NULL) { +giottoToAnnData <- function( + gobject = NULL, + spat_unit = NULL, + feat_type = NULL, + python_path = NULL, + env_name = "giotto_env", + save_directory = NULL) { # Check gobject invalid_obj <- !("giotto" %in% class(gobject)) if (is.null(gobject) || invalid_obj) { @@ -636,10 +622,10 @@ giottoToAnnData <- function(gobject = NULL, } else if (!dir.exists(save_directory)) { warning(wrap_msg("Provided save directory not found. Creating save directory at location:")) - wrap_msg(save_directory) + cat(save_directory) dir.create(save_directory, recursive = TRUE) if (dir.exists(save_directory)) { - wrap_msg("Created directory", save_directory) + cat("Created directory", save_directory) } else { stop(wrap_msg("Unable to create directory. Please change the provided path and try again.")) @@ -658,9 +644,9 @@ giottoToAnnData <- function(gobject = NULL, if (is.null(spat_unit) && is.null(feat_type)) { spat_unit <- unique(expr_dt$spat_unit) feat_type <- unique(expr_dt$feat_type) - } else if (is.null(spat_unit && !is.null(feat_type))) { + } else if (is.null(spat_unit) && !is.null(feat_type)) { spat_unit <- unique(expr_dt$spat_unit) - } else if (!is.null(spat_unit && is.null(feat_type))) { + } else if (!is.null(spat_unit) && is.null(feat_type)) { feat_type <- unique(expr_dt$feat_type) } @@ -774,7 +760,7 @@ giottoToAnnData <- function(gobject = NULL, # Feat Metadata for (su in spat_unit) { for (ft in names(gobject@expression[[su]])) { - cmeta <- getCellMetadata( + cmeta <- get_cell_metadata( gobject = gobject, spat_unit = su, feat_type = ft, @@ -782,7 +768,7 @@ giottoToAnnData <- function(gobject = NULL, set_defaults = FALSE ) - fm <- getFeatureMetadata( + fm <- get_feature_metadata( gobject = gobject, spat_unit = su, feat_type = ft, @@ -805,14 +791,15 @@ giottoToAnnData <- function(gobject = NULL, # Dimension Reductions # error hanldling wrapper to get_dimReduction - try_get_dimReduction <- function(gobject, - spat_unit, - feat_type, - reduction, - reduction_method, - name, - output, - set_defaults) { + try_get_dimReduction <- function( + gobject, + spat_unit, + feat_type, + reduction, + reduction_method, + name, + output, + set_defaults) { tryCatch( { dim_red <- get_dimReduction( @@ -960,13 +947,14 @@ giottoToAnnData <- function(gobject = NULL, # Nearest Neighbor Network # error hanldling wrapper to get_NearestNetwork - try_get_NN <- function(gobject, - spat_unit, - feat_type, - nn_network_to_use, - network_name, - output, - set_defaults) { + try_get_NN <- function( + gobject, + spat_unit, + feat_type, + nn_network_to_use, + network_name, + output, + set_defaults) { tryCatch( { nearest_net <- get_NearestNetwork( @@ -1046,12 +1034,13 @@ giottoToAnnData <- function(gobject = NULL, # Reset indexing variable adata_pos <- 1 - try_get_SN <- function(gobject, - spat_unit, - name, - output, - set_defaults, - verbose) { + try_get_SN <- function( + gobject, + spat_unit, + name, + output, + set_defaults, + verbose) { tryCatch( { spatial_net <- get_spatialNetwork( @@ -1174,10 +1163,11 @@ giottoToAnnData <- function(gobject = NULL, #' @param ... additional params to pass to \code{\link{get_spatial_locations}} #' @returns Seurat object #' @export -giottoToSeurat <- function(gobject, - spat_unit = NULL, - obj_use = NULL, - ...) { +giottoToSeurat <- function( + gobject, + spat_unit = NULL, + obj_use = NULL, + ...) { stop(wrap_txt( "Deprecated. Please use either giottoToSeuratV4() or giottoToSeuratV5()" )) @@ -1197,9 +1187,10 @@ giottoToSeurat <- function(gobject, #' @returns Seurat object #' @keywords seurat interoperability #' @export -giottoToSeuratV4 <- function(gobject, - spat_unit = NULL, - ...) { +giottoToSeuratV4 <- function( + gobject, + spat_unit = NULL, + ...) { # data.table vars feat_type <- name <- dim_type <- nn_type <- NULL # set default spat_unit and feat_type to be extracted as a Seurat assay @@ -1296,7 +1287,7 @@ giottoToSeuratV4 <- function(gobject, } # add cell metadata meta_cells <- data.table::setDF( - getCellMetadata( + get_cell_metadata( gobject = gobject, spat_unit = spat_unit, feat_type = assay_use, @@ -1314,7 +1305,7 @@ giottoToSeuratV4 <- function(gobject, ) # add feature metadata meta_genes <- data.table::setDF( - getFeatureMetadata( + get_feature_metadata( gobject = gobject, spat_unit = spat_unit, feat_type = assay_use, @@ -1465,23 +1456,16 @@ giottoToSeuratV4 <- function(gobject, #' The default values are 'cell' and 'rna' respectively. #' @param gobject Giotto object #' @param spat_unit spatial unit (e.g. 'cell') -#' @param res_type type of 10x image output resolution #' @param ... additional params to pass to \code{\link{get_spatial_locations}} #' @returns Seurat object #' @keywords seurat interoperability -#' @examples -#' g <- GiottoData::loadGiottoMini("visium") -#' -#' giottoToSeuratV5(g) #' @export -giottoToSeuratV5 <- function(gobject, - spat_unit = NULL, - res_type = c("hires", "lowres", "fullres"), - ...) { +giottoToSeuratV5 <- function( + gobject, + spat_unit = NULL, + ...) { # data.table vars feat_type <- name <- dim_type <- nn_type <- NULL - - res_type <- match.arg(res_type, choices = c("hires", "lowres", "fullres")) # set default spat_unit and feat_type to be extracted as a Seurat assay spat_unit <- set_default_spat_unit( @@ -1581,7 +1565,7 @@ giottoToSeuratV5 <- function(gobject, # add cell metadata meta_cells <- data.table::setDF( - getCellMetadata( + get_cell_metadata( gobject = gobject, spat_unit = spat_unit, feat_type = assay_use, @@ -1603,7 +1587,7 @@ giottoToSeuratV5 <- function(gobject, # add feature metadata meta_genes <- data.table::setDF( - getFeatureMetadata( + get_feature_metadata( gobject = gobject, spat_unit = spat_unit, feat_type = assay_use, @@ -1702,18 +1686,15 @@ giottoToSeuratV5 <- function(gobject, } # spatial coordinates - loc_use <- getSpatialLocations( - gobject = gobject, - spat_unit = spat_unit, - output = "spatLocsObj", - copy_obj = TRUE, - ... # allow setting of spat_loc_name through additional params + loc_use <- data.table::setDF( + get_spatial_locations( + gobject = gobject, + spat_unit = spat_unit, + output = "data.table", + copy_obj = TRUE, + ... # allow setting of spat_loc_name through additional params + ) ) - - # flip y vals - loc_use <- flip(loc_use)[] %>% - data.table::setDF() - rownames(loc_use) <- loc_use$cell_ID sobj <- Seurat::AddMetaData(sobj, metadata = loc_use) # add spatial coordinates as new dim reduct object @@ -1758,53 +1739,35 @@ giottoToSeuratV5 <- function(gobject, all_x <- NULL all_y <- NULL - - gimgs <- getGiottoImage(gobject, name = ":all:") - - if (length(gimgs) > 0) { - for (i in seq_along(gimgs)) { - gimg <- gimgs[[i]] - key <- objName(gimg) - imagerow <- loc_use$sdimy - imagecol <- loc_use$sdimx - img_array <- as(gimg, "array") - img_array[, , seq_len(3)] <- img_array[, , seq_len(3)] / 255 - coord <- data.frame( - imagerow = imagerow, imagecol = imagecol, - row.names = loc_use$cell_ID - ) - - scalef <- .estimate_scalefactors( - gimg, - res_type = res_type, - spatlocs = loc_use - ) - - # There does not seem to be a way to tell seurat which image type - # you are using. The lowres scalefactor seems to be the important - # one in mapping the image + + if (length(gobject@largeImages) > 0) { + for (i in seq_along(gobject@largeImages)) { + # spatVec <- terra::as.points( + # gobject@largeImages[[i]]@raster_object) + # geomSpatVec <- terra::geom(spatVec) + # x <- geomSpatVec[,"x"] + # y <- geomSpatVec[,"y"] + imagerow <- gobject@spatial_locs$cell$raw$sdimy + imagecol <- gobject@spatial_locs$cell$raw$sdimx + img <- terra::as.array(gobject@largeImages[[i]]@raster_object) + img[, , seq_len(3)] <- img[, , seq_len(3)] / 255 + coord <- data.frame(imagerow = imagerow, imagecol = imagecol) + scalefactors <- Seurat::scalefactors( - spot = scalef$spot, - fiducial = scalef$fiducial, - hires = scalef$hires, - lowres = scalef[res_type] # this looks like the main one - # so instead of strictly supplying lowres scalef, we use - # the scalef belonging to whichever image was used in Giotto - # since we allow use non-lowres images + spot = gobject@largeImages[[i]]@scale_factor, + fiducial = gobject@largeImages[[i]]@resolution, + hires = max(img), + lowres = min(img) ) - # see https://github.com/satijalab/seurat/issues/3595 newV1 <- new( Class = "VisiumV1", - image = img_array, + image = img, scale.factors = scalefactors, - coordinates = coord, - spot.radius = - scalef$fiducial * scalef$lowres / max(dim(img_array)), - key = paste0(key, "_") + coordinates = coord ) - sobj@images[[key]] <- newV1 + sobj@images[[gobject@largeImages[[i]]@name]] <- newV1 } } @@ -1814,80 +1777,6 @@ giottoToSeuratV5 <- function(gobject, } -#' @param x image object -#' @param res_type type of 10x image output resolution -#' @param spatlocs a data.frame of spatial locations coordinates -#' @noRd -.estimate_scalefactors <- function( - x, res_type = c("hires", "lowres", "fullres"), spatlocs -) { - res_type <- match.arg(res_type, choices = c("hires", "lowres", "fullres")) - - pxdims <- dim(x)[1:2] - edims <- range(ext(x)) - - scalef <- mean(pxdims / edims) - - # assume that lowres and hires follow a general ratio - # may not be that important since the scalefactor should theoretically - # only matter for the image res that we are using - - # this ratio is roughly 3.333334 based on Visium BreastCancerA1 dataset - res_ratio <- 3.333334 - - # fullres should have a scalef of roughly 1. - # No way to guess hires or lowres scalefs so use arbitrary values. - - hres_scalef <- switch(res_type, - "hires" = scalef, - "lowres" = scalef * res_ratio, - "fullres" = 0.08250825 # arbitrary - ) - - lres_scalef <- switch(res_type, - "hires" = scalef / res_ratio, - "lowres" = scalef, - "fullres" = 0.02475247 # arbitrary - ) - - # spot diameter and fid diameter are variable based on how spatial info was - # mapped to the image. Estimate this by getting the center to center - # px distance vs fullsize px dims ratio. - # ! fullsize px dims is the same as edims ! - - coords <- data.table::as.data.table(spatlocs) - # create a delaunay - dnet <- createNetwork( - as.matrix(coords[, c("sdimx", "sdimy")]), - type = "delaunay", - method = "geometry", - include_distance = TRUE, - as.igraph = FALSE, - include_weight = TRUE, - verbose = FALSE - ) - - # expect center to center be most common edge distance - # this gives CC dist as fullres px distance - distances <- sort(unique(dnet$distance)) - cc_px <- distances[which.max(table(dnet$distance))] - - # assume constant ratios between diameters and cc_px - fid_cc_ratio <- 1.045909 - fid_diam <- cc_px * fid_cc_ratio - - spot_cc_ratio <- 0.6474675 - spot_diam <- cc_px * spot_cc_ratio - - scalef_list <- list( - spot = spot_diam, - fiducial = fid_diam, - hires = hres_scalef, - lowres = lres_scalef - ) - return(scalef_list) -} - #' @title Deprecated @@ -1903,10 +1792,11 @@ giottoToSeuratV5 <- function(gobject, #' object. Default is \code{"Vizgen"}. #' @returns giotto object #' @export -seuratToGiotto <- function(sobject, - spatial_assay = "Spatial", - dim_reduction = c("pca", "umap"), - subcellular_assay = "Vizgen") { +seuratToGiotto <- function( + sobject, + spatial_assay = "Spatial", + dim_reduction = c("pca", "umap"), + subcellular_assay = "Vizgen") { stop(wrap_txt( "Deprecated. Please use either seuratToGiottoV4() or seuratToGiottoV5()" )) @@ -1932,13 +1822,14 @@ seuratToGiotto <- function(sobject, #' stored in it. #' @keywords seurat interoperability #' @export -seuratToGiottoV4 <- function(sobject, - spatial_assay = "Spatial", - dim_reduction = c("pca", "umap"), - subcellular_assay = "Vizgen", - sp_network = NULL, - nn_network = NULL, - verbose = TRUE) { +seuratToGiottoV4 <- function( + sobject, + spatial_assay = "Spatial", + dim_reduction = c("pca", "umap"), + subcellular_assay = "Vizgen", + sp_network = NULL, + nn_network = NULL, + verbose = TRUE) { package_check("Seurat") if (is.null(Seurat::GetAssayData( object = sobject, slot = "counts", @@ -1966,10 +1857,9 @@ seuratToGiottoV4 <- function(sobject, # Cell Metadata cell_metadata <- sobject@meta.data # Dimension Reduction - if (sum(vapply( + if (sum(sapply( dim_reduction, - function(x) length(sobject@reductions[[x]]), - FUN.VALUE = integer(1L) + function(x) length(sobject@reductions[[x]]) )) == 0) { dimReduc_list <- NULL } else { @@ -2196,23 +2086,16 @@ seuratToGiottoV4 <- function(sobject, #' @returns A Giotto object converted from Seurat object with all computations #' stored in it. #' @keywords seurat interoperability -#' @examples -#' m_expression <- Matrix::Matrix(rnorm(100), nrow = 10, sparse = TRUE) -#' s <- Seurat::CreateSeuratObject(counts = m_expression) -#' -#' seuratToGiottoV5(s, spatial_assay = "RNA") #' @export -seuratToGiottoV5 <- function(sobject, - spatial_assay = "Spatial", - dim_reduction = c("pca", "umap"), - subcellular_assay = "Vizgen", - sp_network = NULL, - nn_network = NULL, - verbose = TRUE) { +seuratToGiottoV5 <- function( + sobject, + spatial_assay = "Spatial", + dim_reduction = c("pca", "umap"), + subcellular_assay = "Vizgen", + sp_network = NULL, + nn_network = NULL, + verbose = TRUE) { package_check("Seurat") - - # NSE vars - sdimy <- NULL if (is.null(Seurat::GetAssayData( object = sobject, slot = "counts", @@ -2251,21 +2134,13 @@ seuratToGiottoV5 <- function(sobject, # Cell Metadata cell_metadata <- sobject@meta.data - cell_metadata <- data.table::as.data.table( - cell_metadata, keep.rownames = TRUE) - # Feat Metadata feat_metadata <- sobject[[]] - feat_metadata <- data.table::as.data.table( - feat_metadata, keep.rownames = TRUE) - - # rownames of both kept as `rn` # Dimension Reduction - if (sum(vapply( + if (sum(sapply( dim_reduction, - function(x) length(sobject@reductions[[x]]), - FUN.VALUE = integer(1L) + function(x) length(sobject@reductions[[x]]) )) == 0) { dimReduc_list <- NULL } else { @@ -2313,28 +2188,35 @@ seuratToGiottoV5 <- function(sobject, object = sobject, assay = spatial_assay ))) { - spat_coord <- Seurat::GetTissueCoordinates(sobject, - scale = NULL, - cols = c( - "imagerow", - "imagecol" - ) - ) + spat_coord <- Seurat::GetTissueCoordinates(sobject) + # spat_coord = cbind(rownames(spat_coord), + # data.frame(spat_coord, row.names=NULL)) if (!("cell" %in% spat_coord)) { spat_coord$cell_ID <- rownames(spat_coord) - colnames(spat_coord) <- c("sdimy", "sdimx", "cell_ID") + colnames(spat_coord) <- c("sdimx", "sdimy", "cell_ID") } else { - colnames(spat_coord) <- c("sdimy", "sdimx", "cell_ID") + colnames(spat_coord) <- c("sdimx", "sdimy", "cell_ID") } - spat_loc <- data.table::as.data.table(spat_coord) - - # seurat has coords following imaging conventions - # flip them for Giotto - spat_loc[, sdimy := -sdimy] - data.table::setcolorder( - spat_loc, neworder = c("sdimx", "sdimy", "cell_ID")) + spat_loc <- spat_coord + length_assay <- length(colnames(sobject)) + + spat_datatable <- data.table( + cell_ID = character(length_assay), + sdimx = rep(NA_real_, length_assay), + sdimy = rep(NA_real_, length_assay) + ) + + spat_datatable$cell_ID <- colnames(sobject) + match_cell_ID <- match(spat_loc$cell_ID, spat_datatable$cell_ID) + matching_indices <- match_cell_ID + matching_indices <- matching_indices[!is.na(matching_indices)] + spat_datatable[ + matching_indices, + c("sdimx", "sdimy") := list(spat_loc$sdimx, spat_loc$sdimy) + ] + spat_loc <- spat_datatable } else { message("Images for RNA assay not found in the data. Skipping image processing.") @@ -2413,20 +2295,30 @@ seuratToGiottoV5 <- function(sobject, } } - # Find SueratImages, extract them, and pass to create image + # Find SueratImages, extract them, and pass to create seuratobj + for (i in names(sobject@images)) { # check if image slot has image in it - simg <- sobject[[i]] - if ("image" %in% slotNames(simg)) { - img_array <- slot(simg, "image") - if (!is.null(img_array)) { - - scalef <- Seurat::ScaleFactors(simg) - + if ("image" %in% slotNames(sobject@images[[i]])) { + if (!is.null(sobject@images[[i]]@image)) { + # Extract the red (r), green (g), and blue (b) channels + 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) + + # Convert channels to rasters + r <- terra::rast(r) + g <- terra::rast(g) + b <- terra::rast(b) + + # Create Giotto LargeImage gImg <- createGiottoLargeImage( - raster_object = terra::rast(img_array) * 255, - name = i, - scale_factor = 1 / scalef$lowres + raster_object = terra::rast(list(r, g, b)), + name = names(sobject@images) ) } } @@ -2544,30 +2436,25 @@ seuratToGiottoV5 <- function(sobject, ) } } - - gobject <- addCellMetadata( - gobject = gobject, new_metadata = cell_metadata, - by_column = TRUE, column_cell_ID = "rn") - gobject <- addFeatMetadata( - gobject = gobject, new_metadata = feat_metadata, - by_column = TRUE, column_feat_ID = "rn") + gobject <- addCellMetadata(gobject = gobject, new_metadata = cell_metadata) + gobject <- addFeatMetadata(gobject = gobject, new_metadata = feat_metadata) - if (exists("gpoints")) { + if (exists("gpoints") == TRUE) { gobject <- addGiottoPoints( gobject = gobject, gpoints = list(gpoints) ) } - if (exists("gpolygon")) { + if (exists("gpolygon") == TRUE) { gobject <- addGiottoPolygons( gobject = gobject, gpolygons = polygon_list ) } - if (exists("gImg")) { + if (exists("gImg") == TRUE) { gobject <- addGiottoLargeImage( gobject = gobject, largeImages = list(gImg) @@ -2591,6 +2478,11 @@ seuratToGiottoV5 <- function(sobject, #' #' @returns A SpatialExperiment object that contains data from the input Giotto #' object. +#' @examples +#' \dontrun{ +#' mini_gobject <- GiottoData::loadGiottoMini("vizgen") +#' giottoToSpatialExperiment(mini_gobject) +#' } #' @export giottoToSpatialExperiment <- function(giottoObj, verbose = TRUE) { spat_unit <- NULL @@ -2726,7 +2618,7 @@ giottoToSpatialExperiment <- function(giottoObj, verbose = TRUE) { ) } SpatialExperiment::spatialCoords(spe) <- data.matrix( - spatialLocs[, c("sdimx", "sdimy")] + spatialLocs[, seq_len(2)] ) } else { if (verbose) { @@ -2918,17 +2810,21 @@ giottoToSpatialExperiment <- function(giottoObj, verbose = TRUE) { #' networks. This can be a vector of multiple network names. #' @param verbose A boolean value specifying if progress messages should #' be displayed or not. Default \code{TRUE}. +#' @import data.table #' @returns Giotto object #' @examples -#' spe <- STexampleData::Visium_humanDLPFC() -#' -#' spatialExperimentToGiotto(spe, python_path = NULL) +#' \dontrun{ +#' library(SpatialExperiment) +#' example(read10xVisium, echo = FALSE) +#' spatialExperimentToGiotto(spe) +#' } #' @export -spatialExperimentToGiotto <- function(spe, - python_path, - nn_network = NULL, - sp_network = NULL, - verbose = TRUE) { +spatialExperimentToGiotto <- function( + spe, + python_path, + nn_network = NULL, + sp_network = NULL, + verbose = TRUE) { # Create giotto instructions and set python path instrs <- createGiottoInstructions(python_path = python_path) @@ -3168,8 +3064,10 @@ if (requireNamespace("SpatialExperiment", quietly = TRUE)) { #' #' @returns A Giotto object compatible with suite version #' @export -giottoMasterToSuite <- function(gobject, - expression_feat = "rna") { +#' +giottoMasterToSuite <- function( + gobject, + expression_feat = "rna") { master_object <- gobject spatial_locs <- cell_metadata <- feat_metadata <- instructions <- NULL @@ -3349,3 +3247,431 @@ giottoMasterToSuite <- function(gobject, return(gobject) } + + + +## SpatialData object to Giotto #### + +#' @title Convert SpatialData to Giotto +#' @name spatialdataToGiotto +#' @description Converts a saved SpatialData object into a Giotto object +#' +#' @param spatialdata_path path to SpatialData object +#' @param n_key_added equivalent of "key_added" argument from scanpy.pp.neighbors(). +#' If multiple spatial networks are in the anndata object, a list of key_added +#' terms may be provided. +#' If converting an anndata object from giottoToAnnData, a .txt file may be +#' provided, which was generated in that function, +#' i.e. {spat_unit}_{feat_type}_nn_network_keys_added.txt +#' Cannot be "spatial". This becomes the name of the nearest network in the gobject. +#' @param spatial_n_key_added equivalent of "key_added" argument from squidpy.gr.spatial_neighbors. +#' If multiple spatial networks are in the anndata object, a list of key_added +#' terms may be provided. +#' If converting an anndata object from giottoToAnnData, a .txt file may be +#' provided, which was generated in that function, +#' i.e. {spat_unit}_{feat_type}_spatial_network_keys_added.txt +#' Cannot be the same as n_key_added. +#' @param delaunay_spat_net binary parameter for spatial network. If TRUE, the spatial network is a delaunay network. +#' @param spat_unit desired spatial unit for conversion, default NULL +#' @param feat_type desired feature type for conversion, default NULL +#' @param python_path path to python executable within a conda/miniconda environment +#' @param env_name name of environment containing python_path executable +#' +#' @return Giotto object +#' @details Function in beta. Converts a structured SpatialData file into a Giotto object. +#' The returned Giotto Object will take default insructions with the +#' exception of the python path, which may be customized. +#' See \code{\link{changeGiottoInstructions}} to modify instructions after creation. +#' @export + +spatialdataToGiotto <- function( + spatialdata_path = NULL, + n_key_added = NULL, + spatial_n_key_added = NULL, + delaunay_spat_net = TRUE, + spat_unit = NULL, + feat_type = NULL, + python_path = NULL, + env_name = NULL) { + + # File check + if (is.null(spatialdata_path)) { + stop("Please provide a path to SpatialData object for conversion.\n") + } + if (!file.exists(spatialdata_path)) { + stop("The provided path to SpatialData object does not exist.\n") + } + + # Initialize reticulate + instrs <- createGiottoInstructions( + python_path = python_path, + show_plot = FALSE, + return_plot = FALSE, + save_plot = TRUE, + save_dir = NULL, + plot_format = NULL, + dpi = NULL, + units = NULL, + height = NULL, + width = NULL, + is_docker = FALSE, + plot_count = 0, + fiji_path = NULL, + no_python_warn = FALSE + ) + + # Check spatialdata dependencies + spatialdata_installed <- checkPythonPackage(package_name = "spatialdata", env_to_use = env_name) + + # Import sd2g, a python module for parsing SpatialData + sd2g_path <- system.file("python", "sd2g.py", package = "GiottoClass") + reticulate::source_python(sd2g_path) + sdata <- read_spatialdata_from_path(spatialdata_path) + + # Extract expression matrix + expr_df <- extract_expression(sdata) + cID <- extract_cell_IDs(sdata) + fID <- extract_feat_IDs(sdata) + + # Extract spatial locations + spatial_df <- extract_spatial(sdata) + sp <- parse_obsm_for_spat_locs(sdata) + + # Set up metadata + cm <- extract_cell_metadata(sdata) + cm <- as.data.table(cm) + if ("leiden" %in% names(cm)) { + cm$leiden <- as.numeric(cm$leiden) + } + + fm <- extract_feat_metadata(sdata) + fm <- as.data.table(fm) + + # Create baseline Giotto object + gobject <- createGiottoObject( + expression = expr_df, + spatial_locs = spatial_df, + instructions = instrs + ) + + # Attach hires image + raster <- terra::rast(extract_image(sdata)) + giotto_image <- createGiottoLargeImage(raster) + gobject <- addGiottoLargeImage(gobject = gobject, largeImages = c(giotto_image)) + + # Attach metadata + cm <- readCellMetadata(cm) + gobject <- setCellMetadata(gobject, x = cm) + fm <- readFeatMetadata(fm) + gobject <- setFeatureMetadata(gobject, x = fm) + + spat_unit <- activeSpatUnit(gobject) + feat_type <- activeFeatType(gobject) + + # Add PCA + p <- extract_pca(sdata) + if (!is.null(p)) { + pca <- p$pca + evs <- p$eigenvalues + loads <- p$loadings + # Add pca to giottoObject + dobj <- createDimObj( + coordinates = pca, + name = "pca", + spat_unit = spat_unit, + feat_type = feat_type, + method = "pca", + reduction = "cells", + provenance = NULL, + misc = list( + eigenvalues = evs, + loadings = loads + ), + my_rownames = colnames(expr_df) + ) + gobject <- set_dimReduction(gobject = gobject, dimObject = dobj) + } + + # Add UMAP + u <- extract_umap(sdata) + if (!is.null(u)) { + # Add UMAP to giottoObject + dobj <- createDimObj( + coordinates = u, + name = "umap", + spat_unit = spat_unit, + feat_type = feat_type, + method = "umap", + reduction = "cells", + provenance = NULL, + misc = NULL, + my_rownames = colnames(expr_df) + ) + gobject <- set_dimReduction(gobject = gobject, dimObject = dobj) + } + + # Add tSNE + t <- extract_tsne(sdata) + if (!is.null(t)) { + # Add TSNE to giottoObject + dobj <- createDimObj( + coordinates = t, + name = "tsne", + spat_unit = spat_unit, + feat_type = feat_type, + method = "tsne", + reduction = "cells", + provenance = NULL, + misc = NULL, + my_rownames = colnames(expr_df) + ) + gobject <- set_dimReduction(gobject = gobject, dimObject = dobj) + } + + ## Nearest Network + + weights_sd <- NULL + num_NN_nets <- length(n_key_added) + + if (is.null(n_key_added) && !is.null(extract_NN_connectivities(sdata, key_added = n_key_added))) { + num_NN_nets <- 1 + } + + for (i in num_NN_nets) { + if (inherits(n_key_added, "list")) { + n_key_added_it <- n_key_added[[i]] + } else { + n_key_added_it <- n_key_added + } + + weights_sd <- extract_NN_connectivities(sdata, key_added = n_key_added_it) + # adw = methods::as(weights_ad, "TsparseMatrix") + if (!is.null(weights_sd)) { + distances_sd <- extract_NN_distances(sdata, key_added = n_key_added_it) + + nn_dt <- align_network_data(distances = weights_sd, weights = distances_sd) + + # pre-allocate DT variables + from <- to <- weight <- distance <- from_cell_ID <- to_cell_ID <- uniq_ID <- NULL + nn_dt <- data.table::data.table(nn_dt) + + nn_dt[, from_cell_ID := cID[from]] + nn_dt[, to_cell_ID := cID[to]] + nn_dt[, uniq_ID := paste0(from, to)] + nn_dt[order(uniq_ID)] + nn_dt[, uniq_ID := NULL] + vert <- unique(x = c(nn_dt$from_cell_ID, nn_dt$to_cell_ID)) + nn_network_igraph <- igraph::graph_from_data_frame(nn_dt[, .(from_cell_ID, to_cell_ID, weight, distance)], directed = TRUE, vertices = vert) + + nn_info <- extract_NN_info(adata = adata, key_added = n_key_added_it) + + net_type <- "kNN" # anndata default + if (("sNN" %in% n_key_added_it) & !is.null(n_key_added_it)) { + net_type <- "sNN" + net_name <- paste0(n_key_added_it, ".", nn_info["method"]) + } else if (!("sNN" %in% n_key_added_it) & !is.null(n_key_added_it)) { + net_name <- paste0(n_key_added_it, ".", nn_info["method"]) + } else { + net_name <- paste0(net_type, ".", nn_info["method"]) + } + + netObj <- createNearestNetObj( + name = net_name, + network = nn_network_igraph, + spat_unit = spat_unit, + feat_type = feat_type + ) + + gobject <- set_NearestNetwork( + gobject = gobject, + nn_network = netObj, + spat_unit = spat_unit, + feat_type = feat_type, + nn_network_to_use = net_type, + network_name = net_name, + set_defaults = FALSE + ) + } + } + + ## Spatial Network + + s_weights_sd <- NULL + num_SN_nets <- length(spatial_n_key_added) + + # Check for the case where NULL is provided, since the + # anndata object takes the default value for SN + + if (is.null(spatial_n_key_added) && !is.null(extract_SN_connectivities(sdata, key_added = spatial_n_key_added))) { + num_SN_nets <- 1 + } + + for (i in 1:num_SN_nets) { + if (inherits(spatial_n_key_added, "list")) { + spatial_n_key_added_it <- spatial_n_key_added[[i]] + } else { + spatial_n_key_added_it <- spatial_n_key_added + } + + s_weights_sd <- extract_SN_connectivities(sdata, key_added = spatial_n_key_added_it) + if (!is.null(s_weights_sd)) { + s_distances_sd <- extract_SN_distances(sdata, key_added = spatial_n_key_added_it) + ij_matrix <- methods::as(s_distances_sd, "TsparseMatrix") + from_idx <- ij_matrix@i + 1 # zero index!!! + to_idx <- ij_matrix@j + 1 # zero index!!! + + # pre-allocate DT variables + from <- to <- weight <- distance <- from_cell_ID <- to_cell_ID <- uniq_ID <- NULL + sn_dt <- data.table::data.table( + from = from_idx, + to = to_idx, + weight = s_weights_sd@x, + distance = s_distances_sd@x + ) + + sn_dt[, from_cell_ID := cID[from]] + sn_dt[, to_cell_ID := cID[to]] + + sdimx <- "sdimx" + sdimy <- "sdimy" + xbegin_name <- paste0(sdimx, "_begin") + ybegin_name <- paste0(sdimy, "_begin") + xend_name <- paste0(sdimx, "_end") + yend_name <- paste0(sdimy, "_end") + + network_DT <- data.table::data.table( + from = sn_dt$from_cell_ID, + to = sn_dt$to_cell_ID, + xbegin_name = sp[sn_dt$from, sdimx], + ybegin_name = sp[sn_dt$from, sdimy], + xend_name = sp[sn_dt$to, sdimx], + yend_name = sp[sn_dt$to, sdimy], + weight = s_weights_sd@x, + distance = s_distances_sd@x + ) + data.table::setnames(network_DT, + old = c("xbegin_name", "ybegin_name", "xend_name", "yend_name"), + new = c(xbegin_name, ybegin_name, xend_name, yend_name) + ) + data.table::setorder(network_DT, from, to) + + dist_mean <- get_distance(network_DT, method = "mean") + dist_median <- get_distance(network_DT, method = "median") + cellShapeObj <- list( + "meanCellDistance" = dist_mean, + "medianCellDistance" = dist_median + ) + + # TODO filter network? + # TODO 3D handling? + if (delaunay_spat_net) { + spatObj <- create_spat_net_obj( + name = "sNN", + method = "delaunay", + networkDT = network_DT, + cellShapeObj = cellShapeObj + ) + } else { + spatObj <- create_spat_net_obj( + name = "sNN", + method = "non-delaunay", + networkDT = network_DT, + cellShapeObj = cellShapeObj + ) + } + + gobject <- set_spatialNetwork( + gobject = gobject, + spatial_network = spatObj, + name = "sNN" + ) + } + } + return(gobject) +} + + + +## Giotto to SpatialData#### + +#' @title Convert Giotto to SpatialData +#' @name giottoToSpatialData +#' @description Converts a Giotto object to SpatialData object +#' +#' @param gobject giotto object to be converted +#' @param spat_unit spatial unit which will be used in conversion +#' @param feat_type feature type which will be used in conversion +#' @param spot_radius radius of the spots +#' @param python_path path to python executable within a conda/miniconda environment +#' @param env_name name of environment containing python_path executable +#' @param save_directory directory in which the SpatialData object will be saved +#' +#' @return SpatialData object saved on disk. +#' @details Function in beta. Converts and saves a Giotto object in SpatialData format on disk. +#' @export + +giottoToSpatialData <- function( + gobject = NULL, + spat_unit = NULL, + feat_type = NULL, + spot_radius = NULL, + python_path = NULL, + env_name = NULL, + save_directory = NULL) { + + # Initialize reticulate + instrs <- createGiottoInstructions(python_path = python_path) + + # Check spatialdata dependencies + spatialdata_installed <- checkPythonPackage(package_name = "spatialdata", env_to_use = env_name) + + # Import sd2g, a python module for parsing SpatialData + g2sd_path <- system.file("python", "g2sd.py", package = "GiottoClass") + reticulate::source_python(g2sd_path) + + # Get metadata + spat_unit <- activeSpatUnit(gobject) + feat_type <- activeFeatType(gobject) + + # Create a temporary folder to hold anndata + temp <- "temp_anndata/" + + # First, convert Giotto object to AnnData using an existing function + giottoToAnnData( + gobject = gobject, + spat_unit = spat_unit, + feat_type = feat_type, + python_path = python_path, + env_name = env_name, + save_directory = temp + ) + + # Extract GiottoImage + gimg <- getGiottoImage(gobject, image_type = "largeImage") + + # Temporarily save the image to disk + writeGiottoLargeImage( + giottoLargeImage = gimg, + gobject = gobject, + largeImage_name = "largeImage", + filename = "temp_image.png", + dataType = NULL, + max_intensity = NULL, + overwrite = TRUE, + verbose = TRUE + ) + + spat_locs <- getSpatialLocations(gobject, output="data.table") + + # Create SpatialData object + createSpatialData(temp, spat_locs, spot_radius, save_directory) + + # Delete temporary files and folders + unlink("temp_image.png") + unlink("temp_image.png.aux.xml") + unlink(temp, recursive = TRUE) + + # Successful Conversion + cat("Giotto object has been converted and saved to SpatialData object at: ", save_directory, "\n") + +} diff --git a/inst/python/g2ad.py b/inst/python/g2ad.py index b46e329a..33614524 100644 --- a/inst/python/g2ad.py +++ b/inst/python/g2ad.py @@ -133,7 +133,7 @@ def set_adg_pca(adata = None, pca_coord = None, loadings = None, eigenv = None, adata.varm["PCs"] = pc_placehold elif loadings is not None: - adata.varm["PCs"] = loadings + adata.varm["PCs"] = loadings.to_numpy(dtype=float) if eigenv is not None: eigenv_shape = len(eigenv) eigenv = np.array(eigenv) diff --git a/inst/python/g2sd.py b/inst/python/g2sd.py new file mode 100644 index 00000000..5f7d2233 --- /dev/null +++ b/inst/python/g2sd.py @@ -0,0 +1,47 @@ +import anndata as ad +import numpy as np +import pandas as pd +from dask_image.imread import imread +from xarray import DataArray +import geopandas as gpd +from shapely.geometry import Point +import os + +from spatialdata import SpatialData +from spatialdata.models import Image2DModel, ShapesModel, TableModel +from spatialdata.transformations.transformations import Identity + +def createImageModel(): + images = {} + hires_image_path = "temp_image.png" + hires_img = imread(hires_image_path).squeeze().transpose(2,0,1) + hires_img = DataArray(hires_img, dims=("c","y","x")) + images["hires_image"] = Image2DModel.parse(hires_img, transformations={"downscaled_hires": Identity()}) + return images + +def createShapeModel(spat_locs, spot_radius): + shapes_df = pd.DataFrame() + shapes_df['x'] = spat_locs.sdimx + shapes_df['y'] = spat_locs.sdimy + shapes_df['geometry'] = shapes_df.apply(lambda row: Point(row['x'], row['y']), axis=1) + shapes_df = shapes_df.drop(['x','y'], axis=1) + shapes_df['radius'] = spot_radius + shapes_df = shapes_df.rename_axis('spot_id') + gdf = gpd.GeoDataFrame(shapes_df, geometry='geometry') + gdf.set_crs(epsg=4326, inplace=True) + shapes = ShapesModel.parse(gdf) + return shapes + +def createTableModel(temp): + alist = os.listdir(temp)[0] + adata = ad.read_h5ad(os.path.join(temp, alist)) + table = TableModel.parse(adata) + return table + +def createSpatialData(temp, spat_locs, spot_radius, save_directory): + images = createImageModel() + table = createTableModel(temp) + shapes = createShapeModel(spat_locs, spot_radius) + sd = SpatialData(table = table, images = images) + sd.shapes["Shapes"] = shapes + sd.write(save_directory) diff --git a/inst/python/sd2g.py b/inst/python/sd2g.py new file mode 100644 index 00000000..c63aee1e --- /dev/null +++ b/inst/python/sd2g.py @@ -0,0 +1,372 @@ +import warnings +import numpy as np +import pandas as pd +import scipy as sc +import matplotlib.pyplot as plt +from scipy.sparse import csr_matrix, csc_matrix +import dask.array as da +from time import perf_counter + +from spatialdata import SpatialData + +__name__ = "sd2g" +__package__ = "sd2g" + +# Imports and safeguards +def read_spatialdata_from_path(sd_path = None): + if sd_path is None: + print("Please provide a path to the SpatialData object folder.") + assert(False) + try: + sdata = SpatialData.read(sd_path) + except (FileNotFoundError): + print(f"File {sd_path} was not found.") + print("Please try again.") + raise(FileNotFoundError) + return sdata + +# Extract gene expression +def extract_expression(sdata = None): + expr = sdata.table.X.transpose().todense() + expr_df = pd.DataFrame(expr, index=sdata.table.var['gene_ids'].index, columns=sdata.table.obs['array_row'].index) + return expr_df + +# Extract cell IDs +def extract_cell_IDs(sdata = None): + cell_IDs = sdata.table.obs['array_row'].index.tolist() + return cell_IDs + +# Extract feature IDs +def extract_feat_IDs(sdata = None): + feat_IDs = sdata.table.var['gene_ids'].index.tolist() + return feat_IDs + +# Metadata +def extract_cell_metadata(sdata = None): + cell_metadata = sdata.table.obs.reset_index() + cell_metadata = cell_metadata.rename(columns={"index": "cell_ID"}) + return cell_metadata + +def extract_feat_metadata(sdata = None): + feat_metadata = sdata.table.var.reset_index() + feat_metadata = feat_metadata.rename(columns={"index": "feat_ID"}) + return feat_metadata + +# Alternative expression data +def extract_layer_names(sdata = None): + layer_names = None + if len(sdata.table.layers) > 0: + layer_names = [i for i in sdata.table.layers] + return layer_names + +def extract_layered_data(sdata = None, layer_name = None): + layer_names = [i for i in sdata.table.layers] + if layer_name not in layer_names: + print(f"Invalid Key, {layer_name}, for sdata.table.layers") + raise(KeyError) + target_layer = sdata.table.layers[layer_name] + if type(target_layer) == sc.sparse.csr_matrix: + target_layer = target_layer.T + elif type(target_layer) == sc.sparse.csr.csr_matrix: + target_layer = target_layer.T + else: + target_layer = pd.DataFrame(target_layer) + return target_layer + +# Extract spatial information +def extract_spatial(sdata = None): + spatial = sdata.table.obsm['spatial'] + spatial_df = pd.DataFrame(spatial) + spatial_df.columns = ['X', 'Y'] + spatial_df['Y'] = spatial_df['Y'] * -1 + return spatial_df + +def parse_obsm_for_spat_locs(sdata = None): + cID = np.array(extract_cell_IDs(sdata)) + spat_locs = None + spat_key = None + + try: + spat_locs = sdata.table.obsm["spatial"] + except (KeyError): + spat_keys = [i for i in sdata.table.obsm if 'spatial' in i] + if len(spat_keys) > 0: + spat_key = spat_keys[0] + spat_locs = sdata.table.obsm[spat_key] + + if spat_locs is None: + err_mess = '''Spatial locations were not found. If spatial locations should have been found, + please modify the anndata table within SpatialData object to include a keyword-value pair within the obsm slot, + in which the keyword contains the phrase "spatial" and the value corresponds to the spatial locations.\n + In the Giotto Object resulting from this conversion, dummy locations will be used.''' + print(err_mess) + spat_locs = None + return spat_locs + else: + print("Spatial locations found.") + + cID = np.array(cID).reshape(len(cID),1) + spat_locs = np.concatenate((spat_locs,cID), axis = 1) + num_col = spat_locs.shape[1] + + colnames = ["sdimx","sdimy","sdimz","cell_ID"] + conv = {"sdimx":float, "sdimy":float,"sdimz":float,"cell_ID":str} + + if num_col > 3: + spat_locs = pd.DataFrame(spat_locs, columns = colnames) + else: + del colnames[2] + del conv['sdimz'] + spat_locs = pd.DataFrame(spat_locs, columns = colnames) + + spat_locs = spat_locs.astype(conv) + # Giotto y axis convention + spat_locs["sdimy"] = -1 * spat_locs["sdimy"] + return spat_locs + +# Extract hires image +def extract_image(sdata = None): + # Find SpatialData image name for hires image + for key in sdata.images.keys(): + if "hires" in key: + hires_image_name = key + + # Extract image from SpatialData and convert it to numpy array + hires_image = sdata.images[hires_image_name] + hires_image_array = np.transpose(hires_image.compute().data, (1, 2, 0)) # Transpose to (y, x, c) + return hires_image_array + +# Extract PCA +def extract_pca(sdata = None): + o_keys = sdata.table.obsm_keys() + v_keys = sdata.table.varm_keys() + u_keys = None + + pca = dict() + + for ok in o_keys: + if "X_pca" in ok: + pca['pca'] = sdata.table.obsm[ok] + u_keys = sdata.table.uns['pca'].keys() + + for vk in v_keys: + if "PCs" in vk: + pca['loadings'] = sdata.table.varm[vk] + + if type(u_keys) is not type(None): + for uk in u_keys: + if "variance" == uk: + pca['eigenvalues'] = sdata.table.uns['pca'][uk] + + if(len(pca)) == 0: + pca = None + + return pca + +# Extract UMAP +def extract_umap(sdata = None): + o_keys = sdata.table.obsm_keys() + + umap = None + + for ok in o_keys: + if "X_umap" in ok: + umap = sdata.table.obsm[ok] + + return umap + +# Extract tSNE +def extract_tsne(sdata = None): + o_keys = sdata.table.obsm_keys() + + tsne = None + + for ok in o_keys: + if "X_tsne" in ok: + tsne = sdata.table.obsm[ok] + + return tsne + +## NN Network + +def find_NN_keys(sdata = None, key_added = None): + nn_key_list = [] + + if key_added is None: + param_keys = list(sdata.table.uns.keys()) + for pk in param_keys: + if "neighbors" in pk and "spatial" not in pk: + try: + tmp_keys = sdata.table.uns[pk].keys() + except KeyError: + tmp_keys = None + return None + for i in tmp_keys: + nn_key_list.append(sdata.table.uns[pk].keys()) + break #Only returns connectivity and distance keys for one network + elif ".txt" in key_added: + line_keys = [] + with open(key_added) as f: + for line in f.readlines(): + line = line.strip() + line_keys.append(line) + + for key in line_keys: + map_keys = sdata.table.uns[key].keys() + for i in map_keys: + nn_key_list.append(sdata.table.uns[key][i]) + elif key_added and key_added.casefold() != "spatial": + map_keys = sdata.table.uns[key].keys() + for i in map_keys: + nn_key_list.append(sdata.table.uns[key][i]) + elif key_added and key_added.casefold() == "spatial": + s1 = "String 'spatial' cannot be used as n_key_added to retrieve a Nearest Neighbor Network. " + s2 = "This results from conflicting keys for nearest neighbor and spatial networks. " + s3 = "\nSee defaults here:\nhttps://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.neighbors.html\nhttps://squidpy.readthedocs.io/en/stable/api/squidpy.gr.spatial_neighbors.html" + msg = s1+ s2 + s3 + warnings.warn(msg) + + if len(nn_key_list) == 0: + nn_key_list = None + return nn_key_list + +def extract_NN_connectivities(sdata = None, key_added = None): + connectivities = None + nn_key_list = find_NN_keys(sdata = sdata, key_added = key_added) + + if type(nn_key_list) is type(None): + return connectivities + + for nk in nn_key_list: + if "connectivities" in nk: + connectivities = sdata.table.obsp[nk] + + return connectivities + +def extract_NN_distances(sdata = None, key_added = None): + distances = None + nn_key_list = find_NN_keys(sdata = sdata, key_added = key_added) + if type(nn_key_list) is type(None): + return distances + for nk in nn_key_list: + if "distances" in nk: + distances = sdata.table.obsp[nk] + return distances + +def extract_NN_info(sdata = None, key_added = None): + nn_keys = find_NN_keys(sdata, key_added=key_added) + nn_info = None + for nk in nn_keys: + if type(nk) is dict: + nn_info = pd.Series(nk) + return nn_info + +def align_network_data(distances = None, weights = None): + idx_dist_not_sparse = distances.nonzero() + blank = [0 for i in range(len(idx_dist_not_sparse[0]))] + df = pd.DataFrame({"distance":blank.copy(), "weight":blank.copy(), "from":blank.copy(), "to":blank.copy()}) + t0 = perf_counter() + + d_nz = distances[idx_dist_not_sparse] + d_nz = np.array(d_nz).reshape(len(d_nz.T),) + w_nz = weights[idx_dist_not_sparse] + w_nz = np.array(w_nz).reshape(len(w_nz.T),) + + df.loc[:,"distance"] = pd.Series(d_nz) + df.loc[:,"weight"] = pd.Series(w_nz) + with warnings.catch_warnings(): + warnings.simplefilter(action='ignore', category=(DeprecationWarning, FutureWarning)) + # Ignoring the warning here because the desired behavior is maintained + df.loc[:,"from"] = pd.Series(idx_dist_not_sparse[0]) + df.loc[:,"to"] = pd.Series(idx_dist_not_sparse[1]) + + df.loc[:,"from"] += 1 + df.loc[:,"to"] += 1 + # for x, i in enumerate(zip(*dist_not_sparse)): + # to = i[-1] + # from_ = i[0] + # dist = distances[from_,to] + # weig = weights[from_,to] + + # # Correct indexing convention for export to R + # from_ += 1 + # to += 1 + + # df.iloc[x,:] = {"distance":dist, "weight":weig, "from":from_, "to":to} + t1 = perf_counter() + print("Network extraction time:",t1-t0) + return df + +def find_SN_keys(sdata = None, key_added = None): + sn_key_list = [] + prefix = "" + suffix = "neighbors" + + if key_added is None: + map_key = prefix + suffix + try: + tmp_keys = sdata.table.uns[map_key].keys() + except KeyError: + tmp_keys = None + return None + + for i in tmp_keys: + #if type(adata.uns[pk][i]) == type(dict()): continue + sn_key_list.append(sdata.table.uns[map_key][i]) + elif ".txt" in key_added: + line_keys = [] + with open(key_added) as f: + for line in f.readlines(): + line = line.strip() + line_key_added = line + suffix + line_keys.append(line_key_added) + for key in line_keys: + map_keys = sdata.table.uns[key].keys() + for i in map_keys: + sn_key_list.append(sdata.table.uns[key][i]) + + elif key_added is not None: + key_added = key_added + suffix + map_keys = sdata.table.uns[key_added].keys() + for i in map_keys: + #if type(adata.uns[key_added][i]) == type(dict()): continue + sn_key_list.append(sdata.table.uns[key_added][i]) + + if len(sn_key_list) == 0: + sn_key_list = None + return sn_key_list + +def extract_SN_connectivities(sdata = None, key_added = None): + connectivities = None + sn_key_list = find_SN_keys(sdata = sdata, key_added = key_added) + + if type(sn_key_list) is type(None): + return connectivities + + for sk in sn_key_list: + if "connectivities" in sk: + connectivities = sdata.table.obsp[sk] + + return connectivities + +def extract_SN_distances(sdata = None, key_added = None): + distances = None + sn_key_list = find_SN_keys(sdata = sdata, key_added = key_added) + + if type(sn_key_list) is type(None): + return distances + + for sk in sn_key_list: + if "distances" in sk: + distances = sdata.table.obsp[sk] + + return distances + +def extract_SN_info(sdata = None, key_added = None): + sn_keys = find_SN_keys(sdata, key_added=key_added) + sn_info = None + for sk in sn_keys: + if type(sk) is dict: + sn_info = pd.Series(sk) + return sn_info + diff --git a/man/anndataToGiotto.Rd b/man/anndataToGiotto.Rd index 9d228a1d..e894dfe2 100644 --- a/man/anndataToGiotto.Rd +++ b/man/anndataToGiotto.Rd @@ -23,7 +23,7 @@ anndataToGiotto( neighbors(). If multiple spatial networks are in the anndata object, a list of key_added terms may be provided. If converting an anndata object from giottoToAnnData, a .txt file may be provided, which was generated in that -function, i.e. \{spat_unit\}_\{feat_type\}_nn_network_keys_added.txt. Cannot +function, i.e. {spat_unit}_{feat_type}_nn_network_keys_added.txt. Cannot be "spatial". This becomes the name of the nearest network in the gobject.} \item{spatial_n_key_added}{equivalent of "key_added" argument from @@ -31,15 +31,15 @@ squidpy.gr.spatial_neighbors. If multiple spatial networks are in the anndata object, a list of key_added terms may be provided. If converting an anndata object from giottoToAnnData, a .txt file may be provided, which was generated in that function, -i.e. \{spat_unit\}_\{feat_type\}_spatial_network_keys_added.txt +i.e. {spat_unit}_{feat_type}_spatial_network_keys_added.txt Cannot be the same as n_key_added.} \item{delaunay_spat_net}{binary parameter for spatial network. If TRUE, the spatial network is a delaunay network.} -\item{spat_unit}{desired spatial unit to use for conversion, default NULL} +\item{spat_unit}{desired spatial unit for conversion, default NULL} -\item{feat_type}{desired feature type to use for conversion, default NULL} +\item{feat_type}{desired feature type for conversion, default NULL} \item{h5_file}{name to create and on-disk HDF5 file} @@ -62,15 +62,3 @@ exception of the python path, which may be customized. See \code{\link{changeGiottoInstructions}} to modify instructions after creation. } -\examples{ -g <- GiottoData::loadGiottoMini("visium") -# create saved object to test -temp_directory <- tempdir() -giottoToAnnData(g, save_directory = paste0(temp_directory, "/")) - -g2 <- anndataToGiotto(anndata_path = file.path( - temp_directory, - "/cell_rna_converted_gobject.h5ad" -)) -force(g2) -} diff --git a/man/giottoToAnnData.Rd b/man/giottoToAnnData.Rd index 2aa5bd5d..57a08213 100644 --- a/man/giottoToAnnData.Rd +++ b/man/giottoToAnnData.Rd @@ -52,8 +52,3 @@ and feature type pair. The save_directory will be created if it does not already exist. The default save_directory is the working directory. } -\examples{ -g <- GiottoData::loadGiottoMini("visium") - -giottoToAnnData(g, save_directory = paste0(tempdir(), "/")) -} diff --git a/man/giottoToSeuratV5.Rd b/man/giottoToSeuratV5.Rd index 83eae66b..d06ceca7 100644 --- a/man/giottoToSeuratV5.Rd +++ b/man/giottoToSeuratV5.Rd @@ -4,20 +4,13 @@ \alias{giottoToSeuratV5} \title{Convert Giotto to Seurat V5} \usage{ -giottoToSeuratV5( - gobject, - spat_unit = NULL, - res_type = c("hires", "lowres", "fullres"), - ... -) +giottoToSeuratV5(gobject, spat_unit = NULL, ...) } \arguments{ \item{gobject}{Giotto object} \item{spat_unit}{spatial unit (e.g. 'cell')} -\item{res_type}{type of 10x image output resolution} - \item{...}{additional params to pass to \code{\link{get_spatial_locations}}} } \value{ @@ -28,10 +21,5 @@ Converts Giotto object into a Seurat object. This functions extracts specific sets of data belonging to specified spatial unit. The default values are 'cell' and 'rna' respectively. } -\examples{ -g <- GiottoData::loadGiottoMini("visium") - -giottoToSeuratV5(g) -} \keyword{interoperability} \keyword{seurat} diff --git a/man/giottoToSpatialData.Rd b/man/giottoToSpatialData.Rd new file mode 100644 index 00000000..25d1b268 --- /dev/null +++ b/man/giottoToSpatialData.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/interoperability.R +\name{giottoToSpatialData} +\alias{giottoToSpatialData} +\title{Convert Giotto to SpatialData} +\usage{ +giottoToSpatialData( + gobject = NULL, + spat_unit = NULL, + feat_type = NULL, + spot_radius = NULL, + python_path = NULL, + env_name = NULL, + save_directory = NULL +) +} +\arguments{ +\item{gobject}{giotto object to be converted} + +\item{spat_unit}{spatial unit which will be used in conversion} + +\item{feat_type}{feature type which will be used in conversion} + +\item{spot_radius}{radius of the spots} + +\item{python_path}{path to python executable within a conda/miniconda environment} + +\item{env_name}{name of environment containing python_path executable} + +\item{save_directory}{directory in which the SpatialData object will be saved} +} +\value{ +SpatialData object in the form of +} +\description{ +Converts a Giotto object to SpatialData object +} +\details{ +Function in beta. Converts a .h5ad file into a Giotto object. +The returned Giotto Object will take default insructions with the +exception of the python path, which may be customized. +See \code{\link{changeGiottoInstructions}} to modify instructions after creation. +} diff --git a/man/giottoToSpatialExperiment.Rd b/man/giottoToSpatialExperiment.Rd index 32af3950..9cc5f806 100644 --- a/man/giottoToSpatialExperiment.Rd +++ b/man/giottoToSpatialExperiment.Rd @@ -19,3 +19,9 @@ object. \description{ Utility function to convert a Giotto object to a SpatialExperiment object. } +\examples{ +\dontrun{ +mini_gobject <- GiottoData::loadGiottoMini("vizgen") +giottoToSpatialExperiment(mini_gobject) +} +} diff --git a/man/seuratToGiottoV5.Rd b/man/seuratToGiottoV5.Rd index 362f7a0c..b72e2f28 100644 --- a/man/seuratToGiottoV5.Rd +++ b/man/seuratToGiottoV5.Rd @@ -39,11 +39,5 @@ stored in it. \description{ Convert a Seurat V5 object to a Giotto object } -\examples{ -m_expression <- Matrix::Matrix(rnorm(100), nrow = 10, sparse = TRUE) -s <- Seurat::CreateSeuratObject(counts = m_expression) - -seuratToGiottoV5(s, spatial_assay = "RNA") -} \keyword{interoperability} \keyword{seurat} diff --git a/man/spatialExperimentToGiotto.Rd b/man/spatialExperimentToGiotto.Rd index 234064de..ee4c49d7 100644 --- a/man/spatialExperimentToGiotto.Rd +++ b/man/spatialExperimentToGiotto.Rd @@ -35,7 +35,9 @@ Giotto object Utility function to convert a SpatialExperiment object to a Giotto object } \examples{ -spe <- STexampleData::Visium_humanDLPFC() - -spatialExperimentToGiotto(spe, python_path = NULL) +\dontrun{ +library(SpatialExperiment) +example(read10xVisium, echo = FALSE) +spatialExperimentToGiotto(spe) +} } diff --git a/man/spatialdataToGiotto.Rd b/man/spatialdataToGiotto.Rd new file mode 100644 index 00000000..7bbdba3a --- /dev/null +++ b/man/spatialdataToGiotto.Rd @@ -0,0 +1,58 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/interoperability.R +\name{spatialdataToGiotto} +\alias{spatialdataToGiotto} +\title{Convert SpatialData to Giotto} +\usage{ +spatialdataToGiotto( + spatialdata_path = NULL, + n_key_added = NULL, + spatial_n_key_added = NULL, + delaunay_spat_net = TRUE, + spat_unit = NULL, + feat_type = NULL, + python_path = NULL, + env_name = NULL +) +} +\arguments{ +\item{spatialdata_path}{path to SpatialData object} + +\item{n_key_added}{equivalent of "key_added" argument from scanpy.pp.neighbors(). +If multiple spatial networks are in the anndata object, a list of key_added +terms may be provided. +If converting an anndata object from giottoToAnnData, a .txt file may be +provided, which was generated in that function, +i.e. {spat_unit}_{feat_type}_nn_network_keys_added.txt +Cannot be "spatial". This becomes the name of the nearest network in the gobject.} + +\item{spatial_n_key_added}{equivalent of "key_added" argument from squidpy.gr.spatial_neighbors. +If multiple spatial networks are in the anndata object, a list of key_added +terms may be provided. +If converting an anndata object from giottoToAnnData, a .txt file may be +provided, which was generated in that function, +i.e. {spat_unit}_{feat_type}_spatial_network_keys_added.txt +Cannot be the same as n_key_added.} + +\item{delaunay_spat_net}{binary parameter for spatial network. If TRUE, the spatial network is a delaunay network.} + +\item{spat_unit}{desired spatial unit for conversion, default NULL} + +\item{feat_type}{desired feature type for conversion, default NULL} + +\item{python_path}{path to python executable within a conda/miniconda environment} + +\item{env_name}{name of environment containing python_path executable} +} +\value{ +Giotto object +} +\description{ +Converts a saved SpatialData object into a Giotto object +} +\details{ +Function in beta. Converts a .h5ad file into a Giotto object. +The returned Giotto Object will take default insructions with the +exception of the python path, which may be customized. +See \code{\link{changeGiottoInstructions}} to modify instructions after creation. +}