📚 课程目录第2章基因表达数据预处理

基因表达数据预处理：从raw count到干净特征

我经常会被问到这样一个问题："RNA-seq数据拿过来，直接丢进模型行不行？"

说实话，不行。至少在我经手的几十个项目里，没有哪个是这么干的。

原始的表达数据，说白了就是一堆counts数。里面有各种噪声、批次效应、还有缺失值。你要是直接扔给神经网络，模型大概率会学偏。我刚开始做这个的时候也踩过坑，拿原始数据跑了一轮，结果调参调了三天，收敛效果依然很差。后来才发现，问题出在数据预处理阶段。

预处理不是可选项，它是模型训练的"地基"。地基不稳，上层建筑再漂亮也是白搭。

缺失值处理：别让"NA"毁了你的模型

RNA-seq数据里出现缺失值，原因很多。可能是测序深度不够，某个基因恰好没被捕获；也可能是低表达基因在某个样本里压根检测不到。你想想看，如果有个基因在80%的样本里表达值都正常，就几个样本是NA，你直接去掉整个基因，信息损失太大了。

我个人的处理原则是：

基因维度的缺失率 > 50%：直接删掉这个基因，因为它提供的信息太少。
样本维度的缺失率 > 30%：删掉这个样本，可能存在实验质量问题。
少量缺失（< 5%）：用均值、中位数或者KNN插补。

我的习惯

对于深度学习的输入，我通常用基因在全部样本中的中位数来填充缺失值。为什么是中位数？因为表达数据偏态严重，均值容易被极端值拉偏。中位数更稳健。

避坑提醒

千万不要在插补之前做标准化！我曾经见过一个错误案例，先用Z-score标准化，发现缺失值被算成了0，插补完整个分布都扭曲了。顺序搞反，后果很严重。

归一化：消除文库大小差异

归一化解决的是"文库大小不同"的问题。同一个基因，在A样本里测了100万条reads，在B样本里测了200万条。如果不做归一化，模型会误以为B样本里这个基因表达量翻倍了。其实只是测序深度不同而已。

常用的方法有几种，我整理了一张表：

方法	原理	适用场景	我推荐的场景
CPM (Counts Per Million)	每个基因的count除以总count，再乘1e6	快速粗略比较	初筛、探索性分析
TPM (Transcripts Per Million)	考虑基因长度，先除以基因长度再归一化	跨基因比较	需要比较不同基因的表达量
FPKM/RPKM	类似TPM，早期常用，现已逐渐淘汰	历史遗留项目	尽量别用，有偏差
DESeq2 median of ratios	基于每个基因与几何均值的比值的中位数	差异表达分析、深度学习	最推荐
TMM (edgeR)	基于 trimmed mean of M-values	差异表达分析	适合与edgeR配合

我个人习惯用DESeq2的median of ratios方法。原因很简单：它对高表达和低表达基因的偏差校正最鲁棒。尤其是当你的数据里存在大量低表达基因时，CPM这种粗暴方法会让方差变得非常大。

实操命令

在R里，用DESeq2做归一化其实就一行代码：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~ condition)
dds <- estimateSizeFactors(dds)
normalized_counts <- counts(dds, normalized = TRUE)

这里的关键是estimateSizeFactors。它会计算出每个样本的"size factor"，相当于一个缩放系数。归一化后的counts，就可以用来做后续的机器学习任务了。

小技巧

拿到normalized_counts之后，我习惯先取log2变换：log2(normalized_counts + 1)。这能让偏态分布变得更接近正态，对模型收敛很有帮助。

标准化：让所有基因站在同一起跑线

标准化和归一化不是一回事。归一化是消除文库大小差异，标准化是让每个基因的表达值分布对齐——通常是零均值、单位方差。

为什么需要标准化？你想想看，有的基因表达量很高（比如管家基因），均值为几千甚至上万；有的基因表达量很低，均值不到10。梯度下降的时候，高表达基因的梯度会主导更新方向，低表达基因的信息就被淹没了。

标准做法是Z-score：

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(normalized_log2_counts)

标准化之后，每个基因的均值为0，标准差为1。模型就不会"偏袒"高表达基因了。

注意分割

标准化参数（均值和标准差）一定要在训练集上计算，然后用同样的参数去变换验证集和测试集。千万别在全体数据上fit_transform。我早年犯过这个错误，导致模型在验证集上的表现虚高，上线后直接崩了。

批次校正：整合多个实验的数据

实际项目中，你很少能在一个批次里收完所有样本。不同时间、不同操作员、不同试剂批次，都会引入系统性的偏差。这种偏差叫"批次效应"。

有一次我做一个项目，数据来自三个不同实验室。不校正批次的话，样本在PCA图上直接按实验室分成了三群。模型学到的不是生物学差异，而是实验室差异。

常用的批次校正方法：

ComBat (sva包)：基于贝叶斯框架，是目前最流行的方法之一。对已知批次效果很好。
limma (removeBatchEffect)：线性回归的方法，简单高效。
Harmony：基于软聚类和迭代校正，适合大规模单细胞数据，但也可以用于bulk RNA-seq。
scVI / scGen：基于变分自编码器的深度学习方法，适合复杂场景。

R代码（ComBat）

library(sva)
adjusted_counts <- ComBat(dat = log2_normalized_counts,
                          batch = batch_labels,
                          mod = model.matrix(~ condition),
                          par.prior = TRUE,
                          prior.plots = FALSE)

注意mod参数里要包含你的生物学变量，这样做是为了避免校正掉真正的生物学差异。

我的经验

批次校正不是万能的。如果批次和生物学分组完全共线（比如所有case在批次1，所有control在批次2），那谁也校正不了。这种情况我建议你找实验设计的问题，别在数据分析上硬扛。

预处理检查清单

每次我拿到一批RNA-seq数据，都会按下面的清单过一遍：

缺失值检查：统计每个基因和每个样本的缺失率，决定删除还是插补。
低表达基因过滤：通常我会去掉在超过80%的样本中count < 10的基因。
归一化：用DESeq2的median of ratios方法，得到normalized counts。
log2变换：log2(normalized_counts + 1)，稳定方差。
标准化：Z-score，每个基因零均值单位方差。
批次校正：如果存在多个批次，用ComBat或者limma处理。
最终检查：PCA可视化，看看样本是否按生物学分组聚集，而不是按批次或技术变量聚集。

一句话总结

预处理的目的，就是让数据干净、可比、无偏。你花在预处理上的每一个小时，都会在模型训练和调参阶段加倍还给你。