📚 课程目录第3章 SPARK算法原理（上）

SPARK算法原理（上）：统计模型基础——高斯过程（GP）与空间相关性的数学表述

聊到空间转录组的数据建模，有一个工具绕不开——高斯过程。说实话，我第一次接触GP的时候也觉得它有点"玄乎"。什么"函数的分布"、"无限维的多元高斯"，听着就头大。但后来在几个项目里反复用过之后，我发现它其实就是一套很自然的框架：用位置之间的相似性，去推断表达量之间的相关性。

你想想看，空间转录组数据最核心的特征是什么？每个spot不仅有基因表达量，还有一个二维坐标(x, y)。相邻的spot，组织微环境相似，基因表达也倾向于相似——这就是空间相关性。高斯过程，恰恰是刻画这种"位置相近→表达相似"关系最优雅的数学模型。

  核心直觉： 高斯过程 ≈ 多元高斯分布在无限维空间的推广。它不对"函数长什么样"做硬性假设，而是通过核函数来定义"多近的点算相似、相似到什么程度"。

1. 为什么是高斯过程

空间转录组数据有一个麻烦：spot之间的空间依赖是非线性的，而且没有固定的模式。有的基因表达呈现斑块状，有的呈现梯度变化，还有的只在特定区域高表达。传统的线性模型很难同时处理这么多形状。

GP的优势在于——它是一个非参数模型。不预设函数形式，而是让数据自己"说话"。建模时，我们只需要定义两个东西：

均值函数：表达量的基线水平，通常设为常数或零
核函数（协方差函数）：定义任意两个spot之间的表达相关性如何随距离变化

我在项目中遇到过一种情况：用普通线性模型拟合空间数据，残差总是有明显的空间结构，说明模型没抓住空间依赖。换成GP之后，残差白噪化了，拟合效果明显提升。这就是GP的"用武之地"——它专门处理数据点之间的依赖结构。

2. 高斯过程的数学定义

数学上，高斯过程定义为一个随机过程，其中任意有限个点的联合分布都是多元正态分布。写成公式就是：

