Data Integration for Omics

"整合"(integration)在组学分析中有两个不同但同样核心的问题,本系列用两个约 35 分钟的动手教程分别覆盖:

本页面为交互版:三个 + 三个交互实验的所有参数组合已预先计算完毕,拖动滑块即时查看,无需运行任何代码。完整讲解与代码见两个 notebook。

I单模态整合(Single-Modality / Horizontal Integration)

0.1 生物学背景:中心法则与细胞类型

中心法则(central dogma):遗传信息从 DNA → RNA → 蛋白质流动。基因表达(gene expression)指基因被转录成 mRNA 的活跃程度——一个细胞"是谁",很大程度上由它正在表达哪些基因决定。生物体由成百上千种细胞类型(cell type)组成:即使基因组完全相同,不同细胞类型表达的基因组合也不同。区分和刻画细胞类型,是理解发育、疾病(如肿瘤中浸润了哪些免疫细胞)和药物响应的基础。本教程的真实案例数据是 PBMC(外周血单个核细胞)——血液中的各类免疫细胞,包括 T 细胞、B 细胞、NK 细胞和单核细胞。

图:中心法则。来源:yourgenome, Genome Research Limited(CC BY 4.0)

0.2 scRNA-seq 实验:从细胞到计数矩阵

单细胞 RNA 测序让我们能够逐个细胞地测量基因表达。以主流的液滴微流控技术(10x Genomics)为例:(1) 每个液滴包裹一个细胞和一个带条形码的微珠;(2) mRNA 被捕获并逆转录,每个分子带上细胞条形码和 UMI(去除 PCR 重复);(3) 测序后按条形码归还细胞、按基因计数 UMI,得到计数矩阵。因此矩阵中那个 "5" 的物理含义是:这个细胞里检测到了 5 个来自该基因的 mRNA 分子——它是离散的计数,不是连续的浓度。

图:scRNA-seq 从样本到计数矩阵的完整流程。来源:Wang et al., Int J Mol Sci 2023(CC BY 4.0)

0.3 数据的三个性质

在学习任何整合方法之前,先观察一份真实 scRNA-seq 数据(5,197 个 PBMC × 1,000 个高变基因,来自 Part II 的 CITE-seq 数据集):

• 稀疏性:矩阵中 94.5% 的元素是 0——测序深度永远不够用,许多真实表达的基因没被测到(dropout);
• 计数性:数值是离散的 UMI 计数,不同基因的方差跨越好几个数量级——所以需要标准化;
• 文库大小差异:每个细胞的总 UMI 数从 501 到 53,577——这是技术噪声而非生物学差异,所以分析前必须先归一化。

0.4 批次效应从哪来?

批次(batch):在同一轮实验流程中一起处理的一组样本。批次效应的来源几乎无法穷尽:不同日期上机、不同测序仪 run、不同批次的试剂(反转录酶、微珠批号)、不同操作人员、不同技术平台(如 10x v2 vs v3)。这不是假想敌:Seurat 官方教程的基准数据 pbmcsca 汇集了 9 个批次、7 种测序平台的 PBMC——不做整合时,细胞几乎完全按批次/平台聚类,而不是按细胞类型。

0.5 单模态整合的定义

单模态整合要解决的问题:同一种组学(如 scRNA-seq)在不同批次、不同样本、不同平台上测了多次,如何合并成一个一致的数据集?核心障碍是批次效应(batch effect):同一类型的细胞,仅因为在不同批次中测量,表达谱就出现系统性差异。数据模型:

其中 x_ij 是细胞 i 在批次 j 中的观测向量,z_i 是其真实生物学状态,b_j 是批次 j 的系统性偏移。整合的目标:去掉 b_j,保留 z_i。难点在于 z 和 b 都不可直接观测——校正太弱,批次仍然分离;校正太强(过度校正),批次间真实的生物学差异也被一并抹掉。

0.6 如何定量评估:一把尺子不够,要两把

• 生物学保守性:整合后同类细胞是否仍聚在一起。指标:silhouette(按细胞类型),越大越好。对每个细胞 i,记 a(i) 为它到同类型其他细胞的平均距离,b(i) 为它到最近的异类型细胞的平均距离:

