SpaGCN:整合基因表达、空间定位和组织学，通过图卷积网络识别空间域和空间可变基因

摘要：

空间解析转录组学(SRT)技术的最新进展使得在组织微环境背景下的基因表达模式的综合表征成为可能。为了阐明空间基因表达变异，我们提出了SpaGCN，这是一种图卷积网络方法，集成了SRT数据分析中的基因表达、空间定位和组织学。通过图卷积，SpaGCN从相邻点聚合每个点的基因表达，这使得具有连贯表达和组织学的空间域识别成为可能。随后的域引导差异表达(DE)分析在识别的域中检测到表达模式丰富的基因。通过使用SpaGCN分析7个SRT数据集，我们表明它可以检测出比其他竞争方法更丰富的空间表达模式的基因。此外，通过SpaGCN检测到的基因是可转移的，可用于研究其他数据集中基因表达的空间变异。SpaGCN计算速度快，独立于平台，是各种SRT研究的理想工具。

最近SRT技术的进步使得利用组织中的空间信息进行基因表达谱分析成为可能。了解组织中不同细胞的相对位置对于理解疾病病理学是至关重要的，因为空间信息有助于理解细胞的基因表达是如何受到周围环境的影响的。

流行的SRT实验方法大致可以分为两类。第一类是基于单细胞分辨率的原位杂交或测序技术，其中包括seqFISH2,3seqFISH +4, MERFISH5,6, STARmap7 和FISSEQ8 用来测量组织环境中细胞中成百上千个基因的表达水平。第二类是基于原位捕获的技术，其中包括空间条形码和测序，其中包括空间转录组学(ST)9, SLIDE-seq10， SLIDE-seqV2(引用; 11), HDST12 和10x Visium，测量捕获位置(称为点)中数千个基因的表达水平。这些不同的SRT技术使得发现异质组织的复杂转录结构成为可能，并增强了我们对疾病中的细胞机制的理解13,14。

在SRT研究中，一个重要的步骤是识别被定义为在基因表达和组织学上都具有空间一致性的区域的空间域。

我们没有将空间域和SVG（空间可变基因）识别视为独立的问题，而是开发了一种基于图卷积网络(GCN)的方法SpaGCN，它联合考虑了这两个问题。SpaGCN首先通过构建代表数据空间依赖性的无向加权图，整合基因表达、空间位置和组织学来识别空间域。对于每个空间域，接下来，SpaGCN进行检测在该领域中富集的svg。通过将搜索空间限制到空间域，保证了SpaGCN检测到的svg具有空间表达模式。空间域和相应的svg提供了组织中基因表达的空间梯度的全面图景。在分析多种SRT数据方面，SpaGCN是万能的，包括ST、10x Visium、SLIDE-seqV2、STARmap和MERFISH。

结果

SpaGCN的概述和评价。我们以使用基于原位捕获的SRT数据为例解释了SpaGCN的工作流程，但该方法可以很容易地修改以分析其他类型的SRT数据。如图1a所示，SpaGCN首先构建一个图来表示所有点之间的关系，同时考虑了空间位置和组织学信息。接下来，SpaGCN利用一个图卷积层来聚合来自相邻点的基因表达信息。然后，SpaGCN利用聚合后的表达矩阵，使用无监督迭代聚类算法对点进行聚类26。每个聚类被认为是一个空间域，

（图1a）SpaGCN工作流程。a, SpaGCN首先使用图卷积网络(GCN)整合基因表达、空间定位和组织学信息，然后使用无监督迭代聚类将点分离到不同的空间域。GCN基于一个无向加权图，其中每两个点之间的边缘权重由两个点之间的欧氏距离决定，由空间坐标(x,y)和三维坐标z定义，由组织学图像中的RGB值获得

然后通过DE分析从该空间域中检测到在某个域中富集的svg（空间可变基因）(图1b)。当单个基因无法标记某一区域的表达模式时，SpaGCN将构建一个由多个基因组合形成的元基因，来表示该区域的表达模式。

（图1b）对于每一个检测到的空间域，通过域引导DE分析，SpaGCN识别SVGs或meta基因。

