|
| 1 | +library(spacexr) |
| 2 | +library(Matrix) |
| 3 | +library(SingleCellExperiment) |
| 4 | +library(anndataR) |
| 5 | +library(SPLIT) |
| 6 | +library(Seurat) |
| 7 | +library(scuttle) |
| 8 | + |
| 9 | +## VIASH START |
| 10 | +par <- list( |
| 11 | + "input_spatial_with_cell_types" = "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/spatial_with_celltypes.h5ad", |
| 12 | + "input_scrnaseq_reference"= "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/scrnaseq_reference.h5ad", |
| 13 | + "output" = "task_ist_preprocessing/tmp/split_corrected.h5ad", |
| 14 | + "keep_all_cells" = FALSE, |
| 15 | +) |
| 16 | + |
| 17 | +meta <- list( |
| 18 | + 'cpus': 4, |
| 19 | +) |
| 20 | + |
| 21 | +## VIASH END |
| 22 | + |
| 23 | +# Read the input h5ad file and convert to SingleCellExperiment and Seurat |
| 24 | +sce <- read_h5ad(par$input_spatial_with_cell_types, as = "SingleCellExperiment") |
| 25 | +xe <- read_h5ad(par$input_spatial_with_cell_types, as = "Seurat") |
| 26 | + |
| 27 | +# filter out 0 cells |
| 28 | +if (!par$keep_all_cells) { |
| 29 | + cat("Filtering cells with 0 counts\n") |
| 30 | + sce <- sce[, colSums(counts(sce)) > 0] |
| 31 | + xe <- subset(xe, subset = nCount_RNA > 0) |
| 32 | +} |
| 33 | + |
| 34 | +# Extract spatial coordinates and counts matrix |
| 35 | +centroid_x <- colData(sce)$centroid_x |
| 36 | +centroid_y <- colData(sce)$centroid_y |
| 37 | +coords <- data.frame(centroid_x, centroid_y) |
| 38 | +counts <- assay(sce, "counts") |
| 39 | +rownames(coords) <- colData(sce)$cell_id |
| 40 | +puck <- SpatialRNA(coords, counts) |
| 41 | + |
| 42 | +# Read reference scrnaseq |
| 43 | +ref <- read_h5ad(par$input_scrnaseq_reference, as = "SingleCellExperiment") |
| 44 | + |
| 45 | +#filter reference cell types to those with >25 cells (minimum for RCTD) |
| 46 | +valid_celltypes <- names(table(colData(ref)$cell_type))[table(colData(ref)$cell_type) >= 25] |
| 47 | +filtered_ref <- ref[,colData(ref)$cell_type %in% valid_celltypes] |
| 48 | + |
| 49 | +ref_counts <- assay(filtered_ref, "counts") |
| 50 | +# factor to drop filtered cell types |
| 51 | +colData(filtered_ref)$cell_type <- factor(colData(filtered_ref)$cell_type) |
| 52 | +cell_types <- colData(filtered_ref)$cell_type |
| 53 | +names(cell_types) <- colnames(ref_counts) |
| 54 | +reference <- Reference(ref_counts, cell_types, min_UMI = 0) |
| 55 | + |
| 56 | +# check cores |
| 57 | +cores <- 1 |
| 58 | +if ("cpus" %in% names(meta) && !is.null(meta$cpus)) cores <- meta$cpus |
| 59 | +cat(sprintf("Number of cores: %s\n", cores)) |
| 60 | + |
| 61 | +# Run the algorithm |
| 62 | +cat("Running RCTD\n") |
| 63 | +myRCTD <- create.RCTD(puck, reference, max_cores = cores) |
| 64 | +myRCTD <- run.RCTD(myRCTD, doublet_mode = "doublet") |
| 65 | + |
| 66 | +# Get the "spot_class" annotation from RCTD |
| 67 | +# cat("Saving RCTD spot_class\n") |
| 68 | +# results <- myRCTD@results |
| 69 | +# rctd_spot_class <- results$results_df$spot_class |
| 70 | +# names(rctd_spot_class) <- rownames(results$results_df) |
| 71 | +# colData(sce)$RCTD_class <- "not_included" |
| 72 | +# colData(sce)[names(rctd_spot_class),"RCTD_class"] <- as.character(rctd_spot_class) |
| 73 | + |
| 74 | +# Post-process RCTD output |
| 75 | +RCTD <- SPLIT::run_post_process_RCTD(myRCTD) |
| 76 | + |
| 77 | +# Run SPLIT purification |
| 78 | +cat("Running SPLIT\n") |
| 79 | +res_split <- SPLIT::purify( |
| 80 | + counts = GetAssayData(xe, assay = 'RNA', layer = 'counts'), # or any gene x cells counts matrix |
| 81 | + rctd = RCTD, |
| 82 | + DO_purify_singlets = TRUE # optional |
| 83 | +) |
| 84 | + |
| 85 | + |
| 86 | +# create corrected counts layer in original SingleCell object |
| 87 | +cat("Normalizing counts\n") |
| 88 | + |
| 89 | +# Preserve original normalized values before overwriting with corrected normalization |
| 90 | +assay(sce, "normalized_uncorrected") <- assay(sce, "normalized") |
| 91 | + |
| 92 | +# First copy in counts |
| 93 | +assay(sce, "corrected_counts") <- assay(sce, "counts") |
| 94 | + |
| 95 | +# Then, replace only the updated cells |
| 96 | +assay(sce, "corrected_counts")[rownames(res_split$purified_counts), colnames(res_split$purified_counts)] <- res_split$purified_counts |
| 97 | + |
| 98 | +# Library size normalization - see note in resolVI |
| 99 | +size_factors <- librarySizeFactors(assay(sce, "corrected_counts")) |
| 100 | +assay(sce, "normalized") <- assay(logNormCounts(sce, size_factors=size_factors, assay.type = "corrected_counts"),"logcounts") |
| 101 | + |
| 102 | +# Write the final object to h5ad format |
| 103 | +cat("Writing to h5ad\n") |
| 104 | +dir.create(dirname(par$output), showWarnings = FALSE, recursive = TRUE) |
| 105 | +write_h5ad(sce, par$output, mode = "w") |
0 commit comments