全体取平均:s ≈ 1 表示类内紧密、类间分离;s ≈ 0 表示完全混杂。

• 批次混合度:整合后不同批次的细胞是否混在一起。指标:batch mixing rate——每个细胞的 k 个最近邻中来自另一批次的平均比例:

两批次等量时,0.5 = 完全混合(理想),0 = 完全分离(失败)。

1基线:合并 + PCA —— 以及它为什么失败

最朴素的做法:把两个批次的矩阵直接摞起来做 PCA。PCA 寻找方差最大的投影方向,但它不知道哪些方差来自生物学、哪些来自批次——批次效应足够强时,方差最大的方向就是批次方向。数学上,PCA 最小化重建误差,且有显式解(对中心化数据的协方差矩阵做特征分解):

特征值 λ_i 即第 i 个主成分解释的方差。模拟场景:3 种细胞类型、2 个批次;50 个特征中只有 10 个携带细胞类型信息;批次 2 带有强度可控的整体偏移。

2逐批次校正(ComBat 思想)

最直接的想法:既然批次效应是"批次 2 整体偏了一块",那就把每个批次各自中心化,把偏移减掉:

其中 x̄_·j 是批次 j 的均值向量,x̄_·· 是全体均值。完整的 ComBat 还会校正尺度,对每个基因单独建模:

同时用经验贝叶斯在基因间共享信息、收缩小样本下的估计。对平移型批次效应,中心化几乎完美;但真实批次差异可能是旋转/缩放——给模拟数据加一个旋转成分:

3从 autoencoder 到 VAE / scVI

自编码器用 encoder 把 x 压缩成低维 z,再用 decoder 重建 x(min‖x − x̂‖²)。一个经典结论:当 encoder 和 decoder 都是线性映射时,自编码器的最优解与 PCA 张成同一个子空间——因此普通(线性)AE 面对批次效应会和 PCA 一样失败。

scVI 等方法的核心思想一句话:把批次标签作为 decoder 的额外输入,让 latent z 不必再编码批次信息:

既然 decoder 已经"知道"批次,批次造成的差异可以由 decoder 的批次项单独解释,z 就可以专心编码生物学信息。下面的交互展示用约 20 行 numpy 实现的极简线性版本(代码见下方折叠块)与普通线性 AE 的对比:

def cond_ae(X, batch, latent=2, epochs=4000, lr=0.001, beta=0.01, seed=0):
    """x̂ = z·Wd + onehot(batch)·Bd;损失 = 重建误差 + beta·‖z‖²"""
    r = np.random.default_rng(seed)
    B = np.eye(2)[batch]
    We = r.normal(0, 0.05, (X.shape[1], latent))
    Wd = r.normal(0, 0.05, (latent, X.shape[1]))
    Bd = np.zeros((2, X.shape[1]))
    Bd[1] = X[batch == 1].mean(0) - X[batch == 0].mean(0)
    for _ in range(epochs):
        Z = X @ We
        E = Z @ Wd + B @ Bd - X
        We -= lr * (X.T @ (E @ Wd.T) / len(X) + beta * X.T @ Z / len(X))
        Wd -= lr * (Z.T @ E) / len(X)
        Bd -= lr * (B.T @ E) / len(X)
    return X @ We

完整的 scVI 在极简实现之上再进一步:(1) 概率化——encoder 输出高斯分布 q_φ(z|x) = N(μ_φ(x), diag σ_φ²(x)),先验取 p(z) = N(0, I),训练目标为 ELBO:

第一项是重建质量;第二项把 latent 分布向先验收缩(两个高斯之间的 KL 有闭式解)。训练时通过重参数化技巧 z = μ_φ(x) + σ_φ(x) ⊙ ε, ε ~ N(0, I) 使梯度可以穿过采样步骤回传。(2) 批次协变量——decoder 变为 p_θ(x|z, b),正是极简实现中批次项的角色。此外,普通线性 AE 的目标 min ‖X − X·W_e·W_d‖² 与 PCA 完全一致(W_eW_d 遍历所有秩 ≤ k 的矩阵,Baldi & Hornik 1989)——这正是"线性 AE ≈ PCA"的代数形式。

