由于数据的高维和复杂性，空间转录组数据的分析一直是一个具有挑战性的问题，而聚类分析则是空间转录组数据分析的核心问题。本文提出了一种基于图注意力网络的深度学习方法，用于空间转录组数据的聚类分析。首先，对空间转录组数据进行增强，然后使用图注意力网络对节点进行特征提取，最后使用莱顿（Leiden）算法进行聚类分析。通过聚类的评价指标发现，与传统的非空间及空间聚类方法相比，本文提出的方法具有更好的数据分析性能。实验结果表明，本文所提方法可以有效地聚类空间转录组数据，从而能够识别不同的空间区域，为研究空间转录组数据提供了新的工具。
引用本文： 吴瀚文, 高洁. 通过图注意力网络识别空间转录组中的空间域. 生物医学工程学杂志, 2024, 41(2): 246-252. doi: 10.7507/1001-5515.202304030 复制
0 引言
空间转录组学是一种结合了组织学和转录组学的新兴技术，它可以研究细胞在组织中的分布和表达方式。空间转录组（spatial transcriptome，ST）数据通常是由高通量成像技术生成，如空间转录组学和空间RNA测序（spatial RNA sequencing，sRNA-seq）。这些技术可以在组织中捕获数百万个单细胞的转录组数据，并提供高分辨率的空间信息。然而，由于数据的高维和复杂性，空间转录组数据的分析一直是一个具有挑战性的问题。识别空间域（即，在基因表达和组织学上都具有空间一致性的区域）是空间转录组学中最重要的课题之一。目前，空间域识别方法主要分为非空间聚类和空间聚类两类。传统的非空间聚类算法，如K均值聚类算法（k-means clustering algorithm）、鲁汶算法（Louvain method），以基因表达数据作为输入，导致聚类难以符合组织切片[1-2]。而空间聚类方法结合了基因表达、空间位置和形态学，可以解释基因表达的空间依赖性，以更好地匹配空间位置。例如，Zhao等[3]提出的空间贝叶斯方法（BayesSpace）即是基于贝叶斯统计方法，通过引入空间邻居结构来鼓励相邻点属于同一聚类，然后将捕获区域划分为更小的亚区域，实现了聚类精度的提高。Xu等[4]提出的用于空间转录组分析的卷积神经网络（convolutional neural networks，CNN）学习策略（CNN learning strategy for spatial transcriptomics analysis，CosTA）使用CNN聚类结构，通过CNN生成特征，再通过高斯混合模型聚类生成软分配[5]，并使用软分配来更新CNN。Yuan等[6]提出的用于基因的图CNN（graph CNN for genes，GCNG）方法将空间信息编码为图形，使用监督训练将其与表达数据组合，并用图卷积网络对图数据进行特征嵌入，从而实现聚类。Hu等[7]提出的整合基因表达、空间位置和组织学，通过图卷积网络识别空间域和空间可变基因（integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network，SpaGCN）方法应用图卷积网络来整合基因表达和空间位置，并进一步与自监督模块相结合来识别空间域。Pham等[8]提出的空间转录组学数据的下游分析工具包（stLearn）提供了一种组织内归一化技术，该技术基于从形态学图像中收集到的特征和空间位置，通过形态学距离对基因表达进行归一化。Xu等[9]提出的无监督空间嵌入深度表示方法（unsupervised spatially embedded deep representation of spatial transcriptomics，SEDR）采用深度自动编码器网络和图形自动编码器来嵌入空间信息。Teng等[10]提出可用于荧光原位杂交技术的迭代细胞分型（fluorescence in situ hybridization iterative cell typing，FICT）方法，能同时利用空间转录组数据中的表达信息和空间信息来识别细胞类型和亚型。尽管这些空间聚类方法可以在一定程度上识别空间域，但往往是通过节点的空间位置计算节点之间的距离，并由此预先定义节点之间的相似性，而这些预先定义的方法没能利用深度学习的灵活性，故存在一定局限。
为了解决这个问题，本文提出了一种基于图注意力网络（graph attention networks，GAT）的深度学习方法，并将其命名为用图注意力自动编码器进行空间转录组聚类分析（spatial transcriptomics clustering analysis with graph attention auto-encoders，STCAGAE），用于空间转录组数据的聚类分析[11]。本文使用的图注意力网络可以自适应地学习节点之间的相似性，从而更好地聚类空间域。本文利用预先训练好的深度神经网络模型从形态图像块中提取特征向量，然后将提取的特征与基因表达和空间位置数据相结合，以表征空间相邻点的相关性，并创建空间增强基因表达矩阵。同时使用图注意力自动编码器和去噪自动编码器联合生成增强空间转录组数据的潜在表示，最后使用莱顿（Leiden）算法最终的嵌入进行聚类[12]。本文旨在识别不同的空间域，并为研究空间转录组数据提供新的方法。
1 方法
1.1 数据来源及预处理
本文研究所用数据来自十倍基因组学（10 × Genomics）生物技术公司利用维希姆（Visium）技术发布的高通量空间转录组测序平台（10 × Genomics Visium）。该平台可以利用Visium技术，从整个组织切片和各种各样的样本类型中获得整个转录组高通量基因表达分析结果。本文选取了该平台的一组可公开下载的数据集，即人类背外侧前额叶皮层（dorsolateral prefrontal cortex dataset，DLPFC）数据集。DLPFC数据集包含了12个不同的切片数据，其中每个切片都包含6个皮质层和1个白质层。首先，从DLPFC数据集中移除主要组织区域之外的区域。然后使用计算机编程语言Python 3.8（Guido van Rossum，荷兰）中分析单细胞数据的软件包Scanpy 1.9.1（Theis Lab，德国），根据文库大小对原始基因表达数据进行过滤、对数转换和标准化[13]。最后，使用主成分分析（principal components analysis，PCA）对增广基因表达数据进行降维，并将降维数据作为下一个模型训练的输入。
1.2 STCAGAE算法
1.2.1 增强空间数据
空间基因表达技术提供了具有额外空间位置信息和组织形态的转录组范围的基因表达谱。本文使用这两个额外的空间数据实现了节点的基因表达数据增强，得到了对节点更真实的表示，从而提高了STCAGAE的性能。计算点i和点j之间的空间基因表达权重GCij，如式(1)所示：
 |
