Stereo-seq空间转录组数据处理全流程：从FASTQ到Seurat对象

拿到一批真实的Stereo-seq原始数据时，很多人的第一反应可能跟我一样：道理我都懂，空间转录组嘛，但到底从哪个文件下手？10x Visium的Space Ranger教程遍地都是，轮到华大的Stereo-seq，官方文档和社区经验都相对零散，尤其很少有文章能把“从FASTQ一路走到Seurat对象”这条路完完整整铺开讲清楚。我前后处理过好几批不同来源、不同版本SAW产出的Stereo-seq数据，这里就把整条链路的实际操作和踩坑经历整理成文，给准备入坑或者正在被数据折腾的同学做一个可以直接照做的参考。

这篇文章的整体思路是尽量贴近真实项目流程，段落也会按实际工作顺序展开：先理解Stereo-seq的数据结构，再搭好分析环境，然后用SAW或备用方案跑出表达矩阵，最后把矩阵转成Seurat对象并完成简单的质量检查。代码和命令行都会给到，版本相关的坑也会单独列出来。

1. Stereo-seq数据处理前必须懂的几个基础概念

1.1 从DNB到CID：空间坐标到底是怎么编码的

Stereo-seq是华大基于DNB（DNA纳米球）技术开发的高分辨率空间转录组平台。它的核心是芯片上密集排布的DNA纳米球阵列，每一个纳米球表面都连接着一条带空间信息的探针。这条探针由三段序列构成：一小段坐标条形码（CID，Coordinate ID）、一段用于去除PCR重复的UMI，以及一段polyT序列。组织切片贴在芯片上后，经过透化处理，细胞释放出的mRNA会被polyT捕获，随后反转录为cDNA，最终测序得到的数据里就同时包含了转录本序列和它的空间位置信息。

这个设计逻辑和10x Visium完全不同。Visium是离散的55微米Spot，每个Spot之间是物理隔离的；而Stereo-seq的DNB直径只有两百多纳米，在芯片上是连续排布的，因此每一个单独DNB代表一个极微小的空间单位。我们通常不直接拿单个DNB做分析，因为捕获效率有限、每个DNB上能捕获到的分子太少，而是要经过一个“聚合”操作，把一定范围内的DNB合并成更大的分析单元，这就是Stereo-seq数据分析里随处可见的“Bin”概念。

CID条码本质上就是坐标的编码映射。芯片制造时，每个DNB的位置和它上面的CID序列是一一对应的；分析软件拿到读段中的CID后，需要查一张映射表才能把CID还原成芯片上的X/Y坐标。这也是为什么Stereo-seq数据处理通常不能脱离芯片信息文件（BC文件或配置文件）来单独完成。

1.2 原始FASTQ的文件结构和常见命名规律

从测序仪下机后拿到的通常是FASTQ文件。Stereo-seq的文库结构决定了R1和R2各自承担不同的任务：R1主要包含CID和UMI序列，用于空间定位和分子去重；R2则是真正的转录本序列，需要比对到参考基因组。

实际操作中要特别注意FASTQ的读长配置和CID长度是否匹配。不同批次的芯片、不同版本SAW对应的CID长度可能不一样，比如有的方案是25bp CID加10bp UMI，有的方案采用更长或更短的CID设计。如果R1读长被截断导致CID信息不完整，那整个空间定位都会出问题，这种错误通常会在比对和定量阶段表现为大量reads无法分配坐标。

另外，FASTQ的文件命名也不是随便看看就完事。分析前务必确认清楚样本名、芯片编号、lane编号这些信息，因为后续SAW流程里很多参数配置都要跟这些信息保持一致。我遇到过一位同学把两个不同芯片的FASTQ混放在同一个目录下，结果SAW直接报错，排查了大半天才发现是样本目录结构不规范导致的。

1.3 “Bin”的聚合逻辑与下游分析尺度的关系

Bin是Stereo-seq文库分析里绕不开的概念。所谓bin 20，就是把20乘以20个相邻DNB的捕获信号合并成一个分析单元，bin 50则是50乘以50个DNB合并。选择合适的bin大小，直接决定了你后续看到的数据分辨率和数据量。

