Feature: Updates for Space Ranger 4.0

R/preprocessing.R

@@ -498,6 +498,7 @@ GetResidual <- function(
 #' @param filename Name of H5 file containing the feature barcode matrix
 #' @param slice Name for the stored image of the tissue slice
 #' @param bin.size Specifies the bin sizes to read in - defaults to c(16, 8)
+#' @param load.segmentations Should the cell segmentation from space ranger be read in?
 #' @param filter.matrix Only keep spots that have been determined to be over
 #' tissue
 #' @param to.upper Converts all feature names to upper case. Can be useful when
@@ -527,6 +528,7 @@ Load10X_Spatial <- function (
   assay = "Spatial",
   slice = "slice1",
   bin.size = NULL,
+  load.segmentations = FALSE,
   filter.matrix = TRUE,
   to.upper = FALSE,
   image = NULL,
@@ -640,6 +642,22 @@ Load10X_Spatial <- function (
     y = object.list[-1]
   )
 
+  # Possibly make this into a new function
+  if(load.segmentations) {
+    segmentation.assay.name <- "Segmentations"
+    seg.data.dir <- paste0(data.dir, "/segmented_outputs/")
+    seg.counts.path <- paste0(seg.data.dir, "filtered_feature_cell_matrix.h5")
+    segmentation.counts <- Read10X_h5(seg.counts.path)
+    segmentation.count.cellIDs <- colnames(segmentation.counts)
+    #segmentation.image <- Read10X_Image(image.dir = paste0(seg.data.dir, "spatial"), assay = "Segmentations", slice = "slice1.segmentations")
+    segmentation.object <- CreateSeuratObject(segmentation.counts, assay = segmentation.assay.name)
+    #segmentation.object[["slice1.segmentation.image"]] <- segmentation.image
+    segmentation.fov <- Read10X_Segmentations(image.dir = paste0(seg.data.dir, "spatial"), data.dir = data.dir, cell.names = segmentation.count.cellIDs)
+    segmentation.object[["slice1.polygons"]] <- segmentation.fov
+    object <- merge(x = object, y = segmentation.object)
+    DefaultAssay(object = object) <- segmentation.assay.name
+  }
+
   return(object)
 }
 
@@ -1263,8 +1281,8 @@ Read10X_Image <- function(
       image = image
     )
 
-    # As of v5.1.0 `Radius.VisiumV1` no longer returns the value of the 
-    # `spot.radius` slot and instead calculates the value on the fly, but we 
+    # As of v5.1.0 `Radius.VisiumV1` no longer returns the value of the
+    # `spot.radius` slot and instead calculates the value on the fly, but we
     # can populate the static slot in case it's depended on.
     visium.v1@spot.radius <- Radius(visium.v1)
 
@@ -1390,6 +1408,97 @@ Read10X_ScaleFactors <- function(filename) {
   return (scale.factors)
 }
 
+#' Load 10X Genomics Visium Cell Segmentations
+#'
+#' @param image.dir Path to directory with 10X Genomics visium image data;
+#' @param data.dir Directory of the base spaceranger outs
+#' @param image.name Name of the tissue image to be plotted. tissue_lowres_image.png or tissue_hires_image.png
+#' @param assay Name of assay to associate segmentations to
+#' @param slice Name of the slice to associate the segmentations to
+#' @param segmentation.type Which segmentations to load, cell or nucleus. If using nucleus the full matrix from cells is still used
+#' @param cell.names Names of the cells from the matrix.h5
+#'
+#'
+#' @return A FOV
+#'
+#' @export
+#' @concept preprocessing
+#'
+Read10X_Segmentations <- function (image.dir,
+                                   data.dir,
+                                   image.name = "tissue_lowres_image.png",
+                                   assay = "Segmentations",
+                                   slice = "slice1.segmentations",
+                                   segmentation.type = "cell",
+                                   cell.names
+)
+{
+
+  image <- png::readPNG(source = file.path(image.dir, image.name))
+  scale.factors <- Read10X_ScaleFactors(filename = file.path(image.dir, "scalefactors_json.json"))
+  coordinates <- Read10X_GeoJson(data.dir = data.dir, segmentation.type = segmentation.type)
+  key <- Key(slice, quiet = TRUE)
+
+  st_coords <- merge(
+    as.data.frame(sf::st_coordinates(coordinates)),
+    data.frame(
+      cell_id = coordinates$cell_id,
+      cell = coordinates$barcodes
+    ),
+    by.x = 'L2',
+    by.y = 'cell_id'
+  )
+  st_coords <- st_coords[st_coords$cell %in% cell.names, ]
+  image.scale <- ifelse(grepl("lowres", image.name, ignore.case = TRUE), "lowres", "highres")
+  fov <- CreateFOV(st_coords,
+                   type = "segmentation",
+                   #radius = scale.factors[[image.scale]],
+                   assay = assay, key = key)
+  return(fov)
+}
+
+#' Format 10X Genomics GeoJson cell IDs
+#'
+#' A helper function to format cell IDs from the segmentation GeoJson to the same type as in the matrix.h5
+#' The GeoJson has cell IDs as integers (eg 1). They need to be in the format cellid_000000001-1
+#'
+#' @return Vector of formatted cell IDs
+Format10X_GeoJson_CellID <- function(ids, prefix = "cellid_", suffix = "-1", digits = 9) {
+  format_string <- paste0("%0", as.integer(digits), "d")
+
+  formatted_ids <- sapply(ids, function(id) {
+    numeric_part <- sprintf(format_string, as.integer(id))
+    paste0(prefix, numeric_part, suffix)
+  })
+
+  return(formatted_ids)
+}
+
+#' Load 10X Genomics GeoJson
+#'
+#' @param filename Path to a \code{cell_segmentations.geojson} file
+#' @param segmentation.type Which segmentations to load, cell or nucleus. If using nucleus the full matrix from cells is still used
+#' @param scale.factor If scaling the segmentations coordinates for the associated tissue image. "lowres" or "highres"
+#'
+#' @return A FOV
+#'
+#' @export
+#' @concept preprocessing
+#'
+Read10X_GeoJson <- function(data.dir, segmentation.type = "cell", scale.factor = NULL) {
+  segmentation_polygons <- read_sf(file.path(data.dir,"segmented_outputs", paste0(segmentation.type, "_segmentations.geojson")))
+  if (!is.null(scale.factor)) {
+    scale.factor.path <- file.path(data.dir,"segmented_outputs/spatial/scalefactors_json.json")
+    scale.factors <- jsonlite::fromJSON(scale.factor.path)
+    segmentation_polygons$geometry <- segmentation_polygons$geometry*scale.factors[[scale.factor]]
+  }
+
+  segmentation_polygons$barcodes <- Format10X_GeoJson_CellID(segmentation_polygons$cell_id)
+  segmentation_polygons
+}
+
+
+
 #' Read and Load Akoya CODEX data
 #'
 #' @param filename Path to matrix generated by upstream processing.
@@ -4557,7 +4666,7 @@ FindSpatiallyVariableFeatures.Seurat <- function(
     verbose = verbose,
     ...
   )
-  
+
   object <- LogSeuratCommand(object)
 
   return(object)