其中，Si代表节点i的基因表达向量，
代表节点i的基因表达向量的均值，Sj代表节点j的基因表达向量，
代表节点j的基因表达向量的均值。对于具有形态学信息的空间转录组数据，本文首先根据每个点的坐标对一个图像进行分割，得到其部分图像。然后，使用深度学习框架Pytorch 1.8.0（Facebook Inc，美国）中的图像预处理模块转换器（transform）对图像进行变换和增强，包括归一化、旋转、调整锐度等[14]。每个点的特征从预先训练的CNN模型得出，该模型可以将每个节点的图像转换为2 048维的潜在变量。为了更好地表示节点的形态，本文进行了PCA，提取了前50个主成分作为潜在特征。最后，利用余弦距离计算点i与相邻点j之间的形态相似性权重MSij，如式(2)所示：
 |
本文使用空间坐标来确定每个点和所有其他点之间的距离，然后对前4个相邻点之间的距离进行排序，以计算半径r。对于一个给定的点i，当且仅当两个点之间的距离小于r时，认为点j是点i的相邻点，则设表示点i与点j位置关系的权重参数SWij = 1，否则SWij = 0。然后，增强每个点i的基因表达，如果是10 × Genomics Visium平台的数据，将数据增强为如式(3)所示；否则，将数据增强为如式(4)所示：
 |
 |
其中，GEi和GEj代表点i和点i的相邻点j的原始基因表达。
1.2.2 构造图
为了结合给定点的相邻点的相似性，本文使用空间坐标来计算点之间的距离，并使用前6个最近邻构建一个无向连接图。对于无向连接图，在每个点上都添加了自循环。A是邻接矩阵，那么如果i和j是相邻点，则在点i和点j处的值Aij = Aji = 1，否则Aij = 0。
1.2.3 去噪自动编码器
使用Pytorch 1.8.0（Facebook Inc，美国）的线性层实现了一个去噪自动编码器，用于基因表达的潜在表示[15]。编码器E由用户设置的多个完全连接的堆叠线性层组成，将集成的基因表达式X转换为低维表示，如式(5)所示：
 |
