📚 课程目录第3章高质量原始数据处理

高质量原始数据处理

数据质量决定了后续分析的成败，这句话我在不同场合说过很多次。你拿到的10X原始数据里，有空的液滴、破损的细胞、还有两个细胞黏在一起的双细胞——这些东西不清理掉，下游聚类和差异分析全都会跑偏。这一章我们就从头捋一遍：怎么把10X数据正确读进来，怎么判断哪些细胞该扔，以及双细胞检测的基本套路。

  核心原则：宁缺毋滥。过滤太松，噪声会淹没真实信号；过滤太严，稀有细胞群可能被误删。关键是根据你的数据分布来调阈值，而不是硬套某个数值。

1. 读取10X数据

读取这一步看着简单，但格式搞错了后面全得重来。10X官方输出主要有三种格式，我列在下面：

格式	文件	读取方式	适用场景
mtx	matrix.mtx / barcodes.tsv / features.tsv	sc.read_10x_mtx()	Cell Ranger < 3.0 或自定义输出
h5	filtered_feature_bc_matrix.h5	sc.read_10x_h5()	Cell Ranger ≥ 3.0 标准输出
h5ad	*.h5ad	sc.read_h5ad()	已预处理的数据交换

我个人习惯用 h5 格式，一个文件搞定，不用管三个文件的路径对应。但如果你拿到的是 mtx 目录，也别慌——两种方式效果完全一样。

# 方式1：读取 mtx 目录
adata = sc.read_10x_mtx('filtered_feature_bc_matrix/')
# 注意目录里要有 matrix.mtx、barcodes.tsv、features.tsv

# 方式2：读取 h5 文件（推荐）
adata = sc.read_10x_h5('filtered_feature_bc_matrix.h5')

一个小坑：老版本 Cell Ranger 的 features.tsv 第三列是基因类型（Gene Expression / CRISPR 等），读取后记得检查 adata.var['feature_types']，把非编码 RNA 或抗体标签过滤掉。我遇到过有人把 CRISPR 导入了主矩阵，结果聚类全乱套。

# 只保留基因表达数据
adata = adata[:, adata.var['feature_types'] == 'Gene Expression'].copy()

读取完后，看一眼 adata 的维度。我一般先跑 print(adata)，确认细胞数和基因数在合理范围——人血样本通常在 3000~10000 个细胞，20000~30000 个基因。

2. 过滤低质量细胞与基因

过滤的核心思路很简单：把空的、破的、异常的细胞扔掉，同时把那些几乎没人用的基因也清理掉。但阈值怎么定？嗯，这里要注意——不要凭感觉，要看分布。

先算三个关键指标：

n_genes_by_counts：每个细胞检测到多少种基因。太低说明可能是空液滴，太高可能是双细胞。
total_counts：每个细胞的 UMI 总数。太低说明 RNA 捕获失败或细胞已死。
pct_counts_mt：线粒体基因 UMI 占比。太高说明细胞膜破损，细胞质里的线粒体 RNA 漏出来了。

# 标记线粒体基因
adata.var['mt'] = adata.var_names.str.startswith('MT-')

# 计算 QC 指标
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False)

算完之后，我习惯先用小提琴图扫一眼分布：

sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             multi_panel=True)

  我的经验：别一上来就设 min_genes=200。先看图，如果大部分细胞集中在 300~600，那 200 就太松了，可以提到 300。我在一个项目里遇到过建库质量极差，大部分细胞只有 100~200 个基因，硬用 200 过滤就几乎把数据全扔了——后来降到 80 才保住主要细胞群。阈值要跟着数据走。

常用的过滤操作：

# 基础过滤（宽松版）
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# 自定义过滤（根据你的数据分布调整）
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 20, :]

# 看看还剩多少细胞
print(f'剩余细胞数：{adata.n_obs}')

对于基因过滤，min_cells=3 是通用标准。意思是某个基因至少在 3 个细胞里有表达，否则直接丢掉。这样能把那些偶然捕获的噪声基因清掉，减少计算量。

注意：过滤顺序很重要。先过滤细胞，再过滤基因。如果先删基因，有些本来有信号的基因可能因为细胞太少被误删。另外每次过滤后记得用 .copy() 保留原始备份，不然调参调错了得重新读数据。

3. 双细胞检测基础

双细胞就是一个液滴里包了两个细胞。10X 上机时细胞悬液浓度没控制好，双细胞率可能飙到 5%~10%。这些"假细胞"在聚类时会自己形成一个中间态群体，把你的细胞类型搞混。

目前最常用的工具是 Scrublet 和 DoubletFinder。我一般用 Scrublet，因为它快，而且对单个样本处理很方便。

方法	原理	优势	局限
Scrublet	模拟双细胞，比较局部密度	速度快，单样本独立运行	对高度相似的细胞类型不敏感
DoubletFinder	基于聚类后的人工双细胞	能利用整合后的多样本信息	计算慢，依赖聚类参数

Scrublet 的基本用法：

import scrublet

scrub = scrublet.Scrublet(adata.X, expected_doublet_rate=0.06)
doublet_scores, predicted_doublets = scrub.scrub_doublets()

# 把结果加到 obs 里
adata.obs['doublet_scores'] = doublet_scores
adata.obs['predicted_doublets'] = predicted_doublets

# 过滤双细胞
adata = adata[~adata.obs['predicted_doublets'], :].copy()

参数建议：expected_doublet_rate 根据你的上机细胞数估算。10X 官方文档给了一个参考：回收 5000 个细胞时双细胞率约 0.8%，10000 个细胞时约 3.2%。如果你不确定，设 0.05~0.08 作为起点就行。

Scrublet 会画一个直方图，显示模拟双细胞和真实细胞的分数分布。你会看到两个峰——左边是单细胞，右边是双细胞。如果分界不明显，说明样本里双细胞很少，或者细胞类型太相似。

我曾经踩过的坑：用了默认阈值没看分布，结果把一个稀有细胞群当成双细胞全删了。后来发现那个"双细胞峰"其实是某种激活的免疫细胞，本身基因数就高。所以跑完 Scrublet 一定先看直方图，确认阈值合理，再动手过滤。别让它"替你决定"。

双细胞检测放在过滤之后做，数据量小了跑得快。但如果你用的是 DoubletFinder，它需要聚类结果，所以要放到标准化和聚类之前——嗯，那又是另一套流程了。单细胞分析里，什么时候做什么操作，顺序本身就有讲究。

数据清理到这一步，你已经得到一份相对干净的表达矩阵了。接下来就可以放心地做标准化、降维、聚类——那些操作我们后面会细讲，但根基就是这一章的内容。预处理做扎实了，下游分析才能站得住。