参考教程链接：
- https://github.com/dmcable/spacexr
- 官方教程： https://github.com/dmcable/spacexr/tree/master/vignettes
- 文献解读：10X单细胞空间联合分析之十（RCTD） https://www.jianshu.com/p/92a6df0dcb1b
RCTD可以将单细胞类型或细胞类型混合分配到空间转录组spots上。RCTD 有三种模式：
- doublet mode：它为每个Spot分配1-2种细胞类型，推荐用于高空间分辨率的技术，如Slide-seq和MERFISH
- full mode：它为每个Spot分配任意数量的细胞类型，推荐用于空间分辨率较低的技术，如100微米分辨率的Visium
- multi mode：doublet mode的扩展，可以在每个点上发现两种以上的细胞类型，作为full mode的替代选项
每种模态都有一个对应的教程
今天先来学习Doublet mode: spatial transcriptomics vignette. Also, most other vignettes use doublet mode.
包安装
# install.packages("devtools")
devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)
library(spacexr)
library(Matrix)
Introduction
稳健细胞类型分解(Robust Cell Type Decomposition，简称RCTD)是一种从空间转录组数据中学习细胞类型的统计方法。
在本次示例中，我们将为小脑Slide-seq数据集反卷积注释细胞类型。教程使用带注释的snRNA-seq小脑数据集定义细胞类型。
Note：参考也可以是单细胞数据集或细胞类型特定的bulkRNA-seq数据集。
数据预处理
RCTD需要两个数据：带有注释的单细胞转录组参考数据，需要反卷积注释的空间转录组数据。
reference数据
首先是单细胞参考数据集。reference使用RCTD包中Reference函数，这个函数需要三个参数：
- counts：为矩阵或dgCmatrix对象。行名是基因，列名代表barcode/细胞名。counts应该是未转换的原始counts数据
- cell_types：为带有细胞命名的(通过细胞barcode)细胞类型因子。因子的“水平”是可能的细胞类型标识。
- nUMI：可选，具有barcode的每个像素中总counts or UMI列表。如果没有提供，nUMI将被假定为出现在每个像素上的总数。
reference可能来自各种数据类型，但需要将其加载到R环境中。
在本教程中，我们的reference作为两个csv文件存储在' reference /Vignette '文件夹中:
- meta_data.csv：CSV文件(有3列，标题为“barcode”、“cluster”和“nUMI”)，包含每个细胞的nUMI和单元格类型分配。
- dge.csv：标准10x格式的基因表达CSV文件。
### Load in/preprocess your data, this might vary based on your file type
# directory for the reference
refdir <- system.file("extdata",'Reference/Vignette',package = 'spacexr')
# load in counts matrix
counts <- read.csv(file.path(refdir,"dge.csv"))
# Move first column to rownames
rownames(counts) <- counts[,1]
counts[,1] <- NULL
counts[1:4,1:4]
# load in meta_data (barcodes, clusters, and nUMI)
meta_data <- read.csv(file.path(refdir,"meta_data.csv"))
head(meta_data)
# create cell_types named list
cell_types <- meta_data$cluster
names(cell_types) <- meta_data$barcode
# convert to factor data type
cell_types <- as.factor(cell_types)
head(cell_types)
# create nUMI named list
nUMI <- meta_data$nUMI
names(nUMI) <- meta_data$barcode
head(nUMI)
### Create the Reference object
reference <- Reference(counts, cell_types, nUMI)
str(reference)
## Examine reference object (optional)
#observe Digital Gene Expression matrix
print(dim(reference@counts))
#> [1] 384 475
#number of occurences for each cell type
table(reference@cell_types)
## Save RDS object (optional)
saveRDS(reference, file.path(refdir,'SCRef.rds'))
至此，reference已经构建好。
空间转录组数据
接下来，加载空间数据为SpatialRNA对象。SpatialRNA函数需要三个参数：
- coords：数据框或者矩阵对象，为空间像素坐标，行名为每个像素的barcode，列为x和y
- counts：矩阵或者dgCmatrix对象，行名是基因，列名是每个像素的barcode，counts值没有进行转换。
- nUMI：可选
构建SpatialRNA有不同的方法，本教程为两个文件：
- BeadLocationsForR.csv：CSV文件，三列，列名为"barcodes”, “xcoord”, and “ycoord”，包含每个像素的空间坐标；
- MappedDGEForR.csv：每个像素的表达矩阵文件
##=============================== 读取空间数据
# directory for sample Slide-seq dataset
datadir <- system.file("extdata",'SpatialRNA/Vignette',package = 'spacexr')
counts <- read.csv(file.path(datadir,"MappedDGEForR.csv"))
rownames(counts) <- counts[,1]; counts[,1] <- NULL
counts[1:4,1:4]
coords <- read.csv(file.path(datadir,"BeadLocationsForR.csv"))
rownames(coords) <- coords$barcodes
coords$barcodes <- NULL
head(coords)
# In this case, total counts per pixel is nUMI
nUMI <- colSums(counts)
head(nUMI)
### Create SpatialRNA object
puck <- SpatialRNA(coords, counts, nUMI)
str(puck)
### Create SpatialRNA object
puck <- SpatialRNA(coords, counts, nUMI)
str(puck)
## Examine SpatialRNA object (optional)
print(dim(puck@counts))
# histogram of log_2 nUMI
hist(log(puck@nUMI,2))
结果图：可以看到空间数据的每个像素中UMI的数值分布。
image-20221111222523662.png
print(head(puck@coords)) # start of coordinate data.frame
barcodes <- colnames(puck@counts) # pixels to be used (a list of barcode names).
# This list can be restricted if you want to crop the puck e.g.
# puck <- restrict_puck(puck, barcodes) provides a basic plot of the nUMI of each pixel
# on the plot:
p <- plot_puck_continuous(puck, barcodes, puck@nUMI, ylimit = c(0,round(quantile(puck@nUMI,0.9))), title ='plot of nUMI')
结果图：使用x，y坐标展示每个像素中的UMI分布
image-20221111222947585.png
创建RCTD对象
使用create.RCTD函数创建。以下为几个有用的参数：
- max_cores：对于并行处理，使用的线程数
- gene_cutoff, fc_cutoff, gene_cutoff_reg, fc_cutoff_reg：差异表达基因筛选阈值
- UMI_min, UMI_max：每个像素点的阈值筛选
运行RCTD
可以使用run.RCTD函数，这个函数等价于：fitBulk,choose_sigma_c, andfitPixels这三个函数。doublet_mode可以选择上面提到的三种反卷积模式。
myRCTD <- create.RCTD(puck, reference, max_cores = 10)
myRCTD <- run.RCTD(myRCTD, doublet_mode = 'doublet')
str(myRCTD)
myRCTD结构：
image-20221112182253019.png
RCTD结果
RCTD结果保存在myRCTD对象中的@results中
- @results$weights：是一个数据框，为full mode模式下每一个spots（即像素）中细胞类型的权重
- @results$results_df：doublet_mode的结果
- spot_class：因子类型，RCTD的分类，singlet表示每个spot有一种细胞类型，doublet_certain表示每个spot有两种细胞类型，doublet_uncertain表示spot上有2种细胞类型，但只确定有1种，reject表示没有预测到结果。
- first_type：因子类型，预测到的第一种细胞类型
- second_type：因子类型，预测到的第一种细胞类型
- @results$weights_doublet：doublet_mode下每种细胞类型的权重
results <- myRCTD@results
# normalize the cell type proportions to sum to 1.
norm_weights = normalize_weights(results$weights)
#list of cell type names
cell_type_names <- myRCTD@cell_type_info$info[[2]]
spatialRNA <- myRCTD@spatialRNA
## you may change this to a more accessible directory on your computer.
resultsdir <- 'RCTD_Plots'
dir.create(resultsdir)
results$results_df结果示例如下：
image-20221112184204705.png
接下来对结果可视化绘图：
# make the plots
# 绘制full_mode模式下每种细胞类型的可信权重 (saved as
# 'results/cell_type_weights_unthreshold.pdf')
plot_weights(cell_type_names, spatialRNA, resultsdir, norm_weights)
# 绘制full_mode模式下每种细胞类型的权重
# 这里每种细胞类型一幅图，点表示空间上的一个像素或者spot，颜色为权重 (saved as
# 'results/cell_type_weights.pdf')
plot_weights_unthreshold(cell_type_names, spatialRNA, resultsdir, norm_weights)
# 绘制full_mode模式下每种细胞类型预测到的spots数 (saved as
# 'results/cell_type_occur.pdf')
plot_cond_occur(cell_type_names, resultsdir, norm_weights, spatialRNA)
# 绘制doublet_mode模式下每种细胞类型的权重 (saved as
# 'results/cell_type_weights_doublets.pdf')
plot_weights_doublet(cell_type_names, spatialRNA, resultsdir, results$weights_doublet, results$results_df)
其中，绘制full_mode模式下每种细胞类型预测到的spots数结果图如下：
从这里可以看得出，空间数据预测出来的主要为细胞类型10，其次是18。
image-20221112192742832.png
# 所有细胞类型的map结果 (saved as
# 'results/all_cell_types.pdf')
plot_all_cell_types(results$results_df, spatialRNA@coords, cell_type_names, resultsdir)
所有细胞类型的map结果图：
image-20221112194400757.png
绘制doublets：
# doublets
# obtain a dataframe of only doublets
doublets <- results$results_df[results$results_df$spot_class == "doublet_certain",]
# 绘制所有的doublets (saved as
# 'results/all_doublets.pdf')
plot_doublets(spatialRNA, doublets, resultsdir, cell_type_names)
# 对每种细胞类型绘制doublets (saved as
# 'results/all_doublets_type.pdf')
plot_doublets_type(spatialRNA, doublets, resultsdir, cell_type_names)
# a table of frequency of doublet pairs
doub_occur <- table(doublets$second_type, doublets$first_type)
# 绘制doublet出现的堆积柱状图 (saved as
# 'results/doublet_stacked_bar.pdf')
plot_doub_occur_stack(doub_occur, resultsdir, cell_type_names)
doublet_stacked_bar结果图：
image-20221112194946431.png
感觉这个软件的结果图不是很好看，可调整度也不高，我觉得跟教程的数据有关系，后面重新找个数据看看。
教程的空间数据是slide-seq的数据：
img
找了个10x Visium的文献来看看：
ref1：Cell2location maps fine-grained cell types in spatial transcriptomics ( https://doi.org/10.1038/s41587-021-01139-4 )
image-20221112222920831.png