f(x) ~ GP( m(x), k(x, x') )

其中：
m(x) = E[f(x)]           — 均值函数
k(x, x') = Cov(f(x), f(x'))  — 核函数，也叫协方差函数

说白了，GP告诉我们的就是：给定任意一组位置x₁, x₂, ..., xₙ，对应的函数值f(x₁), f(x₂), ..., f(xₙ)服从一个n维正态分布。这个分布的均值由m(x)决定，协方差结构由k(x, x')决定。

为什么会这样？因为GP的本质就是在函数空间上放了一个先验。这个先验的"形状"由核函数控制。如果我们认为真实的空间表达模式是平滑的，那就选一个平滑的核函数；如果认为可能有突变，那就选一个允许跳跃的核函数。

我的个人习惯： 在SPARK建模中，我会先把均值函数设为一个常数（用所有spot的平均表达量）。然后通过核函数去捕捉空间变异。这样分解清晰，调参也方便。

3. 核函数：空间相关性的"编码器"

核函数是整个GP的灵魂。它编码了我们对空间相关性的先验知识。不同的核函数对应不同的空间模式假设。

下表总结了SPARK中常用的几种核函数：

核函数	公式	空间假设	适用场景
平方指数核 (SE)	k(d) = σ² exp(-d² / 2l²)	无限平滑，各向同性	缓变的表达梯度
指数核	k(d) = σ² exp(-d / l)	连续但不平滑	边缘锐利的斑块
马特恩核 (Matérn)	k(d) = σ² (2√ν d/l)ᵛ Kᵥ(2√ν d/l) / Γ(ν)	可调节平滑度	真实组织数据（推荐）
周期核	k(d) = σ² exp(-2 sin²(πd/p)/l²)	周期性重复模式	分层组织或有周期性结构

这里d是两个spot之间的欧氏距离，l是长度尺度参数，控制相关性随距离衰减的速度。l越大，相关性随距离下降得越慢，空间模式越"大尺度"。

我曾经踩过的坑： 有一回处理小脑数据，我默认用了平方指数核，结果模型死活不收敛。后来发现小脑皮层有非常锐利的层间边界，平滑性假设根本不成立。换成马特恩核（ν=1.5）之后，问题迎刃而解。核函数的选择直接影响模型的拟合能力，别偷懒、得试。

4. 从核函数到协方差矩阵

有了核函数，我们就可以构建协方差矩阵Σ。假设我们有n个spot，位置分别是s₁, s₂, ..., sₙ。那么Σ就是一个n×n的矩阵，其中第i行第j列的元素就是k(sᵢ, sⱼ)。

这个矩阵捕捉了所有spot两两之间的空间相关性。对角线上是k(sᵢ, sᵢ) = σ²，表示每个spot自身的方差；非对角线上则是spot i和spot j之间的协方差。

构建过程很简单，但有几个关键细节：

对称性： Σ一定是对称正定矩阵。核函数必须保证这一点
数值稳定性： 当n很大（比如上万个spot）时，Σ可能接近奇异。通常会在对角线上加一个小噪声项（nugget），模型变为 GP + 噪声
计算复杂度： n×n矩阵的求逆是O(n³)。对于空间转录组数据（n可达数千到数万），直接做GP推断是不现实的。SPARK用了加权核函数技巧来绕过这个瓶颈

  SPARK的创新之一： 它不直接对整个协方差矩阵求逆，而是将GP模型转化为一个加权广义线性模型。通过核矩阵的谱分解和随机效应逼近，把计算复杂度降到了O(n)。这也是SPARK能处理大规模空间转录组数据的根本原因。

5. SPARK中的GP建模流程

在SPARK中，GP被用作零假设的生成模型。具体来说，对于每一个基因，SPARK都做以下事情：

标准化表达量： 对基因的表达计数做log变换和尺度化
构建核矩阵： 基于所有spot的坐标，用选定的核函数构建协方差矩阵
拟合GP模型： 用残差极大似然（REML）估计核函数的参数（l, σ²等）
计算得分统计量： 构造一个检验统计量，衡量"加入空间相关性能多大程度上改善拟合"
置换检验/P值校正： 通过置换空间标签来得到经验P值，最后做多重假设检验校正

核心的R代码逻辑大致是这个样子：

# 构建核矩阵（以马特恩核为例）
library(SPARK)

# 输入: locations是n×2的坐标矩阵
#       counts是n×1的表达向量
#       kernel_type = "matern"

# 内部核心步骤（简化版）
dist_matrix <- as.matrix(dist(locations))
sigma <- matern_kernel(dist_matrix, 
                       l = length_scale, 
                       nu = 1.5, 
                       sigma_sq = 1.0)

# 拟合GP模型（REML估计）
fit <- fit_gp(counts, sigma, method = "REML")

# 提取空间变异参数
spatial_variance <- fit$sigma_sq
noise_variance <- fit$tau_sq

# 检验空间变异是否显著
statistic <- compute_score(fit)

当然，SPARK的实际代码要复杂得多。它用了一组加权核函数来捕捉不同尺度的空间模式，而不是只用一个核。这个我们下一章再详细展开。

一个实用的调试建议： 如果你在自己的数据上跑SPARK，发现某些基因的P值全都不显著，先检查一下核函数的长度尺度参数l是否合适。如果l设得太小，相关性只限于最近的几个spot，全局模式检测不到；如果设得太大，所有spot都"强相关"，又区分不出局部变异。我一般会尝试l取spot之间平均距离的0.5倍、1倍、2倍，看看哪个结果更合理。

好，这一章的核心内容就是这些。GP提供了一个优雅的框架来刻画空间相关性——从核函数的定义到协方差矩阵的构建，再到统计检验的实现。理解了这个基础，下一步再看SPARK中加权核函数的巧妙设计，就会顺畅得多。

用一句话总结这一章的核心：高斯过程通过核函数将空间位置映射为相关性结构，再通过协方差矩阵将这种结构嵌入统计模型，最终为SPARK检验空间可变基因提供了坚实的数学基础。你把这个逻辑捋顺了，后面的内容就都是"术"层面的细节了。