模拟数据集生成：用 NumPy 搭建多组学数据

做 MOFA 模型调优之前，我最头疼的一件事就是找不到合适的数据集。公开数据集要么样本量太大，要么预处理不干净，跑一次实验光清理数据就得搭进去半天。所以我自己动手，用 NumPy 生成一套模拟多组学数据。听起来像造假？其实不是，模拟数据能帮你精准控制信号强度、噪声水平、组学间的共享结构，非常适合验证算法效果。

转录组：模拟基因表达矩阵

转录组数据通常是连续值，比如 RNA-seq 的 TPM 或 FPKM。我一般用负二项分布来模拟，因为真实表达数据通常是过离散的。不过为了简单，这里先用正态分布替代，加一点随机噪声就行。

import numpy as np

np.random.seed(42)
n_samples = 100
n_genes = 500
# 假设有 3 个潜在因子（比如处理条件）
latent = np.random.randn(n_samples, 3)
# 基因载荷矩阵
loadings = np.random.randn(3, n_genes) * 0.5
# 生成表达矩阵
expr = latent @ loadings + np.random.randn(n_samples, n_genes) * 0.1
# 转成正数（模拟表达值）
expr = np.abs(expr)
print(f"转录组数据形状: {expr.shape}")

我习惯把噪声标准差设成 0.1，这样信号强度大概占方差的三分之二，不至于一眼就看穿模式。你在实际调参时可以试试调到 0.3 或 0.5，观察 MOFA 对噪声的鲁棒性。

甲基化：用 Beta 分布模拟 β 值

DNA 甲基化数据通常是 0 到 1 之间的 β 值。Beta 分布天然适合这种区间，而且它的形状参数能控制尖峰和拖尾。我早年做表观遗传项目时，发现直接硬截断正态分布会引入假阳性，后来就一直用 Beta 分布来模拟甲基化矩阵。

n_cpgs = 300
# 用潜在因子生成甲基化倾向
methyl_prop = np.random.beta(2, 2, (n_samples, n_cpgs))
# 加入因子结构
shared_factor = latent[:, :2] @ np.random.randn(2, n_cpgs) * 0.3
methyl_beta = np.clip(methyl_prop + shared_factor, 0.01, 0.99)
print(f"甲基化β值范围: {methyl_beta.min():.2f} ~ {methyl_beta.max():.2f}")

注意 Beta 分布的参数 a=2, b=2 时比较均匀，如果改成 a=1, b=1 就变成均匀分布，更适合极端情况测试。

蛋白质组：模拟对数强度值

蛋白质组学数据通常来自质谱检测，表达值一般是对数转换后的强度，近似正态分布。我直接复用转录组的结构，但把噪声加大一点，模拟技术差异。

n_proteins = 200
protein = latent[:, :2] @ np.random.randn(2, n_proteins) * 0.4
protein += np.random.randn(n_samples, n_proteins) * 0.3
protein = np.abs(protein) + 1  # 避免零值
print(f"蛋白质组数据形状: {protein.shape}")

这里特意只用了前两个潜在因子，和转录组的三个因子不完全一样。为什么要这样？因为真实的多组学数据中，不同组学受共享因子影响的程度不同，留一点私有因子才能更真实地模拟 MOFA 需要分解的结构。

合并成多组学列表

MOFA 要求输入是一个列表，每个元素是样本 × 特征的矩阵。我把三个组学压成一个 list，再分别标准化。

data_list = [expr, methyl_beta, protein]
# 标准化（按组学独立做，不能跨组学混在一起）
data_scaled = []
for mat in data_list:
    mean = mat.mean(axis=0)
    std = mat.std(axis=0)
    std[std == 0] = 1  # 防止除零
    data_scaled.append((mat - mean) / std)

# 保存为 NumPy 文件备用
np.savez('multi_omics_sim.npz', *data_scaled)

标准化这一步很多人会忽略，但 MOFA 假设每个组学的特征服从平均 0 方差 1 的分布。如果不标准化，数值大的组学（比如蛋白质组的强度）会主导因子分解，结果就可能偏差。

参数对比一览

组学类型分布假设样本数特征数共享因子数噪声水平
转录组正态（绝对值）10050030.1
甲基化Beta (2,2)+因子偏移10030020.3
蛋白质组对数正态（简化）10020020.3

上面这组参数只是起点，我个人建议你至少跑三次不同噪声级别的实验，看看 MOFA 的因子一致性是否稳定。我曾经因为噪声设太低，导致模型完美拟合垃圾信号，浪费了两周时间。嗯，避坑指南：模拟数据一定要包含噪声和缺失因子结构，别做得太完美。

流程可视化

我画了一张示意图，帮你快速理解整个模拟流水线。

多组学模拟数据生成流程 共享潜在因子（3个） 转录组 500 基因 × 100 样本 DNA甲基化 300 CpG × 100 样本 蛋白质组 200 蛋白 × 100 样本 Z-score 标准化（按组学） Z-score 标准化（按组学） Z-score 标准化（按组学） MOFA 输入列表 3 个矩阵，每个 100 样本

关键要点：三个组学用了不同分布和不同共享因子数目，目的就是让模拟数据既存在共享结构，又保留各自的私有噪声。这样 MOFA 才能展现出它分解“共享因子+私有因子”的核心能力。

常见踩坑与我的解决经验

  • 因子数目不匹配：我在生成甲基化数据时只用了两个共享因子，转录组用了三个。如果后面做 MOFA 训练时统一设成 K=3，甲基化那一层的第三个因子就会退化或变噪声。我建议你分别设不同因子数，看看 MOFA 是否能自动识别哪些因子共享度高。
  • 缺失值处理：真实多组学里总有缺失值，但模拟数据里我没加。你可以主动删除 5% 的甲基化数据，观察 MOFA 的 imputation 效果。我之前做过一次 30% 缺失率的测试，结果因子结构和完整数据几乎一致，说明 MOFA 对缺失挺稳健。
  • 样本对齐：三个组学的样本 ID 必须一一对应。我遇到过有人故意打乱甲基化行顺序，直接用 .values 拼接导致对齐错误。切记保存时用同样的索引顺序。

小技巧：保存模拟数据时，可以顺便把真实的因子矩阵（latent）和载荷（loadings）一起存下来。后面做因子解读时，直接跟真实因子对比，一眼就能看出模型学得准不准。

慎重：不要用过于简单的噪声模型。比如只用加法高斯噪声，MOFA 很容易就完美恢复因子结构。现实中数据噪声是异质的（不同组学、不同样本差异大），建议在蛋白质组和甲基化里加入 10% 的极端值（outlier）模拟技术假象。

好了，现在你手里已经有了一个结构清晰、噪声可控的多组学模拟数据集。后面几章我们会拿它来调 MOFA 的参数、评估因子分解效果。这段代码在我自己的实验里跑了不下五十遍，每次改一点参数就能看到模型行为的变化，值不值得你花十分钟亲手敲一遍？

← 上一章 📖 返回目录 下一章 →