IF10+空转文献复现（二）：空转聚类结果与空转切片结合

生信技能树

发布于 2025-07-29 10:33:01

11000

代码可运行

运行总次数：0

代码可运行

最近我们的空转单细胞专辑会更新一篇文章复现学习，这篇文章于2023年发表在高分杂志Cancer Res上，文献标题为《Spatial Transcriptomics Depict Ligand–Receptor Cross-talk Heterogeneity at the Tumor-Stroma Interface in Long-Term Ovarian Cancer Survivors》。

今天的学习内容为，将空转的聚类结果与空转切片结合起来，绘制带有位置的聚类图。空转的四张片子如下：

读取前面保存的结果

这里用的就是文章提供的代码：step2_map_seurat_clusters.R，里面进行了一些修改~

这里的对象 seurat_object.RData 在上一篇代码中我没有放保存的代码，保存代码如下：

# save the seurat object
fnm = paste0(save_dir,'/seurat_object.RData')
fnm
save(sample4,sample5,sample10,sample12, file = fnm) # Save the object

这个 seurat_object.RData 在文章提供的代码中也可以下载到。

然后读取上一次保存的数据：

### This code performs the following operations:
### 1. Map all ST clusters on a white background
### 2. Map all ST clusters on the H&E image

####################################################################################
library(Seurat)
library(scales)
library(imager)
####################################################################################

save_dir = '../results'
image_dir = '../images'
fnm = paste0(save_dir,'/seurat_object.RData')
fnm
load(fnm)

# ####### Vinplot of signature genes #######
VlnPlot(sample12, features = c("WFDC2", "COL3A1", "CCL5", "ACTA2", "CD79A", "CD68", "VCAN", "KRT18"),
        pt.size = 0.2, ncol = 4)

并查看了这几个基因在sample12中不同亚群里面的表达情况：

基因名称	功能	主要表达细胞类型
WFDC2	WFDC2编码的蛋白质是一种小的分泌蛋白，具有抗蛋白酶活性和潜在的抗菌功能。它可能参与精子成熟，并且在肺上皮细胞中表达。此外，WFDC2在卵巢癌组织中高表达，并且与肿瘤细胞的迁移和增殖相关。	肺上皮细胞、卵巢癌细胞
COL3A1	COL3A1编码III型胶原蛋白，是细胞外基质的重要组成部分，有助于可伸展结缔组织（如肺、皮肤和血管系统）中纤维化的进展。	成纤维细胞、肿瘤相关成纤维细胞、肿瘤浸润性免疫细胞（如巨噬细胞、T细胞等）
CCL5	CCL5是一种趋化因子，主要通过结合受体（如CCR1、CCR3、CCR5）招募免疫细胞（如T细胞、单核细胞、嗜酸性粒细胞等）到炎症或感染部位。	T细胞、巨噬细胞、内皮细胞、血小板
ACTA2	ACTA2编码平滑肌肌动蛋白α2，是平滑肌细胞和某些成纤维细胞的标记基因。	平滑肌细胞、肌成纤维细胞
CD79A	CD79A是B细胞受体复合物的一部分，参与B细胞的激活和免疫反应。	B细胞
CD68	CD68是巨噬细胞和单核细胞的标记基因。	巨噬细胞、单核细胞
VCAN	VCAN编码软骨糖蛋白，是一种细胞外基质蛋白，参与细胞间信号传导和组织的结构维持。	成纤维细胞、基质细胞
KRT18	KRT18编码角蛋白18，是中间丝蛋白家族的成员，主要在上皮细胞中表达，参与细胞结构的维持和细胞应激反应。	上皮细胞

这里貌似只有一个亚群的特征很明显，毕竟每个spot还是有多个细胞混在一起。

定义两个函数来绘制空转位置切片图：

输出目录：

############# save results to the folder ###############
folder = paste0(save_dir,'/map_seurat_clusters')
folder
if (!file.exists(folder)){
  dir.create(folder)
}

定义函数：绘制空转位置切片图，没有HE背景：