4Part I 小结与常见误区

合并 + PCA(不知道什么是批次) → 逐批次校正(假设批次 = 平移) → 条件 AE / scVI(批次 = decoder 协变量)

• 批次效应:整合的目标是去 b_j 保 z_i,二者不可兼得时即为过度校正;
• 双面评估:生物学保守性(silhouette)× 批次混合度(mixing rate),任何单一指标都不足以判断整合质量;
• 逐批次校正只适用于平移型批次效应,对旋转等方向性差异无能为力;
• 条件 AE / scVI 把批次作为 decoder 协变量;其线性情形退化为逐批次校正,神经网络化后才能处理非线性批次效应。

• ❌ 合并后不校正就直接聚类、注释(聚类可能完全由批次驱动);
• ❌ 只用单一指标评估整合——mixing 差时 silhouette 也可以很高,反之亦然;
• ❌ 默认"校正越强越好"——过度校正会抹掉真实的组间差异(如疾病 vs 对照);
• ❌ 不看原始数据就校正——先画未校正的图,确认批次效应确实存在;
• ❌ 把深度模型当作默认选项——小数据、简单批次结构下,线性方法往往足够且更稳健。

延伸阅读

Part II · 多模态整合:Machine Learning for Omics Data Integration

本教程围绕三个核心问题展开:

1. 为什么需要整合多个组学(omics)数据,而不是只分析其中一种?
2. 最朴素的整合策略(直接拼接 + PCA)在什么情况下会失效?
3. CCA 等多视图(multi-view)方法为什么更稳健?在真实数据上表现如何?

学习目标(Learning Objectives)

• 用"样本 × 特征矩阵"的语言描述组学数据,并解释 view(视图)与 latent(潜在状态)两个概念;
• 说明多视图整合为何能产生"1 + 1 > 2"的效果;
• 独立实现两条整合流程:concat + PCA、per-view PCA + CCA;
• 使用 silhouette 系数与 kNN 准确率对整合结果进行定量评估,而非仅凭可视化图主观判断。

受众与先修要求

本教程不涉及的内容(为避免预期错位,特此说明):深度学习方法(scVI 等)的实操、跨批次整合(horizontal integration)的细节、UMAP/t-SNE 可视化、以及除 CITE-seq 之外的数据类型。相关内容见第 4 章的延伸阅读。

时间安排(约 35 分钟)

环境要求:仅依赖 numpy / matplotlib / scikit-learn(pip install -r requirements.txt),全程 CPU。本页面为交互版:所有参数组合已预先计算,拖动滑块即时查看,无需运行代码;动手代码练习见 exercise_pbmc10k.ipynb。若本地环境配置遇到问题,可使用 Google Colab:上传 .ipynb 与 .npz 数据文件即可零修改运行。

0.1 组学数据是一个矩阵

一份组学数据本质上是一个矩阵:每一行是一个样本(本教程中为单细胞),每一列是一个特征(例如某个基因的表达量)。矩阵中的数值表示"该样本中该特征的测量强度"。

            gene_1  gene_2  ...  gene_20000
cell_1      [  0      5            2      ]
cell_2      [  3      0            0      ]      ← 这是一个 view(RNA 转录组)
 ...        [       ...                   ]
cell_300    [  1      0            7      ]

0.2 什么是 view(视图)?

对同一批细胞,可以使用不同的分子测量技术各测一次,得到多个矩阵。每一个矩阵称为一个 view(视图)。例如 CITE-seq 技术可以在测转录组的同时,用抗体测量细胞表面蛋白的丰度,得到第二个矩阵(列数通常只有几十):

            CD3   CD4   CD8 ...  (共 29 种蛋白)
cell_1      [ 12    0     5      ]
 ...        [        ...         ]              ← 这是另一个 view(表面蛋白)
cell_300    [  0    9     0      ]

0.3 整合的目标:恢复潜在状态(latent)

两个 view 描述的是同一批细胞,因此它们背后存在一个共同的、不可直接观测的"细胞状态"(cell state),记作潜在变量 z。每个 view 都可以看作 z 经过某种变换、再叠加噪声后的观测:

