作者，Evil Genius

关于RCTD，本来都不打算更新了，R版本的联合分析我觉得大家自己写写就完了，现在看来，还是需要整理一下。

文章在Robust decomposition of cell type mixtures in spatial transcriptomics | Nature Biotechnology

参考链接在GitHub - dmcable/spacexr: Spatial-eXpression-R: Cell type identification (including cell type mixtures) and cell type-specific differential expression for spatial transcriptomics

至于引用该方法的文章，那就很多了。

RCTD有三种模式:

（1）doublet mode:每个spot分配1-2种细胞类型，推荐用于具有高空间分辨率的技术，如Stereo-seq、HD等（注意bin的大小）;

（2）full mode:每个spot分配任意数量的细胞类型，推荐用于具有低空间分辨率的技术，如Visium;

（3）multi mode : doublet mode的扩展，可以每个spot发现两个以上的细胞类型，作为全模式的替代选项。

很多博主、公司推文写了非常多的介绍和代码示例，大家可以借鉴一下。

不过对于这些网络写手，我更喜欢网络侠客这个称呼。

官网的教程写了很多，针对不同精度的平台都有，列举一下：

• Doublet mode: spatial transcriptomics vignette. Also, most other vignettes use doublet mode.
• Doublet mode on MERFISH: MERFISH nonparametric vignette.
• Full mode: full mode on Visium hippocampus
• Multi mode: multi mode on Visium hippocampus

不过官网分析的结果是真的丑

Seurat官网也借鉴了RCTD的方法，在Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurat

示例代码如下，采用了Doublet mode模式

library(spacexr)

# set up reference
ref <- readRDS("../data/mouse_hippocampus_reference.rds")
ref <- UpdateSeuratObject(ref)
Idents(ref) <- "celltype"

# extract information to pass to the RCTD Reference function
counts <- ref[["RNA"]]$counts
cluster <- as.factor(ref$celltype)
names(cluster) <- colnames(ref)
nUMI <- ref$nCount_RNA
names(nUMI) <- colnames(ref)
reference <- Reference(counts, cluster, nUMI)

# set up query with the RCTD function SpatialRNA
slide.seq <- SeuratData::LoadData("ssHippo")
counts <- slide.seq[["Spatial"]]$counts
coords <- GetTissueCoordinates(slide.seq)
colnames(coords) <- c("x", "y")
coords[is.na(colnames(coords))] <- NULL
query <- SpatialRNA(coords, counts, colSums(counts))


RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD <- run.RCTD(RCTD, doublet_mode = "doublet")
slide.seq <- AddMetaData(slide.seq, metadata = RCTD@results$results_df)

RCTD目前大多数用在高精度的平台，比如Stereo-seq的bin20、bin30的情况，HD的8um情况等。

但是对于真正的文章，一般都要根据自己的课题进行修改。

下方是一个文章分析代码示例

# Load required libraries
library(spacexr)
library(Seurat)
library(ggplot2)
library(dplyr)
library(tidyr)
library(pheatmap)
library(progeny)

# Data Preprocessing

# Load the Spatial data
pdac <- readRDS("/path/to/pdac_most_updated.rds")

# Load the Single cell data
sc <- readRDS("/path/to/sc_rctd.rds")

# Prepare data for RCTD
counts <- sc@assays$RNA@counts
Idents(sc) <- "SCT"
cluster <- as.factor(sc$celltype_nicheDE)
names(cluster) <- colnames(sc)
nUMI <- sc$nCount_RNA
names(nUMI) <- colnames(sc)

# Create the reference object
reference <- Reference(counts, cluster, nUMI)

# Prepare spatial transcriptomics data
counts <- pdac@assays$Spatial@counts

coordinates_list <- lapply(image_names, function(image_name) {
  pos <- GetTissueCoordinates(pdac, image = image_name)
  colnames(pos) <- c('x','y')
  return(pos)
})

coords <- do.call(rbind, coordinates_list)
rownames(coords) <- gsub("^.+\\.", "", rownames(coords))

# Create SpatialRNA object
query <- SpatialRNA(coords, counts, colSums(counts))

# Run RCTD
RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD.full <- run.RCTD(RCTD, doublet_mode = "full")

# Add RCTD results to Seurat object
pdac <- AddMetaData(pdac, metadata = RCTD.full@results$results_df)

# Normalize weights
weights <- RCTD.full@results$weights
norm_weights <- normalize_weights(weights)

# Add RCTD results as a new assay
pdac[["rctd_full"]] <- CreateAssayObject(data = t(as.matrix(norm_weights)))
if (length(pdac@assays$rctd_full@key) == 0) {
    pdac@assays$rctd_full@key <- "rctd_full_"
}

这样就获取了RCTD的解卷积空间细胞矩阵，当然了，通常分析到这里还没结束，我们需要继续分析空间细胞聚类，共定位等分析内容。

我们把RCTD、空间细胞聚类、共定位分析一起封装起来，注意要适用各种空间平台，写好参数的设定。跟cell2location一样的，最好匹配的样本联合分析。

#! usr/bin/R
### zhaoyunfei
### 20241111

suppressMessages({
library(Seurat)
library(compositions)
library(tidyverse)
library(clustree)
library(patchwork)
library(argparse)
library(robustbase)
library(ISCHIA)
library(factoextra)
library(dplyr)
library(scran)
library(ggplot2)
library(spacexr)
library(cluster)
library(showtext)
library(gridExtra)
library(pdftools)
})

parser = ArgumentParser()
parser$add_argument("--sc_rds", help="the sc data",required = T)
parser$add_argument("--spatial_rds", help="the sp data",required = T)
parser$add_argument("--sample", help="the sample name",required = T)
parser$add_argument("--outdir", help="the outdir",default = './')
parser$add_argument("--celltype", help="the annotation for celltype")
parser$add_argument("--normalization_method",default = 'SCT')
parser$add_argument("--reduction",help='Dimensional reduction to perform when finding anchors',choices = c('pcaproject','cca'),default = 'pcaproject')
parser$add_argument("--mode"，help='mode for RCTD',default = 'doublet',choices = c('doublet','full','multi'))
parser$add_argument("--cell_interest", help="the interest of cell type,eg : 'T,B:Fobro,epi'")
args <- parser$parse_args()
str(args)

sc_rds = args$sc_rds
spatial_rds = args$spatial_rds
outdir = args$outdir
celltype = args$celltype
normalization_method = args$normalization_method
reduction = args$reduction
sample = args$sample
mode = args$mode
cell_interest = args$cell_interest

接下来就是脚本主体