############## plot clusters without background function ##############
plot_clu_nbg <- function(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im){
  tmp = SOE@meta.data$seurat_cluster
  Nclu = length(levels(tmp))
  nspots = length(tmp)
  colorcode = hue_pal()(Nclu)
  colorcode = c('orange','green','blue','red','cyan','magenta','yellow','brown')
  mycols <- c(rep(NA,length(tmp)))
for (i in seq(1,nspots)){
    mycols[i] = colorcode[tmp[i]]
  }

  filenm = paste0(folder,'/A',SID,'noBG.png')
  png(filenm, width = 512, height = 512)
  par(mar = c(0,0,0,0)) # set zero margins on all 4 sides
  plot(x = NULL, y = NULL, xlim = c(0,dim(im)[1]), ylim = c(0,dim(im)[2]), pch = '',
       xaxt = 'n', yaxt = 'n', xlab = '', ylab = '', xaxs = 'i', yaxs = 'i',
       bty = 'n') # plot empty figure
# rasterImage(im, xleft = 0, ybottom = 0, xright = dim(im)[1], ytop = dim(im)[2])
  points(topleft[1]+(d[,1]-1)*delta_x+apply(d,1,function(x){return(adjust_x(x[2]))}),
         topleft[2]-(d[,2]-1)*delta_y+apply(d,1,function(x){return(adjust_y(x[1]))}),
         col=mycols,pch=20,cex=2.8)
  dev.off()
}

定义函数：绘制空转位置切片图，有HE背景：

############## plot clusters with background function ##############
plot_clu <- function(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im){
  tmp = SOE@meta.data$seurat_cluster
  Nclu = length(levels(tmp))
  nspots = length(tmp)
  colorcode = hue_pal()(Nclu)
  colorcode = c('orange','green','blue','red','cyan','magenta','yellow','brown')
  mycols <- c(rep(NA,length(tmp)))
for (i in seq(1,nspots)){
    mycols[i] = colorcode[tmp[i]]
  }

  filenm = paste0(folder,'/A',SID,'.png')
  png(filenm, width = 512, height = 512)
  par(mar = c(0,0,0,0)) # set zero margins on all 4 sides
  plot(x = NULL, y = NULL, xlim = c(0,dim(im)[1]), ylim = c(0,dim(im)[2]), pch = '',
       xaxt = 'n', yaxt = 'n', xlab = '', ylab = '', xaxs = 'i', yaxs = 'i',
       bty = 'n') # plot empty figure
  rasterImage(im, xleft = 0, ybottom = 0, xright = dim(im)[1], ytop = dim(im)[2])
  points(topleft[1]+(d[,1]-1)*delta_x+apply(d,1,function(x){return(adjust_x(x[2]))}),
         topleft[2]-(d[,2]-1)*delta_y+apply(d,1,function(x){return(adjust_y(x[1]))}),
         col=mycols,pch=20,cex=1.7)
  dev.off()
}

绘制每个样本的聚类空间位置图

样本A4：

这里绘图之所以能够将空转的聚类结果与HE切片对应起来是因为空转的表达矩阵的列名就是空间的坐标信息。

############## A4 ###############
SID = '4'
SOE = sample4

fnm = paste0(image_dir,'/Sample_',SID,'_ds.jpg')  # downsample to 0.25x
fnm
im = load.image(fnm)
fnmask = paste0(image_dir,'/Mask_sample_',SID,'.jpg')
fnmask
msk = load.image(fnmask)
plot(im)
# xd和yd：获取图像的宽度和高度。
xd = dim(im)[1]
yd = dim(im)[2]
# 定义图像的四个角点坐标（左上、右上、左下、右下），这些坐标用于后续的空间校正。
topleft = c(58,yd-50)
topright = c(1412,yd-64)
bottomleft = c(48,yd-1512)
bottomright = c(1400,yd-1524)
# 计算空间校正参数
# delta_x和delta_y：计算每个单位步长在x和y方向上的像素变化量。
# adjust_x和adjust_y：定义两个函数，用于根据步长调整x和y坐标，以匹配图像的空间位置。
delta_x = mean(abs(topright[1]-topleft[1]),abs(bottomright[1]-bottomleft[1]))/32
delta_y = mean(abs(topright[2]-bottomright[2]),abs(topleft[2]-bottomleft[2]))/34
adjust_x = function(step){return((step-1)*(bottomleft[1]-topleft[1])/34)}
adjust_y = function(step){return((step-1)*(topright[2]-topleft[2])/32)}
# 提取单细胞数据的空间位置信息
tmp = as.matrix(GetAssayData(object = SOE, assay = "SCT", slot = "data"))
tmp1 = colnames(tmp)
tmp1 = substr(tmp1,4,100)
colnames(tmp) = tmp1
x = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][1]})))
y = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][2]})))
d = data.frame(x,y)
# 绘制聚类结果
## Plot all clusters on top of image ##
plot_clu(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)
plot_clu_nbg(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)