为了展示SpaGCN的强度，我们将其应用于7个公开可用的数据集(补充表1)。与Louvain、stLearn和BayesSpace相比，由SpaGCN识别的空间域与已知组织结构的一致性更好。我们还将SpaGCN检测到的svg与SpatialDE和SPARK检测到的svg进行了比较，发现与其他两种方法相比，SpaGCN检测到的svg具有更连贯的表达模式和更好的生物学可解释性。通过Moran’s I和Geary’s C统计，进一步证实了spagcn检测到的SVGs所揭示的空间表达模式的特异性27，两个常用的量化基因表达空间自相关性的指标28,29。

(论文使用不同SRT技术生成的不同物种、区域和组织的数据集上，已经对SpaGCN进行了广泛的测试，但本文未对这些实验进行详细分析。主要注重网络的作用和spaGCN实现的方法）

应用于人类原发性胰腺癌ST数据。为了证明结合组织学信息的重要性，我们分析了使用ST技术生成的人类原发性胰腺癌数据集13。这个数据集包括224个点和16448个基因，以及三个人工标注的组织区域。Louvain仅基于基因表达检测到的肿瘤区域与病理学家注释的肿瘤区域并不紧密匹配(图2a)。

 

stLearn和BayesSpace等空间聚类方法也没有检测到癌症区域。在使用默认参数时，SpaGCN揭示了类似的模式。由于组织学图像显示了癌症和非癌症区域之间的明显差异，这表明组织学对聚类有信息。SpaGCN具有用缩放参数s对组织学建模的灵活性，该参数在为每个点检测邻居时控制给组织学的权重。通过将s的值从1增加到2,SpaGCN检测到一个与人工标注的癌症区域很吻合的集群。值得注意的是，当s设置为默认值1时，SpaGCN能很好地检测到非癌症区域。当s增加到2时，SpaGCN不仅保持了检测非癌症区域的能力，还检测到了癌症区域。这个例子表明，在聚类过程中，SpaGCN能够灵活地融合组织学信息。虽然stLearn可以合并组织学数据，但在定义相邻点时，它对组织学信息的使用是由半径预先固定的。在调整组织学权重方面缺乏灵活性，导致它们的聚类与病理学家的手工注释之间存在差异。

接下来，我们使用SpaGCN、SPARK和SpatialDE检测svg。总共，SpaGCN检测到12个svg，域0、域1和域2分别有3个、8个和1个svg(图2b;补充图1)。

 

应用于人背外侧前额叶皮层10x Visium数据。为了定量显示SpaGCN在空间域检测方面优于Louvain、stLearn和BayesSpace，

为了证明由spagcn检测到的svg对下游分析是有用的，我们在切片151507上进行了K-means聚类，其中来自不同的大脑，使用了由SpaGCN从151673片检测到的所有67个svg。与手动策划的层分配相比，这种聚类分析的调整兰德指数(ARI)为0.23(图3e,f)。我们使用SpatialDE和SPARK检测到的svg进行了类似的分析。当从SpatialDE/SPARK检测的基因中随机选取67个P或Q值为0的svg时，SpatialDE的ARI仅为0.13,SPARK为0.14。为SpatialDE的ari即使svg数量增加，SPARK也没有改善(图3e)。这些结果进一步证实了SPARK和SpatialDE检测到的基因缺乏空间模式。

 尽管很难识别单个基因来标记某些神经元层，但SpaGCN能够找到域特异性的元基因。如图3g所示，SpaGCN检测到区域1、2、4和6的元基因。

讨论

在本文中，我们提出了一种整合基因表达、空间定位和组织学的方法SpaGCN，以建模基因表达的空间依赖性，用于识别空间域和域丰富svg。在使用不同SRT技术生成的不同物种、区域和组织的数据集上，已经对SpaGCN进行了广泛的测试。ST的附加分析9， SLIDE-seqV2(参考文献。 11)和MERFISH5 资料载于补充注释1-3。我们的结果一致表明，SpaGCN可以识别具有一致基因表达和组织学的空间域，并且检测svg和元基因，它们比SpatialDE和SPARK检测的基因具有更清晰的空间表达模式和生物学解释。此外，spangn检测到的SVGs是可转移的，可以用于独立组织切片的下游分析。与SPARK和SpatialDE相比，SpaGCN计算速度快，内存效率高(补充说明4)。

在SpaGCN中的空间域检测步骤是灵活的。首先，在基因表达平滑中，SpaGCN可以调节组织学的权重。对于组织学中组织结构清晰的数据集，更高的权重导致了癌症与非癌区域更清晰的分离。其次，在GCN拟合过程中，图权值被更新，这使得SpaGCN能够学习到一种有效的方法，从每个基因的邻近点聚合基因表达。对于来自不同平台生成的数据，随着捕获组织面积的大小不同，点/细胞之间的空间依赖性也不同。建模空间依赖性的灵活性使得SpaGCN对于不同类型的SRT数据具有通用性。

SpaGCN的一个局限性是空间域检测主要由基因表达驱动，这可能导致被检测区域与底层组织解剖结构的差异。这是基于基因表达的聚类方法普遍存在的问题。SpaGCN的另一个局限性是，在检测到的svg的基因表达模式中，缺乏对空间变异和细胞类型变异的分离。为了解决这些局限性，需要在聚类中联合考虑基因表达和组织学特征的方法。此外，还需要估计细胞类型特异性基因表达，以梳理出细胞类型和空间位置对基因表达变异的贡献。我们预计，在未来的研究中，这些方向的方法发展是有保障的。

方法

数据预处理。

SpaGCN将空间基因表达和组织学图像数据（如果有的话）作为输入。为了便于解释，我们将使用基于原位捕获的SRT数据来说明该方法。空间基因表达数据存储在一个N×D的独特分子标识符（UMI）矩阵中，其中有N个斑点和D个基因，以及每个斑点的（x，y）二维（2D）空间坐标。在少于三个点上表达的基因被剔除。每个点的基因表达值被归一化，即每个基因的UMI计数除以给定点的所有基因的UMI总计数，再乘以10,000，然后转换为自然对数尺度。

将SRT数据转换为图结构数据。

经过预处理，SpaGCN将基因表达和组织学图像数据转换为加权无向图，G(V,E)。在这个图中，每个顶点v∈V代表一个斑点，V中的每两个顶点都通过一条具有特定权重的边连接。我们注意到，spage2vec35也采用了基于图的方法，但目标是对信使RNA分子进行聚类，其中每个节点代表一个mRNA。这种无分割的方法可能比需要分割细胞的方法有优势，因为细胞分割仍然是单细胞分析中最难的问题之一。

两个顶点之间的距离计算。图中任何两个顶点u和v之间的距离反映了两个相应点的相对相似性。这个距离是由两个因素决定的。(1）斑点u和v在组织切片中的物理位置，以及（2）这两个斑点的相应组织学信息。尽管有些斑点在组织中的物理位置很接近，但组织学图像可能显示它们属于不同的组织层。

因此，SpaGCN认为两个斑点是接近的，当且仅当（1）这两个斑点在物理上接近，以及（2）它们在组织学图像中具有相似的组织学特征。为了定义一个考虑到这两个方面的距离指标，SpaGCN将组织切片中的二维空间扩展为一个包含组织学信息的三维空间。对于斑点v，其在组织切片中的物理位置由二维坐标（xv, yv）表示。为了确定斑点v在组织学图像中的对应像素，SpaGCN根据其像素坐标（xpv, ypv）将斑点v映射到组织学图像上。SpaGCN不使用位于（xpv, ypv）的像素的颜色，而是以（xpv, ypv）为中心画一个包含50×50个像素的正方形，并计算所有落在正方形的像素的RGB通道（rv, gv, bv）的平均颜色值。这一步平滑了颜色值，确保颜色不被单个像素所支配。为了得出一个代表组织学图像特征的单一数值，SpaGCN使用RGB值的加权和，如下所示

 其中，Vr=方差（rv），Vg=方差（gv），Vb=方差（bv），对于这个变换中的所有v∈V，对具有较大方差的通道给予较高的权重，以便这个组合值zv能准确地代表组织学图像的模式。

接下来，SpaGCN将zv的规模重新调整为

其中μz是zv的平均值，σx、σy、σz分别是v∈V的xv、yv和zv的标准差，s是一个比例因子。默认情况下，s被设置为1，以确保z∗v具有与xv和yy相同的尺度方差，当目标是增加组织学的权重时，我们将s设置为大于1的值。在扩展的三维空间中，斑点v的坐标被设定为（xv, yv, z∗v ）。最后，每两个点u和v之间的欧几里得距离被计算为

 当组织学信息不可用时，每两个点之间的欧氏距离将只根据空间位置信息来计算。

 计算每条边的权重并构建图。每条边（u,v）的权重衡量点u和点v之间的相关程度，并与它们的距离呈负相关。图形结构G存储在一个N×N的邻接矩阵A=[w (u, v)]中，其中点u和点v之间的边缘权重定义为

超参数l，也被称为特征长度尺度，决定了权重作为距离的函数衰减的速度。SpatialDE24中也采用了类似的函数。让I表示身份矩阵。对于斑点v，A-I的相应行和，用av表示，可以解释为其他斑点对其基因表达的相对贡献。默认情况下，我们选择l的值，使所有斑点的平均数等于一个预先指定的值，例如0.5。对于组织捕获面积较小的SRT平台产生的数据，例如SLIDE-seqV2、STARmap和MERFISH，我们建议选择l的值，使相邻的点/细胞在基因表达聚集中贡献更多的信息。

图形卷积层。

SpaGCN使用主成分分析（PCA）降低了预处理的基因表达矩阵的维度。前50个主成分被用作输入，这对本文分析的所有数据集都很有效。接下来，利用图卷积网络的力量，SpaGCN将基因表达信息和G中的边缘权重连接起来，对节点进行聚类。根据Kipf和Welling36，图卷积层可以写成

 

其中X是由PCA得到的N×50的嵌入矩阵，B是代表卷积层过滤参数的50×50矩阵，δ(-)是一个非线性激活函数，如ReLU。图形卷积层确保B中相应的一行参数将控制X中每个特征的邻域信息的聚合，从而提供了对邻域点提供的信息进行特定特征聚合的灵活性。B中的过滤器参数在图中的所有顶点上共享，并在迭代训练的过程中自动更新。通过图卷积，SpaGCN根据G中指定的边缘权重聚合了基因表达信息。这一层的输出是一个聚合矩阵，包括基因表达、空间位置和组织学信息。图卷积层是基于Kipf和Welling36实现的，其中反向传播是通过谱图卷积的局部一阶近似来操作的。

 

通过聚类识别空间域。接下来，基于上述图卷积层的输出，SpaGCN采用无监督的聚类算法，反复将斑点聚类到不同的空间域。从这个分析中确定的每个聚类被认为是一个空间域，它包含在基因表达和组织学上一致的点。为了初始化聚类中心点，我们在图卷积层的聚集输出矩阵上使用Louvain方法15。如果组织中的域的数量是已知的，Louvain的分辨率参数将被设置为产生相同数量的空间域。否则，我们将分辨率参数从0.2到1.0变化，并选择给出最高Silhouette分数的分辨率37。

为了迭代更新聚类的分配，我们定义了一个指标来衡量从一个点到聚类中心点的距离，使用学生的t分布作为内核。嵌入点i与聚类j的中心点μj之间的距离是指

可以被解释为将单元i分配到群组j的概率。

接下来，我们通过定义一个基于qij的辅助目标分布P来迭代细化群组

 

它提高了分配给高置信度的斑点的权重，并将每个中心点对整个损失函数的贡献归一化，以防止大型集群扭曲了隐藏的特征空间。现在我们有了软分配qij和辅助分布pij，我们可以将目标函数定义为Kullback-Leibler（KL）分歧损失。 

网络参数和聚类中心点同时通过使用随机梯度下降的动量最小化L来优化。这种无监督的迭代聚类算法以前曾被用于scRNA-seq分析，并显示出比Louvain的方法更优越的性能。

聚类后，SpaGCN还为聚类结果提供了一个可选的细化步骤。在这个步骤中，SpaGCN检查了每个斑点及其周围斑点的领域分配。对于一个给定的斑点，如果它周围的斑点有一半以上被分配到不同的领域，这个斑点将被重新标记为与它周围斑点的主要标签相同的领域。由于这一细化步骤只对少数斑点进行重新标记，因此对下游的SVG检测影响不大。

我们只对人类背外侧前额叶皮层10倍Visium数据和STARmap数据进行了聚类细化，并与它们的人工注释进行了明确的领域界限比较。

 

svg的检测。我们感兴趣的是检测在每个空间域都丰富的svg（空间可变基因）。我们注意到，一些基因可能在多个但不相连的区域中表达。虽然它们在某一特定区域内没有唯一表达，但这些基因对于理解基因表达的空间变异仍然是有用的，可以用来形成在某一特定区域内唯一表达的元基因。因此，我们不是使用来自目标域的点与所有其他点进行DE（域引导差异表达）分析，而是首先选择点来形成目标域的邻接集合。目标是检测在目标区域中高度表达但在邻近点中不表达或表达水平较低的基因。补充图25说明了相邻区域是如何被识别的。简单地说，我们在目标域中的每个点周围画一个半径预先指定的圆，来自非目标域但驻留在圆内的点被认为是它的邻居。半径设置为这样，目标域中的每个点平均大约有十个邻居。接下来，收集目标域中所有点的邻居，形成一个邻居集。对于每个非目标域，如果它的点超过50%(默认)在相邻集中，那么这个域就会被选为相邻域。设置这个准则是为了避免某个域被选为邻接域，但它的点中只有一小部分与目标域相邻的情况。

相邻域确定后，接下来，SpaGCN使用Wilcoxon秩和检验在目标域和相邻域的点之间进行DE分析。选择FDR-adjusted P值<0.05的基因作为svg。为了确保只选择目标区域表达模式丰富的基因，我们进一步要求一个基因满足以下三个标准:(1)在目标区域表达该基因的斑点的百分比，即in-fraction，为>80%;(2)对于每个邻近区域，目标区域内表达该基因的斑点的百分比与邻近区域(s)的比值，即in/out分数比，为>1;以及(3)靶域和邻近域(s)之间的表达倍数变化为>1.5。如果用户对寻找特定空间域组合的svg感兴趣，SpaGCN提供了这样做的选项。

空间变量元基因的检测。上面描述的空间域特异性DE分析通常检测到大多数域表达丰富的SVGs。对于未检测到此类SVGs的结构域，我们的目标是识别一组基因，当组合成一个元基因时，在给定的结构域显示出丰富的表达模式。为了识别形成元基因的基因，我们采用了多步骤的方法。首先，我们降低SVG筛选的阈值，例如，将最小折叠变化阈值从1.5更改为1.2，以识别出在目标域中表现出较弱的富集表达模式的基因。在存在多个这样的较弱svg的情况下，我们随机选择其中一个作为基基因，并将其标记为基因0。其次，我们的目标是聚合其他基因到基基因的表达，以增强目标区域的空间格局。为了实现这一目标，我们首先计算基因的平均表达水平0 对于目标区域内的点为e0。那么，所有表达水平高于e的非目标区域的点0 的基因0 被提取出来组成对照组。接下来，我们使用来自目标域的点对对照组中的点使用Wilcoxon秩和检验进行DE分析。选择fdr调整后P值最小、在目标区域表达量较高的基因作为基因0+。同样，我们将对照组的斑点与目标区域的斑点进行DE分析，并选择在对照组中fdr调整P值最小、表达量较高的基因作为基因0−。元基因的表达计算为

其中C0是一个常数，使log（meta_gene1）非负值。对数转换是用来重定表达量，使不同基因的表达水平具有可比性。我们发现，包括阴性基因可以加强没有富集阳性标记基因的领域的空间表达模式。这种算法可以反复用来寻找额外的基因，以形成一个更新的元基因，使目标域的空间模式更加清晰。对于第（t+1）次迭代，元基因的表达量计算为

 

在第（t+1）次迭代中，在加上genet+和减去genet-后，SpaGCN将根据meta_genet+1选择第（t+1）个控制组。新的控制组的大小，即不在目标域但meta_genet+1的表达量高于目标域的斑点的数量，应该小于第t个控制组的大小，以确保meta_genet+1比meta_genet有更清晰的空间模式。同时，预计meta_genet+1在目标组和对照组之间的平均表达量差异会比meta_genet大。因此，在每次迭代时，SpaGCN都会检查是否满足这两个标准，否则将停止搜索额外的基因。这种迭代式元基因搜索的图示见补充图26。

使用Moran's I和Geary's C统计学对SVG进行评估。不同位置的基因表达可能不是独立的。例如，一个基因在附近位置的表达水平可能比相距较远的位置的表达水平在数值上更接近。这种现象被称为空间自相关，它衡量的是一个变量通过空间与自身的关联性。为了评估检测到的SVG是否表现出有组织的空间表达模式，我们使用了Moran's I和Geary's C，这两个常用的统计量来量化基因表达的空间自相关程度。空间自相关可以是正的或负的。当相似的值出现在彼此附近时，就会出现正的空间自相关。负的空间自相关发生在不相似的值相互靠近时。Moran's I指标27是一个相关系数，用于衡量数据集的总体空间自相关性。直观地说，对于一个给定的基因，它衡量一个点与它周围的其他点的相似程度。如果这些斑点相互吸引（或排斥），这意味着这些斑点不是独立的。因此，自相关的存在表明基因表达的空间模式。Moran's I值的范围是-1到1，其中接近1的值表示明确的空间模式，接近0的值表示随机的空间表达，而接近-1的值表示类似棋盘的模式。为了评估一个给定基因的空间变异性，我们用以下公式计算莫兰氏I。

 

其中xi和xi是i和j点的基因表达量，¯x是基因的平均表达量，N是点的总数，wij是用点的二维空间坐标计算的i和j点之间的空间权重，W是wij的总和。

对于每个斑点，我们用空间坐标选择k个最近的邻居。Moran's I统计量对k的选择是稳健的，在我们的分析中被设定为4。如果j点在㼿点的最近邻居中，我们就把wij=1，否则wij=0。

Geary's C是另一个常用于测量空间自相关的统计学。它的计算方法是

 

Geary's C的数值范围是0到2。为了使它与Moran's I在同一刻度上，我们把它转换成[-1,1]的刻度，即

 

其中1表示完全正相关，0表示没有自相关，-1表示完全负相关。对于每个基因，Geary's C的值与Moran's I相似，但不完全相同。

SpaGCN检测的SVG的可转移性。

一个真正的SVG应该表现出可转移性，也就是说，空间表达模式在从同一组织类型收集的不同数据集中应该是相似的。为了显示SpaGCN检测到的SVG的可转移性，首先，我们显示在一个数据集中检测到的SVG在一个独立的数据集中也显示出类似的空间表达模式。其次，我们表明在一个数据集中检测到的SVG可以用来对一个独立数据集中的斑点进行聚类，并达到相对较高的聚类精度。

与其他方法的比较。为了评估SpaGCN在识别空间域方面的表现，我们与Louvain, BaysSpace, stLearn和HMRF进行了比较。我们对所有的方法都使用默认的参数设置，并且在聚类时使用相同的聚类数量。为了评估SpaGCN在检测SVG方面的性能，我们与SPARK和SpatialDE进行比较。我们最初使用它们的默认参数设置来检测SVG。默认情况下，SPARK和SpatialDE都会过滤掉总UMI计数小于10的点。在SPARK中，在少于10%的点中表达的基因被过滤掉。在SpatialDE中，在少于三个点中表达的基因被过滤掉。为了评估过滤标准对SPARK和SpatialDE的影响，我们进一步剔除在少于20%的斑点中表达的基因。检测到的SVG的空间表达模式的特异性由Moran's I和Geary's C统计来评估。