|
| 1 | +#!/usr/bin/env Rscript |
| 2 | +# R/run_rctd.R |
| 3 | +# RCTD cell type deconvolution for 10x Visium |
| 4 | +# Called from Python via rpy2 or standalone |
| 5 | +# Author: Shivani Patel |
| 6 | + |
| 7 | +suppressPackageStartupMessages({ |
| 8 | + library(spacexr) # RCTD |
| 9 | + library(Matrix) |
| 10 | + library(dplyr) |
| 11 | + library(ggplot2) |
| 12 | +}) |
| 13 | + |
| 14 | + |
| 15 | +run_rctd_pipeline <- function(spaceranger_path, reference_path, |
| 16 | + sample_id, outdir) { |
| 17 | + dir.create(outdir, showWarnings = FALSE, recursive = TRUE) |
| 18 | + |
| 19 | + # ── Load Visium spatial data ─────────────────────────────────────────────── |
| 20 | + counts <- load_visium_counts(spaceranger_path) |
| 21 | + coords <- load_visium_coords(spaceranger_path) |
| 22 | + |
| 23 | + # ── Build SpatialRNA object ──────────────────────────────────────────────── |
| 24 | + puck <- SpatialRNA(coords, counts, require.int = FALSE) |
| 25 | + |
| 26 | + # ── Load single-cell reference ───────────────────────────────────────────── |
| 27 | + reference <- load_sc_reference(reference_path) |
| 28 | + |
| 29 | + # ── Run RCTD ────────────────────────────────────────────────────────────── |
| 30 | + message(" Running RCTD deconvolution...") |
| 31 | + myRCTD <- create.RCTD(puck, reference, max_cores = 4) |
| 32 | + myRCTD <- run.RCTD(myRCTD, doublet_mode = "full") |
| 33 | + |
| 34 | + # ── Extract proportions ──────────────────────────────────────────────────── |
| 35 | + results <- myRCTD@results |
| 36 | + norm_weights <- normalize_weights(results$weights) |
| 37 | + proportions <- as.data.frame(as.matrix(norm_weights)) |
| 38 | + |
| 39 | + # ── Save ────────────────────────────────────────────────────────────────── |
| 40 | + write.table(proportions, |
| 41 | + file.path(outdir, paste0(sample_id, "_rctd_proportions.tsv")), |
| 42 | + sep = "\t", quote = FALSE) |
| 43 | + |
| 44 | + # ── Plot spatial proportions ─────────────────────────────────────────────── |
| 45 | + plot_rctd_spatial(proportions, coords, sample_id, outdir) |
| 46 | + |
| 47 | + return(proportions) |
| 48 | +} |
| 49 | + |
| 50 | + |
| 51 | +load_visium_counts <- function(path) { |
| 52 | + h5_file <- file.path(path, "filtered_feature_bc_matrix.h5") |
| 53 | + if (file.exists(h5_file)) { |
| 54 | + return(Seurat::Read10X_h5(h5_file)) |
| 55 | + } |
| 56 | + # Try MEX format |
| 57 | + mex_path <- file.path(path, "filtered_feature_bc_matrix") |
| 58 | + if (dir.exists(mex_path)) { |
| 59 | + return(Seurat::Read10X(mex_path)) |
| 60 | + } |
| 61 | + stop("Cannot find count matrix in: ", path) |
| 62 | +} |
| 63 | + |
| 64 | + |
| 65 | +load_visium_coords <- function(path) { |
| 66 | + coords_file <- file.path(path, "spatial", "tissue_positions_list.csv") |
| 67 | + if (!file.exists(coords_file)) { |
| 68 | + coords_file <- file.path(path, "spatial", "tissue_positions.csv") |
| 69 | + } |
| 70 | + coords <- read.csv(coords_file, header = FALSE) |
| 71 | + colnames(coords) <- c("barcode", "in_tissue", "row", "col", "y", "x") |
| 72 | + coords <- coords[coords$in_tissue == 1, ] |
| 73 | + rownames(coords) <- coords$barcode |
| 74 | + return(coords[, c("x", "y")]) |
| 75 | +} |
| 76 | + |
| 77 | + |
| 78 | +load_sc_reference <- function(ref_path) { |
| 79 | + if (!file.exists(ref_path)) { |
| 80 | + message(" Reference not found — using demo reference") |
| 81 | + return(make_demo_reference()) |
| 82 | + } |
| 83 | + # Load h5ad via reticulate/SeuratDisk |
| 84 | + tryCatch({ |
| 85 | + library(SeuratDisk) |
| 86 | + Convert(ref_path, dest = "h5seurat", overwrite = TRUE) |
| 87 | + seurat_obj <- LoadH5Seurat(gsub("\\.h5ad$", ".h5seurat", ref_path)) |
| 88 | + cell_types <- as.factor(seurat_obj$cell_type) |
| 89 | + counts <- seurat_obj@assays$RNA@counts |
| 90 | + return(Reference(counts, cell_types)) |
| 91 | + }, error = function(e) { |
| 92 | + message(" Reference load failed: ", e$message, " — using demo reference") |
| 93 | + return(make_demo_reference()) |
| 94 | + }) |
| 95 | +} |
| 96 | + |
| 97 | + |
| 98 | +make_demo_reference <- function() { |
| 99 | + set.seed(42) |
| 100 | + n_cells <- 500 |
| 101 | + n_genes <- 2000 |
| 102 | + cell_types <- factor(sample( |
| 103 | + c("Tumor", "T_cell_CD8", "T_cell_CD4", "Macrophage_M1", |
| 104 | + "Macrophage_M2", "CAF", "Endothelial", "B_cell"), |
| 105 | + n_cells, replace = TRUE |
| 106 | + )) |
| 107 | + counts <- matrix( |
| 108 | + rnbinom(n_cells * n_genes, mu = 2, size = 0.5), |
| 109 | + nrow = n_genes, ncol = n_cells, |
| 110 | + dimnames = list(paste0("GENE_", seq_len(n_genes)), |
| 111 | + paste0("CELL_", seq_len(n_cells))) |
| 112 | + ) |
| 113 | + Reference(Matrix::Matrix(counts, sparse = TRUE), cell_types) |
| 114 | +} |
| 115 | + |
| 116 | + |
| 117 | +plot_rctd_spatial <- function(proportions, coords, sample_id, outdir) { |
| 118 | + cell_types <- colnames(proportions) |
| 119 | + shared <- intersect(rownames(proportions), rownames(coords)) |
| 120 | + if (length(shared) == 0) return(invisible(NULL)) |
| 121 | + |
| 122 | + plot_df <- cbind( |
| 123 | + coords[shared, ], |
| 124 | + proportions[shared, ] |
| 125 | + ) |
| 126 | + |
| 127 | + pdf(file.path(outdir, paste0(sample_id, "_rctd_spatial.pdf")), |
| 128 | + width = 4 * min(4, length(cell_types)), |
| 129 | + height = 4 * ceiling(length(cell_types) / 4)) |
| 130 | + |
| 131 | + par(mfrow = c(ceiling(length(cell_types) / 4), min(4, length(cell_types))), |
| 132 | + mar = c(2, 2, 2, 2)) |
| 133 | + |
| 134 | + for (ct in cell_types) { |
| 135 | + vals <- plot_df[[ct]] |
| 136 | + col_scale <- colorRampPalette(c("#f8f9fa", "#e76f51"))(100) |
| 137 | + breaks <- seq(0, max(vals, na.rm = TRUE), length.out = 101) |
| 138 | + col_idx <- findInterval(vals, breaks, rightmost.closed = TRUE) |
| 139 | + |
| 140 | + plot(plot_df$x, plot_df$y, |
| 141 | + col = col_scale[pmax(1, pmin(100, col_idx))], |
| 142 | + pch = 16, cex = 0.5, |
| 143 | + main = gsub("_", " ", ct), |
| 144 | + xlab = "", ylab = "", axes = FALSE) |
| 145 | + } |
| 146 | + dev.off() |
| 147 | +} |
0 commit comments