MOFA数据准备：不同组学数据的标准化与缺失值处理
数据准备这一步，说实话，经常被初学者忽视。大家拿到多组学数据就急着跑MOFA模型，结果出来的因子要么解释不了，要么模型根本不收敛。我早期踩过这个坑，后来才发现——数据准备比建模本身更决定成败。
这章我们聊两件事：标准化和缺失值处理。说白了，就是让不同组学的数据在同一个“坐标系”里说话，同时把空缺的信息补上或合理忽略。我个人习惯把这步称为“整容级预处理”，整得好，模型里子面子都漂亮。
为什么标准化如此关键？
MOFA的核心假设是：不同组学数据背后共享同一个低维潜在空间。但你想想看，基因表达量的数量级是几万到几十万，蛋白质丰度是几到几千，甲基化beta值固定在0~1之间。如果不做标准化，高方差的数据会直接主导因子分解，低方差的组学就被淹没了。嗯，这可不是我们想要的。
| 组学类型 | 典型数值范围 | 常用标准化方法 |
|---|---|---|
| 转录组（RNA-seq） | 0 ~ 10⁶ (counts) | FPKM / TPM + log2转换 |
| 蛋白质组（质谱） | 0 ~ 10² (intensity) | log2转换 + 中位数中心化 |
| 甲基化（450K/850K） | 0 ~ 1 (beta值) | M值转换 logit(beta) |
| 代谢组（LC-MS） | 0 ~ 10⁸ (peak area) | log2转换 + Pareto缩放 |
| 微生物组（16S） | 0 ~ 10⁵ (counts) | CLR (中心化对数比) 转换 |
我个人强烈建议：先用log2转换把偏态分布拉近正态，再对每个特征做Z-score（均值为0，方差为1）。这样不同组学的特征贡献权重才能均衡。我见过太多人只做log2就扔进MOFA，结果第一个因子全是高表达基因，模型解释一塌糊涂。
缺失值处理：三种主流策略
多组学数据里缺失值几乎是必然的。为什么？因为不同检测平台的覆盖度不一样，有的基因在RNA-seq里有表达，在质谱里就是测不到。MOFA要求每个样本在每个组学里至少有一部分特征有值，但允许某些特征完全缺失（只要不是全部）。
处理思路分三条路，我按推荐度排个序：
- 直接建模（推荐）：MOFA的贝叶斯框架本身能处理缺失值——它只会用观测到的数据推断，缺失的地方自动跳过。但注意，每个样本至少要有10%的非缺失特征，否则模型不稳定。
- 插补（谨慎使用）：用kNN、SVD或mice等方法插补缺失值。但插补会引入偏差，尤其是高缺失率（>50%）的特征。我曾在项目中试过kNN插补，结果模型因子解读出了伪关联……后来排查了三天才发现是插补导致的虚假信号。
- 特征过滤（最保守）：删除缺失率超过某个阈值（比如80%）的特征或样本。简单暴力，但可能丢掉稀有但重要的信息。
千万别对log2后的负无穷值直接用0填充！比如log2(0) = -Inf，如果填0，等于把“未表达”误标为“中等表达”。我曾经见过一篇论文的补充数据里就这样处理，然后漂亮的因子载荷图下面藏着逻辑错误。正确的做法是：先标记缺失（NA），再用MOFA的固有处理机制。
实战流程：一份标准化与缺失值处理模板
下面是我常用的R语言处理流程。注意：MOFA官方推荐用R包MOFA2，数据格式是嵌套列表。
# 假设你已经有了三个矩阵：rna_expr, prot_expr, meth_beta
# 样本名必须一致，特征名可以是不同的
library(MOFA2)
# 步骤1：log2转换（对count数据先做CPM）
rna_cpm <- apply(rna_expr, 2, function(x) (x / sum(x)) * 1e6)
rna_log2 <- log2(rna_cpm + 1)
# 步骤2：Z-score标准化（按特征）
rna_scaled <- t(scale(t(rna_log2)))
# 蛋白质组：log2 + 中位数中心化
prot_log2 <- log2(prot_expr + 1)
prot_centered <- prot_log2 - apply(prot_log2, 1, median, na.rm = TRUE)
# 甲基化：转为M值
meth_M <- log2(meth_beta / (1 - meth_beta + 1e-8))
# 步骤3：合并成列表
data_list <- list(
RNA = rna_scaled,
Protein = prot_centered,
Methylation = meth_M
)
# 步骤4：创建MOFA对象（自动处理缺失值）
MOFAobject <- create_mofa(data_list)核心要点：标准化必须在组内按特征做，不能跨组学一起scale。因为不同组学方差结构不同，混合scale会扭曲组间差异。
缺失率检查与特征过滤
实战中，我习惯先检查每个特征的缺失率：
# 检查每个组学特征缺失率
check_missing <- function(mat) {
row_na <- rowSums(is.na(mat)) / ncol(mat)
return(row_na)
}
rna_miss <- check_missing(rna_scaled)
prot_miss <- check_missing(prot_centered)
meth_miss <- check_missing(meth_M)
# 去除缺失率 > 80% 的特征
rna_clean <- rna_scaled[rna_miss <= 0.8, ]
prot_clean <- prot_centered[prot_miss <= 0.8, ]
meth_clean <- meth_M[meth_miss <= 0.8, ]阈值怎么定？我一般用80%作为硬上限——如果某个特征在超过八成的样本里都没有测量值，那它带来的信息太少，反而增加计算噪音。要是你的样本量很大（比如500+），可以放宽到90%；样本少（<50）就收紧到70%。
数据准备流程图
图上能看到两条并行的处理路径：左边是缺失值处理（直接建模/插补/过滤），右边是缺失率检查（确保特征和样本的幸存率）。最终汇合到清洁数据，再交给
create_mofa。我个人建议：先做标准化，再做缺失值处理。因为有些标准化（比如Z-score）在缺失值存在时也能计算，但插补算法通常要求数据尺度一致。几个容易翻车的细节
- 样本顺序必须一致：不同组学的样本名要完全匹配，而且顺序在列表中要严格对应。MOFA不会自动帮你重排。
- 批次效应要不要去除？：如果不同组学是在不同批次测的，我建议先用limma或ComBat去批次，再做标准化。我遇到过某肿瘤项目，RNA-seq和蛋白质组是分开两年测的，结果第一个因子完美分离了年份——而不是生物学差异。
- 不要对二进制特征做Z-score：比如突变状态（0/1）或拷贝数变异（-1/0/1），直接用原始值，或者做中心化但不除以标准差。因为方差本身携带有用信息。
特别提醒：如果你用的是MOFA2的
create_mofa，它会自动检查数据格式并给出警告。但不会检查缺失率是否合理。我曾经跑了一个模型，有个组学里90%的特征缺失（因为测序深度太低），模型照样收敛了，但因子方差解释全是0——因为MOFA以为那些特征不存在。所以一定要在create_mofa之前，自己动手做缺失率筛查。这一套流程下来，你的数据就准备好给MOFA“吃”了。记住：数据准备不是锦上添花，而是地基。地基歪了，楼再高也得塌。嗯，别偷懒，按步骤来，后面你会感激现在的自己。