SAW处理Stereo-seq的流程整理

最新推荐文章于 2025-12-12 16:47:26 发布

原创最新推荐文章于 2025-12-12 16:47:26 发布 · 677 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#服务器 #运维

需求1：需要生成cellbin.gef和tissue.gef

本流程基于华大官方 SAW (Stereo-seq Analysis Workflow) 实现，涵盖测序数据处理→坐标配准→细胞分割→多尺度 GEF 生成全环节，明确标注每一步的输入文件、核心操作、输出文件，适配植物学（如水稻胚乳）空间转录组数据分析场景。

前提：准备必备文件清单

文件类型	具体文件	作用
原始测序数据	`sample_R1.fq.gz`、`sample_R2.fq.gz`	双端测序原始数据，R1 含空间 barcode+UMI，R2 含基因序列
参考基因组数据	`genome.fa`（基因组序列）、`genes.gtf`/`genes.gff3`（基因注释）	用于序列比对和基因表达定量
图像数据	`ssDNA_tifs.tar.gz`（分割的 ssDNA 染色 TIFF 小图）	用于空间坐标配准和细胞边界识别
芯片参考文件	`mask.h5`（芯片点阵 mask 文件）	提供芯片位点分布信息，关联测序 barcode 与物理坐标

前置准备：定义全局变量（简化后续命令）

# 替换为你的样本名称（如 Rice_Endosperm_01）
SAMPLE="your_sample_name"
# 替换为你的数据根目录（所有输入文件存放路径）
ROOT_DIR="/path/to/your/data"
# 替换为参考基因组目录（genome.fa + genes.gtf 所在路径）
REF_DIR="/path/to/reference/genome"
# 替换为芯片 mask.h5 文件路径（STOmics 官方提供对应芯片型号的 mask）
MASK_FILE="/path/to/chip_mask/mask.h5"
# 替换为 ssDNA 图像压缩包路径
IMG_TAR="${ROOT_DIR}/ssDNA_tifs.tar.gz"

# 输出目录（自动创建，避免文件混乱）
OUT_DIR="${ROOT_DIR}/saw_output"
mkdir -p ${OUT_DIR}/mapping ${OUT_DIR}/image ${OUT_DIR}/register ${OUT_DIR}/tissue ${OUT_DIR}/cellbin

步骤 1：序列比对 + 基因表达定量（mapping + count）

# 1.1 序列比对（mapping）
saw mapping \
  -i ${ROOT_DIR}/${SAMPLE}_R1.fq.gz,${ROOT_DIR}/${SAMPLE}_R2.fq.gz \  # 双端测序文件（R1在前，R2在后，逗号分隔）
  -g ${REF_DIR}/genome.fa \                                          # 参考基因组序列
  -a ${REF_DIR}/genes.gtf \                                          # 基因注释文件
  -m ${MASK_FILE} \                                                  # 芯片 mask 文件
  -o ${OUT_DIR}/mapping \                                            # 比对结果输出目录
  -p 16 \                                                             # 线程数（根据服务器配置调整，建议8-32）
  --species rice \                                                    # 物种（如 human/mouse/rice，按需替换）
  --barcode-length 26 \                                              # Stereo-seq 标准 barcode 长度
  --umi-length 10                                                    # Stereo-seq 标准 UMI 长度

# 1.2 基因表达定量（count）
saw count \
  -i ${OUT_DIR}/mapping/align.bam \                                  # mapping 输出的 bam 文件
  -g ${REF_DIR}/genes.gtf \                                          # 基因注释文件（与 mapping 一致）
  -m ${MASK_FILE} \                                                  # 芯片 mask 文件
  -o ${OUT_DIR}/mapping/${SAMPLE}.raw.gef \                          # 输出未配准的 raw.gef
  -p 16 \
  --min-umi 1 \                                                      # 过滤低表达位点（至少1个UMI）
  --min-gene 1                                                       # 过滤无基因表达的位点（至少1个基因）

输出文件：${OUT_DIR}/mapping/${SAMPLE}.raw.gef、${OUT_DIR}/mapping/mapping_stat.txt

步骤 2：图像拼接与预处理（imageTools）

