SPARK空间转录组SVG分析：从统计原理到R包实操全解

空间转录组这两年有多火，不用我多说。但拿到一张带空间坐标的表达矩阵，除了常规的降维聚类、细胞类型注释之外，还有一个绕不开的核心分析——识别空间可变基因（spatially variable genes，SVG）。这正是SPARK这个R包的主场。如果你正在做10x Visium、Slide-seq、Stereo-seq这类空间转录组数据，想知道“哪些基因的表达量和空间位置显著相关”，那SPARK就是绕不开的工具之一。这篇文章我会从原理到实操，完整拆解SPARK的用法、参数逻辑和我在实际项目中踩过的坑，希望给你一份能直接照着跑的参考。

先说清楚一件事：这里的SPARK是一个R包，全称是Spatial Pattern Recognition via Kernels，不是大数据圈那个Apache Spark集群计算框架。两者没有任何关系，但网上搜索时经常会被混在一起，光这一点就够劝退不少新手。所以这篇文章里所有出现SPARK的地方，默认都在说R语言里的这个空间转录组分析包。

我自己的习惯是，接触一个新工具先不急着跑代码，而是先搞清楚它到底在解决什么问题、和同类方法有什么本质区别，这样后面调参、排错才不会抓瞎。这篇博文就按照这个思路来：先讲清楚SPARK的定位和原理，再给一份从安装到可视化的完整实操流程，最后把我在真实数据上遇到的坑集中整理出来。

1. 先搞清楚SPARK在解决什么问题

1.1 高变基因分析在空间转录组里的新含义

做单细胞转录组的朋友对“高变基因”（highly variable genes，HVG）这个概念应该很熟。常规流程里，我们会用Seurat的FindVariableFeatures跑一遍，筛选出那些在不同细胞之间表达量波动最明显的基因，然后用这些基因做下游的PCA、聚类。这个思路背后的逻辑是：表达差异大的基因，通常携带更多生物学信息。

但空间转录组数据多了一个维度——空间位置。这个时候，“高变”的定义就不够用了。一个基因可能在全局表达量上很平稳，方差很小，传统方法根本筛不出它；可如果你把表达量映射到组织切片的空间坐标上，却发现它在某个区域高表达、在另一个区域几乎不表达，形成非常清晰的空间分布模式。这样的基因，从统计角度看方差不大，从生物学角度看却非常关键——它可能标记了一个重要的组织结构或功能区域。

SPARK要识别的就是这一类基因，术语叫空间可变基因（SVG）或空间高变基因。它跟普通高变基因的核心区别在于：HVG只看表达量的离散程度，SVG看的是表达量和空间位置之间的关联性。一个基因只有在表达模式呈现出空间分布结构时，才会被SPARK判定为显著。

这也决定了分析目标不同。用Seurat找HVG，目的是为下游降维聚类准备特征；用SPARK找SVG，目的往往是探索组织区域 marker、发育梯度、疾病区域的空间分子特征，或者为后续的空间域识别、细胞互作分析提供生物学候选基因。

1.2 和其他SVG识别方法比，SPARK强在哪

目前识别SVG的方法大概可以分为几类：基于空间自相关统计量的（比如Moran's I、Geary's C）、基于高斯过程回归的（SPARK就是这一类）、基于神经网络的（如SpatialDE也有类似思路）、以及近年一些基于深度学习的工具。SPARK能在其中站稳脚跟，主要是它在统计建模上做得比较完整。

Moran's I这类方法操作简单，速度也快，但有一个天然短板：它本质上是在检验基因表达和空间邻域均值之间的相关性，对空间模式的假设比较单一，容易漏掉那些非单调、非线性的空间分布模式。SPARK的做法则是用一组核函数去拟合各种可能的空间模式，再通过统计检验判断哪个基因和空间位置的关联显著，能捕获的模式更丰富。

还有一个经常被拿来对比的工具是SpatialDE。SPARK和SpatialDE思想有相似之处，但SPARK在速度上有明显优势，尤其是在基因数比较多的时候。我自己的体感是，SPARK跑一万多个基因在几十个核心下是能接受的，而SpatialDE同样规模下等待时间会长不少。当然SPARK-X出来之后速度又提升了一个量级，这个我们后面专门说。