从我的使用经验来看：

• bin 20通常用于单细胞级别的精细分析，理论上更接近细胞大小，但数据稀疏性很高，很多基因的表达量会变成0，导致下游聚类噪声偏大。
• bin 50是很多项目的默认选择，信号聚合后比较均衡，既能看出组织结构，又不会把细胞混在一起太多。
• bin 100甚至更大则偏向组织区域级别的粗略分析，适合做组织分区的初步探索。

如果后续要用CellBin方案做细胞分割，那就会得到真正意义上的“细胞级”表达矩阵，而不是固定大小的网格单元。SAW流程中CellBin是可选的，但它对组织成像质量、细胞密度和算法参数都很敏感，后面章节里会专门讲。

2. 分析环境与工具选型：到底该用哪条路

2.1 SAW、Stereopy和自写流程怎么选

目前将Stereo-seq原始数据处理成表达矩阵，主要有三条路径：使用华大官方发布的SAW分析流程，使用Stereopy等Python工具包，以及自己写脚本拼接完整流程。

我的建议是，如果数据量不是特别庞大，且目标只是得到一份可用的表达矩阵，直接使用SAW是最稳妥的选择。SAW是华大官方为其平台定制的一站式流程，封装了从FASTQ比对、CID解码、UMI去重，到生成GEM/GEF矩阵文件的完整模块。它虽然不是开源软件，但针对自家数据做了大量适配，跑起来最不容易出幺蛾子。

Stereopy更适合在SAW完成粗加工之后做二次处理，比如读取GEF文件、计算bin表达矩阵、做初步质控、转换格式等。它本身也能读取FASTQ，但如果你指望完全用Stereopy替代SAW从零处理原始数据，会遇到不少文档缺失的问题，毕竟它的核心定位更接近“处理和分析Stere-seq表达数据的Python库”，而不像SAW是一个完整的pipeline。

自己写流程的好处是灵活、可控，坏处是容易出错。原理解密看起来不复杂：用STAR比对R2得到基因注释、从R1解析CID和UMI、再做一个“基因-坐标-UMI”的聚合。但实际执行时会遇到大量细节问题，比如比对后reads的标签管理、参考基因组版本与GTF的匹配、CID映射表的格式解析等。我建议新手先走通SAW再考虑DIY，除非你有充分理由自定义分析。

2.2 必备软件清单和版本核对

根据我自己跑过的流程，以下软件是必须准备或者大概率能用上的：

版本筛选是这个阶段最容易被忽略的问题。SAW本身在持续更新，不同版本产出的文件格式有差异，比如早期版本输出GEM文本，后续版本输出GEF二进制，而Stereopy对GEF格式的支持也依赖特定版本，所以“SAW跑完结果Stereopy读不了”这种问题真不算罕见。

我的经验是，在开跑之前先把软件版本写在一个固定的环境配置文件里，包括SAW版本、STAR版本、参考基因组的Ensembl版本，以及Stereopy和Seurat的版本。这能避免很多莫名其妙的问题。如果数据是由测序服务方提供的，建议直接查看服务说明书里的推荐环境版本，尽量保持一致。

2.3 参考基因组与注释文件准备的细节

大部分转录组分析的第一步就是准备参考基因组索引，Stereo-seq也不例外。需要注意，STAR索引的构建结果跟参考基因组的版本强相关，如果上游数据的比对参数是用Ensembl版本做的，下游基因注释也必须对应同一版本，否则会出现基因ID对不上、基因名大量缺失等问题。

实战中最常见的一个坑：分析人员拿到的参考基因组是Ensembl的GTF，基因注释是ENSG开头的ID，但后期做差异分析或可视化时想用symbol名称，于是打开了另一个版本号的GPL或GFF文件。结果就是部分基因ID匹配不上，下游结果莫名其妙少了一大堆基因。我现在的习惯是，在启动分析前就确定注释来源，并且在转换Seurat对象后单独保存一份基因ID与symbol的对照表，以备后期随时查。

如果是人类或小鼠数据，建议直接用对应版本的Ensembl primary assembly加上相同版本的GTF，NCBI的RefSeq注释也可以，但一定不要混着用。其它物种则要优先参考官方推荐，或者自行用Cell Ranger类似流程的基因组构建方式。

