From 289e6a9b0ae0aff37933d77ffc27cba11b18f8f6 Mon Sep 17 00:00:00 2001 From: iqraAmin Date: Wed, 1 May 2024 20:30:55 -0400 Subject: [PATCH] Added Support for Seurat and Giotto Interoperability for Sequencencing-based Images --- Giotto_Seurat_Interoperability.Rmd | 960 ++++++++++++++++++++++++++ NAMESPACE | 3 - R/interoperability.R | 10 +- man/aggregateStacksExpression.Rd | 4 +- man/aggregateStacksLocations.Rd | 4 +- man/aggregateStacksPolygonOverlaps.Rd | 4 +- man/aggregateStacksPolygons.Rd | 4 +- 7 files changed, 976 insertions(+), 13 deletions(-) create mode 100644 Giotto_Seurat_Interoperability.Rmd diff --git a/Giotto_Seurat_Interoperability.Rmd b/Giotto_Seurat_Interoperability.Rmd new file mode 100644 index 00000000..b55071ca --- /dev/null +++ b/Giotto_Seurat_Interoperability.Rmd @@ -0,0 +1,960 @@ +--- +title: "Giotto_Seurat_interoperability" +author: "Iqra Amin" +date: "2024-03-06" +output: + html_document: default + pdf_document: default +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +# Introduction + +Giotto facilitates seamless interoperability with various tools, including Seurat. The conversion between Giotto and Seurat relies on four primary functions. 'giottoToSeurat_v4' and 'SeuratToGiotto_v4' cater to Seurat version 4 , while 'giottoToSeurat_v5' and 'SeuratToGiotto_v5' are specifically for Seurat version 5. These functions play a vital role in ensuring effortless data transitions between Giotto and Seurat across different versions, ensuring compatibility and ease of use. + +## conversion of Seurat V5 to Giotto + +```{r message=FALSE, warning=FALSE, include=FALSE} +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") + + if (is.null(Seurat::GetAssayData(object = sobject, slot = "counts", assay = spatial_assay))) { + wrap_msg("No raw expression values are provided in spatial_assay") + return(sobject) + } else { + exp <- Seurat::GetAssayData(object = sobject, slot = "counts", assay = spatial_assay) + if (!is.null(sobject@assays$SCT)) { + normexp <- Seurat::GetAssayData(object = sobject, slot = "counts", assay = "SCT") + } + + if ("data" %in% slotNames(sobject@assays[[spatial_assay]])) { + if (!is.null(slot(sobject, "assays")[[spatial_assay]]@data)) { + normexp <- Seurat::GetAssayData(object = sobject, slot = "data", assay = spatial_assay) + } + } + if ("layers" %in% slotNames(sobject@assays[[spatial_assay]])) { + if (!is.null(slot(sobject, "assays")[[spatial_assay]]@layers)) { + normexp <- SeuratObject::LayerData(object = sobject, assay = spatial_assay) + } + } + + # Cell Metadata + cell_metadata <- sobject@meta.data + # Feat Metadata + feat_metadata <- sobject[[]] + + # Dimension Reduction + if (sum(sapply(dim_reduction, function(x) length(sobject@reductions[[x]]))) == 0) { + dimReduc_list <- NULL + } else { + dimReduc_list <- lapply(dim_reduction, function(x) { + dim_coord <- as.matrix(Seurat::Embeddings(object = sobject, reduction = x)) + dim_load <- as.matrix(Seurat::Loadings(object = sobject, reduction = x)) + dim_eig <- Seurat::Stdev(sobject, reduction = x) + if (length(dim_eig) > 0) { + dim_eig <- sapply(dim_eig, function(x) x^2) + } + colnames(dim_coord) <- paste0("Dim.", 1:ncol(dim_coord)) + if (length(dim_load) > 0) { + colnames(dim_load) <- paste0("Dim.", 1:ncol(dim_load)) + } + dimReduc <- create_dim_obj( + name = x, + reduction = "cells", + reduction_method = x, + spat_unit = "cell", + feat_type = "rna", + provenance = NULL, + coordinates = dim_coord, + misc = list(eigenvalues = dim_eig, loadings = dim_load) + ) + + return(dimReduc) + }) + # names(dimReduc_list) <- dim_reduction + } + + # Spatial Locations + if (length(sobject@assays[[spatial_assay]]) == 0) { + spat_loc <- NULL + } else { + # !requires image objects! + if (!is.null(Seurat::Images(object = sobject, assay = spatial_assay))) { + spat_coord <- Seurat::GetTissueCoordinates(sobject, + scale = NULL, + cols = c("imagerow", "imagecol")) + # 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("sdimx", "sdimy", "cell_ID") + } else { + colnames(spat_coord) <- 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.") + spat_loc <- NULL + } + } + # Subcellular + name <- names(sobject@images) + #if (!is.null(subcellular_assay)){ + if (length(sobject@assays[[subcellular_assay]]) == 1){ + spat_coord <- Seurat::GetTissueCoordinates(sobject) + colnames(spat_coord) <- c("sdimx", "sdimy", "cell_ID") + exp <- exp[, c(intersect(spat_coord$cell_ID, colnames(exp)))] + spat_loc <- spat_coord + } + #} + if (!length(sobject@images) == 0) { + for (i in names(sobject@images)) { + if ("molecules" %in% names(sobject@images[[i]]) == TRUE) { + if (!length(sobject@images[[i]][["molecules"]]) == 0) { + assay <- names(sobject@assays) + featnames <- rownames(sobject) + mol_spatlocs <- data.table::data.table() + + for (x in featnames) { + df <- (Seurat::FetchData(sobject[[i]][["molecules"]], vars = x)) + mol_spatlocs <- rbind(mol_spatlocs, df) + } + gpoints <- createGiottoPoints(mol_spatlocs, feat_type = "rna") + if ("centroids" %in% names(sobject@images[[i]])) { + centroids_coords <- sobject@images[[i]]$centroids@coords + centroids_coords <- vect(centroids_coords) + gpolygon <- create_giotto_polygon_object(name = "cell", spatVector = centroids_coords) + } + if ("segmentation" %in% names(sobject@images[[i]])) { + polygon_list <- list() + + for (j in seq(sobject@images$hippo@boundaries$segmentation@polygons)) { + polygon_info <- sobject@images$hippo@boundaries$segmentation@polygons[[j]] + + # Get coordinates from segmentation + seg_coords <- polygon_info@Polygons[[1]]@coords + + # Fetch cell_Id from polygon information + cell_ID <- polygon_info@ID + + # Convert it to SpatVector + seg_coords <- vect(seg_coords) + + # Create giotto_polygon_object + gpolygon <- create_giotto_polygon_object(name = "cell", + spatVector = centroids_coords, + spatVectorCentroids = seg_coords) + + # Add the cell_ID to the list of polygon names + polygon_list[[cell_ID]] <- gpolygon + } + } + } + } + } + } + } + + # Find SueratImages, extract them, and pass to create seuratobj + + for (i in names(sobject@images)) { + # check if image slot has image in it + 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(list(r, g, b)), + name = names(sobject@images), + scale_factor = sobject@images[[i]]@scale.factors$lowres ) + } + } + } + + gobject <- createGiottoObject(exp, + spatial_locs = spat_loc, + dimension_reduction = dimReduc_list + ) + if (exists("normexp") == TRUE) { + exprObj <- create_expr_obj( + name = "normalized", + exprMat = normexp, + spat_unit = "cell", + feat_type = "rna", + provenance = "cell" + ) + 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) + stempgraph <- as(stempgraph, "sparseMatrix") + sobjIgraph <- igraph::graph.adjacency(stempgraph, 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::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) + } + DT <- data.table() + + edges <- igraph::ends(sobjIgraph, es = igraph::E(sobjIgraph), names = TRUE) + DT$from <- edges[, 1] + DT$to <- edges[, 2] + ed_attr <- igraph::edge.attributes(sobjIgraph) + DT$weight <- ed_attr[1] + DT$distance <- ed_attr[2] + spatNetObj <- create_spat_net_obj( + networkDT = DT + ) + gobject <- set_spatialNetwork( + gobject = gobject, + spatial_network = spatNetObj, + name = sp_network[i] + ) + } + } + + ##nn_network + if (!is.null(nn_network)) { + nnNetObj_list <- lapply(seq_along(nn_network), function(i) { + stempgraph <- Seurat::as.Graph(sobject@graphs[[nn_network[i]]]) + #m <- as.matrix(stempgraph) + stempgraph <- as(stempgraph, "sparseMatrix") + sobjIgraph <- igraph::graph.adjacency(stempgraph, 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::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]]) + } + } + gobject <- addCellMetadata(gobject = gobject, new_metadata = cell_metadata) + gobject <- addFeatMetadata(gobject = gobject, new_metadata = feat_metadata) + + + if (exists("gpoints") == TRUE) { + gobject <- addGiottoPoints( + gobject = gobject, + gpoints = list(gpoints) + ) + } + + if (exists("gpolygon") == TRUE) { + gobject <- addGiottoPolygons(gobject = gobject, gpolygons = polygon_list) + } + + if (exists("gImg") == TRUE) { + gobject <- addGiottoLargeImage(gobject = gobject, largeImages = list(gImg)) + } + return(gobject) +} +``` + +### Load required Libraries + +To start, we load the required libraries. + +```{r message=FALSE, warning=FALSE} +library(data.table) +library(Giotto) +library(GiottoData) +library(Seurat) +library(SeuratData) +``` + +### Load Seurat object + +We begin by loading a sample Seurat object named "stxBrain" with the type specified as "anterior1". + +```{r message=FALSE, warning=FALSE} +brain <- LoadData("stxBrain", type = "anterior1") +``` + +### Convert Seurat object to Giotto + +To convert Seurat object to Giotto object, we use the function seuratToGiottoV5(). + +```{r message=FALSE, warning=FALSE} +brainG <- seuratToGiottoV5(sobject = brain, spatial_assay = "Spatial") +``` + +### Print Giotto object + +```{r message=FALSE, warning=FALSE} +brainG +``` + +### Downstream Analysis using Giotto Object + +The downstream analysis of the giotto object being created by seurat object is given below. + +#### High Efficiency Data Processing + +```{r processgiotto, message=FALSE, warning=FALSE} +brainG <- processGiotto( + brainG, + filter_params = list(expression_threshold = 1, + feat_det_in_min_cells = 100, + min_det_feats_per_cell = 10), + norm_params = list(norm_methods = 'standard', + scale_feats = TRUE, + scalefactor = 6000), + stat_params = list(expression_values = 'normalized'), + adjust_params = list(expression_values = c('normalized'), + covariate_columns = 'nr_feats') +) +``` + +#### Dimension Reduction and PCA + +First, we find HVF (highly varialble features) using the loess regression prediction model. + +```{r giotto-hvf, message=FALSE} +brainG <- calculateHVF(gobject = brainG, method = 'cov_loess') +``` + +then use HVF to compute the reduced dimensions using PCA. + +```{r giotto-pca, message=FALSE} +## Select genes highly variable genes that fit specified statistics +# These are both found within feature metadata +feature_metadata = getFeatureMetadata(brainG)[] +featgenes = feature_metadata[hvf == 'yes' & perc_cells > 4 & mean_expr_det > 0.5]$feat_ID + +## run PCA on expression values (default) +brainG <- Giotto::runPCA(gobject = brainG, feats_to_use = featgenes, scale_unit = F, center = F) +``` + +Here, we'll proceed with running UMAP based on PCA dimension reduction and visualize the pre-clustering UMAP. + +```{r giotto-umap, message=FALSE, warning=FALSE} +brainG <- Giotto::runUMAP(brainG, dimensions_to_use = 1:15) +graph1 <- Giotto::plotUMAP(gobject = brainG) + +``` + +#### Clustering + +Before proceeding with clustering cells based on their expression profiles, ensure that the Giotto Object has undergone PCA dimension reduction and includes either t-SNE or UMAP dimension reduction, as well as a defined neighbor network. + +To initialize the clustering process, create a shared nearest neighbor network (sNN) with a specified number of nearest neighbors ('k'): + +```{r giotto-nn} +brainG <- createNearestNetwork(gobject = brainG, type = "sNN", dimensions_to_use = 1:15, k = 15) +``` + +In Giotto, various clustering algorithms are available, such as k-means, Leiden, or Louvain. These algorithms automatically store cluster information within cell metadata, typically following default naming conventions. However, you can easily customize the cluster name by specifying the 'name' argument. + +```{r giotto-kmeans} +## k-means clustering +brainG <- doKmeans(gobject = brainG, dim_reduction_to_use = 'pca') + +#Plot UMAP post-clustering to visualize kmeans +graph2 <- Giotto::plotUMAP( + gobject = brainG, + cell_color = 'kmeans', + show_NN_network = T, + point_size = 2.5 +) +``` + +#### Spatial co-expression + +Using the spatial network previously established, we can employ the binSpect method to identify spatially variable genes. + +```{r giotto-binspect, message=FALSE} +brainG <- Giotto::createSpatialDelaunayNetwork(gobject = brainG) +showGiottoSpatNetworks(brainG) + +ranktest = binSpect( + brainG, bin_method = 'rank', + calc_hub = T, hub_min_int = 5, + spatial_network_name = 'Delaunay_network' +) +``` + +Here, we'll narrow down our analysis by selecting the top 300 spatial genes identified through binSpect. Additionally, we'll illustrate how to identify the genes that exhibit the highest spatial correlation. + +```{r spatial-cor-genes} +# 3.1 cluster the top 500 spatial genes into 20 clusters +ext_spatial_genes = ranktest[1:300,]$feats + +# here we use existing detectSpatialCorGenes function to calculate pairwise distances between genes (but set network_smoothing=0 to use default clustering) +spat_cor_netw_DT = detectSpatialCorFeats( + brainG, + method = 'network', + spatial_network_name = 'Delaunay_network', + subset_feats = ext_spatial_genes +) +# 3.2 identify most similar spatially correlated genes for one gene +top10_genes = showSpatialCorFeats(spat_cor_netw_DT, feats = 'Dsp', show_top_feats = 10) +``` + +Employing the pheatmap function, we can visualize spatial co-expression modules and iteratively adjust the number of clusters (k) for optimal representation. + +```{r coexpression-modules} +# 3.3 identify potenial spatial co-expression +spat_cor_netw_DT = clusterSpatialCorFeats(spat_cor_netw_DT, name = 'spat_netw_clus', k = 7) + +# visualize clusters +graph3 <- heatmSpatialCorFeats( + brainG, + spatCorObject = spat_cor_netw_DT, + use_clus_name = 'spat_netw_clus', + heatmap_legend_param = list(title = NULL), + save_param = list(base_height = 6, base_width = 8, units = 'cm'), + show_plot = TRUE, + return_plot = FALSE, + save_plot = FALSE +) +``` + +We'll proceed by isolating genes from individual co-expression modules and then consolidate them into meta-genes through aggregation. + +```{r balanced-spat-coexpression-feats, message=FALSE, include=FALSE} +# 3.4 create metagenes / co-expression modules +cluster_genes = getBalancedSpatCoexpressionFeats(spat_cor_netw_DT, maximum = 30) +``` + +```{r create-meta-feats, fig.width=10, fig.height=8} +brainG = createMetafeats(brainG, feat_clusters = cluster_genes, name = 'cluster_metagene') + + +graph4 <- spatCellPlot( + brainG, + spat_enr_names = 'cluster_metagene', + cell_annotation_values = as.character(c(1:7)), + point_size = 1, cow_n_col = 3, gradient_style = 'sequential' +) +``` + + +## Conversion of Giotto to Seurat V5 + +With Seurat's recent upgrade to version 5, ensuring compatibility is essential. The 'giottoToSeuratV5()' function simplifies the process by seamlessly converting Giotto objects to the latest Seurat object. + +```{r message=FALSE, warning=FALSE, include=FALSE} +giottoToSeuratV5 <- 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 + spat_unit <- set_default_spat_unit( + gobject = gobject, + spat_unit = spat_unit + ) + + # verify if optional package is installed + package_check(pkg_name = "Seurat", repository = "CRAN") + requireNamespace("Seurat") + + # check whether any raw data exist -- required for Seurat + avail_expr <- list_expression(gobject = gobject, spat_unit = spat_unit) + raw_exist <- avail_expr[, "raw" %in% name, by = feat_type] + # raw_exist <- sapply(gobject@expression_feat,function(x) + # 'raw' %in% names(gobject@expression[[x]])) + if (nrow(raw_exist) > 0) { + assays_all <- raw_exist[, feat_type] + # assays_all <- names(raw_exist[1]) + # assays_all <- union(assays_all,gobject@expression_feat) + } else { + stop("Raw count data not found. Required for Seurat object.") + } + + # create Seurat object when at least one raw data is available + for (i in seq_along(assays_all)) { + assay_use <- assays_all[i] + expr_use <- lapply( + avail_expr[feat_type == assay_use, name], + function(x) { + get_expression_values( + gobject = gobject, + spat_unit = spat_unit, + feat_type = assay_use, + values = x, + output = "exprObj" + ) + } + ) + # expr_use <- gobject@expression[[assay_use]] + names(expr_use) <- unlist(lapply(expr_use, objName)) + slot_use <- names(expr_use) + if (i == 1) { + data_raw <- expr_use[["raw"]][] + sobj <- Seurat::CreateSeuratObject( + counts = data_raw, + assay = assay_use + ) + if ("normalized" %in% slot_use) { + sobj <- Seurat::SetAssayData(sobj, + slot = "data", + new.data = expr_use[["normalized"]][], + assay = assay_use + ) + } + if ("scaled" %in% slot_use) { + sobj <- Seurat::SetAssayData(sobj, + slot = "scale.data", + # does not accept 'dgeMatrix' + new.data = as.matrix(expr_use[["scaled"]][]), + assay = assay_use + ) + } + } else { + if ("raw" %in% slot_use) { + data_raw <- expr_use[["raw"]][] + flag_raw <- 1 + } else { + flag_raw <- 0 + } + if ("normalized" %in% slot_use) { + data_norm <- expr_use[["normalized"]][] + flag_norm <- 1 + } else { + flag_norm <- 0 + } + if (flag_raw == 1) { + assay_obj <- Seurat::CreateAssayObject(counts = data_raw) + } else if (flag_raw == 0 & flag_norm == 1) { + assay_obj <- Seurat::CreateAssayObject(data = data_norm) + } else { + stop(paste0("Raw and normalized data not found for assay ", assay_use)) + } + sobj[[assay_use]] <- assay_obj + if ("scaled" %in% slot_use) { + data_scale <- as.matrix(expr_use[["scaled"]][]) + sobj <- Seurat::SetAssayData(sobj, + slot = "scale.data", + new.data = data_scale, + assay = assay_use + ) + } + } + + # add cell metadata + meta_cells <- data.table::setDF( + getCellMetadata( + gobject = gobject, + spat_unit = spat_unit, + feat_type = assay_use, + output = "data.table", + copy_obj = TRUE + ) + ) + rownames(meta_cells) <- meta_cells$cell_ID + meta_cells <- meta_cells[, -which(colnames(meta_cells) == "cell_ID")] + if (ncol(meta_cells) > 0) { + colnames(meta_cells) <- paste0(assay_use, "_", colnames(meta_cells)) + } + sobj <- Seurat::AddMetaData(sobj, metadata = meta_cells[Seurat::Cells(sobj), ]) + + # add feature metadata + meta_genes <- data.table::setDF( + get_feature_metadata( + gobject = gobject, + spat_unit = spat_unit, + feat_type = assay_use, + output = "data.table", + copy_obj = TRUE + ) + ) + rownames(meta_genes) <- meta_genes$feat_ID + 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 + avail_dr <- list_dim_reductions(gobject = gobject, spat_unit = spat_unit, feat_type = assay_use, ) + if (!is.null(avail_dr)) { + if (nrow(avail_dr) > 0) { + dr_use <- avail_dr[, name] + for (i in seq(nrow(avail_dr))) { + dr_name <- avail_dr[i, name] + dr_type <- avail_dr[i, dim_type] + dr_obj <- get_dimReduction( + gobject = gobject, + output = "dimObj", + spat_unit = spat_unit, + feat_type = assay_use, + reduction_method = dr_type, + name = dr_name + ) + emb_use <- dr_obj[][Seurat::Cells(sobj), ] + if (sum(c("loadings", "eigenvalues") %in% names(slot(dr_obj, "misc"))) == 2) { + loadings_use <- slot(dr_obj, "misc")$loadings + stdev_use <- slot(dr_obj, "misc")$eigenvalues + sobj[[dr_name]] <- Seurat::CreateDimReducObject( + embeddings = as.matrix(emb_use), + loadings = loadings_use, + key = paste0(dr_name, "_"), + stdev = stdev_use, + assay = assay_use + ) + } else { + sobj[[dr_name]] <- Seurat::CreateDimReducObject( + embeddings = as.matrix(emb_use), + key = paste0(dr_name, "_"), + assay = assay_use + ) + } + } + } + } + + + + # network objects + # expression network + avail_nn <- list_nearest_networks(gobject = gobject, spat_unit = spat_unit, feat_type = assay_use) + if (!is.null(avail_nn)) { + if (nrow(avail_nn) > 0) { + for (i in seq(nrow(avail_nn))) { + nn_name <- avail_nn[i, name] + nn_type <- avail_nn[i, nn_type] + nn_use <- get_NearestNetwork( + gobject = gobject, + spat_unit = spat_unit, + feat_type = assay_use, + nn_network_to_use = nn_type, + network_name = nn_name, + output = "data.table" + ) + idx1 <- match(nn_use$from, Seurat::Cells(sobj)) + idx2 <- match(nn_use$to, Seurat::Cells(sobj)) + edge_weight <- nn_use$weight + edge_dist <- nn_use$distance + nn_mtx <- Matrix::sparseMatrix(i = idx1, j = idx2, x = edge_dist, dims = c(ncol(sobj), ncol(sobj))) + rownames(nn_mtx) <- colnames(nn_mtx) <- Seurat::Cells(sobj) + nn_name <- paste0( nn_name) + sGraph <- Seurat::as.Graph(nn_mtx) + sGraph@assay.used <- assay_use + sobj[[nn_name]] <- sGraph + } + } + } + } + + # spatial coordinates + 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 + ) + ) + rownames(loc_use) <- loc_use$cell_ID + sobj <- Seurat::AddMetaData(sobj, metadata = loc_use) + # add spatial coordinates as new dim reduct object + loc_2 <- loc_use[, c("sdimx", "sdimy")] + colnames(loc_2) <- c("spatial_1", "spatial_2") + sobj[["spatial"]] <- Seurat::CreateDimReducObject( + embeddings = as.matrix(loc_2), + assay = names(sobj@assays)[1], + key = "spatial_" + ) + + + + # spatial network + avail_sn <- list_spatial_networks(gobject = gobject, spat_unit = spat_unit) + if (!is.null(avail_nn)) { + if (nrow(avail_sn) > 0) { + sn_all <- avail_sn[, name] + for (i in sn_all) { + snt_use <- get_spatialNetwork( + gobject = gobject, + spat_unit = spat_unit, + name = i, + output = "networkDT" + ) + idx1 <- match(snt_use$from, Seurat::Cells(sobj)) + idx2 <- match(snt_use$to, Seurat::Cells(sobj)) + edge_weight <- snt_use$weight + edge_dist <- snt_use$distance + nn_mtx <- Matrix::sparseMatrix(i = idx1, j = idx2, x = edge_dist, dims = c(ncol(sobj), ncol(sobj))) + rownames(nn_mtx) <- colnames(nn_mtx) <- Seurat::Cells(sobj) + nn_name <- paste0(i) + spatGraph <- Seurat::as.Graph(nn_mtx) + spatGraph@assay.used <- names(sobj@assays)[1] + sobj[[nn_name]] <- spatGraph + } + } + } + + all_x <- NULL + all_y <- NULL + + if(length(gobject@largeImages) > 0){ + for (i in seq_along(gobject@largeImages)){ + imagerow <- gobject@spatial_locs$cell$raw$sdimx + imagecol <- gobject@spatial_locs$cell$raw$sdimy + img <- terra::as.array(gobject@largeImages[[i]]@raster_object) + img[, , 1:3] <- img[, , 1:3] / 255 + coord <- data.frame(imagerow = imagerow, imagecol = imagecol) + rownames(coord) <- gobject@cell_ID$cell + + scalefactors <- Seurat::scalefactors(spot = gobject@largeImages[[i]]@scale_factor[[1]] , + fiducial = gobject@largeImages[[i]]@resolution[[1]], + hires = max(img), + lowres = gobject@largeImages[[i]]@resolution[[1]] + ) + + senv <- rlang::ns_env("Seurat") + + newV1 <- do.call( + new, + args = list( + Class = "VisiumV1", + image = img, + scale.factors = scalefactors, + coordinates = coord + ), + envir = senv + ) + + sobj@images[[gobject@largeImages[[i]]@name]] <- newV1 + ikey <- paste0(gobject@largeImages[[i]]@name, "_") + sobj@images[[i]]@key <- ikey + sobj@images[[i]]@assay <- GiottoClass::activeFeatType(gobject) + sobj@images[[i]]@spot.radius <- 0.012439 + + } + + } + + return(sobj) +} +``` + +### Reverting GiottoObj to Seurat + +In this step, we revert the Giotto object, previously converted from Seurat, back to its original Seurat format: + +```{r echo=TRUE, message=FALSE, warning=FALSE} +brain2 <- giottoToSeuratV5(brainG) +``` + +### Print Seurat Object + +```{r} +brain2 +``` + +### Downstream Analysis of Giotto to Seurat V5: + +Following the conversion of Giotto to Seurat V5, downstream analysis can be conducted as outlined below: + +#### Load libraries + +Load Required Libraries + +```{r echo=FALSE, message=FALSE, warning=FALSE} +library(Seurat) +library(SeuratData) +library(ggplot2) +library(patchwork) +library(dplyr) +``` + +#### Visualization of Cell Counts Distribution + +Visualize the distribution of cell counts using violin plots and spatial feature plots: + +```{r} +plot1 <- VlnPlot(brain2, features = "nCount_rna", pt.size = 0.1) + NoLegend() + +plot2 <- SpatialFeaturePlot(brain2, features = "nCount_rna") + theme(legend.position = "right") +wrap_plots(plot1, plot2) + +``` + +#### Data Transformation + +To enhance the analysis, we apply SCTransform to perform data transformation on the RNA assay: + +```{r, warning=FALSE} +brain2 <- SCTransform(brain2, assay = "rna", verbose = FALSE) +``` + +#### Dimensionality Reduction and Clustering + +Following data transformation, we proceed with dimensionality reduction and clustering. This involves conducting PCA, identifying neighbors, performing clustering, and running UMAP for dimensionality reduction: + +```{r message=FALSE, warning=FALSE} +brain2 <- RunPCA(brain2, assay = "SCT", verbose = FALSE) +brain2 <- FindNeighbors(brain2, reduction = "pca", dims = 1:30) +brain2 <- FindClusters(brain2, verbose = FALSE) +brain2 <- RunUMAP(brain2, reduction = "pca", dims = 1:30) +``` + +#### Visualization of Clusters on UMAP + +To visualize clusters on UMAP, we generate dimensional plots with cluster labels for better interpretation: + +```{r message=FALSE, warning=FALSE} +p1 <- DimPlot(brain2, reduction = "umap", label = TRUE) +p2 <- SpatialDimPlot(brain2, label = TRUE, label.size = 3) +p1 + p2 +``` + +Additionally, we can highlight specific cell populations within the spatial dimensional plot: + +```{r message=FALSE, warning=FALSE} +SpatialDimPlot(brain2, cells.highlight = CellsByIdentities(object = brain2, idents = c(2, 1, 4, 3, + 5, 8)), facet.highlight = TRUE, ncol = 3) +``` + +#### Identification of Spatially Variable Features + +To identify spatially variable features between clusters, the FindMarkers function is employed. Subsequently, we visualize the top spatially variable features using SpatialFeaturePlot: + +```{r message=FALSE, warning=FALSE} +de_markers <- FindMarkers(brain2, ident.1 = 5, ident.2 = 6) +SpatialFeaturePlot(object = brain2, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3) +``` + +we identify spatially variable features within the spatial data. This is achieved through the `FindSpatiallyVariableFeatures` function, which calculates spatial statistics to determine the variability of features across spatial locations. In the following code chunk, we focus on the top 1000 variable features selected from the SCT assay using the Moransi selection method. + +```{r message=FALSE, warning=FALSE} +brain2 <- FindSpatiallyVariableFeatures(brain2, assay = "SCT", features = VariableFeatures(brain2)[1:1000], + selection.method = "moransi") +``` +#### Subset out anatomical regions + +This section focuses on subsetting anatomical regions from the dataset and visualizing them using spatial dimensional plots. + +```{r message=FALSE, warning=FALSE} +cortex <- subset(brain2, idents = c(1, 2, 3, 4, 6, 7)) +cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert = TRUE) +cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert = TRUE) +cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert = TRUE) +p1 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE) +p2 <- SpatialDimPlot(cortex, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3) +p1 + p2 +``` + +#### Working with multiple slices in Seurat + +To incorporate data from multiple slices in our analysis, we first load and transform the posterior slice: + +```{r message=FALSE, warning=FALSE} +brainP <- LoadData("stxBrain", type = "posterior1") +brainP <- SCTransform(brainP, assay = "Spatial", verbose = FALSE) +``` + +Next, we merge the datasets obtained from the anterior and posterior slices: + +```{r message=FALSE, warning=FALSE} + +brain.merge <- merge(brain2, brainP) +``` + +We set the default assay to "SCT" and combine variable features from both datasets: + +```{r message=FALSE, warning=FALSE} +DefaultAssay(brain.merge) <- "SCT" +VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2)) +brain.merge <- RunPCA(brain.merge, verbose = FALSE) +brain.merge <- FindNeighbors(brain.merge, dims = 1:30) +brain.merge <- FindClusters(brain.merge, verbose = FALSE) +brain.merge <- RunUMAP(brain.merge, dims = 1:30) +``` + +After performing PCA, identifying neighbors, clustering, and running UMAP, we visualize the results: + +```{r message=FALSE, warning=FALSE} +DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident")) +``` + +```{r message=FALSE, warning=FALSE} +SpatialDimPlot(brain.merge) +``` + +```{r message=FALSE, warning=FALSE} +SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1")) +``` + +```{r} +sessionInfo() +``` + diff --git a/NAMESPACE b/NAMESPACE index 41d345ea..599b9fe6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -372,9 +372,6 @@ import(sp) import(utils) importClassesFrom(terra,SpatExtent) importClassesFrom(terra,SpatVector) -importFrom(GiottoUtils,getDistinctColors) -importFrom(GiottoUtils,getMonochromeColors) -importFrom(GiottoUtils,getRainbowColors) importFrom(checkmate,assert_character) importFrom(grDevices,dev.size) importFrom(graphics,legend) diff --git a/R/interoperability.R b/R/interoperability.R index 497966a4..6d5363a3 100644 --- a/R/interoperability.R +++ b/R/interoperability.R @@ -2206,7 +2206,10 @@ seuratToGiottoV5 <- function( object = sobject, assay = spatial_assay ))) { - spat_coord <- Seurat::GetTissueCoordinates(sobject) + spat_coord <- Seurat::GetTissueCoordinates(sobject, + scale = NULL, + cols = c("imagerow", + "imagecol")) # spat_coord = cbind(rownames(spat_coord), # data.frame(spat_coord, row.names=NULL)) @@ -2336,7 +2339,10 @@ seuratToGiottoV5 <- function( # Create Giotto LargeImage gImg <- createGiottoLargeImage( raster_object = terra::rast(list(r, g, b)), - name = names(sobject@images) + name = names(sobject@images), + scale_factor = sobject@images[[ + i + ]]@scale.factors$lowres ) } } diff --git a/man/aggregateStacksExpression.Rd b/man/aggregateStacksExpression.Rd index 7d7ab64a..eed57c81 100644 --- a/man/aggregateStacksExpression.Rd +++ b/man/aggregateStacksExpression.Rd @@ -42,9 +42,9 @@ aggregateStacksExpression(g, spat_units = c("z0", "z1"), feat_type = "rna") } \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 7fcad4a2..84d0dd78 100644 --- a/man/aggregateStacksLocations.Rd +++ b/man/aggregateStacksLocations.Rd @@ -36,9 +36,9 @@ aggregateStacksLocations(g, spat_units = c("z0", "z1")) } \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 a2848c0a..df3900c5 100644 --- a/man/aggregateStacksPolygonOverlaps.Rd +++ b/man/aggregateStacksPolygonOverlaps.Rd @@ -35,9 +35,9 @@ aggregateStacksPolygonOverlaps(g, } \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 1b9e972d..0332292d 100644 --- a/man/aggregateStacksPolygons.Rd +++ b/man/aggregateStacksPolygons.Rd @@ -36,9 +36,9 @@ aggregateStacksPolygons(g, spat_units = c("z0", "z1")) } \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}