另外还要提一点：SPARK的输入是原始count矩阵，不需要预先标准化。这点和很多只接受标准化后表达矩阵的工具不同。因为SPARK本身是用广义线性模型带随机效应去拟合count数据的分布，把测序深度的差异作为协变量处理了。这是一个非常贴心的设计，省掉了一步数据处理，也避免了标准化步骤可能带来的偏差。

2. SPARK的核心原理，尽量用大白话讲清楚

2.1 它是怎么捕捉“空间模式”的

SPARK的核心思想，可以理解成把每个基因的表达量当成一个带空间位置信息的观测值，然后看这个观测值能不能被空间位置“解释”。

具体实现上，SPARK引入了核函数（kernel function）的概念。你可以把核函数理解为一种衡量两个spot之间“空间相似度”的标尺。如果两个spot离得近，它们之间就有更高的空间相似度；离得远，相似度就低。SPARK默认提供了几组不同的核参数，相当于用不同尺度去“扫描”空间模式——有的核适合捕捉那种精细的、小范围的局部表达结构，有的核适合捕捉那种大范围的、平滑的梯度变化。

那为什么需要好几组核呢？因为真实组织里的空间模式太复杂了。比如脑组织里，有些基因是沿着皮层深度梯度变化的，跨度很大；有些基因则只在一个很小解剖结构里强表达，范围很窄。一种固定尺度很难同时适配这两种情况。SPARK的策略就是让多种尺度的核一起上，通过统计模型自动判断哪种模式更能解释这个基因的表达变化。

这一步听起来复杂，但实际操作中你不需要操心核函数的细节。SPARK包内部把这一整套逻辑封装好了，用户要做的只是提供表达矩阵和坐标信息。

2.2 显著性检验是怎么做的

光看核函数拟合效果还不够，SPARK还需要回答一个问题：这个基因的空间模式是真的，还是随机碰巧出现的？

它的做法是构建一个统计检验。零假设是这个基因的表达量和空间位置没有任何关系，也就是表达量在空间上是随机分布的。在这个零假设之下，SPARK会计算一个统计量的理论分布，然后把实际观察到的统计量放进去，算出相应的P值。之后再用Benjamini-Hochberg等方法做多重假设检验校正，得到FDR。我们最终筛选显著基因，看的都是这个校正后的P值。

这里有一个细节值得注意：SPARK在拟合每个基因的时候，用的不是简单的线性回归，而是带随机效应的广义线性模型。随机效应的协方差结构正是由前面说的核函数定义的。这意味着模型在估计一个基因是否显著的时候，已经把空间相关性这个因素纳入进去了，而不是把一个一个spot当成完全独立的样本。这一点在统计上是更严谨的——因为空间数据天然存在“邻近spot更相似”的依赖性，如果忽略这种依赖性，很容易产生假阳性。

2.3 SPARK-X为什么能快那么多

SPARK好用是真好用，但最初版本在数据量大的时候，速度真让人着急。后来实验室推出了SPARK-X，从原理上做了一个根本性的简化。

这么说吧，原版SPARK对每个基因都要做一次带随机效应的模型拟合，这个计算量非常大。SPARK-X则绕开了这个步骤，改用一个基于“两次核变换”的思路：先用核函数把空间坐标变换成一个特征矩阵，然后用这个特征矩阵和表达量做协方差计算，最后构造一个近似卡方分布的检验统计量。整个过程中，每个基因的检验都是通过矩阵运算完成的，不需要做复杂的迭代优化，速度自然快了几十倍甚至上百倍。

SPARK-X的代价是，它对空间模式的建模精度会略逊于原版SPARK，在复杂空间结构上可能不如SPARK敏感。但大部分实际分析场景下，SPARK-X的结果已经足够可靠了。我的建议是：如果组织切片spot数不多（比如几千个），直接跑SPARK没问题；如果是Slide-seq或者Stereo-seq这种几万甚至十几万spot的数据，优先用SPARK-X，否则等待时间会让你怀疑人生。

3. 实操：从安装到出结果，照着跑就行

3.1 安装SPARK，注意这几个坑