其中，Zg代表经过编码器得到的低维表示，R代表实数集，N代表节点数量，M代表输入基因数量，d代表最后一个编码器层的维数。然后，向Zg加入噪声得到重构的低维表示，如式(6)所示。接下来，解码器D将重构后的低维表示转换为重构的基因表达矩阵，如式(7)所示：
 '/> |
 '/> |
其中，
代表重构的低维表示，Z代表噪声向量，D代表解码器，X'代表重构的基因表达矩阵，d''代表最终层的维数，N、M和d与上述相同。输入基因和重建表达式的类似性用均方误差来计算，如式(8)所示：
 |
其中，Ll代表去噪自编码器的损失函数，Xi代表节点i的基因表达向量，E（Xi）代表节点i的基因表达向量经过编码器层得到的低维表示，D（∙）代表经过解码器层重构后的表示。
1.2.4 变分图自动编码器
变分图自动编码器使用了基于Pytorch 1.8.0（Facebook Inc，美国）构建的用于编写和训练图神经网络的库Pytorch Geometric 2.2.0（Technische University Dortmund，德国） [16]，其中编码器由图注意力网络组成。图注意力网络将特征矩阵X作为输入，首先通过一个线性变换W做维度转换，如式(9)所示：
 |
其中，
代表特征矩阵经过线性变换之后的列数。接着使用自注意力（self-attention）为每个节点分配注意力，对于点i和点j，它们的初始特征为hi和hj，于是它们经过线性变换W后变为Whi和Whj，为它们分配注意力，如式(10)所示：
 |
其中 ，eij 代表节点 i 对节点 j 的影响力系数，a (·)代表计算节点相关度的函数，如式(11)所示：
 '/> |
在本文中，对于节点i，只计算其邻域内节点j对节点i的影响力系数。为了更好地在不同节点之间分配权重，本文将目标节点与所有相邻点计算出来的相关度进行统一的归一化处理，使用线性整流函数（rectified linear unit，ReLU），如式(12)所示：
 |
其中，
代表归一化后的相关度权重系数，Ni代表节点i的邻域，eik代表节点 i 对节点 i邻域内的节点 k 的影响力系数。每个节点的最终输出特征向量如式(13)所示：
 '/> |
其中，hi'代表节点的最终输出特征，σ代表激活函数。损失函数包括生成图与原始图之间的重构损失，节点表示的向量分布和正态分布的相对熵（relative entropy），如式(14)所示：
 |
