数据预处理:不同组学数据格式化、标准化与归一化

说实话,多组学整合项目里,80%的时间都花在数据预处理上。这话一点都不夸张。我刚开始做MOFA的时候,以为算法调参才是重头戏,结果数据格式不统一直接报错,反复调试了两天。从此以后,我养成了一个习惯:先花半天把数据理顺,再谈建模。

说白了,MOFA要求的输入就是三个核心矩阵——样本×特征、样本×特征、样本×特征。但实际拿到的数据呢?基因表达是矩阵、甲基化是beta值、蛋白组是强度值、代谢组是峰面积。你怎么把它们对齐?这就是本章要解决的核心问题。

数据预处理核心流程 原始多组学数据 RNA/甲基化/蛋白/代谢 维度统一 基因/探针/特征对齐 标准化 Z-score / 均值方差 归一化 Min-Max / 分位数 / Log MOFA输入矩阵 N×M₁, N×M₂, N×M₃ 核心原则 • 样本维度必须一致 • 特征维度可不同 • 缺失值不能留NaN • 数值范围要接近 • 保持数据结构完整

第一步:维度的统一

不同组学数据能不能直接丢进MOFA?答案是不能。你想想看,基因表达矩阵是20000行×100列,甲基化beta值是450000行×100列,蛋白质组是5000行×95列——样本数都不一致。MOFA要求所有组学的样本维度必须完全对齐。我遇到过最坑的情况是,样本命名格式不统一,有的用“TCGA-XX-1234”,有的用“XX1234”,肉眼看着一样,程序一跑就崩。

维度统一的标准动作:

  • 样本对齐:取所有组学样本的交集,去掉单组学独有的样本
  • 行标识统一:样本名全部转为大写或全部使用统一分隔符
  • 特征命名规范:基因名用ENSEMBL ID或Symbol,探针用厂商编号
  • 缺失样本处理:某个组学样本缺失 → 要么删样本,要么插补

第二步:不同组学的具体格式要求

每种组学数据都有它的脾气。做格式转换时,你得分门别类去处理。

组学类型 常见原始格式 MOFA推荐格式 注意事项
基因表达 Counts / FPKM / TPM log2(TPM+1) 矩阵 标准化前先做log转化
DNA甲基化 Beta值 / M值 Beta值矩阵 探针到基因的映射需要提前处理
蛋白组学 强度值 / iBAQ log2(强度) 矩阵 缺失值比例高,注意过滤
代谢组学 峰面积 / 浓度 log转化后矩阵 批次效应明显,需要额外校正
微生物组 OTU表 / 相对丰度 CLR转化后矩阵 稀疏数据,零值处理要慎重

我个人习惯的做法是:先写一个统一的格式转换函数,每个组学单独调,最后输出三个规整的CSV文件。这样后续标准化和归一化就省事了。

第三步:标准化——让不同组学站在同一起跑线

标准化的目的是消除特征之间的量纲差异。说白了,你不能让基因表达量动辄几千、上万,而甲基化beta值只在0到1之间浮动。MOFA的因子载荷对数值尺度非常敏感,不标准化的话,载荷全被高变异特征带跑了。

最常用的三种标准化方法:

  1. Z-score标准化(推荐):

    将每个特征减去均值,再除以标准差。结果服从均值为0、方差为1的分布。这是我默认使用的方法,适用于大部分连续型组学数据。

  2. 最大最小归一化:

    将数据缩放到[0, 1]区间。当数据没有异常值时,效果还可以。但我遇到过极端值把整个压缩空间占满的情况,结果其他特征的信息全丢了。

  3. 分位数标准化:

    强制让所有样本的分布一致。适用于批次效应明显的场景。不过要注意,这样做会破坏原始数据的绝对数值关系。

# Python实现:三种标准化对比
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler, MinMaxScaler

# 读取基因表达矩阵
expr = pd.read_csv('expression.csv', index_col=0)

# Z-score标准化(按特征,即按基因)
scaler_z = StandardScaler()
expr_z = scaler_z.fit_transform(expr.T).T  # 转置再转回来

# Min-Max归一化
scaler_mm = MinMaxScaler()
expr_mm = scaler_mm.fit_transform(expr.T).T

# 分位数标准化(按样本)
from scipy.stats import rankdata
def quantile_normalize(df):
    rank_mean = df.stack().groupby(
        df.rank(method='first').stack().astype(int)
    ).mean()
    return df.rank(method='min').stack().astype(int).map(rank_mean).unstack()

expr_qn = quantile_normalize(expr)

我的经验:如果组学之间动态范围差异极大(比如基因表达范围0-10000,甲基化0-1),建议先做log转化,再做Z-score。这样既能控制方差,又能保留生物学差异。

