空间转录组Stereo-seq分析¶

Visium:   55μm分辨率 → 每个spot含1-10个细胞 → 约5,000 spots
Stereo-seq: 500nm分辨率 → 真正亚细胞级 → 可达数十亿DNB

# SAW V8.0+已经是自包含的tar.gz，解压即可使用
# 从STOmics官方GitHub下载：https://github.com/STOmics/SAW

# 下载SAW（以V8.1.3为例）
wget https://github.com/STOmics/SAW/releases/download/v8.1.3/SAW_v8.1.3.tar.gz  # 下载SAW

# 解压安装
tar -xzf SAW_v8.1.3.tar.gz              # 解压，无需额外配置
cd SAW_v8.1.3                            # 进入目录
export PATH=$PWD:$PATH                   # 添加到环境变量

# ========== SAW分析五大步骤 ==========

# 第一步：Reads分析（比对和基因表达矩阵生成）
# mapping：将FASTQ中的reads与芯片空间信息对应
saw mapping \                            # reads比对
  --fq1 sample_R1.fastq.gz \            # Read1文件（含空间barcode）
  --fq2 sample_R2.fastq.gz \            # Read2文件（含cDNA序列）
  --genomeDir /ref/STAR_index/ \         # STAR基因组索引目录
  --mask /data/chip_mask.h5 \            # 芯片mask文件（含DNB位置信息）
  --outSAMtype BAM \                     # 输出BAM格式
  --threads 16                           # 线程数

# 第二步：组织区域识别与基因表达提取
# tissueCut：自动识别组织覆盖区域
saw tissueCut \                          # 组织分割
  --input mapping_output.gef \           # 上一步生成的GEF文件
  --output tissue.gef \                  # 输出只含组织区域的GEF
  --binSize 200                          # bin200粒度分析

# 第三步：空间聚类分析（bin200级别）
saw spatialCluster \                     # 空间聚类
  --input tissue.gef \                   # 输入组织区域GEF
  --binSize 200 \                        # bin200粒度
  --resolution 0.8 \                     # 聚类分辨率
  --output cluster_result                # 输出目录

# 第四步：CellBin单细胞分割（需要染色图像）
# cellCut：基于显微图像进行细胞核分割
saw cellCut \                            # 细胞分割
  --input tissue.gef \                   # 输入GEF文件
  --image tissue_staining.tif \          # 组织染色图像（ssDNA/DAPI/HE）
  --output cellbin.gef                   # 输出单细胞级GEF

# 第五步：单细胞级聚类
saw cellCluster \                        # 细胞级聚类
  --input cellbin.gef \                  # 单细胞GEF
  --resolution 1.0 \                     # 聚类分辨率
  --output cell_cluster_result           # 输出目录

bin1    → 500nm × 500nm    → 最精细，数据量巨大（适合局部分析）
bin50   → 25μm × 25μm     → 接近亚细胞（适合精细空间分析）
bin100  → 50μm × 50μm     → 类似Visium的spot大小
bin200  → 100μm × 100μm   → 常用起始分析粒度（速度快）
CellBin → 单细胞级别       → 基于图像分割，最有生物学意义

# 安装Stereopy
# pip install stereopy                   # pip安装

import stereo as st                      # 导入Stereopy
import numpy as np                       # 数组操作

# 读取GEF文件
data = st.io.read_gef(
    file_path="tissue.gef",             # GEF文件路径
    bin_size=100                          # bin100粒度
)

# 也可以读取为AnnData格式（兼容Scanpy）
# data = st.io.read_gef_to_anndata("tissue.gef", bin_size=100)

# 质控过滤
data.tl.cal_qc()                         # 计算质控指标
data.plt.violin(gene_names=["total_counts", "n_genes_by_counts"])  # 小提琴图

# 归一化
data.tl.normalize_total(target_sum=1e4)  # 总量归一化到10000
data.tl.log1p()                          # log(x+1)转换

# 高变基因选择
data.tl.highly_variable_genes(
    n_top_genes=3000,                    # 选择top 3000高变基因
    min_mean=0.0125,                     # 平均表达量下限
    max_mean=3                           # 平均表达量上限
)

# PCA降维
data.tl.pca(
    use_highly_genes=True,               # 只用高变基因
    n_pcs=50                             # 计算50个主成分
)

# 邻近图和聚类
data.tl.neighbors(n_pcs=30)              # 构建30维的近邻图
data.tl.leiden(resolution=1.0)           # Leiden聚类

# UMAP降维可视化
data.tl.umap()                           # 运行UMAP
data.plt.umap(color="leiden")            # 按聚类着色的UMAP图

# 空间可视化（核心！）
data.plt.spatial_scatter(
    color="leiden"                        # 按聚类在空间上着色
)

# 差异表达分析
data.tl.rank_genes_groups(
    groupby="leiden",                    # 按聚类分组
    method="wilcoxon"                    # Wilcoxon检验
)

import stereo as st                      # 导入Stereopy
import scanpy as sc                      # 导入Scanpy

# GEF → AnnData（Scanpy格式）
adata = st.io.read_gef_to_anndata(
    file_path="tissue.gef",             # GEF文件
    bin_size=100                          # bin大小
)

# 保存为h5ad（Scanpy标准格式）
adata.write("stereo_seq_data.h5ad")      # 保存

# 之后可以用Scanpy的全部功能分析
sc.pp.normalize_total(adata, target_sum=1e4)  # Scanpy归一化
sc.pp.log1p(adata)                       # log转换
sc.pp.highly_variable_genes(adata)       # 高变基因
sc.tl.pca(adata)                         # PCA
sc.pp.neighbors(adata)                   # 构建图
sc.tl.leiden(adata)                      # 聚类
sc.pl.spatial(adata, color="leiden")     # 空间可视化

# ========== Stereo-seq分析速查 ==========

# SAW核心命令
saw mapping --fq1 R1.fq --fq2 R2.fq --genomeDir REF --mask MASK   # 比对
saw tissueCut --input map.gef --output tissue.gef --binSize 200     # 组织分割
saw spatialCluster --input tissue.gef --binSize 200                 # 聚类
saw cellCut --input tissue.gef --image stain.tif                    # 细胞分割

# Stereopy (Python) 核心流程
data = st.io.read_gef("tissue.gef", bin_size=100)  # 读取
data.tl.normalize_total(target_sum=1e4)             # 归一化
data.tl.log1p()                                     # log转换
data.tl.pca(use_highly_genes=True)                  # PCA
data.tl.neighbors(n_pcs=30)                         # 建图
data.tl.leiden()                                    # 聚类
data.plt.spatial_scatter(color="leiden")             # 空间可视化

# 格式转换
adata = st.io.read_gef_to_anndata("file.gef", bin_size=100)  # → AnnData
adata.write("output.h5ad")                                     # → h5ad