saw imageTools \
  -i ${IMG_TAR} \                                                    # 输入 ssDNA 图像压缩包
  -o ${OUT_DIR}/image \                                              # 图像预处理输出目录
  --stitch \                                                         # 开启图像拼接
  --preprocess \                                                     # 开启标准化预处理（降噪+对比度增强）
  --mask \                                                           # 生成组织区域二值化掩膜
  --pixel-size 0.25 \                                                # 芯片像素尺寸（Stereo-seq 标准 0.25μm/像素）
  --filter-size 3                                                    # 高斯滤波核大小（默认3，无需修改）

输出文件：${OUT_DIR}/image/${SAMPLE}_ssDNA_stitch.tif、${OUT_DIR}/image/${SAMPLE}_tissue_mask.tif、${OUT_DIR}/image/image_preprocess.json

步骤 3：空间坐标配准（register）

saw register \
  -i ${OUT_DIR}/mapping/${SAMPLE}.raw.gef \                          # 输入 raw.gef
  -img ${OUT_DIR}/image/${SAMPLE}_ssDNA_stitch.tif \                 # 拼接后的 ssDNA 全景图
  -mask ${OUT_DIR}/image/${SAMPLE}_tissue_mask.tif \                 # 组织掩膜图
  -m ${MASK_FILE} \                                                  # 芯片 mask 文件
  -o ${OUT_DIR}/register/${SAMPLE}.register.gef \                    # 输出配准后的 register.gef
  --method morphology \                                              # 配准方法（形态学匹配，默认最优）
  --align-threshold 0.8 \                                            # 特征点匹配阈值（0.8为通用值）
  --output-image ${OUT_DIR}/register/${SAMPLE}_registered.tif        # 输出配准后的图像

输出文件：${OUT_DIR}/register/${SAMPLE}.register.gef、${OUT_DIR}/register/registration_param.json、${OUT_DIR}/register/${SAMPLE}_registered.tif

步骤 4：组织筛选 + 多 bin 层级聚合（filter + bin）

# 4.1 组织区域筛选
saw filter \
  -i ${OUT_DIR}/register/${SAMPLE}.register.gef \                    # 输入 register.gef
  -mask ${OUT_DIR}/image/${SAMPLE}_tissue_mask.tif \                 # 组织掩膜图
  -o ${OUT_DIR}/tissue/${SAMPLE}.filtered.gef \                      # 输出筛选后的 GEF
  --stat ${OUT_DIR}/tissue/tissue_filter_stat.txt                    # 筛选统计报告

# 4.2 多 bin 层级聚合
saw bin \
  -i ${OUT_DIR}/tissue/${SAMPLE}.filtered.gef \                      # 输入筛选后的 GEF
  -o ${OUT_DIR}/tissue/${SAMPLE}.tissue.gef \                        # 输出最终 tissue.gef
  --bin-sizes 1,5,10,50,100,200 \                                    # 需聚合的 bin 尺度（按需增减）
  --stat ${OUT_DIR}/tissue/bin_agg_stat.txt                          # 聚合统计报告

输出文件：${OUT_DIR}/tissue/${SAMPLE}.tissue.gef、${OUT_DIR}/tissue/tissue_filter_stat.txt

步骤 5：细胞分割 + cellbin 生成（cellCut + cellCorrect）

# 5.1 细胞分割（cellCut）
saw cellCut \
  -i ${OUT_DIR}/register/${SAMPLE}.register.gef \                    # 输入 register.gef
  -img ${OUT_DIR}/image/${SAMPLE}_ssDNA_stitch.tif \                 # 预处理后的 ssDNA 全景图
  -o ${OUT_DIR}/cellbin/${SAMPLE}.cellbin.raw.gef \                  # 输出初始 cellbin GEF
  --algorithm cellpose \                                             # 分割算法（Cellpose 适配 ssDNA/DAPI 最优）
  --cellpose-model nucleus \                                          # Cellpose 模型（细胞核分割）
  --mask ${OUT_DIR}/image/${SAMPLE}_tissue_mask.tif \                 # 组织掩膜（仅分割组织内细胞）
  -p 16