安装SPARK本身不算复杂，但有几个环境问题很容易卡住人。第一个是R版本不要太老，建议R 4.0以上；第二个是SPARK依赖的RcppArmadillo编译，在Linux服务器上经常需要装一下系统依赖。

# 安装依赖包 install.packages(c("fields", "Rcpp", "RcppArmadillo", "RcppEigen")) install.packages("spatstat") # 从GitHub安装SPARK if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools") devtools::install_github("xzhoulab/SPARK")

我在服务器上安装的时候，遇到最多的报错是RcppArmadillo编译失败。这种情况一般不是R的问题，而是系统层面缺少GSL或者编译工具链不完整。Debian/Ubuntu系统可以用apt-get install gsl-bin libgsl-dev解决，CentOS系列则需要yum install gsl-devel。另外在编译前把Sys.setenv(GSL_CFLAGS="-I/usr/include/gsl")和Sys.setenv(GSL_LIBS="-lgsl -lgslcblas")设一下，能省掉不少麻烦。

装好之后，用library(SPARK)验证一下，如果能正常加载且没有报错提示，安装这一步就算过了。

3.2 准备输入数据：表达矩阵和坐标一个都不能少

SPARK的输入数据是两部分：一个是基因表达矩阵，行是基因，列是spot，值必须是原始count，也就是UMI计数，千万不能先做log标准化；另一个是spot的空间坐标，两列分别对应x轴和y轴。

在实际项目中，数据一般来自Seurat对象或StLearn的h5ad文件。这里给一个从Seurat对象提取数据的常见写法：

library(Seurat) # 假设你的Seurat对象叫obj sp_count <- GetAssayData(obj, assay = "Spatial", slot = "counts") sp_count <- as.matrix(sp_count) # 提取坐标信息：不同平台的位置不同，Visium一般在obj@images[[1]]@coordinates里 coordinate <- as.matrix(obj@images[[1]]@coordinates[, c("col", "row")]) colnames(coordinate) <- c("x", "y")

这里有一个非常关键的点：表达矩阵的列名顺序必须和坐标矩阵的行名顺序完全一致。因为SPARK内部是按位置一一对应的，如果顺序乱了，后面跑出来的结果全无意义。我一般会在构建数据后做一个检查：

stopifnot(all(colnames(sp_count) == rownames(coordinate)))

这句代码如果报错，就说明表达矩阵的spot名和坐标的spot名顺序不一致，需要先对齐再往下走。这个检查我每次都做，治好了我无数次数据错位的毛病。

另外坐标的单位无所谓，是像素还是微米都不影响，SPARK内部会自己做标准化处理，但建议还是统一一下量纲，避免因为坐标数值过大导致计算溢出。

3.3 创建SPARK对象并运行分析

数据准备好之后，创建SPARK对象和跑分析就是核心几步了。

library(SPARK) # 创建SPARK对象 spark_obj <- CreateSPARKObject( counts = sp_count, location = coordinate, percentage = 0.1, # 过滤基因：至少在10%的spot中有表达 min_total_counts = 10 # 过滤基因：总表达量至少为10 ) # 计算总UMI数，作为测序深度的协变量 spark_obj@lib_size <- apply(spark_obj@counts, 2, sum) # 运行SPARK spark_obj <- SPARK( spark_obj, num_core = 8, # 并行核心数，根据服务器情况调整 covariates = "lib_size", # 把测序深度作为协变量 num_kernel = 10, # 核函数数量，一般默认即可 kernel_para = c(seq(1, 5, length.out = 5), seq(1, 10, length.out = 5)) )

这里解释几个参数。percentage和min_total_counts是过滤条件，避免那些在极少spot里表达、信息量太低的基因进入分析，也能减少计算量。num_core就是并行核心数，这个值别一味调大，得看服务器的CPU核数和内存，我自己试过在96核机器上设64反而比32慢，因为进程间通信开销变大了。num_kernel和kernel_para控制核函数的数量和尺度范围，默认参数在绝大多数场景里都能正常工作，如果没有特别理由，用默认值就好。

完成这一步之后，SPARK会在后台对所有基因逐一做模型拟合和显著性检验。跑完之后的检验结果需要单独提取：