数据整合(data integration)的目标:从多个带噪声的 view 中,把共享的 z 尽可能准确地恢复出来。在本教程中,z 可以理解为"每个细胞属于哪种细胞类型"这一信息的低维表示。

第 1、2 章使用模拟数据:z 由我们人为生成,因此拥有标准答案,可以直接检验每种方法恢复 z 的准确程度;第 3 章再将同样的流程应用于真实数据。

两个贯穿本教程的定量指标

• silhouette 系数:对每个细胞 i,记 a(i) 为它到同类型其他细胞的平均距离,b(i) 为它到最近的异类型细胞的平均距离:
  全体取平均:s ≈ 1 表示类内紧密、类间分离,s ≈ 0 表示完全混杂。
• kNN 准确率:用每个细胞最近的 10 个邻居投票预测其类型,以 5 折交叉验证评估(数据分 5 份,轮流 4 份训练、1 份测试取平均),避免用同一批数据既训练又测试造成的高估。
• 标准化:拼接或降维之前,每个特征(列)缩放到均值 0、方差 1,即 x' = (x − μ) / σ,消除量纲差异。

1模拟实验 A:信息互补 —— 为什么需要整合?

首先构造一个经过简化、但在真实生物学中普遍存在的场景:每个 view 只携带一部分信息。

• 在真实世界中:某些细胞类型(如 B 细胞与单核细胞)在转录组上差异明显;而另一些细胞类型(如 CD4 与 CD8 T 细胞)转录组几乎相同,必须依靠表面蛋白才能区分。没有任何单一组学能够区分所有细胞类型。
• 在模拟数据中:4 种细胞排列在 2×2 网格的四个角上。RNA view 只能"观测到"横坐标的信息,protein view 只能"观测到"纵坐标的信息。

结果解读(noise = 0.5 的默认档):每个 view 单独只能将 4 类细胞区分为 2 类(silhouette ≈ 0,kNN 准确率 ≈ 0.5,相当于随机猜两类的水平);将两个 view 拼接后,4 类细胞完全分开(silhouette ≈ 0.74,kNN 准确率 = 1.0)。

整合流程中有一个细节值得注意:拼接之前,两个 view 都必须先做标准化(每列缩放到均值 0、方差 1)。真实数据中不同 view 的量纲差异极大(RNA 原始计数可达数千,蛋白信号通常仅为个位数);若不进行标准化,量纲较大的 view 将在拼接后的矩阵中占据主导地位,使整合退化为单 view 分析。这是多组学整合中最常见的预处理错误。

2模拟实验 B:CCA —— 提取视图间"共享"的信号

上一章中的两个 view 都只含有用信号和随机噪声。现在引入一种更接近真实数据的情况:每个 view 除了共享信号 z 之外,还含有仅属于该 view 的"私有变异"(private variation)。

私有变异在真实数据中非常常见,最典型的例子是批次效应(batch effect):同一批细胞分两次测序,两次结果之间的系统差异与生物学无关,只存在于各自的测量之中。数学上可以写成:

X₁ = z·W₁ 【共享】 + a₁·s₁·V₁ 【view 1 私有】 + ε₁

X₂ = z·W₂ 【共享】 + a₂·s₂·V₂ 【view 2 私有】 + ε₂

当私有变异较强时,直接拼接 + PCA 会失效:PCA 只认"方差最大的方向",它并不知道哪些方向是两个 view 共有的,因此会把私有变异当作数据结构的主要成分。

CCA(Canonical Correlation Analysis,典型相关分析)的思路可以概括为一句话:分别为两个 view 各寻找一个投影方向 a、b,使两个 view 投影后的相关性最大。设两个 view 的协方差矩阵为 Σ₁₁、Σ₂₂,跨 view 协方差为 Σ₁₂,则:

该问题有显式解——对矩阵 M 做奇异值分解:

奇异值 λ_i 即第 i 对投影的相关性。取前 k 对方向即得 k 维整合表示。

由于私有变异在另一个 view 中不存在对应成分,无法产生跨 view 的相关,因此在优化过程中被自动排除;只有共享信号 z 能够同时驱动两个 view 的投影。这就是 CCA 能够"找回" z 的原因。