3. 主流程实战：从FASTQ到表达矩阵是怎么一步步走通的

3.1 第一步：拿到原始数据后先做文件清点

这个步骤看起来没什么技术含量，但能省下很多后面排查的时间。拿到原始数据时，我通常会做以下几件事：

• 确认FASTQ文件是否完整，是否包含R1和R2两个文件对。
• 用seqkit stats或fastqc快速检查测序质量、GC含量、接头污染情况。
• 确认R1读长是否覆盖完整的CID+UMI，R2读长是否足够用于比对。
• 检查FASTQ中reads数量与预期是否一致，排除因为下机异常导致的文件截断问题。
• 确认芯片配置信息（芯片编号、Bin分割规则等）是否与SAW要求的输入一致。

这个阶段如果发现问题，还能及时联系数据生产方重新提供数据或补充说明。一旦进入SAW流程后才发现数据有问题，那才是真正的浪费时间。

3.2 第二步：SAW流程的核心模块和各阶段意义

以SAW v4.x版本为例，它的执行路径大致分为几个阶段。首先是数据清洗和格式转换，把原始FASTQ整理成流程内部需要的中间格式，同时会做基础的质量过滤。然后是序列比对，将转录本序列比对到参考基因组，这一步通常使用STAR完成，产出BAM文件。这一步本身和普通RNA-seq的比对差别不大，但关键在于后续处理。

接下来是CID解码和坐标映射，这也是Stereo-seq和普通转录组分析分道扬镳的地方。流程会把R1中的CID序列与芯片映射表做比对，给每一条read赋予一个精确的空间坐标；再结合比对到基因组的基因信息，形成一个“基因-坐标-UMI”的初始计数列表。最后通过UMI去重得到唯一的分子计数，并按照不同bin大小聚合成标准化的表达矩阵。

SAW跑完会输出多种格式的文件。如果配置了CellBin，还会输出细胞级别的矩阵和分割结果。这一阶段重点检查几个指标：比对率、基因检出数量、测序饱和度、以及有效read占比。举个例子，如果比对率低于60%，可能是参考基因组或建库有问题；如果用来解析CID的read比例很低，那大概率是CID长度和读长配置不匹配。

3.3 第三步：认识GEM和GEF这两种核心文件格式

GEM和GEF是SAW产出的两种基因表达矩阵格式。GEM本质上是文本表格，每一行记录一个基因在某个坐标上的UMI计数，主要字段包括基因名称、X坐标、Y坐标和UMI数量。GEF则是GEM的二进制封装，读取速度更快，还可直接支持多层bin（bin 1、bin 10、bin 20等）数据的快速切片，Stereopy就是基于GEF格式设计的。

这里要特别提醒：GEM的数据量和坐标范围可能非常巨大。像一张组织切片如果铺了1厘米乘1厘米的芯片，DNB数量能达到数亿甚至更多，GEM文件的文本行数会相当恐怖，用read.table这类函数直接读取极易导致内存爆掉。所以实际项目中我通常会优先转成GEF或中间富集后的稀疏矩阵再做后续处理，尽量避免直接读原始GEM。

如果拿到的是GEM文件且没有SAW环境，可以通过一些简单的命令行工具先做排序和去重，再用Python的pandas分块读取，处理完再保存为h5ad或稀疏矩阵格式，这样可以绕开内存瓶颈。

3.4 备用路径：用Stereopy手动构建表达矩阵

如果你手头没有完整的SAW环境，但已经拿到了比对后带有CID和UMI标签的BAM文件，也可以尝试用Stereopy做一部分重建工作。常见做法是先用samtools提取比对到基因组的read对，把R1的CID、UMI和R2的基因注释整合成类似GEM的表格，然后通过Stereopy加载为StereoExpData对象，再导出为AnnData或Seurat可读的格式。

这种方式的优点是可以更加灵活地控制过滤标准，比如自定义去除低质量reads、自定义UMI去重算法等。缺点是需要你对BAM的格式和标签有比较深入的理解，而且不同版本的SAW在BAM中记录的tag字段可能存在差异。对新手而言，我建议先用SAW跑通，再考虑这种自定义操作。

