Merge pull request #133 from BIMSBbioinfo/main

main -> dev

NAMESPACE

@@ -160,6 +160,7 @@ export(importCosMx)
 export(importGenePS)
 export(importGeoMx)
 export(importImageData)
+export(importPhenoCycler)
 export(importSTOmics)
 export(importVisium)
 export(importXenium)

R/import.R

@@ -1358,6 +1358,318 @@ importSTOmics <- function(h5ad.path, assay_name = "STOmics", sample_name = NULL,
   formVoltRon(rawdata, metadata = metadata, coords = coords, main.assay = assay_name, params = params, assay.type = "spot", image_name = image_name, main_channel = channel_name, sample_name = sample_name, ...)
 }
 
+####
+# Akoya ####
+####
+
+####
+## PhenoCycler ####
+####
+
+#' readPhenoCyclerMat
+#' 
+#' Read and Load Akoya CODEX data, adapted from \code{ReadAkoya} function \code{Seurat} package
+#'
+#' @param filename Path to matrix generated by upstream processing.
+#' @param type Specify which type matrix is being provided.
+#' \itemize{
+#'  \item \dQuote{\code{processor}}: matrix generated by CODEX Processor
+#'  \item \dQuote{\code{inform}}: matrix generated by inForm
+#'  \item \dQuote{\code{qupath}}: matrix generated by QuPath
+#' }
+#' @param filter A pattern to filter features by; pass \code{NA} to skip feature filtering
+#' @param inform.quant When \code{type} is \dQuote{\code{inform}}, the quantification level to read in
+#' 
+#' @noRd
+readPhenoCyclerMat <- function(
+    filename,
+    type = c('inform', 'processor', 'qupath'),
+    filter = 'DAPI|Blank|Empty',
+    inform.quant = c('mean', 'total', 'min', 'max', 'std')
+) {
+  if (!requireNamespace("data.table", quietly = TRUE)) {
+    stop("Please install 'data.table' for this function")
+  }
+  # Check arguments
+  if (!file.exists(filename)) {
+    stop(paste("Can't file file:", filename))
+  }
+  type <- tolower(x = type[1L])
+  type <- match.arg(arg = type)
+  ratio <- getOption(x = 'Seurat.input.sparse_ratio', default = 0.4)
+  # Preload matrix
+  sep <- switch(EXPR = type, 'inform' = '\t', ',')
+  mtx <- data.table::fread(
+    file = filename,
+    sep = sep,
+    data.table = FALSE,
+    verbose = FALSE
+  )
+  # Assemble outputs
+  outs <- switch(
+    EXPR = type,
+    'processor' = {
+      # Create centroids data frame
+      centroids <- data.frame(
+        x = mtx[['x:x']],
+        y = mtx[['y:y']],
+        cell = as.character(x = mtx[['cell_id:cell_id']]),
+        stringsAsFactors = FALSE
+      )
+      rownames(x = mtx) <- as.character(x = mtx[['cell_id:cell_id']])
+      # Create metadata data frame
+      md <- mtx[, !grepl(pattern = '^cyc', x = colnames(x = mtx)), drop = FALSE]
+      colnames(x = md) <- vapply(
+        X = strsplit(x = colnames(x = md), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        2L
+      )
+      # Create expression matrix
+      mtx <- mtx[, grepl(pattern = '^cyc', x = colnames(x = mtx)), drop = FALSE]
+      colnames(x = mtx) <- vapply(
+        X = strsplit(x = colnames(x = mtx), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        2L
+      )
+      if (!is.na(x = filter)) {
+        mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)), drop = FALSE]
+      }
+      mtx <- t(x = mtx)
+      if ((sum(mtx == 0) / length(x = mtx)) > ratio) {
+        mtx <- as.sparse(x = mtx)
+      }
+      list(matrix = mtx, centroids = centroids, metadata = md)
+    },
+    'inform' = {
+      inform.quant <- tolower(x = inform.quant[1L])
+      inform.quant <- match.arg(arg = inform.quant)
+      expr.key <- c(
+        mean = 'Mean',
+        total = 'Total',
+        min = 'Min',
+        max = 'Max',
+        std = 'Std Dev'
+      )[inform.quant]
+      expr.pattern <- '\\(Normalized Counts, Total Weighting\\)'
+      rownames(x = mtx) <- mtx[['Cell ID']]
+      mtx <- mtx[, setdiff(x = colnames(x = mtx), y = 'Cell ID'), drop = FALSE]
+      # Create centroids
+      centroids <- data.frame(
+        x = mtx[['Cell X Position']],
+        y = mtx[['Cell Y Position']],
+        cell  = rownames(x = mtx),
+        stringsAsFactors = FALSE
+      )
+      # Create metadata
+      cols <- setdiff(
+        x = grep(
+          pattern = expr.pattern,
+          x = colnames(x = mtx),
+          value = TRUE,
+          invert = TRUE
+        ),
+        y = paste('Cell', c('X', 'Y'), 'Position')
+      )
+      md <- mtx[, cols, drop = FALSE]
+      # Create expression matrices
+      exprs <- data.frame(
+        cols = grep(
+          pattern = paste(expr.key, expr.pattern),
+          x = colnames(x = mtx),
+          value = TRUE
+        )
+      )
+      exprs$feature <- vapply(
+        X = trimws(x = gsub(
+          pattern = paste(expr.key, expr.pattern),
+          replacement = '',
+          x = exprs$cols
+        )),
+        FUN = function(x) {
+          x <- unlist(x = strsplit(x = x, split = ' '))
+          x <- x[length(x = x)]
+          return(gsub(pattern = '\\(|\\)', replacement = '', x = x))
+        },
+        FUN.VALUE = character(length = 1L)
+      )
+      exprs$class <- tolower(x = vapply(
+        X = strsplit(x = exprs$cols, split = ' '),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        1L
+      ))
+      classes <- unique(x = exprs$class)
+      outs <- vector(
+        mode = 'list',
+        length = length(x = classes) + 2L
+      )
+      names(x = outs) <- c(
+        'matrix',
+        'centroids',
+        'metadata',
+        setdiff(x = classes, y = 'entire')
+      )
+      outs$centroids <- centroids
+      outs$metadata <- md
+      # browser()
+      for (i in classes) {
+        df <- exprs[exprs$class == i, , drop = FALSE]
+        expr <- mtx[, df$cols]
+        colnames(x = expr) <- df$feature
+        if (!is.na(x = filter)) {
+          expr <- expr[, !grepl(pattern = filter, x = colnames(x = expr)), drop = FALSE]
+        }
+        expr <- t(x = expr)
+        if ((sum(expr == 0, na.rm = TRUE) / length(x = expr)) > ratio) {
+          expr <- as.sparse(x = expr)
+        }
+        outs[[switch(EXPR = i, 'entire' = 'matrix', i)]] <- expr
+      }
+      outs
+    },
+    'qupath' = {
+      rownames(x = mtx) <- as.character(x = seq_len(length.out = nrow(x = mtx)))
+      # Create centroids
+      xpos <- sort(
+        x = grep(pattern = 'Centroid X', x = colnames(x = mtx), value = TRUE),
+        decreasing = TRUE
+      )[1L]
+      ypos <- sort(
+        x = grep(pattern = 'Centroid Y', x = colnames(x = mtx), value = TRUE),
+        decreasing = TRUE
+      )[1L]
+      centroids <- data.frame(
+        x = mtx[[xpos]],
+        y = mtx[[ypos]],
+        cell = rownames(x = mtx),
+        stringsAsFactors = FALSE
+      )
+      # Create metadata
+      cols <- setdiff(
+        x = grep(
+          pattern = 'Cell: Mean',
+          x = colnames(x = mtx),
+          ignore.case = TRUE,
+          value = TRUE,
+          invert = TRUE
+        ),
+        y = c(xpos, ypos)
+      )
+      md <- mtx[, cols, drop = FALSE]
+      # Create expression matrix
+      idx <- which(x = grepl(
+        pattern = 'Cell: Mean',
+        x = colnames(x = mtx),
+        ignore.case = TRUE
+      ))
+      mtx <- mtx[, idx, drop = FALSE]
+      colnames(x = mtx) <- vapply(
+        X = strsplit(x = colnames(x = mtx), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        1L
+      )
+      if (!is.na(x = filter)) {
+        mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)), drop = FALSE]
+      }
+      mtx <- t(x = mtx)
+      if ((sum(mtx == 0) / length(x = mtx)) > ratio) {
+        mtx <- as.sparse(x = mtx)
+      }
+      list(matrix = mtx, centroids = centroids, metadata = md)
+    },
+    stop("Unknown matrix type: ", type)
+  )
+  return(outs)
+}
+
+
+#' Title
+#'
+#' @param dir.path path to PhenoCycler output folder
+#' @param assay_name the assay name of the SR object
+#' @param sample_name the name of the sample
+#' @param image_name the image name of the Xenium assay, Default: main
+#' @param filename Path to matrix generated by upstream processing.
+#' @param type Specify which type matrix is being provided.
+#' \itemize{
+#'  \item \dQuote{\code{processor}}: matrix generated by CODEX Processor
+#'  \item \dQuote{\code{inform}}: matrix generated by inForm
+#'  \item \dQuote{\code{qupath}}: matrix generated by QuPath
+#' }
+#' @param filter A pattern to filter features by; pass \code{NA} to skip feature filtering
+#' @param inform.quant When \code{type} is \dQuote{\code{inform}}, the quantification level to read in 
+#' @param ... additional parameters passed to \link{formVoltRon} 
+#' 
+#' @importFrom magick image_info image_read
+#'
+#' @export
+importPhenoCycler <- function(dir.path, assay_name = "PhenoCycler", sample_name = NULL, image_name = "main", 
+                              type = c('inform', 'processor', 'qupath'), filter = 'DAPI|Blank|Empty', inform.quant = c('mean', 'total', 'min', 'max', 'std'), ...){
+
+  # raw counts, metadata and coordinates 
+  listoffiles <- list.files(paste0(dir.path, "/processed/segm/segm-1/fcs/compensated/"), full.names = TRUE)
+  datafile <- listoffiles[grepl("_compensated.csv$", listoffiles)][1]
+  if(file.exists(datafile)){
+    rawdata <- readPhenoCyclerMat(filename = datafile, type = type, filter = filter, inform.quant = inform.quant)
+  } else {
+    stop("There are no files named ending with '_compensated.csv' in the processed/segm/segm-1/fcs/compensated/ subfolder")
+  }
+
+  # cell id 
+  cellid <- paste0("cell", colnames(rawdata$matrix))
+
+  # coordinates
+  coords <- rawdata$centroids[,c("x", "y")]
+  rownames(coords) <- cellid
+  coords <- as.matrix(coords)
+
+  # metadata
+  metadata <- rawdata$metadata
+  rownames(metadata) <- cellid
+
+  # data
+  rawdata <- rawdata$matrix
+  colnames(rawdata) <- cellid
+  rownames(rawdata) <- gsub("(\t|\r|\n)", "", rownames(rawdata))
+
+  # images 
+  image_dir <- paste0(dir.path, "/processed/stitched/reg001/")
+  list_files <- list.files(image_dir)
+  if(!dir.exists(image_dir)){
+    message("There are no images of channels!")
+    image_list <- NULL
+  } else {
+    if(!any(grepl(".tif$", list_files))){
+      stop("The folder doesnt have any images associated with channels!")
+    } else{
+      image_channel_names <- sapply(list_files, function(x) {
+        name <- strsplit(x, split = "_")[[1]]
+        name <- gsub(".tif$", "", name[length(name)])
+        name <- make.unique(name)
+        return(name)
+      })
+      image_channel_names <- c(image_channel_names[grepl("DAPI", image_channel_names)][1], 
+                               image_channel_names[image_channel_names %in% rownames(rawdata)])
+      list_files <- names(image_channel_names)
+      image_list <- lapply(list_files, function(x){
+        magick::image_read(paste0(image_dir, "/", x))
+      })
+      names(image_list) <- image_channel_names   
+      coords[,2] <- magick::image_info(image_list[[1]])$height - coords[,2]
+    }
+  }
+
+  # voltron object
+  object <- formVoltRon(data = rawdata, metadata = metadata, image = image_list, coords = coords, assay.type = "cell", 
+                        sample_name = sample_name, main.assay = assay_name, image_name = image_name, ...)
+
+  # return
+  object
+}
+
 ####
 # Image Data ####
 ####

R/visualization.R

@@ -1253,9 +1253,11 @@ vrEmbeddingFeaturePlot <- function(object, embedding = "pca", features = NULL, n
 
     # get data
     if(feat %in% vrFeatures(object)){
-      datax$score <- normdata[feat,]
+      # datax$score <- normdata[feat,]
+      datax$score <- normdata[feat, rownames(datax)]
     } else {
-      datax$score <- metadata[,feat]
+      # datax$score <- metadata[,feat]
+      datax$score <- metadata[rownames(datax),feat]
     }
 
     # get image information and plotting features