📚 课程目录第2章数据预处理流程

数据预处理流程：从原始数据到干净矩阵

好，咱们进入第二章。数据预处理，说白了就是给原始测序数据“洗澡”。空间转录组的下机数据，说白了是一堆FASTQ文件和图像文件，乱得很。你想想看，如果这步没做好，后面的特征学习、聚类分析，全是垃圾。

我个人习惯把预处理分成四个阶段：原始数据解析、Barcode/UMI校正、质量控制和过滤。下面咱们一个一个来。

原始数据解析：H5与loupe文件

原始数据一般有两类来源。10x Visium平台的官方输出有两种格式：H5文件和loupe文件。我刚开始做空间转录组时，第一反应是直接读loupe文件，后来发现H5文件更灵活。

它们有什么区别？看表：

格式	内容	适用场景
H5 (h5ad)	稀疏表达矩阵 + 空间坐标 + 元数据	编程分析（Python/R）
Loupe (cloupe)	10x官方可视化格式	Loupe Browser交互浏览

我个人习惯在解析完H5文件后，先校验一下数据完整性。为什么？我在一个项目里遇到过，H5文件里空间坐标的索引顺序和表达矩阵不一致，导致映射错误。你想想看，如果坐标对不上，后面所有的空间聚类都是乱的。

import scanpy as sc
adata = sc.read_visium('path/to/spatial/', count_file='filtered_feature_bc_matrix.h5')
# 校验空间坐标与barcode是否对齐
print(adata.obs_names[:5])
print(adata.obsm['spatial'][:5])

小技巧： 解析完成后，建议用 scanpy.pp.calculate_qc_metrics 快速看一眼总UMI数和基因数。如果总UMI数异常低（比如不到100），可能是解析时过滤条件太严，回去检查参数。

Barcode与UMI校正：别让噪音污染数据

为什么要校正Barcode？空间转录组的每个spot都有一个独特的barcode序列。但测序过程中会有错误，比如碱基误读。这时我们需要通过白名单（whitelist）来匹配最接近的barcode。

你可能会问，UMI到底是什么？简单说，UMI是一段随机序列，用来标记每一个原始的mRNA分子。它的作用就是区分PCR重复和生物学重复。

  核心观点： Barcode校正解决“这是哪个spot”的问题。UMI校正解决“这到底是一个真实分子还是PCR伪造的”问题。两者缺一不可。

我在一个CD8+ T细胞空间项目中，发现某些spot的UMI计数异常高。后来排查发现，是因为UMI校正时没有处理“编辑距离为1”的序列，把同一个分子的不同测序版本当成了多个分子。嗯，这里要注意，常用的校正工具会默认只保留与白名单完全匹配的barcode，我建议用umi_tools的whitelist参数，跑一遍“编辑距离1”的模糊匹配。

# 命令行用umi_tools校正barcode
umi_tools whitelist --stdin=reads.fastq \
                    --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN \
                    --plot-prefix=read_whitelist \
                    --log=whitelist.log

# 然后提取校正后的barcode与UMI
umi_tools extract --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN \
                  --stdin=reads.fastq \
                  --stdout=extracted.fastq \
                  --read2-in=read2.fastq \
                  --read2-out=read2_extracted.fastq

避坑指南： 我曾经试过跳过Barcode校正，直接去聚合UMI计数。结果一个spot里多出好几倍的表达量。说白了就是PCR重复没有被正确合并。一定要把Barcode和UMI当两个独立的维度去处理。

质量控制（QC）指标：四个金标准

QC指标说白了就是检查这些点好不好。我习惯用四个核心指标：

总UMI数：每个spot捕获到的mRNA分子总数。太少说明spot可能坏死或捕获失败。
检测到的基因数：每个spot有非零表达的基因数量。太少说明细胞密度过低。
线粒体基因比例：线粒体来源的UMI占总UMI的比例。过高说明细胞破损。
核糖体基因比例：ribosomal RNA的占比，过高可能影响下游分析。

你可能会问，这些指标的阈值怎么定？说实话，没有固定标准。我一般先画一个violin plot和scatter plot，观察分布。

import scanpy as sc
import matplotlib.pyplot as plt

# 计算QC指标
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None,
                            log1p=False, inplace=True)

# 看分布
sc.pl.violin(adata, ['n_counts_by_counts', 'n_genes_by_counts', 'pct_counts_mt'],
             multi_panel=True)
plt.show()

指标	常见阈值	说明
总UMI数	500 - 50000	低于500的spot基本是空点
检测基因数	200 - 5000	低于200可能是捕获失败
线粒体比例	10% - 20%	超过20%的spot大概率已破损
核糖体比例	无严格阈值	极端值要留意，但不算硬性剔除条件

  个人经验： 我在一个肿瘤样本项目里，刚开始用默认阈值（UMI>200）过滤，结果把重要的肿瘤边缘spots全删了。后来发现那些spots的UMI本来就不高，但基因多样性很高。所以，我建议你先把数据可视化再看阈值，不要无脑套公式。

过滤低质量spots及基因：果断但不要过度

过滤，说白了就是砍掉坏点。但砍多了心疼，砍少了噪音大。我一般分两步：

过滤低质量spot：基于上面的QC指标，用布尔索引删除不合格的观察点。
过滤低表达基因：那些在很少spot里表达的基因，留着只会增加多重检验负担。

# 过滤spot
adata = adata[adata.obs['n_counts_by_counts'] > 500, :]
adata = adata[adata.obs['n_genes_by_counts'] > 200, :]
adata = adata[adata.obs['pct_counts_mt'] < 20, :]

# 过滤基因：至少在3个spot里表达
sc.pp.filter_genes(adata, min_cells=3)

注意： 过滤顺序很重要。必须先过滤spot，再过滤基因。因为我过滤掉了一些低质量spot后，有些基因原本在少数几个坏点里表达，被误判为“在3个spot里表达”。顺序不对，结果会偏差。

我习惯在过滤完后再跑一次QC指标，对比前后差异。如果过滤后剩余的spot占比低于60%，那可能阈值设得太严格了，回头调整一下。

小技巧： 如果想保留更多数据，可以用“soft filtering”方法。不直接删除，而是给每个spot赋予一个权重，低质量点的权重低一些。但说实话，对于大多数项目，硬过滤就够用了。

对了，还有一个坑。我在一个脑组织数据中，发现过滤后仍有一些spot的UMI分布峰度异常。后来发现是因为组织切片边缘的spots天生就少细胞。嗯，这种边缘效应不能当作低质量处理，需要保留。

过滤完，我们会得到一个干净的表达矩阵。这个矩阵就是后续所有分析的基础。你想想看，如果这步做对了，后面就顺了。