4. 核心目标：把表达矩阵成功转换为Seurat对象

4.1 转换前必须想清楚的几个数据结构问题

Seurat对象的底层实际上是一个基因乘以细胞的稀疏矩阵，加上一个记录细胞注释信息的meta.data表。所以，无论你的数据来自GEM、GEF还是h5ad，转换前都要先明确两件事：表达矩阵的行是基因还是细胞？坐标信息有没有完整保存下来？

根据Stereo-seq数据自动化处理后的常规格式，如果是bin矩阵，通常行是基因、列是bin编号，每个bin对应一个空间坐标；如果是CellBin矩阵，则行是基因、列是细胞，每个细胞有分割后的质心坐标。无论哪种，构造Seurat对象时矩阵都应该是基因在行、细胞在列，这一点和普通单细胞转录组的要求完全一致。

另一个关键问题是，坐标怎么存储。Seurat的Visium对象通过images槽和spatial坐标来关联组织图像，但Stereo-seq如果没有完整的组织图像和对应的缩放信息，不建议强行去构建SpatialImage对象。一个更通用的做法是把坐标存放在meta.data的x、y两列中，后续用DimPlot、FeaturePlot或者自绘散点图时，根据x和y做映射即可。我几乎所有的Stereo-seq项目都是这么处理的，简单、灵活且不会引入额外依赖。

4.2 三种常用导入方式的详细代码实现

第一种方式是通过AnnData中转，这也是目前最顺滑的路径。先在Python里用Stereopy读出GEF文件，再转成h5ad，最后在R里读取并构建Seurat对象。

import stereopy as st # 读取GEF文件，bin_type可以是bin20、bin50、cellbin data = st.io.read_gef("sample.bin50.gef", bin_type="bin50") # 转为Scanpy的AnnData对象 adata = data.to_scanpy() # 保存为h5ad，方便R端读取 adata.write_h5ad("sample_bin50.h5ad")

library(anndata) library(Seurat) ad <- read_h5ad("sample_bin50.h5ad") # 确认矩阵维度：行是基因，列是细胞/bin counts <- t(as.matrix(ad$X)) # 构建Seurat对象，assay名称设为Spatial obj <- CreateSeuratObject(counts = counts, assay = "Spatial", min.cells = 3, min.features = 10) # 把坐标信息加入meta.data obj$coord_x <- ad$obs$X obj$coord_y <- ad$obs$Y # 保存一下对象 saveRDS(obj, "sample_bin50_seurat.rds")

第二种方式是直接读取CSV或TSV格式的表达矩阵。如果SAW流程产出的矩阵已经整理成了行列清晰的文本文件，可以直接在R里读取。

library(Matrix) library(Seurat) # 文本矩阵通常是基因行、细胞列 counts <- readMM("matrix.mtx") rownames(counts) <- read.csv("genes.tsv", header = FALSE)$V1 colnames(counts) <- read.csv("barcodes.tsv", header = FALSE)$V1 obj <- CreateSeuratObject(counts = counts, assay = "Spatial") # 读取坐标文件 coords <- read.csv("positions.csv", row.names = 1) obj$coord_x <- coords$x obj$coord_y <- coords$y

第三种方式是在R里直接读取CSV，但要求矩阵行列不乱。这种方法适合矩阵不是特别大的场景。如果矩阵很大，我会先把文本转成稀疏矩阵格式，再用ReadMtx函数读取。

4.3 Seurat对象构建后的坐标保存与可视化配置

Seurat对象的meta.data是保存坐标最常见的位置。需要注意坐标列名不要和Seurat默认的列名冲突，否则后续操作可能出错。比如不要命名为x或y这种过于简洁的名字，虽然技术上可以，但容易在合并多个样本时跟其他元数据混淆，我一般用coord_x、coord_y。

可视化的时候，可以直接用ggplot2画空间散点图，也可以用Seurat自己的DimPlot配合自定义坐标轴。这里给出一个简单的示例：