# 5.2 细胞分割结果校准（cellCorrect）
saw cellCorrect \
  -i ${OUT_DIR}/cellbin/${SAMPLE}.cellbin.raw.gef \                  # 输入初始 cellbin GEF
  -img ${OUT_DIR}/image/${SAMPLE}_ssDNA_stitch.tif \                 # ssDNA 全景图
  -o ${OUT_DIR}/cellbin/${SAMPLE}.cellbin.gef \                      # 输出最终校准后的 cellbin.gef
  --stat ${OUT_DIR}/cellbin/cell_correct_stat.txt \                  # 校准统计报告
  --correct-threshold 0.1                                            # 校准阈值（默认0.1，无需修改）

输出文件：${OUT_DIR}/cellbin/${SAMPLE}.cellbin.gef、${OUT_DIR}/cellbin/cell_segmentation.json、${OUT_DIR}/cellbin/cell_correct_stat.txt

关键补充说明

环境依赖：需提前安装 SAW 并配置环境（推荐通过 conda 安装：conda create -n saw -c stomics saw）；
芯片 mask 文件：需从 STOmics 官网下载对应芯片型号的 mask.h5（如 S2/S4 芯片）；
线程数调整：-p 16 为线程数，若服务器资源有限，可改为 -p 8；
物种参数：--species 支持 human/mouse/rice/arabidopsis 等，非模式物种可省略该参数；
报错排查：若某步骤报错，优先检查输入文件路径是否正确、文件格式是否为标准格式（如 GTF 需为版本 2/3，TIFF 需为 8/16 位）

简化版本需求2：仅需要生成registered.gef

前置变量定义（同之前，仅保留核心路径）

SAMPLE="your_sample_name"
ROOT_DIR="/path/to/your/data"
REF_DIR="/path/to/reference/genome"
MASK_FILE="/path/to/chip_mask/mask.h5"
IMG_TAR="${ROOT_DIR}/ssDNA_tifs.tar.gz"

# 输出目录（仅保留核心环节）
OUT_DIR="${ROOT_DIR}/saw_output"
mkdir -p ${OUT_DIR}/mapping ${OUT_DIR}/image ${OUT_DIR}/register

步骤 1：mapping + count（生成 raw.gef）

# 1.1 序列比对
saw mapping \
  -i ${ROOT_DIR}/${SAMPLE}_R1.fq.gz,${ROOT_DIR}/${SAMPLE}_R2.fq.gz \
  -g ${REF_DIR}/genome.fa \
  -a ${REF_DIR}/genes.gtf \
  -m ${MASK_FILE} \
  -o ${OUT_DIR}/mapping \
  -p 16 \
  --species rice \  # 按需替换为 human/mouse 等
  --barcode-length 26 \
  --umi-length 10

# 1.2 基因表达定量
saw count \
  -i ${OUT_DIR}/mapping/align.bam \
  -g ${REF_DIR}/genes.gtf \
  -m ${MASK_FILE} \
  -o ${OUT_DIR}/mapping/${SAMPLE}.raw.gef \
  -p 16 \
  --min-umi 1 \
  --min-gene 1

步骤 2：imageTools（仅做拼接 + 预处理，无需生成组织掩膜）

saw imageTools \
  -i ${IMG_TAR} \
  -o ${OUT_DIR}/image \
  --stitch \          # 仅开启拼接（核心）
  --preprocess \      # 仅开启标准化预处理（核心）
  --pixel-size 0.25 \ # Stereo-seq 标准像素尺寸
  --filter-size 3
# 注：省略 --mask 参数，无需生成组织掩膜（仅配准不需要）

步骤 3：register（生成配准后的 register.gef）

saw register \
  -i ${OUT_DIR}/mapping/${SAMPLE}.raw.gef \
  -img ${OUT_DIR}/image/${SAMPLE}_ssDNA_stitch.tif \  # imageTools 输出的拼接全景图
  -m ${MASK_FILE} \                                    # 无需传入组织掩膜（-mask 参数可省略）
  -o ${OUT_DIR}/register/${SAMPLE}.register.gef \
  --method morphology \
  --align-threshold 0.8 \
  --output-image ${OUT_DIR}/register/${SAMPLE}_registered.tif

四、关键精简点说明

省略步骤：完全跳过 filter + bin（生成 tissue.gef）、cellCut + cellCorrect（生成 cellbin.gef），仅保留到 register 步骤；
imageTools 简化：省略 --mask 参数，无需生成组织掩膜（组织掩膜仅用于筛选 tissue.gef 的背景位点，配准本身不需要）；
register 简化：省略 -mask 输入参数，仅依赖拼接后的 ssDNA 全景图完成坐标对齐。

