一篇新鲜出炉IF10+Q1的空间转录组数据挖掘文章

今天来学习一篇空间转录组数据的挖掘类文章，于2025年3月28号发表在杂志Cell Reports Medicine上，标题为《A spatially resolved transcriptome landscape during thyroid cancer progression》，这是一篇比较典型的空间转录组数据挖掘类的文章。来看看~

实验设计

文章对四种类型的甲状腺组织进行了空间转录组和单细胞转录组测序：

• Para-tumor thyroid (PT) tissue（癌旁甲状腺组织）：与肿瘤组织相邻的正常或非肿瘤性甲状腺组织。这种组织通常用于研究对照；
• Papillary thyroid cancer (PTC)（甲状腺乳头状癌）：甲状腺乳头状癌是最常见的甲状腺恶性肿瘤，占分化型甲状腺癌（DTC）的80%以上。它起源于甲状腺滤泡上皮细胞，具有一定的侵袭性和转移能力；
• Locally advanced PTC (LPTC)（局部晚期甲状腺乳头状癌）：肿瘤已经超出甲状腺包膜，侵犯到周围组织（如气管、食管、颈部软组织等），但尚未发生远处转移；
• Anaplastic thyroid carcinoma (ATC)（甲状腺未分化癌）：一种高度恶性的甲状腺肿瘤，起病迅速，进展快，预后极差。ATC的细胞形态极不规则，细胞核大而多形，核分裂象多见。肿瘤组织常呈浸润性生长，可快速侵犯周围组织和器官。

分析流程

空转数据分析：

• 定量：hg38参考基因组，Space Ranger (version 1.3.1)定量，Seurat降维聚类分群
• 数据过滤：基因 count>=10，至少在3个spots中表达；spots至少有200个基因表达
• 降维聚类分群：SCT标准化，PCA分析，RunPCA，FindNeighbors，FindClusters，and RunUMAP，SpatialDimPlot可视化

Note：这个地方文章中其实还缺少一块，没有提到数据是否进行了整合分析以及是否去批次。文章提供的代码是没有合并在一起分析的，但是文章的空转的图又都是使用同一个cluster编号表示聚类结果（先不管了）

单细胞转录组数据分析：

• 文章中的数据 + GSE184362 + GSE148673 统一进行Cell Ranger定量，Harmony进行整合去批次；
• 细胞过滤：基因数>=200，count数>=1000，线粒体基因比例<=20%；
• DoubletFinder双细胞过滤；
• 细胞注释：已知细胞markers

单细胞与空转整合

这里采用的比较经典的反卷积方法：The RCTD module within the spaceXR software，我们刚刚介绍过RCTD的用法：空间转录组数据注释分析：RCTD反卷积(Nature Biotechnology IF: 33.1)

文章中的反卷积大饼子饼图结果如下：

高级分析：

除了上面的基础分析，作者还对空转数据做了两个高级分析。

第一个是 stage-specific tumor leading edge remodeling，肿瘤边界分析，但是这个分析是借用 有经验的肿瘤病理学家在HE染色切片上手动圈出来的：

作者分析了不同的疾病阶段中肿瘤边界区域中的分子和细胞组成异质性变化。

第二个高级分析问空转细胞通讯：CellphoneDB (v4.1.0)43 and NicheNet。在空转中鉴定到了高可信的互作对：COL8A1-ITHB1 in PTC, LAMB2-ITGB4 in LPTC, and SERPINE1-PLAUR in ATC

此外作者还筛选到了一个特殊的亚群：SERPINE1+ fibroblast，这个亚群与疾病的发展以及预后相关。

看下数据

作者将数据上传到了GEO：.cgi?acc=GSE250521

下载：GSE250521_RAW.tar，使用bash进行整理

# bash命令
# 解压
tar -xvf GSE250521_RAW.tar
# 穿件每个样本的文件夹
ls *gz | perl -ne 'chomp;/(.*)_visium_(.*)/;print"mkdir $1\n";' |uniq
ls *gz | perl -ne 'chomp;/(.*)_visium_(.*)/;print"mkdir $1\n";' |uniq |sh
# 将数据移动到每个文件夹下
ls *gz | perl -ne 'chomp;/(.*)_visium_(.*)/;print"mv $_ $1/$2\n";' |less
ls *gz | perl -ne 'chomp;/(.*)_visium_(.*)/;print"mv $_ $1/$2\n";' |sh
# 解压图片
gunzip  */*tissue_hires_image.png.gz
gunzip  */*scalefactors_json.json.gz

整理后数据格式如下：

.
├── GSM7980860_N-1
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   ├── matrix.mtx.gz
│   ├── scalefactors_json.json
│   ├── tissue_hires_image.png
│   └── tissue_positions_list.csv.gz
├── GSM7980861_N-2
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   ├── matrix.mtx.gz
│   ├── scalefactors_json.json
│   ├── tissue_hires_image.png
│   └── tissue_positions_list.csv.gz
......

读取：

###
### Create: juan zhang
### Blog: /
### Forum:  .html
### Update Log: 2025-04-06   by juan zhang (492482942@qq)
### 

rm(list=ls())
library(ggsci)
library(dplyr) 
library(future)
library(Seurat)
library(clustree)
library(cowplot)
library(data.table)
library(ggplot2)
library(patchwork)
library(stringr)
library(qs)
library(Matrix)

# 创建目录
getwd()
gse <- "GSE250521"
dir.create(gse)

################## 读取数据，有三个文件与tissue_hires_image.png,scalefactors_json.json.gz
samples <- list.dirs("GSE250521/visium/", recursive = F, full.names = F)
samples
scRNAlist <- lapply(samples, function(pro){
#pro <- samples[1]
print(pro)

# 先读取 h5
  paste0("GSE250521/visium/",pro,"/")
  data <- Read10X(data.dir = paste0("GSE250521/visium/",pro,"/"))
  dim(data)
  data[1:5,1:5]
  object <- CreateSeuratObject(counts = data, assay = "Spatial", min.cells = 3, project = pro)
  object

# 再读取
  image <- Read10X_Image(image.dir = paste0("GSE250521/visium/",pro), 
                         image.name = "tissue_hires_image.png",
                         filter.matrix = TRUE)
  image
  dim(image)
  image <- image[Cells(x = object)]
  DefaultAssay(object = image) <- "Spatial"

# 添加图片到object中
  object[[pro]] <- image
  object@images[[1]]@scale.factors$lowres <- object@images[[1]]@scale.factors$hires

# 标准化降维聚类
  object <- SCTransform(object, assay = "Spatial")
  object <- RunPCA(object, assay = "SCT") %>% 
    FindNeighbors(., reduction = "pca", dims = 1:30) %>%
    FindClusters() %>%
    RunUMAP(., reduction = "pca", dims = 1:30)
#SpatialDimPlot(object,pt.size.factor = 1.6,image.alpha = 0.6)

return(object)
})
names(scRNAlist) <-  samples
scRNAlist

# merge
sce.all <- merge(scRNAlist[[1]], y=scRNAlist[-1], add.cell.ids=samples)
sce.all

总共 57997个spots：

An object of class Seurat 
45906 features across 57997 samples within 2 assays 
Active assay: SCT (22149 features, 0 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: Spatial

画个基础spots聚类图看看：

p <- SpatialDimPlot(sce.all,pt.size.factor = 1.9,ncol = 4)
ggsave(filename = "spatial.png", width = 20, height = 18, plot = p)