library(ggplot2) plot_df <- data.frame( coord_x = obj$coord_x, coord_y = obj$coord_y, cluster = obj$seurat_clusters ) ggplot(plot_df, aes(x = coord_x, y = coord_y, color = cluster)) + geom_point(size = 0.1) + theme_minimal()

如果后续要用到SpatialDimPlot这类原生空间函数，需要额外构建SpatialImage对象，并配置images槽位。但经过实测，对于Stereo-seq的纯坐标可视化，不建SpatialImage反而更方便，因为原生Spatial函数默认期望的是基于组织图像的坐标系，并不是我们这里简单的X/Y坐标。强行套用反而会得到被压缩、翻转甚至丢失图像信息的图。

4.4 导入后的快速质量检查清单

每次构造完Seurat对象，我第一件事不是急着跑聚类，而是先检查以下几个指标：

• 总UMI数分布是否合理。加上nCount_RNA和nFeature_RNA之后看直方图，能快速判断是否存在低质量cell/bin。
• 基因数过少的bin可能来自组织边缘或者透化不足的区域，基因数过多的bin则需要留意是否是多个细胞被合并到同一个bin中。
• 坐标范围是否与实际组织形状吻合。在空间散点图上，点的分布应该能看到组织轮廓；如果看到点全部挤在一条线上，说明坐标顺序或缩放出了问题。
• 基因名称中是否包含大量NA或未知编号，如果出现这种情况，一般是注释版本不一致导致的。

在Seurat里可以直接这样快速计算：

obj <- NormalizeData(obj) obj <- FindVariableFeatures(obj) obj <- ScaleData(obj) obj <- RunPCA(obj) obj <- FindNeighbors(obj, dims = 1:20) obj <- FindClusters(obj, resolution = 0.5) obj <- RunUMAP(obj, dims = 1:20)

这些步骤和普通单细胞分析完全一致，唯一区别是不要跑RunTSNE，因为空间转录组的降维和聚类基本都用UMAP或直接在空间坐标上观察聚类结果。

5. 常见问题与排查技巧实录

5.1 坐标信息缺失或坐标系混乱

坐标问题在Stereo-seq数据分析中几乎是必踩的坑。比较常见的情况是坐标列全部为0，或者坐标范围异常。这通常和读取GEF时bin_type没有对应上有关。比如数据本身是bin 50聚合的，但你在Stereopy里按bin 20读取，就会导致坐标范围和矩阵维度对不上。

另一种情况是生成Seurat对象时，坐标列被错误地当成了表达矩阵的一部分，导致矩阵维度异常。这类问题的排查思路很简单：先打印表达矩阵的维度，再打印meta.data的行数，两者必须完全一致；如果坐标表里有重复的行名或缺失的行名，也会导致合并错位。我习惯的做法是，在读入坐标文件时把行名设置为和表达矩阵列名一样的标准格式，然后用match函数做一次严格的顺序匹配，而不是直接按顺序合并。

5.2 基因注释版本不一致导致基因名大量缺失

这个问题在转换Seurat对象时表现得非常隐蔽。矩阵里保存的基因名可能是Ensembl ID，但是GTF文件的版本不对，导致输出时所有基因都变成了NA；或者矩阵里是symbol，而Seurat内部某些函数默认需要的是Ensembl ID，进而导致后续FeaturePlot无法显示。

我踩过一次比较深的坑：参考基因组用的GRCh38.p13，GTF却是从另一个网站下载的GRCh38.p13.patch14版本，前几千个基因没问题，但到某条染色体后续区域就开始出现基因名对不上的情况。后来我在转换前加了一步基因名标准化，先把所有基因ID都统一成symbol，再构建Seurat对象，问题就解决了。建议在构建对象之前先用identical函数对比表达矩阵和注释文件中的基因集，把不一致的基因先处理掉。

5.3 内存爆炸与大矩阵压缩技巧

Stereo-seq的空间数据量很大，一张切片在bin 20级别可能就会产生几十万甚至上百万个bin。如果直接全部加载到Seurat中，8G内存的机器基本是撑不住的，更不用说后续还要做标准化和聚类。