五、最终输出

仅保留核心文件：

${OUT_DIR}/register/${SAMPLE}.register.gef：配准后的目标 GEF 文件（含 bin1 分辨率 + 精准物理坐标）；
${OUT_DIR}/register/${SAMPLE}_registered.tif：配准后的 ssDNA 图像（用于验证配准效果）。

六、验证配准效果（可选但建议）

生成 register.gef 后，可通过 Stereopy 快速验证配准是否成功：

Registered.tif的作用

registered.tif 在后续的 filter/bin/cellCut 等核心分析步骤中，确实不是必需的输入文件，它的核心价值并非 “驱动下游分析”，而是作为「配准有效性的可视化验证工具」和「结果解读的辅助参考」。

一、`registered.tif` 不是下游分析必需文件的原因

Stereo-seq 下游分析的核心是基因表达数据的空间坐标信息，而非图像本身：

下游步骤的核心输入是 register.gefregister.gef 已经包含了校准后的物理坐标（与 registered.tif 完全对齐），后续的 filter（组织筛选）、bin（多尺度聚合）、cellCut（细胞分割）都是直接读取 register.gef 中的坐标和表达数据，无需再调用 registered.tif。
组织筛选依赖的是 mask.tif，而非 registered.tiffilter 步骤区分组织 / 背景的依据是上游 imageTools 生成的 mask.tif（二值化掩膜），而非配准后的图像；cellCut 步骤的细胞分割也是基于原始 stitch.tif（或手动修正的分割图），因为 stitch.tif 保留了更完整的荧光信号细节。

二、`registered.tif` 的真正作用（不可替代的 3 个场景）

虽然它不是分析必需文件，但却是保障分析结果可靠性的关键参考，核心用途集中在验证和解读：

验证配准是否成功（最核心用途）这是 registered.tif 的首要价值。你可以通过两个方式验证：
- 视觉对比：用 ImageJ 同时打开 stitch.tif 和 registered.tif，观察两者的组织形态是否一致 —— 若配准成功，registered.tif 应与 stitch.tif 完全重合（无偏移、旋转）；若配准失败，registered.tif 会出现明显的形态扭曲或位移。
- 数据 - 图像关联验证：用 Stereopy 读取 register.gef，绘制高表达基因的空间热图，将热图叠加到 registered.tif 上 —— 若配准有效，基因表达的热点区域应与 registered.tif 中的细胞核密集区域完全对应。
  
  没有 registered.tif，你无法直观判断 register.gef 的坐标是否准确，后续分析的可靠性就无从保障。
辅助解读空间分析结果当你得到下游分析结果（如 cellbin 的细胞类型注释、区域差异基因）时，registered.tif 是最佳的可视化背景图：
- 例如：将细胞类型的空间分布图叠加在 registered.tif 上，可直接观察 “某种细胞类型是否富集在组织的特定区域”；
- 例如：将差异基因的表达热图与 registered.tif 叠加，可判断基因表达是否与细胞核的分布特征相关。
  
  相比原始 stitch.tif，用 registered.tif 作为背景图能确保图像与数据的坐标 100% 对齐，避免叠加时出现错位。
用于实验记录和论文配图科研论文中展示空间转录组结果时，通常需要提供组织形态图 + 基因表达热图的叠加图，registered.tif 是最适合的组织形态底图 —— 它已经过坐标校准，与数据完全匹配，无需额外调整就能用于作图。

三、总结：`registered.tif` 是 “验证工具”，而非 “分析原料”

角色定位	是否必需	核心用途
下游分析输入文件	❌ 非必需	无（下游步骤仅需 `register.gef`/`mask.tif`）
配准有效性验证工具	✅ 建议保留	直观判断坐标校准是否准确，避免后续分析出错
结果可视化背景图	✅ 推荐使用	辅助解读空间表达模式，用于论文配图

简单来说：没有 registered.tif，你也能跑完整个流程，但你无法确定流程中最关键的「配准步骤」是否真的成功—— 这也是为什么 saw register 步骤默认会输出这个文件。

“merge 真实图像” 的实操逻辑（以 Stereopy 为例）

import stereopy as st
import matplotlib.pyplot as plt

# 1. 读取配准后的GEF（含表达+校准坐标）
data = st.io.read_gef("sample.register.gef")

