Datasets
Models
Diff of
/R/import.R
[000000]
..
[413088]
Switch to side-by-side view
--- a +++ b/R/import.R @@ -0,0 +1,2174 @@ +#### +# 10X Genomics #### +#### + +#### +## Xenium #### +#### + +#' importXenium +#' +#' Importing Xenium data +#' +#' @param dir.path path to Xenium output folder +#' @param selected_assay selected assay from Xenium +#' @param assay_name the assay name of the SR object +#' @param sample_name the name of the sample +#' @param use_image if TRUE, the DAPI image will be used. +#' @param morphology_image the name of the lowred morphology image. Default: morphology_lowres.tif +#' @param resolution_level the level of resolution within Xenium OME-TIFF image, see \link{generateXeniumImage}. Default: 7 (553x402) +#' @param overwrite_resolution if TRUE, the image "file.name" will be generated again although it exists at "dir.path" +#' @param image_name the image name of the Xenium assay, Default: main +#' @param channel_name the channel name of the image of the Xenium assay, Default: DAPI +#' @param import_molecules if TRUE, molecule assay will be created along with cell assay. +#' @param verbose verbose +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom magick image_read image_info +#' @importFrom utils read.csv +#' @importFrom data.table fread +#' @importFrom ids random_id +#' +#' @export +#' +importXenium <- function (dir.path, selected_assay = "Gene Expression", assay_name = "Xenium", sample_name = NULL, use_image = TRUE, + morphology_image = "morphology_lowres.tif", resolution_level = 7, overwrite_resolution = FALSE, + image_name = "main", channel_name = "DAPI", import_molecules = FALSE, verbose = TRUE, ...) +{ + # cell assay + if(verbose) + message("Creating cell level assay ...") + + # raw counts + datafile <- paste0(dir.path, "/cell_feature_matrix.h5") + if(file.exists(datafile)){ + rawdata <- import10Xh5(filename = datafile) + if(any(names(rawdata) %in% selected_assay)){ + rawdata <- as.matrix(rawdata[[selected_assay]]) + } else { + stop("There are no assays called ", selected_assay, " in the h5 file!") + } + } else { + stop("There are no files named 'cell_feature_matrix.h5' in the path") + } + + # image + if(use_image){ + generateXeniumImage(dir.path, file.name = morphology_image, resolution_level = resolution_level, overwrite_resolution = overwrite_resolution, + verbose = verbose) + image_file <- paste0(dir.path, "/", morphology_image) + if(file.exists(image_file)){ + image <- image_read(image_file) + } else { + stop("There are no spatial image files in the path") + } + # scale the xenium image instructed by 10x Genomics help page + scaleparam <- 0.2125*(2^(resolution_level-1)) + } else { + image <- NULL + } + + # coordinates + coord_file <- paste0(dir.path, "/cells.csv.gz") + if(file.exists(coord_file)){ + Xenium_coords <- utils::read.csv(file = coord_file) + coords <- as.matrix(Xenium_coords[,c("x_centroid", "y_centroid")]) + colnames(coords) <- c("x","y") + if(use_image) { + coords <- coords/scaleparam + imageinfo <- unlist(magick::image_info(image)[c("height")]) + range_coords <- c(0,imageinfo) + } else { + range_coords <- range(coords[,2]) + } + coords[,2] <- range_coords[2] - coords[,2] + range_coords[1] + } else { + stop("There are no file named 'cells.csv.gz' in the path") + } + + # segments + segments_file <- paste0(dir.path, "/cell_boundaries.csv.gz") + if(file.exists(segments_file)){ + segments <- as.data.frame(data.table::fread(segments_file)) + segments <- segments[,c("cell_id", "vertex_x", "vertex_y")] + colnames(segments) <- c("cell_id", "x", "y") + if(use_image){ + segments[,c("x","y")] <- segments[,c("x","y")]/scaleparam + } + segments[,"y"] <- range_coords[2] - segments[,"y"] + range_coords[1] + segments <- segments %>% dplyr::group_split(cell_id) + segments <- as.list(segments) + names(segments) <- rownames(coords) + } else { + stop("There are no file named 'cell_boundaries.csv.gz' in the path") + } + + # create VoltRon object for cells + cell_object <- formVoltRon(rawdata, metadata = NULL, image = image, coords, segments = segments, main.assay = assay_name, + assay.type = "cell", image_name = image_name, main_channel = channel_name, sample_name = sample_name, + feature_name = ifelse(selected_assay == "Gene Expression", "RNA", "main"), ...) + + # molecule assay + if(!import_molecules){ + return(cell_object) + } else { + if(verbose) + message("Creating molecule level assay ...") + # transcripts + transcripts_file <- paste0(dir.path, "/transcripts.csv.gz") + if(!file.exists(transcripts_file)){ + stop("There are no file named 'transcripts.csv.gz' in the path") + } else { + # get subcellur data components + subcellular_data <- data.table::fread(transcripts_file) + subcellular_data <- subcellular_data[,c("transcript_id", colnames(subcellular_data)[!colnames(subcellular_data) %in% "transcript_id"]), with = FALSE] + colnames(subcellular_data)[colnames(subcellular_data)=="transcript_id"] <- "id" + colnames(subcellular_data)[colnames(subcellular_data)=="feature_name"] <- "gene" + subcellular_data <- subcellular_data[subcellular_data$qv >= 20, ] + + # coordinates + coords <- as.matrix(subcellular_data[,c("x_location", "y_location", "z_location")]) + colnames(coords) <- c("x", "y", "z") + if(use_image){ + coords[,c("x", "y")] <- coords[,c("x", "y")]/scaleparam + } + coords[,"y"] <- range_coords[2] - coords[,"y"] + range_coords[1] + + # metadata + # mol_metadata <- subcellular_data[,colnames(subcellular_data)[!colnames(subcellular_data) %in% c("cell_id", "transcript_id", "x_location", "y_location")], with = FALSE] + mol_metadata <- subcellular_data[,colnames(subcellular_data)[!colnames(subcellular_data) %in% c("cell_id", "transcript_id", "x_location", "y_location", "z_location")], with = FALSE] + set.seed(nrow(mol_metadata$id)) + mol_metadata$postfix <- paste0("_", ids::random_id(bytes = 3, use_openssl = FALSE)) + mol_metadata$assay_id <- "Assay1" + mol_metadata[, id:=do.call(paste0,.SD), .SDcols=c("id", "postfix")] + + # coord names + rownames(coords) <- mol_metadata$id + + # create VoltRon assay for molecules + mol_assay <- formAssay(coords = coords, image = image, type = "molecule", main_image = image_name, main_channel = channel_name) + + # merge assays in one section + if(verbose) + message("Merging assays ...") + sample.metadata <- SampleMetadata(cell_object) + object <- addAssay(cell_object, + assay = mol_assay, + metadata = mol_metadata, + assay_name = paste0(assay_name, "_mol"), + sample = sample.metadata["Assay1", "Sample"], + layer = sample.metadata["Assay1", "Layer"]) + + # connectivity + connectivity <- data.table::data.table(transcript_id = mol_metadata$id, cell_id = subcellular_data[["cell_id"]]) + if(is.numeric(connectivity$cell_id)){ + connectivity <- subset(connectivity, cell_id != -1) + connectivity[["cell_id"]] <- vrSpatialPoints(cell_object)[connectivity[["cell_id"]]] + } else { + connectivity <- subset(connectivity, cell_id != "UNASSIGNED") + connectivity$cell_assay_id <- "_Assay1" + connectivity[, cell_id:=do.call(paste0,.SD), .SDcols=c("cell_id", "cell_assay_id")] + connectivity$cell_assay_id <- NULL + } + + # add connectivity of spatial points across assays + object <- addLayerConnectivity(object, + connectivity = connectivity, + sample = sample.metadata["Assay1", "Sample"], + layer = sample.metadata["Assay1", "Layer"]) + + # return + return(object) + } + } +} + +#' generateXeniumImage +#' +#' Generate a low resolution DAPI image of the Xenium experiment +#' +#' @param dir.path Xenium output folder +#' @param increase.contrast increase the contrast of the image before writing +#' @param resolution_level the level of resolution within Xenium OME-TIFF image. Default: 7 (553x402) +#' @param overwrite_resolution if TRUE, the image "file.name" will be generated again although it exists at "dir.path" +#' @param output.path The path to the new morphology image created if the image should be saved to a location other than Xenium output folder. +#' @param file.name the name of the lowred morphology image. Default: morphology_lowres.tif +#' @param verbose verbose +#' @param ... additional parameters passed to the \link{writeImage} function +#' +#' @importFrom EBImage writeImage +#' +#' @details +#' The Xenium morphology_mip.ome.tif file that is found under the outs folder comes is an hyperstack of different resolutions of the DAPI image. +#' \link{generateXeniumImage} allows extracting only one of these layers by specifying the \code{resolution} parameter (Default: 7 for 553x402) among 1 to 8. +#' Lower incides of resolutions have higher higher resolutions, e.g. 1 for 35416x25778. Note that you may need to allocate larger memory of Java to import +#' higher resolution images. +#' +#' @export +#' +generateXeniumImage <- function(dir.path, increase.contrast = TRUE, resolution_level = 7, overwrite_resolution = FALSE, + output.path = NULL, file.name = "morphology_lowres.tif", verbose = TRUE, ...) { + + # file path to either Xenium output folder or specified folder + file.path <- paste0(dir.path, "/", file.name) + output.file <- paste0(output.path, "/", file.name) + + # check if the file exists in either Xenium output folder, or the specified location + if((file.exists(file.path) | file.exists(paste0(output.file))) & !overwrite_resolution){ + if(verbose) + message(file.name, " already exists!") + } else { + if (!requireNamespace('RBioFormats')) + stop("Please install RBioFormats package to extract xml from the ome.tiff file!: BiocManager::install('RBioFormats')") + if(dir.exists(paste0(dir.path, "/morphology_focus"))){ + if(verbose) + message("Loading morphology_focus_0000.ome.tif ...") + morphology_image_lowres <- RBioFormats::read.image(paste0(dir.path, "/morphology_focus/morphology_focus_0000.ome.tif"), + resolution = resolution_level, + subset=list(C=1)) + } else if(file.exists(paste0(dir.path, "/morphology_mip.ome.tif"))) { + if(verbose) + message("Loading morphology_mip.ome.tif ...") + morphology_image_lowres <- RBioFormats::read.image(paste0(dir.path, "/morphology_mip.ome.tif"), resolution = resolution_level) + } + + # pick a resolution level + image_info <- morphology_image_lowres@metadata$coreMetadata + if(verbose) + message(" Image Resolution (X:", image_info$sizeX, " Y:", image_info$sizeY, ") ...") + + # increase contrast using EBImage + if(increase.contrast) { + if(verbose) + message(" Increasing Contrast ...") + morphology_image_lowres <- (morphology_image_lowres/max(morphology_image_lowres)) + } + + # write to the same folder + if(verbose) + message(" Writing Tiff File ...") + if(is.null(output.path)){ + EBImage::writeImage(morphology_image_lowres, file = file.path, ...) + } else { + EBImage::writeImage(morphology_image_lowres, file = output.file, ...) + } + } + invisible() +} + +#### +## Visium #### +#### + +#' importVisium +#' +#' Importing Visium data +#' +#' @param dir.path path to Visium output folder +#' @param selected_assay selected assay from Visium +#' @param assay_name the assay name +#' @param sample_name the name of the sample +#' @param image_name the image name of the Visium assay, Default: main +#' @param channel_name the channel name of the image of the Visium assay, Default: H&E +#' @param inTissue if TRUE, only barcodes that are in the tissue will be kept (default: TRUE) +#' @param resolution_level the level of resolution of Visium image: "lowres" (default) or "hires" +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom magick image_read image_info +#' @importFrom rjson fromJSON +#' @importFrom utils read.csv +#' +#' @export +#' +importVisium <- function(dir.path, selected_assay = "Gene Expression", assay_name = "Visium", sample_name = NULL, image_name = "main", channel_name = "H&E", inTissue = TRUE, resolution_level = "lowres", ...) +{ + # raw counts + listoffiles <- list.files(dir.path) + datafile <- listoffiles[grepl("filtered_feature_bc_matrix.h5", listoffiles)][1] + datafile <- paste0(dir.path, "/", datafile) + if(file.exists(datafile)){ + rawdata <- import10Xh5(filename = datafile) + if(any(names(rawdata) %in% selected_assay)){ + rawdata <- as.matrix(rawdata[[selected_assay]]) + } else { + stop("There are no assays called ", selected_assay, " in the h5 file!") + } + } else { + stop("There are no files named 'filtered_feature_bc_matrix.h5' in the path") + } + + # resolution + if(!resolution_level %in% c("lowres","hires")) + stop("resolution_level should be either 'lowres' or 'hires'!") + + # image + image_file <- paste0(dir.path, paste0("/spatial/tissue_", resolution_level, "_image.png")) + if(file.exists(image_file)){ + image <- magick::image_read(image_file) + info <- image_info(image) + } else { + stop("There are no spatial image files in the path") + } + + # coordinates + coords_file <- list.files(paste0(dir.path, "/spatial/"), full.names = TRUE) + coords_file <- coords_file[grepl("tissue_positions",coords_file)] + if(length(coords_file) == 1){ + if(grepl("tissue_positions_list.csv", coords_file)) { + coords <- utils::read.csv(file = coords_file, header = FALSE) + colnames(coords) <- c("barcode", "in_tissue", "array_row", "array_col", "pxl_row_in_fullres", "pxl_col_in_fullres") + } else { + coords <- utils::read.csv(file = coords_file, header = TRUE) + } + } else if(length(coords_file) > 1) { + stop("There are more than 1 position files in the path") + } else { + stop("There are no files named 'tissue_positions.csv' in the path") + } + + if(inTissue){ + coords <- coords[coords$in_tissue==1,] + rawdata <- rawdata[,colnames(rawdata) %in% coords$barcode] + } + coords <- coords[match(colnames(rawdata), coords$barcode),] + spotID <- coords$barcode + coords <- as.matrix(coords[,c("pxl_col_in_fullres", "pxl_row_in_fullres")], ) + colnames(coords) <- c("x", "y") + rownames(coords) <- spotID + + # scale coordinates + scale_file <- paste0(dir.path, "/spatial/scalefactors_json.json") + if(file.exists(scale_file)){ + scalefactors <- rjson::fromJSON(file = scale_file) + scales <- scalefactors[[paste0("tissue_", resolution_level, "_scalef")]] + params <- list( + nearestpost.distance = 200*scales, # distance to nearest spot + spot.radius = scalefactors$spot_diameter_fullres*scales, + vis.spot.radius = scalefactors$spot_diameter_fullres*scales*2, + spot.type = "circle") + coords <- coords*scales + coords[,2] <- info$height - coords[,2] + } else { + stop("There are no files named 'scalefactors_json.json' in the path") + } + + # create VoltRon + formVoltRon(rawdata, metadata = NULL, image, coords, main.assay = assay_name, params = params, assay.type = "spot", + image_name = image_name, main_channel = channel_name, sample_name = sample_name, + feature_name = ifelse(selected_assay == "Gene Expression", "RNA", "main"), ...) +} + +#' importVisiumHD +#' +#' Importing VisiumHD data +#' +#' @param dir.path path to Visium output folder +#' @param bin.size bin size of the VisiumHD output (Exp: "2", "8" and "16") +#' @param selected_assay selected assay from Visium +#' @param assay_name the assay name +#' @param sample_name the name of the sample +#' @param image_name the image name of the Visium assay, Default: main +#' @param channel_name the channel name of the image of the Visium assay, Default: H&E +#' @param inTissue if TRUE, only barcodes that are in the tissue will be kept (default: TRUE) +#' @param resolution_level the level of resolution of Visium image: "lowres" (default) or "hires" +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom magick image_read image_info +#' @importFrom rjson fromJSON +#' @importFrom utils read.csv +#' +#' @export +#' +importVisiumHD <- function(dir.path, bin.size = "8", selected_assay = "Gene Expression", assay_name = "VisiumHD", sample_name = NULL, image_name = "main", channel_name = "H&E", inTissue = TRUE, resolution_level = "lowres", ...){ + + # raw counts + bin.size <- formatC(as.numeric(bin.size), width = 3, format = "d", flag = "0") + dir.path <- paste0(dir.path, "binned_outputs/square_", bin.size, "um/") + listoffiles <- list.files(dir.path) + datafile <- listoffiles[grepl("filtered_feature_bc_matrix.h5", listoffiles)][1] + datafile <- paste0(dir.path, "/", datafile) + if(file.exists(datafile)){ + rawdata <- import10Xh5(filename = datafile) + if(any(names(rawdata) %in% selected_assay)){ + rawdata <- as(rawdata[[selected_assay]], "CsparseMatrix") + } else { + stop("There are no assays called ", selected_assay, " in the h5 file!") + } + } else { + stop("There are no files named 'filtered_feature_bc_matrix.h5' in the path") + } + + # resolution + if(!resolution_level %in% c("lowres","hires")) + stop("resolution_level should be either 'lowres' or 'hires'!") + + # image + image_file <- paste0(dir.path, paste0("/spatial/tissue_", resolution_level, "_image.png")) + if(file.exists(image_file)){ + image <- magick::image_read(image_file) + info <- magick::image_info(image) + } else { + stop("There are no spatial image files in the path") + } + + # coordinates + if(!requireNamespace("arrow")) + stop("Please install arrow package to import VisiumHD!: install.packages('arrow')") + coords_file <- list.files(paste0(dir.path, "/spatial/"), full.names = TRUE) + coords_file <- coords_file[grepl("tissue_positions",coords_file)] + if(length(coords_file) == 1){ + coords <- data.table::as.data.table(arrow::read_parquet(coords_file, as_data_frame = FALSE)) + } else if(length(coords_file) > 1) { + stop("There are more than 1 position files in the path") + } else { + stop("There are no files named 'tissue_positions.parquet' in the path") + } + if(inTissue){ + coords <- subset(coords, in_tissue == 1) + rawdata <- rawdata[,colnames(rawdata) %in% coords$barcode] + } + coords <- coords[match(colnames(rawdata), coords$barcode),] + spotID <- coords$barcode + coords <- as.matrix(coords[,c("pxl_col_in_fullres", "pxl_row_in_fullres")], ) + colnames(coords) <- c("x", "y") + rownames(coords) <- spotID + + # scale coordinates + scale_file <- paste0(dir.path, "/spatial/scalefactors_json.json") + if(file.exists(scale_file)){ + scalefactors <- rjson::fromJSON(file = scale_file) + scales <- scalefactors[[paste0("tissue_", resolution_level, "_scalef")]] + params <- list( + nearestpost.distance = scalefactors$spot_diameter_fullres*scales*(3/sqrt(2)), + spot.radius = scalefactors$spot_diameter_fullres*scales, + vis.spot.radius = scalefactors$spot_diameter_fullres*scales, + spot.type = "square") + coords <- coords*scales + coords[,2] <- info$height - coords[,2] + } else { + stop("There are no files named 'scalefactors_json.json' in the path") + } + + # create VoltRon + formVoltRon(rawdata, metadata = NULL, image, coords, main.assay = assay_name, params = params, assay.type = "spot", + image_name = image_name, main_channel = channel_name, sample_name = sample_name, + feature_name = ifelse(selected_assay == "Gene Expression", "RNA", "main"), ...) +} + +#### +## Auxiliary #### +#### + +#' import10Xh5 +#' +#' import the sparse matrix from the H5 file +#' +#' @param filename the path to h5 file +#' +#' @importFrom Matrix sparseMatrix +#' +#' @export +import10Xh5 <- function(filename){ + + # check package + if(!requireNamespace("rhdf5")){ + stop("You have to install the rhdf5 package!: BiocManager::install('rhdf5')") + } + + # check file + if(file.exists(filename)){ + matrix.10X <- rhdf5::h5read(file = filename, name = "matrix") + } else { + stop("There are no files named ", filename," in the path") + } + + # get data, barcodes and feature + features <- matrix.10X[["features"]][["name"]] + feature_type <- matrix.10X[["features"]][["feature_type"]] + cells <- matrix.10X[["barcodes"]] + mat <- matrix.10X[["data"]] + indices <- matrix.10X[["indices"]] + indptr <- matrix.10X[["indptr"]] + sparse.mat <- sparseMatrix(i = indices + 1, p = indptr, + x = as.numeric(mat), dims = c(length(features), length(cells)), + repr = "T") + colnames(sparse.mat) <- cells + rownames(sparse.mat) <- make.unique(features) + + # separate feature types + matrix.10X <- list() + for(feat in unique(feature_type)){ + cur_features <- features[feature_type %in% feat] + cur_mat <- sparse.mat[features %in% cur_features,] + matrix.10X[[feat]] <- cur_mat + } + + return(matrix.10X) +} + +#### +# Nanostring #### +#### + +#### +## GeoMx #### +#### + +#' importGeoMx +#' +#' Import GeoMx data +#' +#' @param dcc.path path to the folder where the dcc files are found +#' @param pkc.file path to the pkc file +#' @param summarySegment the metadata csv (sep = ";") or excel file, if the file is an excel file, \code{summarySegmentSheetName} should be provided as well. +#' @param summarySegmentSheetName the sheet name of the excel file, \code{summarySegment} +#' @param assay_name the assay name, default: GeoMx +#' @param image the reference morphology image of the GeoMx assay +#' @param segment_polygons if TRUE, the ROI polygons are parsed from the OME.TIFF file +#' @param ome.tiff the OME.TIFF file of the GeoMx experiment if exists +#' @param resolution_level the level of resolution within GeoMx OME-TIFF image, Default: 3 +#' @param image_name the image name of the Visium assay, Default: main +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom dplyr %>% full_join +#' @importFrom utils read.csv +#' @importFrom magick image_info image_read +#' +#' @export +#' +importGeoMx <- function(dcc.path, pkc.file, summarySegment, summarySegmentSheetName, assay_name = "GeoMx", + image = NULL, segment_polygons = FALSE, ome.tiff = NULL, resolution_level = 3, image_name = "main", ...) +{ + # Get pkc file + if(file.exists(pkc.file)){ + pkcdata <- readPKC(pkc.file) + } else { + stop("pkc file is not found!") + } + + # Get dcc file + if(file.exists(dcc.path)){ + dcc_files <- dir(dcc.path, pattern = ".dcc$", full.names = TRUE) + if(length(dcc_files) == 0){ + stop("no dcc files are found under ", dcc.path) + } else { + dcc_files <- dcc_files[!grepl("A01.dcc$", dcc_files)] + dcc_filenames <- dir(dcc.path, pattern = ".dcc$", full.names = FALSE) + dcc_filenames <- dcc_filenames[!grepl("A01.dcc$", dcc_filenames)] + dcc_filenames <- gsub(".dcc$", "", dcc_filenames) + dcc_filenames <- gsub("-", "_", dcc_filenames) + dccData <- sapply(dcc_files, readDCC, simplify = FALSE, USE.NAMES = FALSE) + names(dccData) <- dcc_filenames + } + } else { + stop("path to dcc files does not exist!") + } + + # merge dcc files + rawdata <- NULL + for(i in seq_len(length(dccData))){ + cur_data <- dccData[[i]]$Code_Summary + colnames(cur_data) <- c("RTS_ID", dcc_filenames[i]) + if(i == 1){ + rawdata <- cur_data + } else { + suppressMessages(rawdata <- rawdata %>% full_join(cur_data)) + } + } + rawdata[is.na(rawdata)] <- 0 + rownames(rawdata) <- rawdata$RTS_ID + rawdata <- rawdata[,!colnames(rawdata) %in% "RTS_ID"] + + # get negative probes and targets + # NegProbes <- pkcdata$RTS_ID[pkcdata$Target == "NegProbe-WTX"] + NegProbes <- pkcdata$RTS_ID[grepl("NegProbe", pkcdata$Target)] + + # negative probes + rawdata_neg <- rawdata[rownames(rawdata) %in% NegProbes, ] + rownames(rawdata_neg) <- paste0(rownames(rawdata_neg), "_", + pkcdata$Target[match(rownames(rawdata_neg), pkcdata$RTS_ID)]) + rawdata_neg <- as.matrix(rawdata_neg) + + # other probes + rawdata <- rawdata[!rownames(rawdata) %in% NegProbes, ] + rownames(rawdata) <- pkcdata$Target[match(rownames(rawdata), pkcdata$RTS_ID)] + rawdata <- as.matrix(rawdata) + + # get segment summary + if(file.exists(summarySegment)){ + if(grepl(".xls$|.xlsx$", summarySegment)){ + if (!requireNamespace('xlsx')) + stop("Please install xlsx package for using the read.xlsx function!: install.packages('xlsx')") + if(!is.null(summarySegmentSheetName)){ + segmentsummary <- xlsx::read.xlsx(summarySegment, sheetName = summarySegmentSheetName) + } else { + stop("Please provide 'summarySegmentSheetName' for the excel sheet name!") + } + } else if(grepl(".csv$", summarySegment)) { + segmentsummary <- utils::read.csv(summarySegment, row.names = NULL, header = T, sep = ";") + } + rownames(segmentsummary) <- gsub(".dcc$", "", segmentsummary$Sample_ID) + rownames(segmentsummary) <- gsub("-", "_", rownames(segmentsummary)) + if(all(dcc_filenames %in% rownames(segmentsummary))){ + segmentsummary <- segmentsummary[dcc_filenames, ] + } else{ + stop("Some GeoMx dcc files are not represented in the segment summary file!") + } + } else { + stop(summarySegment, " is not found!") + } + + # get image + if(!is.null(image)){ + if(file.exists(image)){ + image <- magick::image_read(image) + } else { + stop(image, " is not found!") + } + } else { + stop("Please provide a image for the GeoMx data!") + } + geomx_image_info <- magick::image_info(image) + + # get coordinates + coords <- segmentsummary[,c("X","Y")] + colnames(coords) <- c("x", "y") + rownames(coords) <- segmentsummary$ROI.name + coords <- rescaleGeoMxPoints(coords, segmentsummary, geomx_image_info) + + # get ROI segments (polygons) + segments <- list() + if(!is.null(ome.tiff)){ + + # get segments + segments <- importGeoMxSegments(ome.tiff, segmentsummary, geomx_image_info) + segments <- segments[rownames(coords)] + + # parse channels if ome.tiff is given + channels <- importGeoMxChannels(ome.tiff, segmentsummary, geomx_image_info, resolution_level = resolution_level) + image <- c(list(scanimage = image), channels) + } + + # create VoltRon for non-negative probes + object <- formVoltRon(rawdata, metadata = segmentsummary, image, coords, segments, main.assay = assay_name, assay.type = "ROI", + feature_name = "RNA", image_name = image_name, ...) + + # add negative probe assay as new feature set + object <- addFeature(object, assay = assay_name, data = rawdata_neg, feature_name = "NegProbe") + + # return + return(object) +} + +#' readPKC +#' +#' Read a NanoString Probe Kit Configuration (PKC) file. Adapted from \code{GeomxTools} package. +#' +#' @param file A character string containing the path to the PKC file. +#' @param default_pkc_vers Optional list of pkc file names to use as default if more than one pkc version of each module is provided. +#' +#' @importFrom rjson fromJSON +#' +#' @noRd +readPKC <- function (file, default_pkc_vers = NULL) +{ + pkc_json_list <- lapply(file, function(pkc_file) { + rjson::fromJSON(file = pkc_file) + }) + pkc_names <- unlist(lapply(file, function(pkc_file) { + base_pkc_name <- gsub(".pkc", "", trimws(basename(pkc_file))) + return(base_pkc_name) + })) + names(pkc_json_list) <- pkc_names + pkc_modules <- basename(unlist(lapply(pkc_names, sub, pattern = "_[^_]+$", + replacement = ""))) + names(pkc_modules) <- pkc_names + header <- list(PKCFileName = sapply(pkc_json_list, function(list) list[["Name"]]), + PKCModule = pkc_modules, + PKCFileVersion = sapply(pkc_json_list, function(list) list[["Version"]]), + PKCFileDate = sapply(pkc_json_list, function(list) list[["Date"]]), + AnalyteType = sapply(pkc_json_list, function(list) list[["AnalyteType"]]), + MinArea = sapply(pkc_json_list, function(list) list[["MinArea"]]), + MinNuclei = sapply(pkc_json_list, function(list) list[["MinNuclei"]])) + + multi_idx <- duplicated(header[["PKCModule"]]) + multi_mods <- unique(header[["PKCModule"]][multi_idx]) + if (length(multi_mods) < 1) { + if (!is.null(default_pkc_vers)) { + warning("Only one version found per PKC module. ", + "No PKCs need to be combined. ", "Therefore, no default PKC versions will be used.") + } + } + else { + use_pkc_names <- lapply(multi_mods, function(mod) { + mod_idx <- header[["PKCModule"]] == mod + max_vers <- as.numeric(as.character(max(as.numeric_version(header[["PKCFileVersion"]][mod_idx])))) + max_name <- names(header[["PKCFileVersion"]][header[["PKCFileVersion"]] == max_vers]) + return(max_name) + }) + names(use_pkc_names) <- multi_mods + if (!is.null(default_pkc_vers)) { + default_names <- unlist(lapply(default_pkc_vers, function(pkc_file) { + base_pkc_name <- gsub(".pkc", "", trimws(basename(pkc_file))) + return(base_pkc_name) + })) + default_mods <- unlist(lapply(default_names, sub, pattern = "_[^_]+$", + replacement = "")) + dup_defaults <- default_names[duplicated(default_mods) | + duplicated(default_mods, fromLast = TRUE)] + if (!all(default_names %in% names(header[["PKCFileName"]]))) { + removed_pkcs <- default_pkc_vers[!default_names %in% + names(header[["PKCFileName"]])] + stop("Could not match all default PKC versions with a PKC file name. ", + "Check default file names match exactly to a PKC file name.\n", + paste0("Unmatched default PKC versions: ", + removed_pkcs)) + } + else if (length(dup_defaults) > 0) { + stop("There should only be one default PKC version per module. ", + "Ensure only one version per module in default PKCs list.\n", + "Multiple default PKC version conflicts: ", + paste(dup_defaults, collapse = ", ")) + } + else { + use_pkc_names[default_mods] <- default_names + } + } + } + rtsid_lookup_df <- generate_pkc_lookup(pkc_json_list) + rtsid_lookup_df$Negative <- grepl("Negative", rtsid_lookup_df$CodeClass) + rtsid_lookup_df$RTS_ID <- gsub("RNA", "RTS00", rtsid_lookup_df[["RTS_ID"]]) + # rtsid_lookup_df <- S4Vectors::DataFrame(rtsid_lookup_df) + rtsid_lookup_df <- data.table::data.table(rtsid_lookup_df) + if (length(multi_mods) > 0) { + for (mod in names(use_pkc_names)) { + mod_vers <- names(header[["PKCModule"]])[header[["PKCModule"]] == + mod] + mod_lookup <- rtsid_lookup_df[rtsid_lookup_df$Module %in% + mod_vers, ] + mod_tab <- table(mod_lookup$RTS_ID) + remove_rts <- names(mod_tab[mod_tab != length(mod_vers)]) + if (length(remove_rts) > 0) { + warning("The following probes were removed from analysis", + " as they were not found in all PKC module versions used.\n", + # paste(utils::capture.output(print(subset(rtsid_lookup_df, + # subset = RTS_ID %in% remove_rts))), collapse = "\n") + ) + rtsid_lookup_df <- subset(rtsid_lookup_df, subset = !RTS_ID %in% + remove_rts) + } + remove_vers <- mod_vers[mod_vers != use_pkc_names[mod]] + rtsid_lookup_df <- subset(rtsid_lookup_df, subset = !Module %in% + remove_vers) + warning("The following PKC versions were removed from analysis", + " as they were overridden by a newer PKC version or", + " were overridden by a user-defined default PKC version.\n", + paste(remove_vers, collapse = ", ")) + header <- lapply(header, function(elem) { + elem[!names(elem) %in% remove_vers] + }) + } + } + # S4Vectors::metadata(rtsid_lookup_df) <- header + return(rtsid_lookup_df) +} + +#' readDCC +#' +#' Read a NanoString GeoMx Digital Count Conversion (DCC) file. +#' +#' @param file A character string containing the path to the DCC file. +#' +#' @noRd +readDCC <- function(file) +{ + lines <- trimws(readLines(file)) + trimGalore <- grep("trimGalore", lines) + if (length(trimGalore) > 0) { + Raw <- grep("Raw", lines) + lines <- lines[-c(trimGalore:(Raw - 1))] + } + lines <- gsub("SoftwareVersion,\"GeoMx_NGS_Pipeline_ ", + "SoftwareVersion,", lines) + lines <- gsub("SoftwareVersion,\"GeoMx_NGS_Pipeline_", "SoftwareVersion,", + lines) + lines <- gsub("SoftwareVersion,DRAGEN_GeoMx_", "SoftwareVersion,", + lines) + lines[grepl("SoftwareVersion", lines)] <- gsub("\"", "", + lines[grepl("SoftwareVersion", lines)]) + tags <- c("Header", "Scan_Attributes", "NGS_Processing_Attributes", "Code_Summary") + output <- sapply(tags, function(tag) { + bounds <- charmatch(sprintf(c("<%s>", "</%s>"), tag), + lines) + if (anyNA(bounds) || bounds[1L] + 1L >= bounds[2L]) + lines[integer(0)] + else lines[(bounds[1L] + 1L):(bounds[2L] - 1L)] + }, simplify = FALSE) + for (tag in c("Header", "Scan_Attributes", "NGS_Processing_Attributes")) { + while (length(bad <- grep(",", output[[tag]], invert = TRUE)) > + 0L) { + bad <- bad[1L] + if (bad == 1L) + stop(sprintf("%s section has malformed first line", + tag)) + fixed <- output[[tag]] + fixed[bad - 1L] <- sprintf("%s %s", fixed[bad - + 1L], fixed[bad]) + output[[tag]] <- fixed[-bad] + } + output[[tag]] <- strsplit(output[[tag]], split = ",") + output[[tag]] <- structure(lapply(output[[tag]], function(x) if (length(x) == + 1L) + "" + else x[2L]), names = lapply(output[[tag]], `[`, 1L), + class = "data.frame", row.names = basename(file)) + } + cols <- c("FileVersion", "SoftwareVersion") + if (!(all(cols %in% colnames(output[["Header"]])))) + stop("Header section must contain \"FileVersion\" and \"SoftwareVersion\"") + output[["Header"]][, cols] <- lapply(output[["Header"]][, + cols], numeric_version) + fileVersion <- output[["Header"]][1L, "FileVersion"] + if (!(numeric_version(fileVersion) %in% numeric_version(c("0.01", + "0.02")))) + stop("\"FileVersion\" in Header section must be 0.01 or 0.02") + for (section in c("Header", "Scan_Attributes", "NGS_Processing_Attributes")) { + valid <- .validDccSchema(output[[section]], fileVersion, + section) + if (!isTRUE(valid)) + stop(valid) + } + output[["Header"]][["Date"]] <- as.Date(output[["Header"]][["Date"]], + format = "%Y-%m-%d") + cols <- c("Raw", "Trimmed", "Stitched", "Aligned", "umiQ30", + "rtsQ30") + output[["NGS_Processing_Attributes"]][, cols] <- lapply(output[["NGS_Processing_Attributes"]][, + cols], as.numeric) + names(output[["Scan_Attributes"]])[names(output[["Scan_Attributes"]]) == + "ID"] <- "SampleID" + output[["Code_Summary"]] <- paste0("RTS_ID,Count\n", paste(output[["Code_Summary"]], + collapse = "\n")) + output[["Code_Summary"]] <- utils::read.csv(textConnection(output[["Code_Summary"]]), + colClasses = c(RTS_ID = "character", Count = "numeric")) + output[["Code_Summary"]][["Count"]] <- as.integer(round(output[["Code_Summary"]][["Count"]])) + rn <- output[["Code_Summary"]][["RTS_ID"]] + if ((ndups <- anyDuplicated(rn)) > 0L) { + warning(sprintf("removed %d rows from \"Code_Summary\" due to duplicate rownames", + ndups)) + ok <- which(!duplicated(rn, fromLast = FALSE) & !duplicated(rn, + fromLast = TRUE)) + rn <- rn[ok] + output[["Code_Summary"]] <- output[["Code_Summary"]][ok, + , drop = FALSE] + } + rownames(output[["Code_Summary"]]) <- rn + output[["NGS_Processing_Attributes"]][, "DeduplicatedReads"] <- sum(output[["Code_Summary"]][["Count"]]) + return(output) +} + + +#' importGeoMxSegments +#' +#' Import ROI polygons from the OME.TIFF file +#' +#' @param ome.tiff the OME.TIFF file of the GeoMx Experiment +#' @param summary segmentation summary data frame +#' @param imageinfo image information +#' +#' @noRd +importGeoMxSegments <- function(ome.tiff, summary, imageinfo){ + + # check file + if(file.exists(ome.tiff)){ + if(grepl(".ome.tiff$|.ome.tif$", ome.tiff)){ + if (!requireNamespace('RBioFormats')) + stop("Please install RBioFormats package to extract xml from the ome.tiff file!: BiocManager::install('RBioFormats')") + if (!requireNamespace('XML')) + stop("Please install XML package to extract xml from the ome.tiff file!") + omexml <- RBioFormats::read.omexml(ome.tiff) + } else if(grepl(".xml$", ome.tiff)){ + omexml <- XML::xmlParse(file = ome.tiff) + } else { + stop("Please provide either an ome.tiff or .xml file!") + } + omexml <- XML::xmlToList(omexml, simplify = TRUE) + } else { + stop("There are no files named ", ome.tiff," in the path") + } + + # get ROIs + ROIs <- omexml[which(names(omexml) == "ROI")] + + # get masks for each ROI + mask_lists <- list() + for(i in seq_len(length(ROIs))){ + cur_ROI <- ROIs[[i]] + + # if the shape is a polygon + if("Polygon" %in% names(cur_ROI$Union)){ + coords <- strsplit(cur_ROI$Union$Polygon, split = "\\n") + coords <- strsplit(coords$Points, split = " ")[[1]] + coords <- sapply(coords, function(x) as.numeric(strsplit(x, split = ",")[[1]]), USE.NAMES = FALSE) + coords <- as.data.frame(t(coords)) + colnames(coords) <- c("x", "y") + coords <- rescaleGeoMxPoints(coords, summary, imageinfo) + coords <- data.frame(id = cur_ROI$Union$Label[["Text"]], coords) + mask_lists[[cur_ROI$Union$Label[["Text"]]]] <- data.frame(coords) + + # if the shape is an ellipse + } else if("Ellipse" %in% names(cur_ROI$Union)){ + coords <- as.numeric(cur_ROI$Union$Ellipse[c("X","Y", "RadiusX", "RadiusY")]) + coords <- as.data.frame(matrix(coords, nrow = 1)) + colnames(coords) <- c("x", "y", "rx", "ry") + coords[,c("x", "y")] <- rescaleGeoMxPoints(coords[,c("x", "y")], summary, imageinfo) + coords$rx <- coords$rx * imageinfo$width/summary$Scan.Width[1] + coords$ry <- coords$ry * imageinfo$height/summary$Scan.Height[1] + coords <- data.frame(id = cur_ROI$Union$Label[["Text"]], coords) + mask_lists[[cur_ROI$Union$Label[["Text"]]]] <- coords + } + } + + mask_lists +} + +#' rescaleGeoMxROIs +#' +#' Rescale GeoMx point (center or polygon corners of ROI) coordinates +#' +#' @param pts coordinates of the cells from the Xenium assays +#' @param summary segmentation summary data frame +#' @param imageinfo image information +#' +#' @noRd +rescaleGeoMxPoints <- function(pts, summary, imageinfo){ + + xRatio = imageinfo$width/summary$Scan.Width[1] + yRatio = imageinfo$height/summary$Scan.Height[1] + pts$x = (pts$x - summary$Scan.Offset.X[1]) * xRatio + pts$y = (pts$y - summary$Scan.Offset.Y[1]) * yRatio + pts$y = imageinfo$height - pts$y + pts <- as.matrix(pts) + + # return + return(pts) +} + +#' importGeoMxChannels +#' +#' Rescale GeoMx channels with respect to the scan image +#' +#' @param ome.tiff the OME.TIFF file of the GeoMx Experiment +#' @param summary segmentation summary data frame +#' @param imageinfo image information +#' @param resolution_level the resolution level (1-7) of the image parsed from the OME.TIFF file +#' +#' @param RBioFormats read.image +#' @param EBImage as.Image +#' @param grDevices as.raster +#' @param magick image_read +#' +#' @noRd +importGeoMxChannels <- function(ome.tiff, summary, imageinfo, resolution_level){ + + # check file + if(file.exists(ome.tiff)){ + if(grepl(".ome.tiff$|.ome.tif$", ome.tiff)){ + if (!requireNamespace('RBioFormats')) + stop("Please install RBioFormats package to extract xml from the ome.tiff file!: BiocManager::install('RBioFormats')") + if (!requireNamespace('XML')) + stop("Please install XML package to extract xml from the ome.tiff file!") + omexml <- RBioFormats::read.omexml(ome.tiff) + omexml <- XML::xmlToList(omexml, simplify = TRUE) + } else { + warning("ome.tiff format not found!") + return(NULL) + } + } else { + stop("There are no files named ", ome.tiff," in the path") + } + + # get all channels + ome.tiff <- RBioFormats::read.image(ome.tiff, resolution = resolution_level) + + # get frame information + nframes <- ome.tiff@metadata$coreMetadata$imageCount + frames <- EBImage::getFrames(ome.tiff) + frames <- lapply(frames, function(x){ + img <- magick::image_read(grDevices::as.raster(EBImage::as.Image(x))) + rescaleGeoMxImage(img, summary, imageinfo, resolution_level = resolution_level) + }) + + # get channel names + omexml <- omexml$StructuredAnnotations[seq(from = 1, by = 2, length.out = nframes)] + channel_names <- lapply(omexml, function(x){ + x$Value$ChannelInfo$BiologicalTarget + }) + names(frames) <- channel_names + + # return frames + return(frames) +} + +#' rescaleGeoMxImage +#' +#' Rescale GeoMx channels with respect to the scan image +#' +#' @param img coordinates of the cells from the Xenium assays +#' @param summary segmentation summary data frame +#' @param imageinfo image information +#' @param resolution_level the resolution level (1-7) of the image parsed from the OME.TIFF file +#' +#' @param magick image_crop +#' +#' @noRd +rescaleGeoMxImage <- function(img, summary, imageinfo, resolution_level){ + + # adjust offset and scan size to the resolution level + offset.x <- summary$Scan.Offset.X[1]/(2*(resolution_level-1)) + offset.y <- summary$Scan.Offset.Y[1]/(2*(resolution_level-1)) + scan.width <- summary$Scan.Width[1]/(2*(resolution_level-1)) + scan.height <- summary$Scan.Height[1]/(2*(resolution_level-1)) + + # crop image given adjusted offsets + new_img <- magick::image_crop(img, geometry = paste0(scan.width, "x", scan.height, "+", offset.x, "+", offset.y)) + + # resize image + new_img <- magick::image_resize(new_img, geometry = paste0(imageinfo$width, "x", imageinfo$height)) + + # return + return(new_img) +} + +#### +## CosMx #### +#### + +#' importCosMx +#' +#' Import CosMx data +#' +#' @param tiledbURI the path to the tiledb folder +#' @param assay_name the assay name, default: CosMx +#' @param image the reference morphology image of the CosMx assay +#' @param image_name the image name of the CosMx assay, Default: main +#' @param ome.tiff the OME.TIFF file of the CosMx experiment if exists +#' @param import_molecules if TRUE, molecule assay will be created along with cell assay. +#' @param verbose verbose +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom data.table data.table +#' @importFrom ids random_id +#' +#' @export +#' +importCosMx <- function(tiledbURI, assay_name = "CosMx", + image = NULL, image_name = "main", ome.tiff = NULL, import_molecules = FALSE, verbose = TRUE, ...) +{ + # check tiledb and tiledbsc + if (!requireNamespace("tiledb", quietly = TRUE)) + stop("Please install the tiledb package: \n + remotes::install_github('TileDB-Inc/TileDB-R', force = TRUE, ref = '0.17.0')") + if (!requireNamespace("tiledbsc", quietly = TRUE)) + stop("Please install the tiledbsc package: \n + remotes::install_github('tiledb-inc/tiledbsc', force = TRUE, ref = '8157b7d54398b1f957832f37fff0b173d355530e')") + + # get tiledb + if(verbose) + message("Scanning TileDB array for cell data ...") + tiledb_scdataset <- tiledbsc::SOMACollection$new(uri = tiledbURI, verbose = FALSE) + + # raw counts + counts <- tiledb_scdataset$somas$RNA$X$members$counts$to_matrix(batch_mode = TRUE) + counts <- as.matrix(counts) + + # cell metadata + metadata <- tiledb_scdataset$somas$RNA$obs$to_dataframe() + + # coordinates + coords <- as.matrix(metadata[,c("x_slide_mm", "y_slide_mm")]) + colnames(coords) <- c("x","y") + + # transcripts + if(import_molecules){ + if(verbose) + message("Scanning TileDB array for molecule data ...") + subcellular <- tiledb::tiledb_array( + tiledb_scdataset$somas$RNA$obsm$members$transcriptCoords$uri, + return_as="data.table")[] + colnames(subcellular)[colnames(subcellular)=="target"] <- "gene" + } + + # get slides and construct VoltRon objects for each slides + slides <- unique(metadata$slide_ID_numeric) + + # for each slide create a VoltRon object with combined layers + vr_list <- list() + for(slide in slides){ + + # cell assay + if(verbose) + message("Creating cell level assay for slide ", slide, " ...") + + # slide info + cur_coords <- coords[metadata$slide_ID_numeric == slide,] + cur_counts <- counts[,rownames(cur_coords)] + cur_metadata <- metadata[rownames(cur_coords),] + + # create VoltRon object + cell_object <- formVoltRon(data = cur_counts, metadata = cur_metadata, image = image, coords = cur_coords, + main.assay = assay_name, assay.type = "cell", image_name = image_name, main_featureset = "RNA", ...) + cell_object$Sample <- paste0("Slide", slide) + + # molecule assay + if(import_molecules){ + + # get slide + if(verbose) + message("Creating molecule level assay for slide ", slide, " ...") + cur_subcellular <- subset(subcellular, slideID == slide) + + # coordinates + mol_coords <- as.matrix(cur_subcellular[,c("x_FOV_px", "y_FOV_px")]) + colnames(mol_coords) <- c("x", "y") + + # get subcellular data components + mol_metadata <- cur_subcellular[,colnames(cur_subcellular)[!colnames(cur_subcellular) %in% c("CellId", "cell_id", "x_FOV_px", "y_FOV_px")], with = FALSE] + set.seed(nrow(mol_metadata)) + mol_metadata$id <- rownames(mol_metadata) + mol_metadata$postfix <- paste0("_", ids::random_id(bytes = 3, use_openssl = FALSE)) + mol_metadata$assay_id <- "Assay1" + mol_metadata[, id:=do.call(paste0,.SD), .SDcols=c("id", "postfix")] + + # coord names + rownames(mol_coords) <- mol_metadata$id + + # create VoltRon assay for molecules + mol_assay <- formAssay(coords = mol_coords, image = image, type = "molecule", main_image = image_name) + + # merge assays in one section + if(verbose) + message("Merging assays for slide ", slide, " ...") + sample.metadata <- SampleMetadata(cell_object) + cell_object <- addAssay(cell_object, + assay = mol_assay, + metadata = mol_metadata, + assay_name = paste0(assay_name, "_mol"), + sample = sample.metadata["Assay1", "Sample"], + layer = sample.metadata["Assay1", "Layer"]) + } + vr_list <- append(vr_list, cell_object) + } + + # return + if(verbose) + message("Merging slides ...") + vr <- merge(vr_list[[1]], vr_list[-1]) +} + +#' generateCosMxImage +#' +#' Generates a low resolution Morphology image of the CosMx experiment +#' +#' @param dir.path CosMx output folder +#' @param increase.contrast increase the contrast of the image before writing +#' @param output.path The path to the new morphology image created if the image should be saved to a location other than Xenium output folder. +#' @param verbose verbose +#' @param ... additional parameters passed to the \link{writeImage} function +#' +#' @importFrom magick image_read image_contrast +#' @importFrom EBImage writeImage +#' +#' @export +generateCosMxImage <- function(dir.path, increase.contrast = TRUE, output.path = NULL, verbose = TRUE, ...) { + + # check package + if(!requireNamespace("reshape2")){ + stop("You have to install the reshape2 package!: install.packages('reshape2')") + } + + # file path to either Xenium output folder or specified folder + file.path <- paste0(dir.path, "/CellComposite_lowres.tif") + output.file <- paste0(output.path, "/CellComposite_lowres.tif") + + # check if the file exists in either Xenium output folder, or the specified location + if(file.exists(file.path) | file.exists(paste0(output.file))){ + if(verbose) + message("CellComposite_lowres.tif already exists! \n") + return(NULL) + } + + # FOV positions of CosMx + list_of_files <- list.files(dir.path) + fov_positions_path <- paste0(dir.path, "/", list_of_files[grepl("fov_positions_file.csv$",list_of_files)][1]) + fov_positions <- read.csv(fov_positions_path) + + # manipulate fov positions matrix + if(verbose) + message("Getting FOV Positions \n") + relative_fov_positions <- fov_positions + x_min <- min(relative_fov_positions$x_global_px) + y_min <- min(relative_fov_positions$y_global_px) + x_gap <- diff(unique(fov_positions$x_global_px))[1] + y_gap <- diff(unique(fov_positions$y_global_px))[1] + relative_fov_positions[,c("x_global_px","y_global_px")] <- t(apply(relative_fov_positions[,c("x_global_px","y_global_px")], 1, function(cell){ + c((cell[1]-x_min)/x_gap,(cell[2]-y_min)/y_gap) + })) + relative_fov_positions <- relative_fov_positions[order(relative_fov_positions$y_global_px, decreasing = TRUE),] + + # Combine Images of the FOV grid + if(verbose) + message("Loading FOV tif files \n") + image.dir.path <- paste0(dir.path,"/CellComposite/") + morphology_image_data <- NULL + for(i in relative_fov_positions$fov){ + image_path <- paste0(image.dir.path, "CellComposite_F", str_pad(as.character(i), 3, pad = 0), ".jpg") + image_data <- magick::image_read(image_path) %>% magick::image_resize("x500") %>% magick::image_raster() + if(is.null(morphology_image_data)) + dim_image <- apply(image_data[,seq_len(2)], 2, max) + scale_dim <- relative_fov_positions[i,2:3]*dim_image + image_data[,seq_len(2)] <- image_data[,seq_len(2)] + + rep(1, nrow(image_data)) %o% as.matrix(scale_dim)[1,] + morphology_image_data <- rbind(morphology_image_data, image_data) + } + morphology_image_data_array <- reshape2::acast(morphology_image_data, y ~ x) + morphology_image <- magick::image_read(morphology_image_data_array) %>% magick::image_resize("x800") + + # pick a resolution level + morphology_image_info <- image_info(morphology_image) + if(verbose) + message("Image Resolution (X:", morphology_image_info$width, " Y:", morphology_image_info$height, ") \n") + + # increase contrast using EBImage + if(increase.contrast) { + if(verbose) + message("Increasing Contrast \n") + morphology_image <- magick::image_contrast(morphology_image, sharpen = 1) + } + + # write to the same folder + if(verbose) + message("Writing Tiff File \n") + if(is.null(output.path)){ + EBImage::writeImage(magick::as_EBImage(morphology_image), file = file.path, ...) + } else { + EBImage::writeImage(magick::as_EBImage(morphology_image), file = output.file, ...) + } + + return(NULL) +} + +#### +# Spatial Genomics #### +#### + +#### +## GenePs #### +#### + +#' importGenePS +#' +#' Importing GenePS data +#' +#' @param dir.path path to Xenium output folder +#' @param assay_name the assay name of the SR object +#' @param sample_name the name of the sample +#' @param use_image if TRUE, the DAPI image will be used. +#' @param resolution_level the level of resolution within TIFF image. Default: 7 (971x638) +#' @param image_name the image name of the Xenium assay, Default: main +#' @param channel_name the channel name of the image of the Xenium assay, Default: DAPI +#' @param import_molecules if TRUE, molecule assay will be created along with cell assay. +#' @param verbose verbose +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom magick image_read image_info +#' @importFrom utils read.csv +#' @importFrom data.table fread +#' @importFrom ids random_id +#' @importFrom grDevices as.raster +#' @importFrom EBImage as.Image +#' +#' @export +#' +importGenePS <- function (dir.path, assay_name = "GenePS", sample_name = NULL, use_image = TRUE, + resolution_level = 7, image_name = "main", channel_name = "DAPI", import_molecules = FALSE, + verbose = TRUE, ...) +{ + # cell assay + if(verbose) + message("Creating cell level assay ...") + + # raw counts + DataOutputfiles <- list.files(paste0(dir.path, "DataOutput/")) + datafile <- DataOutputfiles[grepl("_cellxgene.csv", DataOutputfiles)] + if(length(datafile) > 0){ + rawdata <- utils::read.csv(file = paste0(dir.path, "DataOutput/", datafile), row.names = NULL) + rawdata <- rawdata[,colnames(rawdata)[!colnames(rawdata) %in% "X"]] + rawdata <- t(rawdata) + } else { + stop("There are no files ending with '_cellxgene.csv' in the path") + } + + # image + if(use_image){ + + # check RBioFormats + if (!requireNamespace('RBioFormats')) + stop("Please install RBioFormats package to extract xml from the ome.tiff file!: BiocManager::install('RBioFormats')") + + # check image + image_file <- paste0(dir.path, "/images/DAPI.tiff") + if(file.exists(image_file)){ + image <- RBioFormats::read.image(image_file, resolution = resolution_level) + image <- EBImage::as.Image(image) + image <- grDevices::as.raster(image) + image <- magick::image_read(image) + } else { + stop("There are no image files in the path") + } + # scale the xenium image instructed by 10x Genomics help page + scaleparam <- ceiling(2^(resolution_level-1)) + } else { + image <- NULL + } + + # coordinates + coordsfile <- DataOutputfiles[grepl("_cellcoordinate.csv", DataOutputfiles)] + if(length(coordsfile) > 0){ + coords <- utils::read.csv(file = paste0(dir.path, "DataOutput/", coordsfile)) + coords <- as.matrix(coords[,c("center_x", "center_y")]) + colnames(coords) <- c("x", "y") + if(use_image) { + coords <- coords/scaleparam + imageinfo <- unlist(magick::image_info(image)[c("height")]) + range_coords <- c(0,imageinfo) + } else { + range_coords <- range(coords[,2]) + } + coords[,2] <- range_coords[2] - coords[,2] + range_coords[1] + } else { + stop("There are no files ending with '_cellcoordinate.csv' in the path") + } + + # create VoltRon object for cells + cell_object <- formVoltRon(rawdata, metadata = NULL, image = image, coords, main.assay = assay_name, + assay.type = "cell", image_name = image_name, main_channel = channel_name, sample_name = sample_name, + feature_name = "RNA", ...) + + # add molecules + if(!import_molecules){ + return(cell_object) + } else { + if(verbose) + message("Creating molecule level assay ...") + transcripts_file <- DataOutputfiles[grepl("_transcriptlocation.csv", DataOutputfiles)] + if(length(transcripts_file) == 0){ + stop("There are no files ending with '_transcriptlocation.csv' in the path") + } else { + + # get subcellur data components + subcellular_data <- data.table::fread(paste0(dir.path, "DataOutput/", transcripts_file)) + subcellular_data$id <- seq_len(nrow(subcellular_data)) + subcellular_data <- subcellular_data[,c("id", colnames(subcellular_data)[!colnames(subcellular_data) %in% "id"]), with = FALSE] + colnames(subcellular_data)[colnames(subcellular_data)=="name"] <- "gene" + + # coordinates + coords <- as.matrix(subcellular_data[,c("x", "y")]) + colnames(coords) <- c("x", "y") + if(use_image){ + coords <- coords/scaleparam + } + coords[,"y"] <- range_coords[2] - coords[,"y"] + range_coords[1] + + # metadata + mol_metadata <- subcellular_data[,colnames(subcellular_data)[!colnames(subcellular_data) %in% c("cell", "x", "y")], with = FALSE] + set.seed(nrow(mol_metadata$id)) + mol_metadata$postfix <- paste0("_", ids::random_id(bytes = 3, use_openssl = FALSE)) + mol_metadata$assay_id <- "Assay1" + mol_metadata[, id:=do.call(paste0,.SD), .SDcols=c("id", "postfix")] + + # coord names + rownames(coords) <- mol_metadata$id + + # create VoltRon assay for molecules + mol_assay <- formAssay(coords = coords, image = image, type = "molecule", main_image = image_name, main_channel = channel_name) + + # merge assays in one section + if(verbose) + message("Merging assays ...") + sample.metadata <- SampleMetadata(cell_object) + object <- addAssay(cell_object, + assay = mol_assay, + metadata = mol_metadata, + assay_name = paste0(assay_name, "_mol"), + sample = sample.metadata["Assay1", "Sample"], + layer = sample.metadata["Assay1", "Layer"]) + + # connectivity + connectivity <- data.table::data.table(transcript_id = mol_metadata$id, cell_id = subcellular_data[["cell"]]) + connectivity <- subset(connectivity, cell_id != 0) + connectivity[["cell_id"]] <- vrSpatialPoints(cell_object)[connectivity[["cell_id"]]] + + # add connectivity of spatial points across assays + object <- addLayerConnectivity(object, + connectivity = connectivity, + sample = sample.metadata["Assay1", "Sample"], + layer = sample.metadata["Assay1", "Layer"]) + + # return + return(object) + } + } +} + +#### +# BGI Genomics #### +#### + +#### +## STOmics #### +#### + +#' importSTOmics +#' +#' Importing STOmics (Stereo-Seq) data +#' +#' @param h5ad.path path to h5ad file of STOmics output +#' @param assay_name the assay name +#' @param sample_name the name of the sample +#' @param image_name the image name of the Visium assay, Default: main +#' @param channel_name the channel name of the image of the Visium assay, Default: H&E +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom methods as +#' @export +#' +importSTOmics <- function(h5ad.path, assay_name = "STOmics", sample_name = NULL, image_name = "main", channel_name = "H&E", ...) +{ + # check package + if(!requireNamespace('anndataR')) + stop("Please install anndataR package!: devtools::install_github('scverse/anndataR')") + + # get h5ad data + stdata <- anndataR::read_h5ad(h5ad.path, to = "HDF5AnnData") + + # observation and feature names + obs_names <- stdata$obs_names + var_names <- stdata$var_names + + # raw counts + rawdata <- Matrix::t(stdata$X) + rownames(rawdata) <- var_names + colnames(rawdata) <- obs_names + rawdata <- methods::as(rawdata, 'CsparseMatrix') + + # metadata + metadata <- stdata$obs + rownames(metadata) <- obs_names + + # coordinates + coords <- stdata$obsm$spatial + rownames(coords) <- obs_names + + # scale coordinates + binsize <- stdata$uns$bin_size + params <- list( + spot.radius = 0.5 + (binsize-1), + vis.spot.radius = 0.5 + (binsize-1)) + + # create VoltRon + 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, feature_name = "RNA", ...) +} + +#### +# Akoya #### +#### + +#### +## PhenoCycler #### +#### + +#' importPhenoCycler +#' +#' Importing PhenoCycler data +#' +#' @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 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 verbose verbose +#' @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'), verbose = TRUE, ...){ + + # 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)){ + if(verbose) + 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, feature_name = "RNA", ...) + + # return + object +} + +#' readPhenoCyclerMat +#' +#' Read and Load Akoya CODEX data, adapted from the \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 +#' +#' @importFrom methods as +#' +#' @noRd +readPhenoCyclerMat <- function( + filename, + type = c('inform', 'processor', 'qupath'), + filter = 'DAPI|Blank|Empty', + inform.quant = c('mean', 'total', 'min', 'max', 'std') +) { + # Check arguments + if (!file.exists(filename)) { + stop("Can't find 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) + mtx <- methods::as(mtx, 'CsparseMatrix') + } + 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 + 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) + expr <- methods::as(expr, 'CsparseMatrix') + } + 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) + mtx <- methods::as(mtx, 'CsparseMatrix') + } + list(matrix = mtx, centroids = centroids, metadata = md) + }, + stop("Unknown matrix type: ", type) + ) + return(outs) +} + +#### +# Non-Commercial #### +#### + +#### +## OpenST #### +#### + +#' importOpenST +#' +#' Importing OpenST data +#' +#' @param h5ad.path path to h5ad file of STOmics output +#' @param assay_name the assay name +#' @param sample_name the name of the sample +#' @param image_name the image name of the Visium assay, Default: main +#' @param channel_name the channel name of the image of the Visium assay, Default: H&E +#' @param verbose verbose +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom methods as +#' @importFrom Matrix t +#' +#' @export +importOpenST <- function(h5ad.path, assay_name = "OpenST", sample_name = NULL, image_name = "main", channel_name = "H&E", verbose = TRUE, ...) +{ + # check package + if(!requireNamespace('anndataR')) + stop("Please install anndataR package!: devtools::install_github('scverse/anndataR')") + + # get h5ad data + stdata <- anndataR::read_h5ad(h5ad.path, to = "HDF5AnnData") + + # observation and feature names + obs_names <- stdata$obs_names + var_names <- stdata$var_names + + # raw counts + rawdata <- Matrix::t(stdata$X) + rownames(rawdata) <- var_names + colnames(rawdata) <- obs_names + rawdata <- methods::as(rawdata, 'CsparseMatrix') + + # metadata + metadata <- stdata$obs + rownames(metadata) <- obs_names + + # coordinates + coords <- stdata$obsm$spatial_3d_aligned + rownames(coords) <- obs_names + zlocation <- unique(coords[,3]) + + # get individual sections as voltron data + sections <- unique(metadata$n_section) + zlocation <- zlocation[order(sections)] + connectivity <- data.frame(Var1 = rep(seq_len(length(sections)), length(sections)), + Var2 = rep(seq_len(length(sections)), each = length(sections))) + sections <- sections[order(sections)] + vr_data_list <- list() + if(verbose) + message("Creating Layers ...") + for(i in seq_len(length(sections))){ + ind <- metadata$n_section == sections[i] + spatialpoints <- rownames(metadata[metadata$n_section == sections[i],]) + cur_data <- rawdata[,spatialpoints] + cur_metadata <- metadata[spatialpoints,] + cur_coords <- coords[ind,c(1,2)] + rownames(cur_coords) <- spatialpoints + vr_data_list[[i]] <- formVoltRon(data = cur_data, metadata = cur_metadata, coords = cur_coords, + main.assay = assay_name, sample_name = paste0("Section", sections[i]), + image_name = image_name, main_channel = channel_name, feature_name = "RNA", ...) + } + + # create VoltRon + sample_name <- ifelse(is.null(sample_name), "Sample", sample_name) + vr_data <- mergeVoltRon(vr_data_list[[1]], vr_data_list[-1], samples = sample_name) + + # set zlocations and adjacency of layer in the vrBlock + vr_data <- addBlockConnectivity(vr_data, connectivity = connectivity, zlocation = zlocation, sample = sample_name) + + # return + vr_data +} + +#### +## DBIT-Seq #### +#### + +#' importDBITSeq +#' +#' Importing DBIT-Seq data +#' +#' @param path.rna path to rna count matrix +#' @param path.prot path to protein count matrix +#' @param size the size of the in situ pixel (defualt is 10 (micron)) +#' @param assay_name the assay name +#' @param sample_name the name of the sample +#' @param image_name the image name of the Visium assay, Default: main +#' @param channel_name the channel name of the image of the Visium assay, Default: H&E +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom utils read.table +#' +#' @export +importDBITSeq <- function(path.rna, path.prot = NULL, size = 10, assay_name = "DBIT-Seq", sample_name = NULL, image_name = "main", channel_name = "H&E", ...) +{ + # count matrix RNA + rnadata <- utils::read.table(path.rna, header = TRUE, sep = "\t", row.names = 1) + rnadata <- t(as.matrix(rnadata)) + + # count matrix Protein + protdata <- utils::read.table(path.prot, header = TRUE, sep = "\t", row.names = 1) + protdata <- t(as.matrix(protdata)) + + # coords + coords <- sapply(colnames(rnadata), function(x) as.numeric(strsplit(x, split = "x")[[1]])) + coords <- t(coords) + colnames(coords) <- c("x", "y") + + # get DBIT-Seq parameters + params <- list( + spot.radius = size/2, + vis.spot.radius = size, + nearestpost.distance = size*2*sqrt(2)) + coords <- coords*size*3/2 + + # make voltron object + object <- formVoltRon(data = rnadata, coords = coords, image = NULL, assay.type = "spot", params = params, image_name = "main", + main.assay = assay_name, sample_name = sample_name, feature_name = "RNA", ...) + + # add protein assay + if(!is.null(path.prot)){ + + # # create protein assay + # new_assay <- formAssay(data = protdata, + # coords = coords, + # type = "spot") + # new_assay@image <- object[["Assay1"]]@image + # sample.metadata <- SampleMetadata(object) + + # # add new assay + # object <- addAssay(object, + # assay = new_assay, + # metadata = Metadata(object), + # assay_name = paste(assay_name, "Prot", sep = "-"), + # sample = sample.metadata["Assay1", "Sample"], + # layer = sample.metadata["Assay1", "Layer"]) + # + # # add connectivity of spatial points across assays + # connectivity <- cbind(vrSpatialPoints(object, assay = "Assay1"), + # vrSpatialPoints(object, assay = "Assay2")) + # object <- addConnectivity(object, + # connectivity = connectivity, + # sample = sample.metadata["Assay1", "Sample"], + # layer = sample.metadata["Assay1", "Layer"]) + + # add negative probe assay as new feature set + object <- addFeature(object, assay = assay_name, data = protdata, feature_name = "Protein") + + } + + # return + return(object) +} + +#### +# Image Data #### +#### + +#' importImageData +#' +#' import an image as VoltRon object +#' +#' @param image a single or a list of image paths or magick-image objects +#' @param tile.size the size of tiles +#' @param segments Either a list of segments or a GeoJSON file. This will result in a second assay in the VoltRon object to be created +#' @param image_name the image name of the Image assay, Default: main +#' @param channel_names the channel names of the images if multiple images are provided +#' @param ... additional parameters passed to \link{formVoltRon} +#' +#' @importFrom magick image_read image_info +#' @importFrom data.table data.table +#' +#' @examples +#' # single image +#' imgfile <- system.file("extdata", "DAPI.tif", package = "VoltRon") +#' vrdata <- importImageData(imgfile, image_name = "main") +#' +#' # multiple images +#' imgfile <- c(system.file("extdata", "DAPI.tif", package = "VoltRon"), +#' system.file("extdata", "DAPI.tif", package = "VoltRon")) +#' vrdata <- importImageData(imgfile, image_name = "main", channel_name = c("DAPI", "DAPI2")) +#' +#' @export +importImageData <- function(image, tile.size = 10, segments = NULL, image_name = "main", channel_names = NULL, ...){ + + # images and channel names + if(!is.null(channel_names)){ + if(length(image) != length(channel_names)) + stop("Provided channel names should of the same length as the images!") + if(any(!is.character(channel_names))) + stop("Invalid channel names!") + } + + # get image + if(!is.list(image)){} + image <- as.list(image) + image <- sapply(image, function(img){ + if(!inherits(img, "magick-image")){ + if(!is.character(img)){ + stop("image should either be a magick-image object or a file.path") + } else{ + if(file.exists(img)){ + img <- magick::image_read(img) + } else { + stop(img, " is not found!") + } + } + } + img + }, USE.NAMES = TRUE, simplify = FALSE) + + # channel names + if(!is.null(channel_names)){ + names(image) <- channel_names + } + + # check image size + imageinfo <- vapply(image, function(img) { + info <- magick::image_info(img) + c(info$width, info$height) + }, numeric(2)) + unique_width <- unique(imageinfo[1,]) + unique_height <- unique(imageinfo[2,]) + if(length(unique_width) == 1 && length(unique_height) == 1){ + imageinfo <- list(width = imageinfo[1,1], height = imageinfo[2,1]) + } + + # coordinates + even_odd_correction <- (!tile.size%%2)*(0.5) + x_coords <- seq((tile.size/2) + even_odd_correction, imageinfo$width, tile.size)[seq_len(imageinfo$width %/% tile.size)] + y_coords <- seq((tile.size/2) + even_odd_correction, imageinfo$height, tile.size)[seq_len(imageinfo$height %/% tile.size)] + y_coords <- rev(y_coords) + coords <- as.matrix(expand.grid(x_coords, y_coords)) + colnames(coords) <- c("x", "y") + rownames(coords) <- paste0("tile", seq_len(nrow(coords))) + + # metadata + metadata <- data.table::data.table(id = rownames(coords)) + + # create voltron object with tiles + object <- formVoltRon(data = NULL, metadata = metadata, image = image, coords, main.assay = "ImageData", assay.type = "tile", params = list(tile.size = tile.size), image_name = image_name, ...) + + # check if segments are defined + if(is.null(segments)){ + return(object) + } else { + + # check if segments are paths + if(inherits(segments, "character")){ + if(grepl(".geojson$", segments)){ + segments <- generateSegments(geojson.file = segments) + } else { + stop("Only lists or GeoJSON files are accepted as segments input!") + } + } + + # make coordinates out of segments + coords <- t(vapply(segments, function(dat){ + apply(dat[,c("x", "y")], 2, mean) + }, numeric(2))) + rownames(coords) <- names(segments) + + # make segment assay + assay <- formAssay(coords = coords, segments = segments, image = image, type = "ROI", main_image = image_name) + + # add segments as assay + sample_metadata <- SampleMetadata(object) + object <- addAssay(object, + assay = assay, + metadata = data.frame(row.names = paste0(names(segments), "_Assay2")), + assay_name = "ROIAnnotation", + sample = sample_metadata$Sample, + layer = sample_metadata$Layer) + + # return + return(object) + } +} + + +#' generateSegments +#' +#' The function to import segments from a geojson file +#' +#' @param geojson.file the GeoJSON file, typically generated by QuPath software +#' +#' @importFrom rjson fromJSON +#' @importFrom dplyr tibble +#' +#' @export +generateSegments <- function(geojson.file){ + + # get segments + if(inherits(geojson.file, "character")){ + if(file.exists(geojson.file)){ + segments <- rjson::fromJSON(file = geojson.file) + } else { + stop("geojson.file doesn't exist!") + } + } else { + stop("geojson.file should be the path to the GeoJSON file!") + } + + # parse polygons as segments/ROIs + segments <- lapply(segments, function(x){ + type <- x$geometry$type + poly <- x$geometry$coordinates + if(grepl("Polygon", type)){ + poly <- as.data.frame(matrix(unlist(poly[[1]]), ncol = 2, byrow = TRUE)) + } + colnames(poly) <- c("x", "y") + dplyr::tibble(poly) + }) + + # attach names to segments + segments <- mapply(function(x, sgt){ + dplyr::tibble(data.frame(id = x, sgt)) + }, seq_len(length(segments)), segments, SIMPLIFY = FALSE) + + # generate ROI names + names(segments) <- paste0("ROI", seq_len(length(segments))) + + # return + return(segments) +} + +#' generateGeoJSON +#' +#' generating geojson files from segments +#' +#' @param segments the segments, typically from \link{vrSegments}. +#' @param file the GeoJSON file, typically to be used by QuPath software. +#' +#' @importFrom rjson fromJSON +#' +#' @export +generateGeoJSON <- function(segments, file){ + + if(!requireNamespace('geojsonR')) + stop("Please install geojsonR package for using geojsonR functions") + + # get segments + if(!inherits(file, "character")){ + stop("file should be the path to the GeoJSON file!") + } + + # reshape segments + segments <- mapply(function(id, sgt){ + poly <- na.omit(as.matrix(sgt[,c("x", "y")])) + poly <- rbind(poly, poly[1,,drop = FALSE]) + poly <- as.list(data.frame(t(poly))) + names(poly) <- NULL + init <- geojsonR::TO_GeoJson$new() + geometry <- init$Polygon(list(poly), stringify = TRUE) + feature <- list(type = "Feature", + id = id, + geometry = geometry[!names(geometry) %in% "json_dump"], + properties = list(objectType = "annotation")) + feature + }, names(segments), segments, SIMPLIFY = FALSE, USE.NAMES = FALSE) + + # save as json + segments <- rjson::toJSON(segments) + write(segments, file = file) +} \ No newline at end of file