📚 课程目录第3章空间转录组数据获取与预处理

3. 空间转录组数据获取与预处理

空间转录组的数据处理，说白了就是把你从测序仪拿到的原始文件，一步步变成可以分析的基因表达矩阵。这一步做不好，后面所有的空间分析都是空中楼阁。我自己早期踩过不少坑，今天把这些经验掰开揉碎讲给你。

3.1 原始数据解析：从BAM/FASTQ开始

拿到空间转录组数据，第一步就是处理原始文件。以10x Visium为例，你通常会得到BAM或者FASTQ文件。我个人习惯先看一眼BAM里的tag信息——嗯，这步很多人会跳过，但我觉得挺关键。

BAM文件里，每个比对片段都带了关键的tag，比如CB（barcode序列）、UB（UMI序列）、GN（基因名）。这些tag是构建表达矩阵的基础。我记得早期做10x数据时，bam里tag含义查了好久文档，后来才发现官方工具已经帮你解析好了大部分。

解析流程速览：

FASTQ → 比对（STAR/cellranger）→ 排序标记 → BAM
BAM → tag提取 → barcode校正 → UMI去重 → 表达计数

如果你用cellranger pipeline，一条命令就能完成大部分工作：

cellranger count --id=sample1 \
                 --transcriptome=/ref/refdata-gex-GRCh38-2020-A \
                 --fastqs=/data/fastq \
                 --sample=mysample \
                 --localcores=16

但如果你需要自己搭建流程，或者处理非10x平台的数据，那就得手动走STAR比对这一步。

STAR --genomeDir /ref/star_genome \
     --readFilesIn R1.fastq R2.fastq \
     --outFileNamePrefix sample1_ \
     --outSAMtype BAM SortedByCoordinate

我的经验：比对参数不要无脑用默认。空间转录组的读长结构和单细胞略有不同，R1通常包含barcode和UMI，R2才是转录本序列。如果搞反了，比对率会直接崩掉。我建议先跑100万条reads测一下比对率，正常应该在85%以上。

3.2 基因表达矩阵构建

从BAM到表达矩阵，核心就两步：barcode校正和UMI去重。说白了，就是告诉你每个spot里每个基因到底检测到了多少唯一的mRNA分子。

为什么要校正barcode？因为测序会有错误。比如实际序列是"ACGTACGT"，但测成了"ACGTACCT"。cellranger会用whitelist来做校正，把错误的barcode匹配到最近的正确barcode上。

UMI去重就更关键了。同一个mRNA分子可能被扩增多次，产生多个测序片段。如果不做去重，表达量会被严重高估。我记得有一次分析一个肺癌样本，没做UMI去重前"看起来"基因表达特别高，去重后直接砍掉了60%的计数——这才是真实情况。

最终构建的矩阵长这样：

基因 \ spot	spot_AAACCCAAGACTGGTA-1	spot_AAACCCAAGAGCTCAA-1	spot_AAACCCAAGATAGGCT-1
TP53	12	0	5
EGFR	0	8	1
CD3E	45	3	0

这是一个稀疏矩阵，大部分值是0。空间转录组的数据稀疏度通常在90%以上，所以存储时强烈建议用稀疏格式（比如Market Exchange Format或HDF5）。

3.3 质控：UMI计数与低质量spots过滤

矩阵到手，别急着分析。先做质控。我遇到过一批胃癌样本，UMI count分布特别奇怪——有的spot有5万UMI，有的只有几百。后来发现是组织透化时间不够，导致边缘区域RNA捕获效率极低。

常用的质控指标有三个：

UMI总数：反映每个spot的RNA捕获效率
检测到的基因数：反映复杂度和组织覆盖
线粒体基因占比：反映细胞状态（高线粒体比例通常意味着细胞受损或死亡）

这里给出一组我常用的阈值（不同组织类型需要调整）：

质控指标	下限	上限	说明
UMI总数	500	30000	低于500视为空spot，高于30000可能是doublet
检测基因数	200	6000	低于200组织覆盖不足，高于6000可能是技术异常
线粒体基因占比	—	20%	超过20%提示细胞损伤，建议剔除

避坑提醒：我曾经遇到过一批数据，线粒体比例整体偏高（30%左右），一开始想全部过滤掉。后来看了组织切片才知道，那片区域是肿瘤缺氧区，线粒体比例高是生物学真实信号，不是技术问题。所以质控一定要结合组织形态学判断，不要硬套阈值。

下面是用Scanpy做质控的代码示例：

import scanpy as sc

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

sc.pp.calculate_qc_metrics(adata, percent_top=None,
                           log1p=False, inplace=True)

# 过滤低质量spots
adata = adata[adata.obs["n_genes_by_counts"] > 200, :]
adata = adata[adata.obs["n_genes_by_counts"] < 6000, :]
adata = adata[adata.obs["total_counts"] > 500, :]
adata = adata[adata.obs["total_counts"] < 30000, :]
adata = adata[adata.obs["pct_counts_mt"] < 20, :]

print(f"质控后保留 {adata.n_obs} 个spots")

3.4 归一化与批次校正

质控完，接下来做归一化。为什么要归一化？因为不同spot的测序深度不同，直接比较原始计数会引入偏差。说白了，就是一个spot测了5万条reads，另一个spot只测了1万条reads，不代表前者基因表达真的高5倍。

常用的归一化方法有两种：

LogNormalize：每个spot的计数除以该spot总计数，乘以缩放因子（默认10000），再log1p变换。简单粗暴，适合快速探索。
SCTransform：基于正则化负二项回归，考虑了基因表达和测序深度的关系。我个人更推荐这个，尤其当数据质量差异较大时。

# LogNormalize (Seurat)
adata = sc.pp.normalize_total(adata, target_sum=1e4)
adata = sc.pp.log1p(adata)

# SCTransform (Seurat)
# 推荐在R中使用，Python可通过sctransform包调用

接下来是批次校正。空间转录组项目中，你可能会有多个样本，来自不同病人、不同切片、不同批次。如果不校正，聚类结果会首先按批次分开，而不是按生物学差异分开。

我习惯用Harmony做空间转录组的批次校正。参数调节比较直观，而且能保留空间结构信息。ComBat也常用，但更适用于bulk数据，对空间数据来说有时会过度校正。

批次校正的核心原则：

先可视化（UMAP/tSNE）看是否有明显批次分离
选择校正方法时，优先考虑能保留数据结构的算法
校正后一定要验证：批次混合了，但生物学差异还在

举个例子，我之前处理过一批结直肠癌样本，来自3个不同批次。用Harmony校正前，UMAP上3个批次完全分开；校正后，样本按细胞类型聚在一起，但肿瘤区域和正常区域仍然分离得很好——这就是我们要的效果。

import harmonypy as hp

# 假设已经做好PCA
ho = hp.run_harmony(
    adata.obsm["X_pca"],
    adata.obs,
    "batch_id",
    max_iter_harmony=10
)
adata.obsm["X_pca_harmony"] = ho.Z_corr.T

个人习惯：归一化和批次校正做完后，我会马上跑一个 quick cluster + marker gene 验证。如果已知的marker基因（比如EPCAM上皮、PTPRC免疫）能够在cluster中正确富集，说明预处理步骤没问题。这一步只要10分钟，但能帮你省掉后面几天走弯路的时间。

到这里，空间转录组的数据预处理就完成了。你手里现在是一个干净、归一化、批次校正后的表达矩阵，可以直接做空间聚类、差异分析、空间可变基因检测等下游任务了。