spark_obj <- SPARK.test( spark_obj, check_positive = TRUE, # 只检验表达量方差为正的基因 output = "all" # 输出所有基因的检验结果 ) # 提取P值和FDR res <- spark_obj@res_mtest head(res)

res_mtest这个数据框里包含了每个基因的检验结果，包括原始P值（combined_pvalue）和校正后的FDR（adjusted_pval）。我们最终筛选显著基因，看的就是adjusted_pval这一列，一般取小于0.05作为显著阈值。

3.4 提取显著基因并可视化

筛选显著基因很简单：

svg_genes <- rownames(res[res$adjusted_pval < 0.05, ]) cat("识别到", length(svg_genes), "个空间可变基因\n")

到这里SPARK的分析就基本完成了。接下来比较自然的操作是把这些基因映射回Seurat对象，和已有的聚类、注释信息结合起来看。

# 把SPARK结果加回Seurat对象 obj <- AddMetaData(obj, metadata = res[colnames(obj)?理论上不是，这里是基因不是spot，所以不应该这样加] # 正确做法：把显著基因列表存下来，然后用FeaturePlot可视化这些基因的空间表达 DefaultAssay(obj) <- "Spatial" for (g in svg_genes[1:6]) { print(SpatialFeaturePlot(obj, features = g, image.alpha = 0.5)) }

注意这里有一个初学者容易踩的坑：res_mtest是基因维度的结果，不能直接AddMetaData到Seurat对象上（Seurat的metadata是spot维度的），要可视化就直接用SpatialFeaturePlot针对基因画，或者自己用ggplot2把坐标和表达量叠加画图。

我自己写过一个快速可视化函数，用ggplot2把目标基因的表达量按spot位置画成散点图，颜色深浅表示表达量高低，和SpatialFeaturePlot等价，但更灵活，感兴趣的朋友可以自己封装一下。核心就三步：取表达量、取坐标、合并成一个数据框，然后geom_point(aes(x=x, y=y, color=expr))一把梭。

3.5 数据量太大？直接上用SPARK-X

如果你的数据spot数超过一两万个，我建议直接跑SPARK-X，体验完全不一样。安装和SPARK在同一个包里面，不需要额外装东西。

# SPARK-X的输入需要先做一次过滤 sparkx_obj <- CreateSPARKObject( counts = sp_count, location = coordinate, percentage = 0.1, min_total_counts = 10 ) sparkx_obj@lib_size <- apply(sparkx_obj@counts, 2, sum) # 运行SPARK-X sparkx_obj <- SPARKX(sparkx_obj, num_core = 8) # 提取结果 res_x <- sparkx_obj@res_mtest svg_genes_x <- rownames(res_x[res_x$adjusted_pval < 0.05, ])

你看，从接口到结果获取，SPARK-X和SPARK保持了高度一致，最大的区别是运行时间。我用一份大约8000个spot、18000个基因的Visium乳腺数据做过对比：SPARK跑完用了大概半小时，SPARK-X只用了不到2分钟。而且两份结果取交集后，重合度相当高，绝大多数核心SVG都能被同时识别出来。所以如果你的数据没有苛刻的空间模式精度要求，SPARK-X完全可以当成主力工具。

4. 实操中的常见问题与排查技巧

4.1 问题速查表

我把实际使用中遇到和听同行反馈过的典型问题整理成了一张表，方便你遇到问题时快速定位。

4.2 关于结果解读的几点心得

SPARK跑出来的显著基因列表，我一般会分成两类去看。第一类是那些传统的marker基因，它们本身在已知组织结构里就有明确意义，比如脑组织里的层特异性基因、肿瘤组织里的增殖区域基因。这类基因出现在SVG列表里，是对分析流程的一个验证——如果已知的marker都没被识别出来，那大概率数据或参数有问题。

第二类是你没预期到的基因。这类基因往往才是SPARK这个工具价值的真正体现。它们可能是之前没有建立空间认知的基因，在表达总量上不起眼，但空间分布极度结构化。按照我的经验，拿到这类基因后值得做两步验证：一是回到原始切片图像上，染色或免疫组化确认一下；二是查一下这个基因在单细胞数据里的细胞类型特异表达情况，如果它在某个空间区域的占优细胞类型里特异表达，就形成了“空间域—细胞类型—功能基因”的完整证据链。

还有一点，P值排序不能完全代表生物学重要性。SPARK给出的P值是统计意义上的显著性，一个P值极小的基因，可能空间分布确实极度非随机，但生物学意义却有限。反过来，一些P值中等偏显著的基因，倒可能是某个过渡区域的有趣线索。所以筛选基因时我会把adjusted_pval < 0.05作为硬阈值，但排序时不光看P值，还会参考表达量在空间上的分化程度，比如高表达区域和低表达区域的表达倍数差。

4.3 我踩过的一个典型坑

有一次我处理一份小鼠大脑的Visium数据，跑完SPARK后惊了：KNN分类时效果很好的几个已知皮层marker竟然不在显著列表里，反而是几个核糖体蛋白基因冲到了最前面。排查了半天，发现问题出在过滤参数上。

我当时的min_total_counts设得太低，导致大量低表达基因进入了后续分析。低表达基因在空间分布上天然比较稀疏，一稀疏就容易被模型识别出“空间结构”——因为它们只在少数spot里有表达，其他spot都是零，这种模式从空间相关性的角度看确实显著，但生物学上大多没有意义。核糖体蛋白基因恰好属于表达广泛但局部波动的类型，捣乱效果尤其明显。

后来我把percentage调到0.1（即基因至少在10%的spot里有表达）、min_total_counts调到10以上，再跑一次，已知的marker基因就都回来了，结果明显合理多了。这个教训让我意识到，SPARK的参数虽然不多，但过滤条件直接决定分析质量，特别是稀疏数据场景下，宁可多过滤掉一些低表达基因，也不要让大量“零膨胀”基因进入模型。

4.4 多切片数据怎么办

单个切片跑SPARK很简单，但很多项目是多切片甚至跨样本比较。SPARK的核心实现是针对单个组织切片设计的，如果你有多个切片，一个惯常做法是对每个切片分别跑SPARK，得到各自的SVG列表，再取交集或并集。

取交集意味着你只保留在所有切片里都显著的空间可变基因，这种基因更稳健，适合作为这个组织类型或疾病状态的保守空间特征；取并集则能保留在不同切片中特异性出现的空间模式基因，适合做异质性探索。没有绝对的好坏之分，取决于你的科学问题。如果切片数量多且来自不同生物学条件下，我的做法是把并集做下游分析，但在生物学解读时，特别注意那些只在一个条件下出现的基因，它们往往是表型相关线索。

如果切片来自同一条件下多个生物学重复，还有一种做法是把所有spot合并成一个大坐标空间后统一跑SPARK，前提是各个切片的坐标能对齐到同一个参考坐标系，且批次效应不强。我自己不太推荐这种方式，因为切片间的测序深度、组织形态差异会产生额外噪声，反而干扰模型判断。

5. 这工具后续还能怎么扩展

SPARK跑出来的SVG列表，本身就是一个很好的基因集。你可以用它做基因集富集分析，看看这些空间可变基因主要集中在哪些通路；也可以把它们当作特征基因，结合空间域（spatial domain）的结果，看每个域特异表达了哪些SVG；还有人是把SVG和转录因子调控网络联合起来分析，找那些空间模式一致、可能处于同一调控模块的基因对，这个思路在发育生物学里特别有潜力。

目前空间转录组分析工具的迭代速度非常快，SPARK之后又有不少新方法出来，但SPARK/SPARK-X凭借其稳健的统计学框架和成熟的R生态，依然是SVG分析里的有力选项。做分析的时候我的建议是不必追求方法最新，关键是理解手里的工具解决的是什么问题、结果的边界在哪里。SPARK就是这样一个工具：如果你能在跑代码之前想清楚它检验的是什么假设，那这套结果你就真正用活了。

最后分享一个小经验：如果只是探索性分析，先跑SPARK-X快速过滤一遍，再用SPARK对筛选出的top基因做精细检验，这个组合方案既省时间又能保证结果质量，我目前用得最多的就是这个策略。