第四步:归一化——样本间的校正

标准化是按特征(基因/蛋白质)做的。归一化是按样本做的。为什么要做样本归一化?因为不同样本的测序深度、捕获效率、检测灵敏度不同。你不校正的话,MOFA可能学到的是技术噪声而不是生物学信号。

我遇到过最典型的案例:一个批次里的样本测序深度是另一个批次的两倍,结果MOFA的第一个因子全在捕获批次差异。后来我做了TMM归一化,才把真正的生物学信号分离出来。

常用的样本归一化方法:

方法 适用场景 说明
CPM/RPKM/TPM RNA-seq / 蛋白组 消除测序深度影响,保留比较性
TMM/边缘R 差异表达分析场景 考虑文库组成效应,比CPM更稳健
分位数归一化 甲基化芯片 强制分布一致,适用于探针级别数据
总强度归一化 代谢组 / 蛋白组 每个样本除以总信号强度

避坑指南:我踩过的数据预处理坑

我曾经犯过以下错误:

  • 标准化顺序搞反:先做样本归一化,再做特征标准化。结果样本归一化把生物学差异都抹掉了,后面标准化又强行拉回来,模型结果毫无意义。
  • 缺失值处理不当:直接删掉有缺失的行,导致样本大量减少。后来改用了kNN插补,效果好了很多。
  • 忽略零值:有些组学数据零值不是缺失,而是真实表达。如果盲目替换成NaN或用插补,会引入严重偏差。
  • 批次效应残留:做了归一化但没做批次校正。MOFA虽然能捕获批次因子,但会影响下游解释。我建议先用ComBat或limma校正批次,再进入预处理流程。
# 完整的预处理流程示例(Python)
def preprocess_for_mofa(expr_df, meth_df, prot_df):
    """
    三组学数据统一预处理
    输入:三个DataFrame,行=样本,列=特征
    输出:三个预处理后的矩阵
    """
    # 1. 样本对齐:取交集
    common_samples = list(set(expr_df.index) & set(meth_df.index) & set(prot_df.index))
    common_samples.sort()
    print(f"对齐后样本数:{len(common_samples)}")
    
    # 2. 缺失值处理(kNN插补)
    from sklearn.impute import KNNImputer
    imputer = KNNImputer(n_neighbors=5)
    
    # 3. 对高动态范围数据做log转化
    expr_log = np.log2(expr_df.loc[common_samples] + 1)
    prot_log = np.log2(prot_df.loc[common_samples] + 1)
    
    # 4. 特征标准化(Z-score)
    from sklearn.preprocessing import StandardScaler
    scaler = StandardScaler()
    expr_scaled = scaler.fit_transform(expr_log)
    meth_scaled = scaler.fit_transform(meth_df.loc[common_samples])
    prot_scaled = scaler.fit_transform(prot_log)
    
    # 5. 转换为DataFrame,保持行列名
    expr_out = pd.DataFrame(expr_scaled, index=common_samples, columns=expr_df.columns)
    meth_out = pd.DataFrame(meth_scaled, index=common_samples, columns=meth_df.columns)
    prot_out = pd.DataFrame(prot_scaled, index=common_samples, columns=prot_df.columns)
    
    return expr_out, meth_out, prot_out

# 调用示例
expr_ready, meth_ready, prot_ready = preprocess_for_mofa(expr, meth, prot)

什么时候需要标准化?什么时候没必要?

这个问题我问过自己很多次。答案其实很简单:看你用什么数据。

  • 基因表达、蛋白组、代谢组 → 必须标准化。因为数值范围和方差差异巨大。
  • 甲基化beta值 → 可标准化也可不标准化。因为beta值天然在0-1之间,特征间量纲一致,直接输入也能收敛。但我个人建议还是做Z-score,能让因子解释更简洁。
  • 突变矩阵(0/1) → 不需要标准化。因为已经是二值数据,标准化反而会破坏结构。
  • 临床指标(年龄、分期等) → 连续变量标准化,分类变量编码后不用。

核心原则总结:

  • MOFA对数据尺度敏感 → 所有连续型组学必须做标准化
  • 标准化选Z-score,基本不会出错
  • 样本归一化在标准化之前做,效果更佳
  • 缺失值不要直接删,尤其是高维组学数据
  • 预处理完成后,务必保存一份快照,方便复现

嗯,这部分其实不难,但确实是整个MOFA流程里最磨人的环节。搞定了数据预处理,后面建模调参就顺风顺水了。我个人建议,每次处理完数据后,跑一个简单的相关性热图检查一下,看看组学内和组学间的结构是否合理。如果发现异常,赶紧回头排查,别等到建模了才发现数据有问题。