样本A5：

############## A5 ###############
SID = '5'
SOE = sample5

fnm = paste0(image_dir,'/Sample_',SID,'_ds.jpg')  # downsample to 0.25x
fnm
im = load.image(fnm)
fnmask = paste0(image_dir,'/Mask_sample_',SID,'.jpg')
fnmask
msk = load.image(fnmask)
plot(im)
xd = dim(im)[1]
yd = dim(im)[2]
topleft = c(30,yd-12)
topright = c(1382,yd-20)
bottomleft = c(22,yd-1458)
bottomright = c(1368,yd-1464)
delta_x = mean(abs(topright[1]-topleft[1]),abs(bottomright[1]-bottomleft[1]))/32
delta_y = mean(abs(topright[2]-bottomright[2]),abs(topleft[2]-bottomleft[2]))/34
adjust_x = function(step){return((step-1)*(bottomleft[1]-topleft[1])/34)}
adjust_y = function(step){return((step-1)*(topright[2]-topleft[2])/32)}
tmp = as.matrix(GetAssayData(object = SOE, assay = "SCT", slot = "data"))
tmp1 = colnames(tmp)
tmp1 = substr(tmp1,4,100)
colnames(tmp) = tmp1
x = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][1]})))
y = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][2]})))
d = data.frame(x,y)
## Plot all clusters on top of image ##
plot_clu(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)
plot_clu_nbg(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)

样本A10：

############## A10 ###############
SID = '10'
SOE = sample10

fnm = paste0(image_dir,'/Sample_',SID,'_clear.jpg')  # downsample to 0.25x
im = load.image(fnm)
fnmask = paste0(image_dir,'/Mask_sample_',SID,'.jpg')
msk = load.image(fnmask)
plot(im)
xd = dim(im)[1]
yd = dim(im)[2]
topleft = c(18,yd-22)
topright = c(1046,yd-18)
bottomleft = c(18,yd-1116)
bottomright = c(1043,yd-1112)
delta_x = mean(abs(topright[1]-topleft[1]),abs(bottomright[1]-bottomleft[1]))/32
delta_y = mean(abs(topright[2]-bottomright[2]),abs(topleft[2]-bottomleft[2]))/34
adjust_x = function(step){return((step-1)*(bottomleft[1]-topleft[1])/34)}
adjust_y = function(step){return((step-1)*(topright[2]-topleft[2])/32)}
tmp = as.matrix(GetAssayData(object = SOE, assay = "SCT", slot = "data"))
tmp1 = colnames(tmp)
tmp1 = substr(tmp1,5,100)
colnames(tmp) = tmp1
x = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][1]})))
y = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][2]})))
d = data.frame(x,y)
## Plot all clusters on top of image ##
plot_clu(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)
plot_clu_nbg(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)

样本A12：

############## A12 ###############
SID = '12'
SOE = sample12

fnm = paste0(image_dir,'/Sample_',SID,'_clear.jpg')  # downsample to 0.25x
im = load.image(fnm)
fnmask = paste0(image_dir,'/Mask_sample_',SID,'.jpg')
msk = load.image(fnmask)
plot(im)
xd = dim(im)[1]
yd = dim(im)[2]
topleft = c(32,yd-30)
topright = c(984,yd-22)
bottomleft = c(40,yd-1114)
bottomright = c(992,yd-1106)
delta_x = mean(abs(topright[1]-topleft[1]),abs(bottomright[1]-bottomleft[1]))/32
delta_y = mean(abs(topright[2]-bottomright[2]),abs(topleft[2]-bottomleft[2]))/34
adjust_x = function(step){return((step-1)*(bottomleft[1]-topleft[1])/34)}
adjust_y = function(step){return((step-1)*(topright[2]-topleft[2])/32)}
tmp = as.matrix(GetAssayData(object = SOE, assay = "SCT", slot = "data"))
tmp1 = colnames(tmp)
tmp1 = substr(tmp1,5,100)
colnames(tmp) = tmp1
x = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][1]})))
y = as.numeric(unlist(lapply(tmp1, function(x){strsplit(x,"x")[[1]][2]})))
d = data.frame(x,y)
## Plot all clusters on top of image ##
plot_clu(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)
plot_clu_nbg(SID,SOE,d,topleft,adjust_x,adjust_y,delta_x,delta_y,im)