2. 质量控制基础：RNA-seq数据质量指标、缺失值处理与过滤策略
做RNA-seq分析，第一步就是质量控制。说白了一句话：垃圾数据进，分析结果必然是垃圾。
我刚开始做单细胞项目那会儿，上来就闷着头跑聚类，结果画出来的UMAP图丑到没法看。后来一查才知道，样本里死细胞比例太高，测序质量一塌糊涂。从那以后，“先做质控”这件事就成了我绕不开的本能。
2.1 核心质量指标
质量控制有哪些指标要盯？我个人习惯分成三类：
- 测序质量指标——读长质量得分、GC含量、接头污染
- 比对相关指标——比对率、唯一比对率、链特异性
- 基因检测水平指标——检测到的基因数量、文库复杂度、ERCC spike-in
这里我列一张表，方便你对照检查：
| 指标类别 | 具体指标 | 正常范围参考 | 常见问题 |
|---|---|---|---|
| 测序质量 | Phred得分 (Q30) | >80% | 低Q30提示测序质量差 |
| 测序质量 | GC含量 | 40%-60% | 偏离可能提示污染或偏好性 |
| 比对 | 总比对率 | >70% | 低比对率表示rRNA污染或样本质量问题 |
| 比对 | 唯一比对率 | >60% | 重复率高可能是PCR过度扩增 |
| 基因检测 | 检测到的基因数 | 因文库不同而异 | 基因数过低提示低质量文库 |
| 基因检测 | 文库复杂度 | 接近1为理想 | 复杂度低表示PCR重复多 |
FastQC 做初始质量检查，再用 MultiQC 汇总所有样本的质量报告。我一般会跳过那些Q30低于70%的样本，省得后面浪费时间。
2.2 缺失值处理
RNA-seq数据中缺失值是家常便饭。基因表达矩阵里经常有大片的0或者NA值。你想想看，一个基因在部分样本里不表达，这到底是真实的生物学特征，还是测序没测到？
我见过最离谱的情况是，有人直接把所有含有NA的基因全部删除，结果损失了超过40%的基因集。嗯，这么做肯定不对。
常用的缺失值处理策略有三种：
- 忽略——在统计检验中自动排除NA，适用于偶尔缺失
- 填充——用均值、中位数、kNN等方法填充缺失值
- 过滤——删除缺失比例过高的基因或样本
我个人的建议是：先过滤，再填充。具体怎么操作？
过滤掉那些在超过80%样本里都没有表达的基因。剩下的缺失值，我习惯用表达矩阵的基因均值填充。为什么？因为单细胞数据里的缺失往往不是因为实验问题，而是因为dropout效应。用均值至少保留了整体分布特征。
R里面操作起来很简单：
# 过滤：保留在至少20%样本中有表达的基因
keep_genes <- rowSums(expr_matrix > 0) >= 0.2 * ncol(expr_matrix)
expr_filtered <- expr_matrix[keep_genes, ]
# 填充缺失值：用基因均值
expr_filled <- apply(expr_filtered, 1, function(x) {
x[is.na(x)] <- mean(x, na.rm = TRUE)
return(x)
})2.3 过滤策略
过滤是质量控制的主战场。你需要决定：哪些样本、哪些基因可以留下？
常见的过滤策略可以归纳为三层：
- 样本级过滤：删除总表达量过低的样本、删除异常高表达的样本、删除Q30不合格的样本
- 基因级过滤：删除在绝大多数样本中不表达的基因、删除mtRNA比例过高的基因、删除ERCC spike-in基因（如果不需要）
- 综合过滤：用PCA检查离群样本、用库大小标准化前的热图观察全局表达量差异
我曾经在一个胃癌转录组项目里，看到某个样本的线粒体基因比例飙到了35%。正常值也就5%左右。一查才知道，那管样本在提取时细胞破裂，线粒体RNA大量释放。我直接在过滤步骤里把它干掉了，后面结果才正常。
下面这段代码展示了我在项目里常用的过滤流程：
# 样本级过滤
sample_libsize <- colSums(expr_filtered)
sample_q30 <- pheno_data$Q30_percent # 假设置表里有Q30数据
good_samples <- sample_libsize > 100000 &
sample_libsize < 10000000 &
sample_q30 > 80
expr_filtered <- expr_filtered[, good_samples]
# 基因级过滤
gene_count <- rowSums(expr_filtered > 0)
keep_genes <- gene_count >= 0.05 * ncol(expr_filtered)
expr_final <- expr_filtered[keep_genes, ]
# 输出过滤后的矩阵
dim(expr_final)这里需要注意过滤阈值的设置。为什么选5%？说白了就是保留那些至少在5%样本中有表达的基因。如果样本有50个，那就是至少要在2.5个样本中有表达（取整就是2个）。这是一个经验值，源自经典的Deseq2过滤推荐。但具体项目里你可以微调。
核心原则：
- 宁可多留，不要少删——过滤太强会丢失生物信号
- 必须有证据——不要因为“看起来异常”就随便删样本，看看PCA图、热图、QC报告
- 记录过滤步骤——你删了多少基因、多少样本，必须记下来！我见过太多人最后复现不出来就是因为忘记自己删了什么
2.4 质量控制流程总览
为了让你对整个质量控制流程有个全局印象，我画了一张流程图。你一看就明白先做什么、后做什么。
流程很清晰吧？从原始数据开始，先做质量评估，然后分两条线走：一条处理缺失值，一条做基因和样本过滤。两条线在最终的表达矩阵汇合。这时候你的数据就可以安心往下跑了。
rna_qc() 函数，输入表达矩阵和样本注释，输出过滤后的数据和QC报告。这样既省事又不容易出错。
2.5 常见避坑指南
最后，分享几个我踩过的坑：
- 不要盲目相信默认阈值——FastQC的“通过/警告/失败”判定非常保守，尤其对于长读长测序数据，很容易标“失败”。我一般自己看Q30和GC含量曲线，而不是直接拿判据删数据。
- 缺失值填充不是万能的——如果某个基因在90%的样本里都是0，你用均值填充，会强行让那些样本的基因表达水平被抬高。我曾经就因为这样搞出了一个假阳性的差异基因集，查了三天才找到原因。
- 过滤的顺序有讲究——先过滤样本，再过滤基因。因为如果先按基因过滤，你可能会保留那些只在异常样本里高表达的基因，后期分析会很头疼。
质量控制说到底就一句话：把假的数据当真的分析，结果就是假的。 认真过一遍质控流程，省下后面排查Bug的时间，这笔账怎么算都值。