02 特征矩阵构建

ML 模型不吃 raw counts。它要"特征 × 样本"的数值矩阵,已归一化、已批次校正、已分割。这一章把 TCGA-LIHC 374 样本从 HTSeq counts 一路处理到一份可以直接喂给 03 章特征选择的 5000 维标准化矩阵。

为什么 ML 不能直接用 counts

raw counts 在 RNA-seq 里有三个性质让经典 ML 不舒服:

1. 均值-方差耦合。基因表达越高,counts 方差也越高。LR / SVM / 神经网络都假设特征方差大致一致(或至少有界),raw counts 直接破这个假设。
2. 跨样本测序深度差异。两个样本的 total counts 差一倍很常见,基因 A 在样本 1 的 100 counts 和样本 2 的 200 counts 不是同一回事。
3. 极端右偏分布。housekeeping 基因 counts 几千上万,大多数基因在百量级以下,直接喂模型权重会被高表达基因主导。

差异分析专栏(DESeq2 / edgeR / limma-voom)用 GLM 显式建模这些性质,所以 raw counts 直接喂没问题。ML 没有这种天然的统计模型护栏,我们要在数据进模型之前就把分布修整好。

三种归一化怎么选

方法	适用场景	输出范围	ML 友好度
vst() (DESeq2)	样本数适中(20-1000)、想保留表达水平差异	大致 [3, 17]	⭐⭐⭐⭐⭐ 默认推荐
rlog() (DESeq2)	样本数 < 30,或表达低	类似 vst	⭐⭐⭐⭐ vst 速度优势更大
log2(TPM+1)	跨数据集对比、单个样本评分	[0, 16]	⭐⭐⭐⭐ 但跨样本归一化弱
log2(FPKM+1)	兼容老论文 / IGV 浏览	类似 TPM	⭐⭐ 别用 — 同基因不同样本不可比
voom() (limma)	limma 流程内部,作为预处理	log-CPM + 权重	⭐⭐⭐ 工具特定

实操约定:默认用 vst(blind=FALSE),既保留 design 信息又能 ML 友好。下面的脚本用的就是这个。

library(DESeq2)
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ tissue_type)
vsd <- vst(dds, blind = FALSE)
expr_vst <- assay(vsd)  # 行=基因, 列=样本, 范围 ~3-17

为什么 blind=FALSE?blind=TRUE 会忽略 design 公式,在批次效应大的数据上反而把生物学信号也磨掉。我们后面要分 train/test,不希望归一化阶段就提前看到 test 信号,所以 design 设成 ~ tissue_type(只看核心生物学变量)是合理折中。

选高变基因(HVG)— 把 60K 基因压到 5000

TCGA-LIHC 经过 Ensembl 注释后是 60660 个基因 × 374 样本。直接喂 ML 不仅慢,还有大量"沉默基因"在拉低信噪比。先选 top 5000 高变基因(HVG):

gene_var <- apply(expr_vst, 1, var)
hvg_top5000 <- names(sort(gene_var, decreasing = TRUE))[1:5000]
expr_hvg <- expr_vst[hvg_top5000, ]

为什么是 5000 不是 1000 或 10000?经验:

• < 1000:容易漏掉中等表达但有判别力的基因
• 1000-5000:大多数 ML 任务的甜蜜区,03 章特征选择会进一步压到 50-200

• 10000:噪声基因占比上升,03 章选择速度慢但收益不明显

5000 也是 TCGA-LIHC 在 03 / 04 / 06 章贯穿的标准维度,后续工具的 demo data 都从这个 5000 维子集里抽取。

批次效应 — ComBat-seq vs 不做

TCGA 数据里"批次"主要是 plate、center、年份。如果你的研究问题恰好和这些混淆(比如 2010-2012 vs 2015-2017 的肿瘤样本),不做批次校正会把"年份效应"误判为"生物学差异"。

判断要不要做批次校正的检查清单:

1. 在 PCA 上按 batch 染色 — PC1 / PC2 是不是把不同 batch 拉开?
2. 在差异分析里 condition 系数被 batch 系数压住?
3. batch 与你的 outcome label 高度相关(用 chisq.test 看)?

任何一条命中就要做。TCGA-LIHC 374 样本恰好不是高度混淆(plate 平衡得不错),所以本专栏默认不做 ComBat-seq。但下面是做的写法:

library(sva)
counts_adj <- ComBat_seq(counts, batch = coldata$plate,
                          group = coldata$tissue_type)
# 然后再走 vst → HVG → split 流程

⚠️ 关键约定:ComBat-seq 一定在 split 之前做(否则 train 上 fit 的模型会"偷窥"test)。如果用 ComBat-seq,后续整个数据集走的都是校正后版本,test 也是。这与下一节"特征预处理只在 train 上 fit"是不同的层次。

Train / Val / Test 三分割

关键约定有四条:

