第3节：SVA（Surrogate Variable Analysis）原理

SVA的数学思想：为什么需要潜在变量？

我第一次接触到SVA这个名词时，说实话有点懵。什么叫“可替代变量”？要理解这个，咱们得先想明白一个问题：你的数据里，到底藏着哪些你不知道的东西？

你跑了一个RNA-seq实验，样本分两组：处理组和对照组。你仔细控制了实验条件，甚至连移液器都校准了三遍。结果跑差异表达分析，出来的p-value分布歪七扭八，假阳性一大堆。为什么？

其实，真实世界的数据里，永远有一些你测量不到的因素在作祟。比如：

• 批次效应：今天做的实验和明天做的，试剂可能微调过
• 操作人员变化：不同人提取RNA，效率不一样
• 样本处理时间：有些样本放了半小时，有些放了俩小时
• 芯片差异：哪怕同一批芯片，不同通道灵敏度也不同

这些因素，你很难全部记录在表格里。但它们确实在影响你的数据。SVA就是用来处理这个问题的——它帮你把这些“看不见的干扰”给找出来。

SVA的核心步骤——一步一步拆解

SVA的全称是Surrogate Variable Analysis，中文叫替代变量分析。它的核心假设很简单：你的数据矩阵Y，可以分解成三部分：

Y = Xβ +   (已知的生物学变量，比如分组)
      + Bh    (已知的批次效应，如果有的话)
      + ΓS    (未知的异质性，用替代变量S来捕获)
      + R     (随机噪声)

这里面，S就是我们要估计的“替代变量”（Surrogate Variables），Γ是它们的系数矩阵。SVA的任务就是从Y中把S和Γ给拆出来。

具体怎么做？我拆成几个步骤讲。

1. 第一步：用已知变量做回归，得到残差

   先把你的已知变量（比如分组、已知的批次）放进模型，拟合一下。然后计算残差矩阵R。这个R里，主要就剩未知的异质性和随机噪声了。

2. 第二步：在残差里找“异常”基因

   残差矩阵R里，有些基因的波动特别大。这些基因很可能跟未知的异质性有关。SVA会用一种统计方法——我简化说——去识别那些表达模式明显偏离预期的基因。这些基因被叫做“异常表达基因”或“潜在相关基因”。

3. 第三步：用这些基因构造替代变量

   有了这些异常基因，就可以用奇异值分解（SVD）之类的办法，从它们的表达谱里提取出几个主成分。这些主成分，就是替代变量S。每个样本会对应一组S值，代表这个样本在“未知因素”上的强度。

4. 第四步：把S加入模型重新分析

   最后，把估计出来的S作为协变量，加回你的差异表达模型里。这样，那些乱七八糟的未知因素就被“控制住”了，真实信号就能露出来。

一张图看懂SVA的数学结构

下面这张图，是我根据实际项目经验画的。它把SVA的数据流和模型结构串起来了：

SVA的三个关键假设——你得心里有数

任何数学模型都有假设。SVA也不例外。我个人习惯在开班前就把这些假设给学生讲清楚，免得后面跑偏。

假设	说明	违反后果
大部分基因不受未知因素影响	SVA依赖“大多数基因是稳定的”这个前提来估计残差和基线	如果几乎所有基因都被干扰，SVA会失效，估计出来的S不可靠
未知因素与已知因素不高度相关	如果S和分组变量强相关，模型无法区分哪个是真实信号	会导致信号混淆：你有可能把真正的处理效应当作批次给去掉了
样本量足够大	SVA需要一定数量的样本来稳定估计替代变量	小样本（比如每组3个）下，SVA的表现不稳定，容易过拟合

如何判断是否需要SVA？

你可能会问：我什么时候需要上SVA？嗯，这里有几个信号：

• 你的p-value分布明显偏向0或1，呈“U型”或“簸箕型”
• PCA图里，样本没有按照你的分组聚在一起，反而按别的方式聚
• 你发现某些基因的表达量跟样本处理批次高度相关，但你没法解释为什么

如果出现这些情况，SVA值得一试。但记住，它不是第一步——它应该在你已经检查过已知批次、做过基础校正之后，再用来兜底。

R语言中的SVA实战接口

最后，给一个快速上手指南。在R里，SVA的主要实现是Bioconductor的sva包。核心函数就两个：

# 安装（如果还没装）
BiocManager::install("sva")

# 加载
library(sva)

# 准备数据：表达矩阵 (基因x样本)
# 准备变量：已知分组和批次
mod <- model.matrix(~ group + batch, data = pheno)
mod0 <- model.matrix(~ 1, data = pheno)  # 只包含截距

# 运行SVA，自动估计替代变量
sv_obj <- sva(dat = expr_matrix, mod = mod, mod0 = mod0)

# 提取替代变量
sv <- sv_obj$sv

# 加入模型
mod_sva <- cbind(mod, sv)

# 然后就可以用mod_sva来做差异分析了