1. 数据准备与Seurat环境搭建大家好我是老张在生物信息分析这个行当里摸爬滚打了十来年尤其喜欢捣鼓各种组学数据。最近几年空间转录组技术也就是咱们常说的“空转”火得不行。它能把基因表达信息精准地定位回组织切片上让我们看到基因在哪个位置“说话”这对于理解肿瘤微环境、胚胎发育这些复杂过程简直是神器。但很多刚入门的朋友拿到数据后看着一堆文件常常会懵这些文件都是啥该先动哪个今天我就用一个真实的单样本数据GSE281978里的GSM8633891手把手带你走一遍用Seurat进行预处理和质量控制的完整流程。我会尽量把每一步的原理和操作都掰开揉碎了讲保证你跟着做一遍就能上手。首先咱们得把数据文件理清楚。从10X Genomics Visium平台下机的数据通常会包含好几个文件每个都有它的使命。我习惯在分析前先建一个独立的项目文件夹比如就叫“0-GSM8633891”然后把所有文件都放进去并且把名字改得清晰易懂。别小看这一步规范的文件管理能让你后续的分析少踩很多坑。这些文件里最核心的是filtered_feature_bc_matrix.h5这是经过初步过滤的基因表达矩阵行是基因列是每个捕获点spot里面的数值就是UMI计数。tissue_positions_list.csv这个文件记录了每个spot的空间坐标信息告诉你它在芯片上的行列位置以及在组织图像上的像素坐标这是空间信息的关键。scalefactors_json.json文件存放了缩放因子它像一把尺子负责把spot的坐标和高分辨率的组织图像对齐。最后tissue_hires_image.png是高分辨率的HE染色组织切片图我们最终要把基因表达信息叠加在这张图上进行可视化。把这几个文件准备好咱们的分析就有了坚实的基础。接下来是R环境的搭建。我强烈建议使用RStudio来操作界面友好调试方便。首先我们需要安装并加载几个核心的R包。Seurat是主角它是处理单细胞和空间转录组数据的瑞士军刀。ggplot2用于绘图hdf5r包是为了读取HDF5格式的表达矩阵文件。有时候数据量大计算慢我们还可以用future包来设置并行计算加快速度。你可以用下面这几行代码来搞定环境准备。我习惯在脚本开头先用rm(listls())清空一下工作环境避免之前运行的变量干扰本次分析。# 清空环境确保分析纯净 rm(listls()) # 安装并加载必要的R包 # install.packages(c(Seurat, ggplot2, hdf5r, future)) library(Seurat) library(ggplot2) library(hdf5r) library(future) # 设置工作目录到你的数据文件夹 setwd(path/to/your/0-GSM8633891) # 可选设置并行计算以提升大样本处理速度 plan(multisession, workers 4) # 根据你的电脑核心数调整环境搭好了数据也备齐了咱们就可以正式开始把数据读进R里创建Seurat对象了。这是万里长征的第一步也是把原始数据转化为可分析对象的关键一步。记住预处理阶段的目标是构建一个干净、可靠的数据对象为后续的降维、聚类和差异分析打下坚实基础。很多后续分析出现奇怪结果追根溯源往往是预处理没做好所以咱们在这一步一定要耐心、仔细。2. 创建Seurat对象与数据初探数据读入是正式分析的起点。对于空间数据我们需要分两步走先读基因表达矩阵再读空间坐标和图像信息。首先使用Read10X_h5()函数读取那个核心的H5文件。读进来之后我习惯先用dim()函数看一眼数据的维度心里有个数知道有多少个基因和多少个spot。然后用CreateSeuratObject()函数创建Seurat对象。这里有个关键参数min.cells我一般设为3意思是只保留在至少3个spot中表达的基因。这个过滤非常温和主要是为了剔除那些在极少数spot里偶然出现的、可能由测序错误产生的基因信号能有效减少数据噪音又不会损失太多真实生物信号。# 1. 读取基因表达矩阵 data - Read10X_h5(./GSM8633891/filtered_feature_bc_matrix.h5) dim(data) # 查看原始矩阵维度例如 [1] 36601 1187 # 2. 创建Seurat对象 object - CreateSeuratObject(counts data, assay Spatial, # 指定为空间分析 min.cells 3, # 过滤在少于3个spot中表达的基因 project GSM8633891) # 给项目起个名 object # 查看对象摘要信息创建好对象后控制台会打印出对象的摘要比如“20050 features across 1187 samples”这意味着经过初步过滤我们还剩下20050个基因和1187个spot。这个对象现在只包含了表达矩阵。接下来我们需要把“空间”属性加进去。这就要用到Read10X_Image()函数来读取空间图像信息。这里我通常使用高分辨率图像tissue_hires_image.png这样最后可视化出来的图更清晰。filter.matrixTRUE这个参数很重要它会自动确保图像信息只包含那些在表达矩阵中存在的spot避免数据错位。# 3. 读取空间图像信息 image - Read10X_Image(image.dir ./GSM8633891/, image.name tissue_hires_image.png, filter.matrix TRUE) image # 查看图像信息 # 4. 确保图像信息与表达矩阵的spot一一对应 image - image[Cells(x object)] # 5. 设置图像的默认assay为“Spatial” DefaultAssay(object image) - Spatial # 6. 将图像信息添加到Seurat对象中 object[[GSM8633891]] - image # 7. 统一缩放因子避免后续可视化报错 # 有时高、低分辨率图像的缩放因子不同这里统一成高分辨率的 objectimages[[GSM8633891]]scale.factors$lowres - objectimages[[GSM8633891]]scale.factors$hires # 给对象起个简单的别名方便后续调用 sce.all - object对象创建完成后我强烈建议花几分钟时间探索一下数据结构做到心中有数。你可以查看一下表达矩阵的前几行前几列看看计数长什么样。更重要的是查看对象的元数据sce.allmeta.data这里面自动计算了两个非常重要的质控指标nCount_Spatial每个spot的总UMI数和nFeature_Spatial每个spot检测到的基因数。这两个指标是后续判断spot质量好坏的核心依据。你还可以用table(sce.all$orig.ident)看看样本信息确认只有一个样本。最后我习惯用qs包的qsave()函数把当前这个完整的对象保存下来。qs格式的读写速度比R原生的rds快很多尤其是对象很大的时候能节省不少等待时间。# 探索数据结构 # 查看表达矩阵片段 as.data.frame(LayerData(sce.all, assaySpatial, layercounts)[1:5, 1:5]) # 查看元数据的前几行 head(sce.allmeta.data) # 查看样本分布 table(sce.all$orig.ident) # 保存对象方便下次直接加载避免重复读取 library(qs) qsave(sce.all, file ./GSM8633891/sce.all.qs)3. 数据质控与可视化质控Quality Control, QC是预处理中最关键的一环目的是把低质量的spot和不可靠的基因信号过滤掉留下高质量的数据进行下游分析。如果质控不严格那些低质量的spot比如只覆盖了很少细胞或处于组织边缘捕获效率低的spot会像噪音一样干扰后续的聚类和差异分析甚至导致完全错误的生物学结论。对于空间数据质控主要围绕两个核心指标展开每个spot的总UMI数nCount_Spatial和每个spot检测到的基因数nFeature_Spatial。首先我们需要直观地看看这些指标的分布情况。这里我推荐两个必看的图。第一个是 violin plot小提琴图用VlnPlot()函数绘制。它能同时展示所有spot的nCount和nFeature的分布范围、中位数以及离散程度。你可以一眼看出数据整体质量如何有没有特别离谱的离群点。通常高质量的数据这两个指标的分布会比较集中有一个明显的峰。# 绘制UMI数和基因数的小提琴图 plot1 - VlnPlot(sce.all, features c(nCount_Spatial, nFeature_Spatial), pt.size 0.1, # 点的大小0.1可以避免点过于密集 ncol 2) # 并排显示两个图 NoLegend() # 不显示图例让图更简洁 print(plot1) ggsave(filename ./1-QC/VlnPlot.pdf, width 12, height 6, plot plot1)光看整体分布还不够我们更需要知道这些“好”的spot和“坏”的spot在组织切片上具体位于什么位置。这就是第二个关键图空间特征图Spatial Feature Plot。使用SpatialFeaturePlot()函数我们可以把nCount_Spatial或nFeature_Spatial的数值映射到组织切片上用颜色深浅来表示数值高低。这张图能直观地揭示空间技术本身可能引入的偏差。比如你可能会发现组织边缘的spot普遍UMI数较低颜色浅这可能是因为切片边缘捕获效率下降或者某些区域因为组织折叠、气泡等原因导致信号异常高或异常低。这些空间上的系统性偏差必须在质控时被识别出来。# 绘制总UMI数的空间分布图 plot2 - SpatialFeaturePlot(sce.all, features nCount_Spatial, pt.size.factor 1.5) # 调整点的大小以适应图像 print(plot2) ggsave(filename ./1-QC/FeaturePlot_nCount.pdf, width 10, height 8, plot plot2) # 绘制检测基因数的空间分布图 plot3 - SpatialFeaturePlot(sce.all, features nFeature_Spatial, pt.size.factor 1.5) print(plot3) ggsave(filename ./1-QC/FeaturePlot_nFeature.pdf, width 10, height 8, plot plot3)看完图心里就有谱了。接下来就是设定阈值进行过滤。过滤没有绝对的金标准需要结合数据分布、图像观察以及生物学常识来定。一个常用的起始阈值是过滤掉nFeature_Spatial小于200的spot。因为一个spot如果只检测到很少的基因很可能它没有覆盖到完整的细胞或者处于无组织的背景区域。我们可以用subset()函数轻松实现过滤。这里要特别注意过滤的松紧程度会影响结果。过滤太狠可能会损失有生物学意义的稀有细胞群体过滤太松则会引入大量噪音。我建议先从保守的阈值开始完成初步分析后如果发现聚类结果很混乱再回头考虑是否要收紧过滤标准。# 基于nFeature_Spatial过滤低质量spot # 这里设定阈值 nFeature_Spatial 200 sce.all - subset(sce.all, subset nFeature_Spatial 200) # 过滤后再次查看对象摘要确认spot数量减少了 sce.all除了这两个核心指标有时我们还会关注线粒体基因MT-开头的占比。在单细胞数据中高线粒体基因占比通常意味着细胞处于应激或凋亡状态。在空间数据中虽然一个spot包含多个细胞但异常高的线粒体占比也可能指示该区域细胞状态不佳或存在技术假象。你可以用PercentageFeatureSet()函数来计算这个比例。不过对于Visium数据由于spot较大线粒体基因的解读需要更谨慎我通常把它作为一个辅助参考指标不会单独用它做硬性过滤。# 可选计算线粒体基因比例 # 首先识别线粒体基因 mito_genes - grep(^MT-, rownames(sce.all), ignore.case TRUE, value TRUE) sce.all - PercentageFeatureSet(sce.all, features mito_genes, col.name percent.mito) # 可以绘制线粒体比例的空间分布图看看 SpatialFeaturePlot(sce.all, features percent.mito, pt.size.factor 1.5)4. 数据标准化与基础降维聚类质控之后我们得到了一个相对干净的数据集。但原始的表达矩阵还不能直接用于下游分析因为其中存在两大干扰一是测序深度差异有的spot捕获的RNA总量就是多有的就是少二是基因表达量的分布严重右偏少数高表达基因占据了大部分计数。数据标准化就是为了消除这些技术噪音让基因表达在不同spot之间具有可比性。在Seurat中目前最推荐的方法是SCTransform。SCTransform是一个强大的标准化方法它基于负二项分布模型同时处理了测序深度归一化、方差稳定化并且能有效识别高变基因。它比传统的LogNormalize方法表现更好尤其是在处理具有大量零值的单细胞或空间数据时。运行SCTransform后数据会被存储在一个名为“SCT”的新assay中。# 使用SCTransform进行标准化和特征选择 sce.all - SCTransform(sce.all, assay Spatial, # 指定对哪个assay进行操作 verbose TRUE) # 显示运行过程信息标准化完成后我们就可以进行降维了。降维的目的是把成千上万个基因的信息压缩到几十个能代表数据主要变异方向的主成分PC中。这大大简化了数据复杂度并且去除了噪音。首先用RunPCA()函数进行主成分分析。这里我们使用标准化后的“SCT”assay。运行后可以用ElbowPlot()函数查看碎石图帮助决定使用多少个主成分PC进行后续分析。通常选择拐点elbow之前的PCs。# 在主成分分析中使用SCT assay sce.all - RunPCA(sce.all, assay SCT, verbose FALSE) # 绘制碎石图辅助决定PC数量 ElbowPlot(sce.all, ndims 50)基于选定的PCs例如前30个我们可以进行细胞邻域图构建和聚类。FindNeighbors()函数基于降维后的空间这里是PCA空间计算spot之间的相似性构建一个K近邻图。FindClusters()函数则在这个图上进行社区发现将相似的spot归为同一个cluster。resolution参数控制聚类的粒度值越大得到的cluster数量越多划分越细。初学者可以从0.4到1.2之间尝试几个值看看聚类结果是否稳定且有生物学意义。# 基于PCA结果构建邻域图并进行聚类 sce.all - FindNeighbors(sce.all, reduction pca, dims 1:30) # 使用前30个主成分 sce.all - FindClusters(sce.all, resolution 0.8, # 聚类分辨率可调整 verbose FALSE) # 使用UMAP进行非线性降维可视化 sce.all - RunUMAP(sce.all, reduction pca, dims 1:30)现在我们可以把降维和聚类的结果可视化出来这是检验我们预处理和标准化效果的第一步。通常我们会并排看两张图一张是UMAP图展示spot在二维空间的分布和聚类情况另一张是空间聚类图展示这些cluster在原始组织切片上的空间位置。两者结合能告诉我们聚类结果是否具有空间特异性。比如如果某个cluster在UMAP上是一团在空间图上也是连续的一片组织区域那这个结果就非常理想很可能对应了某种特定的组织学结构或细胞类型。# 绘制UMAP聚类图 p1 - DimPlot(sce.all, reduction umap, label TRUE, # 在cluster中心显示标签 label.size 6) ggtitle(UMAP of spots) # 绘制空间聚类图 p2 - SpatialDimPlot(sce.all, label TRUE, label.size 3, pt.size.factor 1.8, # 调整点大小 image.alpha 0.6) # 设置组织图像透明度 ggtitle(Spatial distribution of clusters) # 将两张图并排显示 library(patchwork) combined_plot - p1 | p2 print(combined_plot) ggsave(filename ./1-QC/Cluster_Visualization.pdf, width 16, height 7, plot combined_plot)5. 空间差异表达与特征识别经过聚类我们得到了spot的分组。接下来一个很自然的问题就是这些不同的空间区域cluster在基因表达上有什么不同这就是差异表达分析。我们可以用FindMarkers()函数来比较任意两个cluster之间哪些基因的表达有显著差异。比如我们想看看cluster 5和cluster 6有什么不同。# 寻找cluster 5 vs cluster 6的差异表达基因 de_markers - FindMarkers(sce.all, ident.1 5, # 第一个cluster ident.2 6) # 第二个cluster如果为NULL则与其余所有cluster比较 # 查看结果按调整后的p值p_val_adj升序排列最重要的基因排在最前 head(de_markers, n 10) # 将排名前三的差异基因在空间上可视化 top_genes - rownames(de_markers)[1:3] SpatialFeaturePlot(sce.all, features top_genes, alpha c(0.1, 1), # 设置颜色透明度范围 ncol 3) # 一行显示三个基因除了这种有监督的基于已知聚类差异分析空间转录组还有一个独特的分析角度寻找那些本身表达模式就具有空间变异性的基因而不依赖于我们先前的聚类结果。这能帮助我们发现一些新的、可能定义未知生物边界的基因。Seurat提供了FindSpatiallyVariableFeatures()函数来实现这一功能它内置了多种算法比如“markvariogram”基于地理统计学变差函数和“moransi”基于莫兰空间自相关指数。我实测下来“moransi”方法对许多生物学数据都表现不错。# 寻找具有空间变异模式的基因 # 这里从高变基因中选取前1000个进行分析以加快计算速度 sce.all - FindSpatiallyVariableFeatures(sce.all, assay SCT, # 使用标准化后的数据 features VariableFeatures(sce.all)[1:1000], selection.method moransi, verbose TRUE) # 提取空间变异性最强的基因 top_spatial_features - head(SpatiallyVariableFeatures(sce.all, selection.method moransi), 6) # 可视化这些基因的空间表达模式 SpatialFeaturePlot(sce.all, features top_spatial_features, ncol 3, alpha c(0.1, 1), pt.size.factor 1.5)运行完这段代码你会得到一批基因。这些基因的表达不是随机分布的而是在组织切片上呈现出明显的梯度、斑块或边界模式。例如你可能发现一个基因在肿瘤核心区域高表达在边缘区域低表达或者沿着某个组织结构呈条带状表达。这些基因往往是理解组织空间功能分化的关键线索。你可以把它们作为候选标志物去查阅文献看看它们是否与已知的生物学过程或细胞类型相关。这一步常常能带来意想不到的发现也是空转数据分析最迷人的地方之一——直接从数据中挖掘空间规律。6. 基于参考数据的细胞类型映射Visium技术的一个特点是每个spot的直径有55微米通常会捕获到多个细胞的RNA混合物。因此一个spot的转录组信号是其中所有细胞类型的“平均值”。为了解析每个spot内具体的细胞类型组成我们可以借助已有的单细胞RNA测序scRNA-seq数据作为参考进行细胞类型映射Cell Type Mapping或去卷积Deconvolution。这就像有一本已知细胞类型的“密码本”我们去解读每个spot的“混合密码”。Seurat提供了一个非常强大的整合与映射流程核心思想是寻找参考数据集单细胞和查询数据集空间之间的“锚点”anchors然后通过锚点将参考数据集的细胞类型标签概率性地转移到空间数据上。首先我们需要准备一个高质量的、经过良好注释的单细胞参考数据集。这个参考数据应该与研究样本的组织类型尽可能匹配。# 假设我们已经有一个准备好的单细胞参考Seurat对象 ref # ref对象的metadata中应包含一个名为‘celltype’的列存储了细胞类型注释 library(qs) ref - qread(./0-reference/sce_celltype.qs) # 加载参考数据 # 可视化参考数据的细胞类型 DimPlot(ref, group.by celltype, label TRUE, repel TRUE)接下来我们需要对参考数据和我们的空间数据进行“对话前的统一语言”处理即使用相同的标准化方法。这里我们双方都使用SCTransform。为了加快参考数据的处理速度可以设置ncells参数进行下采样这通常不会损失太多信息。# 对参考数据进行SCTransform处理可下采样加速 ref - SCTransform(ref, ncells 3000, verbose FALSE) %% RunPCA(verbose FALSE) %% RunUMAP(dims 1:30) # 确保空间数据也使用SCT assay如果之前没做需重新运行SCTransform # sce.all - SCTransform(sce.all, assay Spatial, verbose FALSE) %% RunPCA(verboseFALSE)现在寻找锚点并进行标签转移。FindTransferAnchors()函数是关键它会在两个数据集间寻找对应的基因表达模式对。TransferData()函数则执行实际的转移为空间数据中的每个spot计算属于各种参考细胞类型的概率分数。# 寻找参考数据和空间数据之间的锚点 anchors - FindTransferAnchors(reference ref, query sce.all, normalization.method SCT, # 使用SCT标准化后的数据 dims 1:30, # 用于寻找锚点的PCA维度 reference.reduction pca) # 指定参考数据的降维方式 # 进行细胞类型预测 predictions - TransferData(anchorset anchors, refdata ref$celltype, # 参考数据的标签列 dims 1:30, weight.reduction sce.all[[pca]]) # 将预测结果添加到空间数据的metadata中 sce.all - AddMetaData(sce.all, metadata predictions)预测完成后我们可以进行可视化。最直接的是查看某个spot最主要的预测细胞类型。我们可以把预测的“主导类型”作为spot的一个新分组用SpatialDimPlot()画出来看看不同的细胞类型在空间上是如何分布的。更精细的做法是我们可以查看特定细胞类型如“T细胞”或“成纤维细胞”的预测概率在空间上的分布。概率高的区域说明该spot中含有较多这类细胞。# 1. 查看每个spot被预测最多的细胞类型first.type sce.all$predicted.id - predictions$predicted.id p1 - SpatialDimPlot(sce.all, group.by predicted.id, label TRUE, label.size 3, pt.size.factor 1.5) # 2. 可视化特定细胞类型的预测概率分布 DefaultAssay(sce.all) - SCT # 切换回SCT assay进行特征绘图 p2 - SpatialFeaturePlot(sce.all, features c(T, fibroblasts), # 假设参考数据中有这些细胞类型 pt.size.factor 1.5, ncol 2) p1 | p2这种映射方法为我们理解空间数据提供了强大的注释工具。你可以看到肿瘤区域是否富集了特定的免疫细胞或者血管周围有哪些类型的基质细胞。这极大地增强了空间数据的解释能力将抽象的基因表达模式与具体的细胞生物学联系了起来。当然这个方法的效果高度依赖于参考数据的质量以及与查询数据的匹配度在实际分析中可能需要尝试不同的参考数据集或调整参数。