3真实数据实战:CITE-seq

下面将前两章学到的流程应用于真实数据。

3.1 生物学背景:PBMC 与免疫细胞类型

PBMC(peripheral blood mononuclear cells,外周血单个核细胞)是从外周血中分离的免疫细胞混合物,是免疫学研究最常用的样本之一。其中的主要细胞类型:

"CD" 是 cluster of differentiation(分化簇) 的缩写——一套国际统一的细胞表面蛋白命名系统,免疫学家正是靠这些表面蛋白的组合来定义和分选细胞类型。

图:PBMC 主要细胞类型与定义性 marker。细胞插图来源:Servier Medical Art(CC BY 4.0),本教程组合加工

3.2 CITE-seq:一次实验,两个 view

实验原理:CITE-seq 在常规液滴单细胞测序的基础上,加入一组抗体-寡核苷酸偶联物(如 TotalSeq):每种抗体识别一种表面蛋白,并带有一段可测序的条形码。同一个液滴内,细胞的 mRNA 和抗体条形码被一起捕获测序——因此两个 view 天然按细胞对齐,这是多模态整合的前提。

图:CITE-seq 的抗体-寡核苷酸偶联物结构(A)与单细胞工作流程(B)。来源:Song et al., Proteomics 2025(CC BY 4.0)

为什么要测两个 view?中心法则并不是 1:1 的管道:mRNA 丰度与蛋白丰度之间隔着翻译调控、蛋白降解等环节,二者经常不一致。一个著名的例子是 CD4/CD8 T 细胞的区分:转录组几乎相同,但表面蛋白差异明确——只靠 RNA 几乎无法分开(3.4 节会用数据验证)。这就是多模态测量的生物学动机。

高变基因(highly variable genes):细胞间表达差异最大的基因。两万个基因里大部分携带的细胞类型信息很少,常规做法是只保留信息最集中的前 1000–3000 个。每个细胞的总 UMI 数从 501 到 53,577(技术噪声,不是生物学差异):

3.3 参考标签从何而来?

免疫学中有一套成熟的表面蛋白标记(marker)规则,例如:CD19 阳性 → B 细胞;CD3 阳性且 CD4 阳性 → CD4 T 细胞。我们据此为每个细胞赋予一个参考类型。

需要说明:标签由蛋白 view 派生而来,因此它是一个"参考答案"而非绝对金标准——在评估时,蛋白 view 会天然占有一定的优势,解读结果时请记住这一点。

数据文件 citeseq_pbmc.npz 已完成预处理:RNA 经过文库大小标准化 + log(1+x);ADT 经过 CLR(centered log-ratio)变换——每个细胞内,各蛋白的 log(1+x) 值减去该细胞的均值,以消除蛋白计数中的组成偏差:

从 10x 原始数据生成该文件的完整流程见 prepare_data.py。

注:真实数据分析的第一步是质量控制(QC)。该数据中存在少量极端离群细胞,会将 PCA 的主要方向拉偏;下列结果已删除落在前后 0.5% 分位之外的 98 个细胞(保留 5,099 个)。

3.4 整合流程对比

流程 1:单 view 分析与朴素拼接。下列三幅图为固定结果(图中为 1,500 个细胞的子采样,指标在全部 5,099 个细胞上计算):

如何解读?在 2 维 silhouette 上,concat + PCA 最高——整合后的整体结构与参考标签的一致性确实优于任一单 view;但 kNN 准确率上蛋白 view 仍然略高(参考标签来自蛋白 marker,蛋白 view 天然占优);此外,2 个主成分无法容纳全部结构,而真实分析从不只使用 2 维。

流程 2:每个 view 先 PCA 降维,再运行 CCA(实际应用中的标准做法)。拖动滑块,调整进入 CCA 前各 view 保留的主成分数:

一个更公平的比较方式是:使用更多组分(RNA 取 20 个主成分、ADT 取 10 个主成分)构建整合表示,然后考察每一种细胞类型各自的预测准确率。整合的价值将在这一比较中真正显现:

这才是真实数据中的信息互补(注意高亮的 CD8 T 一行):