其中，Lg代表变分图自动编码器的损失函数，
代表二进制交叉熵函数，
代表相对熵，
，代表正态分布。
同时，本文利用多头注意力以稳定学习过程和增加模型容量。多头注意力就是将同一个输入分别送入不同的注意力层，分别计算隐藏状态，然后将不同的特征表示连接起来（或相加）。经过不断调整参数并实验，本文发现模型中的参数取8个注意力头和特征相加时，DLPFC数据集的不同切片聚类得到的调整兰德系数（adjusted Rand index，ARI）值都最大，因此模型达到最佳效果。
1.2.5 聚类与可视化
在得到嵌入之后，使用Leiden算法来识别空间域，以识别相似的细胞和模式。为了获得最佳分辨率，在0.1～2.5之间进行网格搜索，步长为0.01，直到达到必要的簇数。为了可视化，使用了统一流形逼近和投影[17]。
综上，本文提出STCAGAE算法，如图1所示。通过对整合基因表达、空间位置和组织形态信息的低维表示进行建模来表征空间域。为了建立形态学特征矩阵，首先使用预先训练好的深度学习网络对苏木精-伊红（hematoxylin-eosin，H&E）染色的组织地形图数据进行处理。结合形态特征和空间相邻点信息，每个点的基因表达均得到增强。然后，利用去噪自编码器学习从集成特征空间到低维表示空间的非线性映射，以减少模型的过拟合。同时，基于节点之间的欧氏距离计算图邻接矩阵。将增强后的基因表达矩阵和图邻接矩阵输入到一个变分图自动编码器中，从而通过与相应的空间相邻点的集成表示来生成空间嵌入。最终的潜在嵌入是通过连接集成表示和空间嵌入来形成的。最后使用Leiden算法进行空间域聚类，识别出从右上到左下依次是皮质层1～6、白质层的不同空间域。
图1
STCAGAE算法的流程示意图
Figure1.
The flow chart of STCAGAE algorithm
2 结果
本文在一个经典空间转录组数据集——DLPFC数据集上进行了实验，验证了STCAGAE算法的有效性。STCAGAE算法可以将细胞聚类成多个簇，并且簇之间具有明显的空间分布和表达差异。与传统的聚类方法相比，STCAGAE算法在聚类性能上更有优势。
本文在DLPFC数据集的编号为151507的切片上通过STCAGAE算法对空间域进行聚类。同时，比较了其他空间算法在该数据集上的聚类性能。通过ARI值计算识别空间域的精度，如图2所示，STCAGAE算法的ARI值达到了0.620，优于通过深度学习识别空间转录组学中的空间域（identifying spatial domains in spatial transcriptomics by deep learning，DeepST）方法的ARI值0.559[18-19]。
图2
DLPFC数据集151507切片上的ARI值对比
Figure2.
The comparison of ARI on section 151507 in the DLPFC dataset
通过STCAGAE算法，又在DLPFC数据集的编号为151673的切片上进行了空间域的聚类。同时也比较了其他空间算法在该数据集上的聚类性能。通过ARI值计算识别空间域的精度，如图3所示ARI值达到了0.620，优于DeepST算法的ARI值0.615。
图3
DLPFC数据集151673切片上ARI值的对比
Figure3.
The comparison of ARI on section 151507 in the DLPFC dataset
对于DLPFC数据集的151676切片，本文通过STCAGAE算法进行了空间域的聚类。通过ARI值计算识别空间域的精度，如图4所示ARI值达到了0.638。然后，本文在该数据集上通过其他空间算法进行实验，并比较各种聚类算法的性能，发现STCAGAE效果优于软件包Scanpy 1.9.1（Theis Lab实验室，德国）（ARI = 0.29）、SEDR（ARI = 0.42）、SpaGCN（ARI = 0.44）、BayesSpace（ARI = 0.44）和用于空间解析转录组学的自适应图注意力自动编码器（deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder，STAGATE）方法[20]（ARI = 0.60）。
图4
DLPFC数据集上151676切片的人工注释以及STCAGAE聚类结果
Figure4.
Human annotation and clustering result of section 151676 in the DLPFC dataset
本文在DLPFC数据集上测试STCAGAE的性能，发现STCAGAE可以有效地将数据划分为6个皮质层和1个白质层[21-22]。其中，皮质层1主要由神经胶质细胞、轴突、树突和少量神经元组成，对神经元活动的调节非常敏感，功能涉及到信息处理、信号传递和空间感知等方面；皮质层2是DLPFC的主要输入层，具有调节和控制信息传递的功能，在决策制定、工作记忆和灵活的认知控制等高级认知过程中扮演着重要的角色；皮质层3主要包含锥形神经元，在视觉加工和感知处理中扮演重要的角色；皮质层4主要包含大锥形神经元，参与计划、决策和行动控制等许多高级认知过程；皮质层5主要由大锥形神经元和星形神经元等神经元组成，对调节神经元活动起到重要的作用，同时也与情绪调节、行动控制和认知过程等方面有关；皮质层6主要由星形神经元和横向神经元等组成，在感知处理、计划和动作控制等方面具有重要的作用；白质层是DLPFC的最内层，主要由神经元轴突组成。白质对于DLPFC与其他脑区的信息交流具有重要作用，同时也与情绪调节、决策制定和计划等各种高级认知过程密切相关。与其他聚类方法相比，STCAGAE实现了显著的改善。STCAGAE清晰地勾画出了层的边界，并获得了最好的聚类精度。而通过聚类得到的空间域中的各种细胞类型的有机结合与相互作用，能够行使特定的生物学功能。
3 结论
本文通过整合空间位置、组织学和基因表达来建模空间嵌入表示，提出了一种基于图注意力网络的深度学习方法，用于空间转录组数据的聚类分析。实验结果表明，本文提出的方法可以有效地聚类空间转录组数据，并且与其他的聚类方法相比具有更好的性能。总之，本文所提方法可以帮助研究人员更好地理解细胞在组织中的分布和表达方式，并为相关领域的研究提供新的工具和思路。
本文通过图注意力网络对空间转录组数据进行了特征提取，并通过Leiden算法进行聚类。值得注意的是，图注意力网络中的注意力是静态的，也就是说对于某个节点，它的相邻点的注意力系数是固定不变的，这一点后续会加以改进。另外，进行特征提取之后，聚类分析的方法有很多种，包括基于传统聚类算法的方法、基于自监督学习的方法、基于端对端训练的方法等。不同聚类方法的效果不同，什么样的方法能在不同的空间转录组数据上得到最好的效果，这一点还有待进一步探索。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：吴瀚文主要负责算法程序设计、数据记录与分析、论文编写；高洁主要负责项目主持、数据分析指导、论文审阅修订。
0 引言
空间转录组学是一种结合了组织学和转录组学的新兴技术，它可以研究细胞在组织中的分布和表达方式。空间转录组（spatial transcriptome，ST）数据通常是由高通量成像技术生成，如空间转录组学和空间RNA测序（spatial RNA sequencing，sRNA-seq）。这些技术可以在组织中捕获数百万个单细胞的转录组数据，并提供高分辨率的空间信息。然而，由于数据的高维和复杂性，空间转录组数据的分析一直是一个具有挑战性的问题。识别空间域（即，在基因表达和组织学上都具有空间一致性的区域）是空间转录组学中最重要的课题之一。目前，空间域识别方法主要分为非空间聚类和空间聚类两类。传统的非空间聚类算法，如K均值聚类算法（k-means clustering algorithm）、鲁汶算法（Louvain method），以基因表达数据作为输入，导致聚类难以符合组织切片[1-2]。而空间聚类方法结合了基因表达、空间位置和形态学，可以解释基因表达的空间依赖性，以更好地匹配空间位置。例如，Zhao等[3]提出的空间贝叶斯方法（BayesSpace）即是基于贝叶斯统计方法，通过引入空间邻居结构来鼓励相邻点属于同一聚类，然后将捕获区域划分为更小的亚区域，实现了聚类精度的提高。Xu等[4]提出的用于空间转录组分析的卷积神经网络（convolutional neural networks，CNN）学习策略（CNN learning strategy for spatial transcriptomics analysis，CosTA）使用CNN聚类结构，通过CNN生成特征，再通过高斯混合模型聚类生成软分配[5]，并使用软分配来更新CNN。Yuan等[6]提出的用于基因的图CNN（graph CNN for genes，GCNG）方法将空间信息编码为图形，使用监督训练将其与表达数据组合，并用图卷积网络对图数据进行特征嵌入，从而实现聚类。Hu等[7]提出的整合基因表达、空间位置和组织学，通过图卷积网络识别空间域和空间可变基因（integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network，SpaGCN）方法应用图卷积网络来整合基因表达和空间位置，并进一步与自监督模块相结合来识别空间域。Pham等[8]提出的空间转录组学数据的下游分析工具包（stLearn）提供了一种组织内归一化技术，该技术基于从形态学图像中收集到的特征和空间位置，通过形态学距离对基因表达进行归一化。Xu等[9]提出的无监督空间嵌入深度表示方法（unsupervised spatially embedded deep representation of spatial transcriptomics，SEDR）采用深度自动编码器网络和图形自动编码器来嵌入空间信息。Teng等[10]提出可用于荧光原位杂交技术的迭代细胞分型（fluorescence in situ hybridization iterative cell typing，FICT）方法，能同时利用空间转录组数据中的表达信息和空间信息来识别细胞类型和亚型。尽管这些空间聚类方法可以在一定程度上识别空间域，但往往是通过节点的空间位置计算节点之间的距离，并由此预先定义节点之间的相似性，而这些预先定义的方法没能利用深度学习的灵活性，故存在一定局限。
为了解决这个问题，本文提出了一种基于图注意力网络（graph attention networks，GAT）的深度学习方法，并将其命名为用图注意力自动编码器进行空间转录组聚类分析（spatial transcriptomics clustering analysis with graph attention auto-encoders，STCAGAE），用于空间转录组数据的聚类分析[11]。本文使用的图注意力网络可以自适应地学习节点之间的相似性，从而更好地聚类空间域。本文利用预先训练好的深度神经网络模型从形态图像块中提取特征向量，然后将提取的特征与基因表达和空间位置数据相结合，以表征空间相邻点的相关性，并创建空间增强基因表达矩阵。同时使用图注意力自动编码器和去噪自动编码器联合生成增强空间转录组数据的潜在表示，最后使用莱顿（Leiden）算法最终的嵌入进行聚类[12]。本文旨在识别不同的空间域，并为研究空间转录组数据提供新的方法。
1 方法
1.1 数据来源及预处理
本文研究所用数据来自十倍基因组学（10 × Genomics）生物技术公司利用维希姆（Visium）技术发布的高通量空间转录组测序平台（10 × Genomics Visium）。该平台可以利用Visium技术，从整个组织切片和各种各样的样本类型中获得整个转录组高通量基因表达分析结果。本文选取了该平台的一组可公开下载的数据集，即人类背外侧前额叶皮层（dorsolateral prefrontal cortex dataset，DLPFC）数据集。DLPFC数据集包含了12个不同的切片数据，其中每个切片都包含6个皮质层和1个白质层。首先，从DLPFC数据集中移除主要组织区域之外的区域。然后使用计算机编程语言Python 3.8（Guido van Rossum，荷兰）中分析单细胞数据的软件包Scanpy 1.9.1（Theis Lab，德国），根据文库大小对原始基因表达数据进行过滤、对数转换和标准化[13]。最后，使用主成分分析（principal components analysis，PCA）对增广基因表达数据进行降维，并将降维数据作为下一个模型训练的输入。
1.2 STCAGAE算法
1.2.1 增强空间数据
空间基因表达技术提供了具有额外空间位置信息和组织形态的转录组范围的基因表达谱。本文使用这两个额外的空间数据实现了节点的基因表达数据增强，得到了对节点更真实的表示，从而提高了STCAGAE的性能。计算点i和点j之间的空间基因表达权重GCij，如式(1)所示：
|
其中，Si代表节点i的基因表达向量， 代表节点i的基因表达向量的均值，Sj代表节点j的基因表达向量， 代表节点j的基因表达向量的均值。对于具有形态学信息的空间转录组数据，本文首先根据每个点的坐标对一个图像进行分割，得到其部分图像。然后，使用深度学习框架Pytorch 1.8.0（Facebook Inc，美国）中的图像预处理模块转换器（transform）对图像进行变换和增强，包括归一化、旋转、调整锐度等[14]。每个点的特征从预先训练的CNN模型得出，该模型可以将每个节点的图像转换为2 048维的潜在变量。为了更好地表示节点的形态，本文进行了PCA，提取了前50个主成分作为潜在特征。最后，利用余弦距离计算点i与相邻点j之间的形态相似性权重MSij，如式(2)所示：
|
本文使用空间坐标来确定每个点和所有其他点之间的距离，然后对前4个相邻点之间的距离进行排序，以计算半径r。对于一个给定的点i，当且仅当两个点之间的距离小于r时，认为点j是点i的相邻点，则设表示点i与点j位置关系的权重参数SWij = 1，否则SWij = 0。然后，增强每个点i的基因表达，如果是10 × Genomics Visium平台的数据，将数据增强为如式(3)所示；否则，将数据增强为如式(4)所示：
|
|
其中，GEi和GEj代表点i和点i的相邻点j的原始基因表达。
1.2.2 构造图
为了结合给定点的相邻点的相似性，本文使用空间坐标来计算点之间的距离，并使用前6个最近邻构建一个无向连接图。对于无向连接图，在每个点上都添加了自循环。A是邻接矩阵，那么如果i和j是相邻点，则在点i和点j处的值Aij = Aji = 1，否则Aij = 0。
1.2.3 去噪自动编码器
使用Pytorch 1.8.0（Facebook Inc，美国）的线性层实现了一个去噪自动编码器，用于基因表达的潜在表示[15]。编码器E由用户设置的多个完全连接的堆叠线性层组成，将集成的基因表达式X转换为低维表示，如式(5)所示：
|
其中，Zg代表经过编码器得到的低维表示，R代表实数集，N代表节点数量，M代表输入基因数量，d代表最后一个编码器层的维数。然后，向Zg加入噪声得到重构的低维表示，如式(6)所示。接下来，解码器D将重构后的低维表示转换为重构的基因表达矩阵，如式(7)所示：
| '/> |
| '/> |
其中，代表重构的低维表示，Z代表噪声向量，D代表解码器，X'代表重构的基因表达矩阵，d''代表最终层的维数，N、M和d与上述相同。输入基因和重建表达式的类似性用均方误差来计算，如式(8)所示：
|
其中，Ll代表去噪自编码器的损失函数，Xi代表节点i的基因表达向量，E（Xi）代表节点i的基因表达向量经过编码器层得到的低维表示，D（∙）代表经过解码器层重构后的表示。
1.2.4 变分图自动编码器
变分图自动编码器使用了基于Pytorch 1.8.0（Facebook Inc，美国）构建的用于编写和训练图神经网络的库Pytorch Geometric 2.2.0（Technische University Dortmund，德国） [16]，其中编码器由图注意力网络组成。图注意力网络将特征矩阵X作为输入，首先通过一个线性变换W做维度转换，如式(9)所示：
|
其中，代表特征矩阵经过线性变换之后的列数。接着使用自注意力（self-attention）为每个节点分配注意力，对于点i和点j，它们的初始特征为hi和hj，于是它们经过线性变换W后变为Whi和Whj，为它们分配注意力，如式(10)所示：
|
其中 ，eij 代表节点 i 对节点 j 的影响力系数，a (·)代表计算节点相关度的函数，如式(11)所示：
| '/> |
在本文中，对于节点i，只计算其邻域内节点j对节点i的影响力系数。为了更好地在不同节点之间分配权重，本文将目标节点与所有相邻点计算出来的相关度进行统一的归一化处理，使用线性整流函数（rectified linear unit，ReLU），如式(12)所示：
|
其中，代表归一化后的相关度权重系数，Ni代表节点i的邻域，eik代表节点 i 对节点 i邻域内的节点 k 的影响力系数。每个节点的最终输出特征向量如式(13)所示：
| '/> |
其中，hi'代表节点的最终输出特征，σ代表激活函数。损失函数包括生成图与原始图之间的重构损失，节点表示的向量分布和正态分布的相对熵（relative entropy），如式(14)所示：
|
其中，Lg代表变分图自动编码器的损失函数，代表二进制交叉熵函数，代表相对熵，，代表正态分布。
同时，本文利用多头注意力以稳定学习过程和增加模型容量。多头注意力就是将同一个输入分别送入不同的注意力层，分别计算隐藏状态，然后将不同的特征表示连接起来（或相加）。经过不断调整参数并实验，本文发现模型中的参数取8个注意力头和特征相加时，DLPFC数据集的不同切片聚类得到的调整兰德系数（adjusted Rand index，ARI）值都最大，因此模型达到最佳效果。
1.2.5 聚类与可视化
在得到嵌入之后，使用Leiden算法来识别空间域，以识别相似的细胞和模式。为了获得最佳分辨率，在0.1～2.5之间进行网格搜索，步长为0.01，直到达到必要的簇数。为了可视化，使用了统一流形逼近和投影[17]。
综上，本文提出STCAGAE算法，如图1所示。通过对整合基因表达、空间位置和组织形态信息的低维表示进行建模来表征空间域。为了建立形态学特征矩阵，首先使用预先训练好的深度学习网络对苏木精-伊红（hematoxylin-eosin，H&E）染色的组织地形图数据进行处理。结合形态特征和空间相邻点信息，每个点的基因表达均得到增强。然后，利用去噪自编码器学习从集成特征空间到低维表示空间的非线性映射，以减少模型的过拟合。同时，基于节点之间的欧氏距离计算图邻接矩阵。将增强后的基因表达矩阵和图邻接矩阵输入到一个变分图自动编码器中，从而通过与相应的空间相邻点的集成表示来生成空间嵌入。最终的潜在嵌入是通过连接集成表示和空间嵌入来形成的。最后使用Leiden算法进行空间域聚类，识别出从右上到左下依次是皮质层1～6、白质层的不同空间域。
图1
STCAGAE算法的流程示意图
Figure1.
The flow chart of STCAGAE algorithm
2 结果
本文在一个经典空间转录组数据集——DLPFC数据集上进行了实验，验证了STCAGAE算法的有效性。STCAGAE算法可以将细胞聚类成多个簇，并且簇之间具有明显的空间分布和表达差异。与传统的聚类方法相比，STCAGAE算法在聚类性能上更有优势。
本文在DLPFC数据集的编号为151507的切片上通过STCAGAE算法对空间域进行聚类。同时，比较了其他空间算法在该数据集上的聚类性能。通过ARI值计算识别空间域的精度，如图2所示，STCAGAE算法的ARI值达到了0.620，优于通过深度学习识别空间转录组学中的空间域（identifying spatial domains in spatial transcriptomics by deep learning，DeepST）方法的ARI值0.559[18-19]。
图2
DLPFC数据集151507切片上的ARI值对比
Figure2.
The comparison of ARI on section 151507 in the DLPFC dataset
通过STCAGAE算法，又在DLPFC数据集的编号为151673的切片上进行了空间域的聚类。同时也比较了其他空间算法在该数据集上的聚类性能。通过ARI值计算识别空间域的精度，如图3所示ARI值达到了0.620，优于DeepST算法的ARI值0.615。
图3
DLPFC数据集151673切片上ARI值的对比
Figure3.
The comparison of ARI on section 151507 in the DLPFC dataset
对于DLPFC数据集的151676切片，本文通过STCAGAE算法进行了空间域的聚类。通过ARI值计算识别空间域的精度，如图4所示ARI值达到了0.638。然后，本文在该数据集上通过其他空间算法进行实验，并比较各种聚类算法的性能，发现STCAGAE效果优于软件包Scanpy 1.9.1（Theis Lab实验室，德国）（ARI = 0.29）、SEDR（ARI = 0.42）、SpaGCN（ARI = 0.44）、BayesSpace（ARI = 0.44）和用于空间解析转录组学的自适应图注意力自动编码器（deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder，STAGATE）方法[20]（ARI = 0.60）。
图4
DLPFC数据集上151676切片的人工注释以及STCAGAE聚类结果
Figure4.
Human annotation and clustering result of section 151676 in the DLPFC dataset
本文在DLPFC数据集上测试STCAGAE的性能，发现STCAGAE可以有效地将数据划分为6个皮质层和1个白质层[21-22]。其中，皮质层1主要由神经胶质细胞、轴突、树突和少量神经元组成，对神经元活动的调节非常敏感，功能涉及到信息处理、信号传递和空间感知等方面；皮质层2是DLPFC的主要输入层，具有调节和控制信息传递的功能，在决策制定、工作记忆和灵活的认知控制等高级认知过程中扮演着重要的角色；皮质层3主要包含锥形神经元，在视觉加工和感知处理中扮演重要的角色；皮质层4主要包含大锥形神经元，参与计划、决策和行动控制等许多高级认知过程；皮质层5主要由大锥形神经元和星形神经元等神经元组成，对调节神经元活动起到重要的作用，同时也与情绪调节、行动控制和认知过程等方面有关；皮质层6主要由星形神经元和横向神经元等组成，在感知处理、计划和动作控制等方面具有重要的作用；白质层是DLPFC的最内层，主要由神经元轴突组成。白质对于DLPFC与其他脑区的信息交流具有重要作用，同时也与情绪调节、决策制定和计划等各种高级认知过程密切相关。与其他聚类方法相比，STCAGAE实现了显著的改善。STCAGAE清晰地勾画出了层的边界，并获得了最好的聚类精度。而通过聚类得到的空间域中的各种细胞类型的有机结合与相互作用，能够行使特定的生物学功能。
3 结论
本文通过整合空间位置、组织学和基因表达来建模空间嵌入表示，提出了一种基于图注意力网络的深度学习方法，用于空间转录组数据的聚类分析。实验结果表明，本文提出的方法可以有效地聚类空间转录组数据，并且与其他的聚类方法相比具有更好的性能。总之，本文所提方法可以帮助研究人员更好地理解细胞在组织中的分布和表达方式，并为相关领域的研究提供新的工具和思路。
本文通过图注意力网络对空间转录组数据进行了特征提取，并通过Leiden算法进行聚类。值得注意的是，图注意力网络中的注意力是静态的，也就是说对于某个节点，它的相邻点的注意力系数是固定不变的，这一点后续会加以改进。另外，进行特征提取之后，聚类分析的方法有很多种，包括基于传统聚类算法的方法、基于自监督学习的方法、基于端对端训练的方法等。不同聚类方法的效果不同，什么样的方法能在不同的空间转录组数据上得到最好的效果，这一点还有待进一步探索。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：吴瀚文主要负责算法程序设计、数据记录与分析、论文编写；高洁主要负责项目主持、数据分析指导、论文审阅修订。
