CellTrek发表于2022年的Nature Biotechnology，题为《Spatial charting of single-cell transcriptomes in tissues》。CellTrek可以结合单细胞和空间转录组数据准确地定位组织内单个细胞的位置，并构建空间细胞图谱。gitHub在https://github.com/navinlabcode/CellTrek
一. CellTrek算法的工作流程
图1:CellTrek工作流程概述。
- CellTrek首先将ST和scRNA-seq的原始矩阵合并，接着将合并后的矩阵进行共降维。然后提取其中的空间转录组数据构建一个多元随机森林模型(RF)，其中空间坐标是结果，潜在特征是预测因子；
- 对ST数据进行二维空间插值以增强ST的spots。然后，将训练好的RF模型应用于共嵌数据（ST插值），以生成RF距离矩阵，该矩阵将使用最近邻对（MNN）转换为稀疏图；
- 最后基于稀疏图，将单细胞转录组数据套入RF模型中，构建Spot-Cell表达相似性矩阵，为每个单细胞结果添加空间坐标信息并实现可视化作图。
众所周知，10X空转数据的一个spots不是单细胞分辨率的，而是包含十余个细胞的混合点。而CellTrek算法可以充分利用单细胞转录组数据，将单细胞信息映射至空间转录组切片上 。值得注意的是，不同于之前的去卷积方法，CellTrek可以实现空转数据的单细胞分辨率。
二. CellTrek代码实现
CellTrek安装比较简单：
代码语言：javascript复制library(devtools)
install_github("navinlabcode/CellTrek")
R包加载：
代码语言：javascript复制library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
library(readr)
library(CellTrek)
library(viridis)
library(ConsensusClusterPlus)
library(SeuratObject)
Step1.加载示例数据
代码语言：javascript复制brain_st_cortex <- read_rds("./Rawdata/brain_st_cortex.rds")
brain_sc <- read_rds("./Rawdata/brain_sc.rds")
## Rename the cells/spots with syntactically valid names
brain_st_cortex <- RenameCells(brain_st_cortex, new.names=make.names(Cells(brain_st_cortex)))
brain_sc <- RenameCells(brain_sc, new.names=make.names(Cells(brain_sc)))
可视化空转数据：
代码语言：javascript复制## Visualize the ST data
SpatialDimPlot(brain_st_cortex)
image-20231105144510038
可视化单细胞数据：
代码语言：javascript复制## Visualize the scRNA-seq data
DimPlot(brain_sc, label = T, label.size = 4.5)
image-20231105144526583
Step2. 使用CellTrek进行细胞映射
代码语言：javascript复制#### 2. 使用 CellTrek 进行细胞映射
# We first co-embed ST and scRNA-seq datasets using traint
brain_traint <- CellTrek::traint(st_data=brain_st_cortex,
sc_data=brain_sc,
sc_assay='RNA',
cell_names='cell_type')
## 我们可以检查共嵌入的结果，以查看这两种数据模态之间是否存在重叠。
DimPlot(brain_traint, group.by = "type")
image-20231105144622926
代码语言：javascript复制# 在共嵌入之后，我们可以将单个细胞映射到它们的空间位置。
# 在这里，我们使用非线性插值（intp = T，intp_lin=F）方法来增强ST的spots
brain_celltrek <- CellTrek::celltrek(st_sc_int=brain_traint, int_assay='traint', sc_data=brain_sc, sc_assay = 'RNA',
reduction='pca', intp=T, intp_pnt=5000, intp_lin=F, nPCs=30, ntree=1000,
dist_thresh=0.55, top_spot=5, spot_n=5, repel_r=20, repel_iter=20, keep_model=T)$celltrek
# 在细胞映射完成后，我们可以使用 celltrek_vis 交互式可视化 CellTrek 的结果
brain_celltrek$cell_type <- factor(brain_celltrek$cell_type,
levels=sort(unique(brain_celltrek$cell_type)))
CellTrek::celltrek_vis(brain_celltrek@meta.data %>%
dplyr::select(coord_x, coord_y, cell_type:id_new),
brain_celltrek@images$anterior1@image,
brain_celltrek@images$anterior1@scale.factors$lowres)
Step3. CellTrek细胞共定位分析
代码语言：javascript复制#### 3. 细胞共定位分析
# 基于 CellTrek 的结果，我们可以使用 SColoc 总结不同细胞类型之间的共定位模式。
# 在这里，我们以谷氨酸能神经元细胞类型作为示例（建议删除一些细胞类型，例如，n<20，细胞数量非常少的细胞类型）。
# 我们首先从我们的映射结果中子集化谷氨酸能神经元细胞类型。
glut_cell <- c('L2/3 IT', 'L4', 'L5 IT', 'L5 PT', 'NP', 'L6 IT', 'L6 CT', 'L6b')
names(glut_cell) <- make.names(glut_cell)
brain_celltrek_glut <- subset(brain_celltrek, subset=cell_type %in% glut_cell)
#然后使用scoloc进行共定位分析：
brain_celltrek_glut$cell_type <- factor(brain_celltrek_glut$cell_type,
levels=glut_cell)
## 我们从图中提取最小生成树（MST）的结果。
brain_sgraph_KL_mst_cons <- brain_sgraph_KL$mst_cons
rownames(brain_sgraph_KL_mst_cons) <- colnames(brain_sgraph_KL_mst_cons) <- glut_cell[colnames(brain_sgraph_KL_mst_cons)]
## 然后，我们提取meta.data数据（包括细胞类型及其频率信息）。
brain_cell_class <- brain_celltrek@meta.data %>% dplyr::select(id=cell_type) %>% unique
brain_celltrek_count <- data.frame(freq = table(brain_celltrek$cell_type))
brain_cell_class_new <- merge(brain_cell_class, brain_celltrek_count, by.x ="id", by.y = "freq.Var1")
brain_sgraph_KL <- CellTrek::scoloc(brain_celltrek_glut,
col_cell='cell_type',
use_method='KL', eps=1e-50)
CellTrek::scoloc_vis(brain_sgraph_KL_mst_cons,
meta_data=brain_cell_class_new )
image-20231105145004137
Step4. 目标细胞类型的空间加权共表达分析
代码语言：javascript复制# 基于CellTrek的结果，我们可以使用SCoexp模块进一步研究感兴趣的细胞类型中的共表达模式。
# 在这里，我们将使用共识聚类(CC)方法以L5 IT单元为例。首先从图表结果中提取L5 IT单元格。
brain_celltrek_l5 <- subset(brain_celltrek, subset=cell_type=='L5 IT')
brain_celltrek_l5@assays$RNA@scale.data <- matrix(NA, 1, 1)
brain_celltrek_l5$cluster <- gsub('L5 IT VISp ', '', brain_celltrek_l5$cluster)
DimPlot(brain_celltrek_l5, group.by = 'cluster')
# 我们选择前2000个可变基因(不包括线粒体、核糖体和高零基因)
brain_celltrek_l5 <- FindVariableFeatures(brain_celltrek_l5)
vst_df <- brain_celltrek_l5@assays$RNA@meta.features %>% data.frame %>% mutate(id=rownames(.))
nz_test <- apply(as.matrix(brain_celltrek_l5[['RNA']]@data), 1, function(x) mean(x!=0)*100)
hz_gene <- names(nz_test)[nz_test<20]
mt_gene <- grep('^Mt-', rownames(brain_celltrek_l5), value=T)
rp_gene <- grep('^Rpl|^Rps', rownames(brain_celltrek_l5), value=T)
vst_df <- vst_df %>% dplyr::filter(!(id %in% c(mt_gene, rp_gene, hz_gene))) %>% arrange(., -vst.variance.standardized)
feature_temp <- vst_df$id[1:2000]
# 我们使用 scoexp 进行空间加权基因共表达分析。
brain_celltrek_l5_scoexp_res_cc <- CellTrek::scoexp(celltrek_inp=brain_celltrek_l5,
assay='RNA',
approach='cc',
gene_select = feature_temp,
sigm=140,
avg_cor_min=.4,
zero_cutoff=3,
min_gen=40, max_gen=400)
# 我们可以使用热图可视化共表达式模块。
brain_celltrek_l5_k <- rbind(data.frame(gene=c(brain_celltrek_l5_scoexp_res_cc$gs[[1]]), G='K1'),
data.frame(gene=c(brain_celltrek_l5_scoexp_res_cc$gs[[2]]), G='K2')) %>%
magrittr::set_rownames(.$gene) %>% dplyr::select(-1)
pheatmap::pheatmap(brain_celltrek_l5_scoexp_res_cc$wcor[rownames(brain_celltrek_l5_k), rownames(brain_celltrek_l5_k)],
clustering_method='ward.D2', annotation_row=brain_celltrek_l5_k, show_rownames=F, show_colnames=F,
treeheight_row=10, treeheight_col=10, annotation_legend = T, fontsize=8,
color=viridis(10), main='L5 IT spatial co-expression')
image-20231105151101908
代码语言：javascript复制# 我们确定了3个不同的基因模块。基于我们识别的共表达模块，可以计算模块的分数。
brain_celltrek_l5 <- AddModuleScore(brain_celltrek_l5,
features=brain_celltrek_l5_scoexp_res_cc$gs,
name='CC_', nbin=10, ctrl=50, seed=42)
## 可视化1
FeaturePlot(brain_celltrek_l5,
grep('CC_', colnames(brain_celltrek_l5@meta.data),
value=T), ncol = 1)
image-20231105151156465
代码语言：javascript复制## 可视化2
SpatialFeaturePlot(brain_celltrek_l5,
grep('CC_', colnames(brain_celltrek_l5@meta.data), value=T))
image-20231105151243508
三. 彩蛋
鉴于有些读者抨击我们只会“抄英文教程”，“照搬英文帖”，这里我基于对CellTrek的理解复现下图E（这是官方教程里没有的）。
image-20231105151822596
图E计算了目标细胞类型（L2/3T）与其余细胞类型在空间上的距离，结果表明L4细胞亚群在空间上距离L2/3T最近（见图C-E，图D在上文复现过了）：
第一步，获取细胞类型和空间坐标信息：
代码语言：javascript复制# coord_x和coord_y是空间坐标：
inp_df <- brain_celltrek_glut@meta.data %>% dplyr::select(cell_names = dplyr::one_of('cell_type'),
coord_x, coord_y)
inp_df$coord_x = 270-inp_df$coord_x
head(inp_df)
image-20231105164435672
感兴趣的小伙伴可以理解一下为啥需要“270-inp_df$coord_x”。
第二步，计算目标细胞与其他细胞的聚类：
代码语言：javascript复制output <- kdist(inp_df = inp_df,
ref = "L2/3 IT", #目标细胞类型
ref_type = 'all',
que = glut_cell, #其余的细胞
k = 10,
new_name = "L23ITvs.Others",
keep_nn = F)
head(output$kdist_df) #最后的结果
# L23ITvs.Others
# F2S4_160422_009_G01 59.72424
# F2S4_160422_011_G01 62.78655
# F2S4_160428_007_D01 58.06266
# F2S4_160516_020_D01 53.70867
# F2S4_180110_072_C01 376.52573
# F2S4_161129_006_F01 757.52189
res = output$kdist_df
res$barcode = row.names(res)
inp_df$barcode = row.names(inp_df)
res = left_join(res, inp_df)
head(res)
# L23ITvs.Others barcode cell_names coord_x coord_y
# 1 59.72424 F2S4_160422_009_G01 L2/3 IT -861.4217 5370.832
# 2 62.78655 F2S4_160422_011_G01 L2/3 IT -903.3036 5301.647
# 3 58.06266 F2S4_160428_007_D01 L2/3 IT -837.6159 5328.688
# 4 53.70867 F2S4_160516_020_D01 L2/3 IT -877.2844 5331.978
# 5 376.52573 F2S4_180110_072_C01 L5 IT -2085.6755 3122.758
# 6 757.52189 F2S4_161129_006_F01 L2/3 IT -4477.5370 2710.916
第三步，可视化：
代码语言：javascript复制library(ggpubr)
ggboxplot(data = res,
x = "cell_names",
y = "L23ITvs.Others",
fill = "cell_names",
title = "K-distance to L2/3 IT cells")
stat_compare_means(method = "kruskal.test")
theme(plot.title = element_text(color="black",hjust = 0.5),
axis.text.x = element_text(angle = 90, hjust = 1,vjust = 0.5), #,vjust = 0.5
legend.position = "none") labs(y = "Distance")
image-20231105165306137
基本和原文一摸一样。
我再额外辩解一下：
- 其实“英文教程”，本身对于一小部分小伙伴来说就有一些困难，阅读上的困难其次，主要还是代码理解上；
- 第二个小问题是，有些“英文教程”有大量的bug，我们生信技能树的小伙伴在写帖子的时候，肯定会测试一遍教程代码，有bug会帮忙解决，然后把解决方案写进教程里，要知道有的bug部分小伙伴花几天甚至几周都没办法解决，所以“抄英文教程”本身也不会是一个简单的事情；
- 第三个就是，绝大部分的教程，我们都会加入自己的理解和知识体系的整理串联，有不对的地方也很正常，我们非常欢迎友好的探讨和质疑，希望受到别人的质疑来提升自己，本身也是我分享笔记的一个初衷；
- 最后，写帖子真的不挣钱...我们基本都是用爱发电的（说实话，我写帖子一般要写大半天，有的甚至一两天，有这功夫，我都能挣不少钱...）。所以还请珍惜我们技能树的小伙伴...
- END -