• CD8 T 细胞:RNA view 单独的预测准确率仅约 0.5(在单细胞转录组中区分 CD4 与 CD8 T 细胞是公认的难题),而蛋白 view 可以轻松达到约 0.9;拼接之后,CD8 T 被恢复至约 0.87——并且无需事先知道"这一任务必须依赖蛋白 view"。
• 由此可见,整合的意义不在于全面超越单 view,而在于:当你无法预知哪个 view 对哪类细胞更重要时,整合表示在任何一个类别上都不会失效。

那么 CCA 表现如何?在这份数据上,CCA 反而逊于 concat——注意 CD8 T 一列:CCA 与 RNA only 一样下降到约 0.5。其机制与模拟实验 B 完全对应:CD8 的判别信号几乎只存在于蛋白 view 中(RNA 难以区分 CD4/CD8),对 CCA 而言这正是需要剔除的"私有变异"。CCA 只保留 view 间共享的信息;当某一类别的关键信息仅存在于单个 view 时,CCA 会将该信息一并丢弃。

因此,在两个 view 高度一致、私有变异较弱的数据中(如这份 PBMC),朴素的 concat 已足够好;CCA(以及 Seurat、scVI 等方法)真正的应用场景是跨批次、跨平台整合——那里的"私有变异"(批次效应)足够强,共享信号的提取才有必要。没有普适最优的方法,只有与数据结构相匹配的方法。

4总结与延伸阅读

本教程的主线

数据整合 = 从多个 view 中恢复共享的细胞状态 z。方法按对数据结构的假设由简到繁排列:

concat + PCA(基线) → CCA(利用配对信息,剔除 view 私有变异)

• 信息互补:整合的价值在于 view 间的互补,而非数据量的简单叠加;
• 预处理:拼接之前每个 view 各自标准化;高维 view 先经 PCA 降维,再进入整合流程;
• 共享 vs 私有:区分共享信号与私有变异是多视图方法的核心思想。真实数据上方法之间的差距通常远小于模拟数据,只有通过定量评估(而非肉眼观察图形)才能确认整合是否带来了真实的提升;
• 整合并非越多越好:过强的整合可能抹掉你真正关心的生物学差异(过度校正,overcorrection)。因此任何整合结果都必须与"未整合"的基线进行对比,确认收益是真实的。

延伸阅读(本教程未涉及,但值得了解的方法)

本教程的 CITE-seq 数据仅包含两个 view、约 5000 个细胞。当数据规模达到数万细胞、或细胞状态呈连续轨迹(如分化、细胞周期)时,线性方法将不再够用,此时需要非线性方法(核技巧、UMAP、autoencoder、scVI 等)。

常见误区清单

• ❌ 将"拼接"等同于整合本身,不考虑各 view 之间的尺度与噪声差异;
• ❌ 拼接前忘记标准化(量纲大的 view 会主导全部结果);
• ❌ 特征数远大于样本数时直接运行 CCA,导致过拟合(应先 PCA 降维);
• ❌ 仅凭"图形美观"下结论,而不计算任何定量指标(silhouette、kNN 准确率等);
• ❌ 默认"越复杂的方法越好",而不先分析数据中实际存在的变异结构(共享?私有?非线性?);
• ❌ 将 ADT(蛋白)计数当作 RNA 处理:ADT 数据存在环境抗体带来的背景噪声,且每个细胞的 ADT 总计数本身与细胞类型相关,不能直接套用 RNA 的标准化流程。最低要求是 CLR 变换(本教程的数据已完成),更标准的做法是 DSB;同时需留意同型对照(IgG control)抗体的质控作用。

配套文件(同目录)

• exercise_pbmc10k.ipynb — 学生练习(在新数据集 PBMC 10k 上独立完成整合流程)
• citeseq_pbmc.npz — 教程数据:PBMC 5k CITE-seq(1.4 MB)
• citeseq_pbmc10k.npz — 练习数据:PBMC 10k CITE-seq(1.7 MB)
• prepare_data.py — 从 10x 原始数据复现上述 npz 的预处理脚本
• requirements.txt — 仅 numpy / matplotlib / scikit-learn / ipywidgets