1. stratified split(按 outcome label 分层)— 防止 test 集类比例失真
2. 样本级别独立(同一病人多个样本不能跨 split — TCGA 一般不会有,但 GEO / 自建数据要小心)
3. 早做、做一次、持久化(写到 derived/ 目录,后续 03/04/05/06 章引用同一份 split)
4. 特征预处理(z-score / 缺失填充)在 train 上 fit,test 上 transform — 先 split 再 scale,顺序反了就泄漏

ML W2 里 ml-classifier 工具自动按这个约定做 split,本章脚本也是同样约定。比例约定:70/15/15 (二分类) 或 60/20/20(多分类需要更多 val 调参):

set.seed(42)
library(caret)
idx_train <- createDataPartition(coldata$tissue_type, p = 0.7, list = FALSE)
train_set <- coldata[idx_train, ]
holdout <- coldata[-idx_train, ]
idx_val <- createDataPartition(holdout$tissue_type, p = 0.5, list = FALSE)
val_set <- holdout[idx_val, ]
test_set <- holdout[-idx_val, ]

写出三份 sample id 表到 derived/,后续章节直接读。BioF3 W2 的 demo data 就是按这个约定生成的。

缺失值处理

蛋白组、甲基化、临床表会有 NA。TCGA-LIHC 表达矩阵几乎没有 NA(测序数据完整),但临床表(stage / grade / OS)会有。约定:

场景	默认做法
整列 NA 比例 > 30%	删掉这列(整个特征不可用)
整列 NA 比例 5-30%	用 train 上的 median / mode 填充,test 用同一个值
整列 NA 比例 < 5%	同上,或者删掉这些 NA 行
同时多个特征 NA 模式相关	用 mice 或 missForest 多重插补(更复杂)

注意:插补值要从 train 上计算,绝不能从全数据集(否则泄漏)。

# 1. fit on train
median_train <- median(train_set$age, na.rm = TRUE)

# 2. apply to all sets
train_set$age[is.na(train_set$age)] <- median_train
val_set$age[is.na(val_set$age)] <- median_train
test_set$age[is.na(test_set$age)] <- median_train

类不平衡处理

TCGA-LIHC 的 Tumor:Normal 大约 4:1 不算极端,但乳腺癌亚型 PAM50 里 Normal-like 只占 5% 就是真不平衡了。处理选项:

方法	适用场景	代价
不做(用 PR-AUC 评估)	比例 ≥ 1:5	最稳,模型可解释性最好
class weight	比例 1:5 - 1:20	几乎所有 sklearn / caret 都内置 class.weights = "balanced"
SMOTE / ADASYN	比例 < 1:20	合成样本可能"造假",降低真实泛化能力
降采样多数类	数据量极大	损失信息但训练快
focal loss	深度学习场景	需要自定义损失函数

实务建议(对应教程 module04):默认不做,在 ml-classifier 工具的 caret 调用里直接传 class.weights。SMOTE 看起来"漂亮"但实际场景下经常带来过度自信的模型,生信论文里 reviewer 见过太多 SMOTE 滥用。

配套脚本

ml02_features_sci.R(本专栏 W3 配套):

Rscript scripts/machine-learning/ml02_features_sci.R

输出到 ~/biof3-data/tcga-lihc/derived/:

• expr_vst_hvg5000.csv(5000 × 374)
• coldata.csv(374 × 临床列)
• train_ids.csv / val_ids.csv / test_ids.csv
• feature_stats.csv(每个 HVG 的方差 / 表达均值,用于 03 章 sanity)

03/04/05/06 章直接 read.csv() 读这几份,无需重跑全流程。

在线工具对接

把这一章产物直接拖到在线工具:

• 把 expr_vst_hvg5000.csv 转置后(行=样本)→ 上传到 ml-feature-select 跑特征选择
• 把 03 章 共识的 top 200 → 上传到 ml-classifier 跑 4 算法对比

常见坑

❌ vst 在 split 之后做 — 一般情况下 vst 在 split 之前做 OK(它是非监督归一化),但更严谨可以在 train 上 vst(dds) 然后 varianceStabilizingTransformation() 把 test 应用同样的 trans → 实务上对 ML 影响 < 1%。BioF3 默认 split 之前做,简化流程。

❌ HVG 选择泄漏 — 在全数据集选 HVG 再 split,这是常见但有泄漏的做法。严格做法是 train 上选 HVG,然后用同一个 gene list 提取 test。本专栏默认做"先选 HVG 后 split",因为 HVG 选择不直接看 outcome label,泄漏风险低于特征选择(03 章会显式区分这两件事)。

❌ train/test 比例 50/50 — 数据量小(< 100)时常见错觉。70/30 或 80/20 更标准,test 集太大其实不会让结果更可信。

❌ k-fold CV 之外再额外做 test 集 — 注意 CV 是训练阶段的调参工具,test 集是终验。两者都要,不互相替代。

❌ 拿 Tumor 样本做 ComBat-seq 但用 Normal 做 group — group 参数要传"想保留的生物学变量",不是分组 1/2。

本章状态

✅ Wave 3 正文完成(2026-05-27)。配套脚本 ml02_features_sci.R 在 W2-A 阶段已生成 derived 数据,TCGA-LIHC 374 样本 × 5000 HVG 已就位。