# 2. 读取registered.tif（含校准后的真实图像）
img = st.io.read_image("sample.registered.tif")

# 3. 融合可视化（merge图像与基因表达数据）
fig, ax = plt.subplots(1,1,figsize=(10,8))
# 先绘制registered.tif作为底图
ax.imshow(img, cmap="gray")
# 再叠加基因表达热图（坐标完全对齐）
st.pl.spatial_scatter(
    data,
    gene="Os01g0100100",  # 目标基因
    ax=ax,
    dot_size=2,
    cmap="Reds",
    alpha=0.7  # 透明化，不遮挡底图
)
plt.savefig("gene_img_merge.png", dpi=300)

核心结论先明确：配准成功后，registered.tif 和 stitch.tif 的视觉形态看起来几乎一致，但坐标体系完全不同 —— 这就是 stitch.tif 不能直接叠加的根本原因。

两者的视觉相似性，很容易让人误以为坐标可以直接对应，但实际上配准的核心就是调整坐标体系，而非改变图像的视觉形态。以下是详细拆解：

一、`stitch.tif` 与 `registered.tif` 的核心差异：坐标体系不同

特征	`stitch.tif`（原始拼接图）	`registered.tif`（配准后图像）
坐标本质	基于图像像素的相对坐标（单位：像素）	基于芯片物理位点的绝对坐标（单位：μm）
坐标原点	图像左上角（默认像素坐标原点）	芯片的物理原点（与 `registered.gef` 位点坐标原点一致）
坐标校准关系	与 `registered.gef` 无关联	与 `registered.gef` 位点坐标一一对应

举个直观例子（假设配准前存在偏移）

stitch.tif 的像素坐标假设 stitch.tif 中一个细胞核的中心在像素坐标 (500, 800)（第 500 列、第 800 行像素），这个坐标只代表它在图像中的相对位置，和芯片上的基因捕获位点没有任何关系。
配准的核心操作：坐标转换saw register 步骤会计算一个坐标转换矩阵（比如纠正 x 方向偏移 50 像素、y 方向偏移 30 像素，同时校准像素与微米的换算比例），然后对 stitch.tif 做两个关键处理：
- 把像素坐标换算为芯片物理坐标（比如 1 像素 = 0.25 μm → 像素 (500,800) 对应物理坐标 (125, 200) μm）；
- 把坐标原点调整为芯片的物理原点（与 registered.gef 中位点的坐标原点对齐）。
registered.tif 的物理坐标经过转换后，这个细胞核的中心坐标变成 (125, 200) μm，和 registered.gef 中对应位点的坐标完全一致—— 此时叠加基因表达数据，才能精准对应。

而 stitch.tif 没有经过这个坐标转换，它的像素坐标和 registered.gef 的物理位点坐标是 “两套独立的坐标系”，直接叠加必然错位。

二、为什么视觉上看起来几乎一致？

你说的 “配准成功后两者完全一致”，指的是图像的像素排列（视觉形态）一致，而非坐标一致。原因如下：

配准的本质是 “坐标校准”，而非 “图像变形”只有当原始图像与基因数据存在旋转、拉伸时，registered.tif 才会在视觉上略有变化；如果配准前仅存在平移偏移（最常见情况），registered.tif 只是坐标原点移动，像素排列完全不变 —— 视觉上和 stitch.tif 几乎一模一样，但坐标体系已经完全不同。
像素尺寸校准不影响视觉形态配准过程会把 stitch.tif 的 “像素单位” 换算为 “微米单位”（比如 0.25 μm / 像素），这个换算只改变坐标的数值和单位，不改变像素的排列方式 —— 视觉上自然看不出差异。

三、直接叠加 `stitch.tif` 的后果：坐标错位

用一个简单的逻辑链就能理解：

registered.gef 中的位点坐标是 物理坐标（μm），比如位点 A 的坐标是 (125, 200) μm；
stitch.tif 的坐标是 像素坐标，比如图像中对应位置的像素坐标是 (500, 800)；
直接叠加时，软件会默认把 registered.gef 的 (125, 200) 当成像素坐标，对应到 stitch.tif 的 (125, 200) 像素位置 —— 这个位置和实际细胞核的位置 (500, 800) 完全不同，最终表现为基因热图和组织形态错位。