我的处理方案通常有两种。第一种是先把低质量的bin过滤掉，比如nFeature_RNA < 100的bin直接剔除，这样能减少很多无效矩阵空间；第二种是利用稀疏矩阵的特性，在读取时就尽量使用Matrix包，而不是将数据全部转成dense矩阵。在R里as.matrix是很容易导致内存翻倍的操作，因此在AnnData转Seurat时不要轻易把稀疏矩阵转成稠密矩阵，直接保留稀疏格式即可。

5.4 CellBin和Bin的选择误区

CellBin和固定大小的Bin是两种完全不同的空间单位。CellBin的核心是借助组织图像和分割算法，把DNB信号聚合成以细胞为边界的单元；而固定Bin只是按照空间网格机械聚合。所以CellBin矩阵的可解释性更强，基因表达也更接近真实单细胞的情况。

但CellBin不是万能的。如果组织切片质量差、细胞密集堆积、或者没有高质量的荧光染色图像，分割结果就会很糟糕。常见标志是：分割出的细胞总数量明显偏离预期、大量细胞面积异常小或异常大、或者空间图上分隔结果和形态学完全对不上。遇到这种情况，我建议退回到bin 50或者bin 100的固定网格矩阵，虽然分辨率低一些，但至少数据是稳定可靠的。

除了质量因素，CellBin的分析还需要考虑它和10x Visium的兼容性。如果你打算用基于离散spot开发的算法来跑CellBin矩阵，比如某些识别空间域的算法，需要先确认算法是否支持非规则空间单元，还是必须依赖规则网格。这个选型决定了后期能不能直接复用很多现成工具。

5.5 空间数据特有的质量控制和空点问题

空间转录组和单细胞测序还有一个明显的不同：空间矩阵里会包含大量位于组织之外的背景bin。这些背景点表达的基因数量很少、总UMI很低，但它们一样会进入下游分析，如果不及时过滤，会在聚类时形成一大群“噪声细胞”，干扰后面的差异分析和空间域识别。

我在做QC时，除了看常规的nFeature_RNA和nCount_RNA，还会特别关注空间位置位于组织轮廓外的低表达点。可以在SAW阶段用组织掩膜文件提前把背景点去掉，也可以在Seurat阶段通过设置subset条件过滤。一个实用技巧是：先画出coord_x和coord_y的密度图，通常能很清楚地看到组织的形态，如果某些点明显偏离主轮廓区域，就可以考虑删除。

6. 最后再分享几个实操中的小经验

先说说我自己的项目流程习惯。拿到GEF文件后，我不会直接一口气把所有bin都转出来，而是会先在bin 100的尺度上快速做一次全流程试探，确认整个数据链路没有问题，再用bin 50甚至CellBin级别去跑正式分析。这样能节省很多时间，尤其是在芯片面积大的时候，bin 100跑的很快，几分钟就能发现大致问题。

其次，关于把坐标保存在meta.data中的注意事项。虽然坐标不是表达矩阵的必需部分，但对空间转录组来说它是灵魂，建议从一开始就把它跟对象深度绑定。每次保存RDS之后，我还会单独导出一份坐标CSV备份，避免后续R版本升级导致对象读取异常时，还能从备份恢复。

还有一个小建议，就是给Seurat对象建立一套规范的命名体系。因为Stereo-seq数据量大、样本多，如果对象的assay名称、样本ID、坐标列名不统一，后面合并多个样本时会非常痛苦。我在构建对象时都会额外添加样本来源信息到meta.data中，这样即使多个样本合并成一个Seurat对象，也能随时通过meta.data进行拆分和对比。

关于Seurat处理Stereo-seq数据，我个人的体会是：不要过于纠结是否需要构建一个跟10x Visium完全一致的Spatial对象，尤其是图像信息这块。Stereo-seq的坐标体系与Visium的组织图像坐标差异很大，强行适配反而会引入新的问题。灵活地把坐标放到meta.data里，才是更通用、更省事的做法。只要表达矩阵和坐标都在，绝大多数下游分析需求都可以顺畅落地。

如果这篇文章帮你少踩了一个坑，或者让你少加了一晚上的班，那说明它写得很值。后面如果遇到具体报错信息，也欢迎拿你实际的数据和日志来沟通，结合具体版本环境来排查会更高效。