DGI算法在乳腺癌空间转录组学分析上的应用

Application of DGI algorithm in spatial transcriptomic analysis of breast cancer

摘要

为了实现在单细胞水平上量化肿瘤空间异质性，选取10×Genomics平台上的乳腺癌空间转录组数据集为研究对象，使用深度图互信息（Deep graph infomax，DGI）模型对乳腺癌细胞进行聚类研究。结果显示，DGI算法展示出较好的聚类性能，调整兰德系数达到0.55，聚类结果接近人工注释分层且边界平滑，能够出色识别出乳腺癌标记基因和簇4与簇8之间的差异表达基因，富集结果表明这些基因与乳腺癌的发生发展有非常密切的关系。分析结果可能为乳腺癌患者找到作为临床诊断和治疗依据的标志物，对乳腺癌诊断和预后产生新的见解。

Abstract

In order to quantify the spatial heterogeneity of tumors at the single cell level, the spatial transcriptomic dataset of breast cancer on the 10×Genomics platform was selected as the study object, and deep graph infomax (DGI) model was used to cluster breast cancer cells. The results showed that the DGI algorithm showed good clustering performance, and the adjusted Rand index reached 0.55. The clustering results were close to manual annotation stratification, the boundary was smooth, and the breast cancer marker genes and differentially expressed genes between cluster 4 and cluster 8 were well identified. The enrichment results showed that these genes were closely related to the occurrence and development of breast cancer. The results of this analysis may provide new insights into the diagnosis and prognosis of breast cancer by identifying markers for clinical diagnosis and treatment of breast cancer patients.

Graphical abstract

关键词

乳腺癌 / 聚类 / 标记基因 / 差异基因

Key words

breast cancer / clustering / marker gene / differential gene

目前，乳腺癌已经成为严重危害妇女健康的主要恶性肿瘤之一，且发病率呈上升趋势。在中国，乳腺癌的发病率居女性恶性肿瘤的首位^[1]。因此，乳腺癌的早期筛查、诊断和治疗已经成为人们日益关注的问题。在乳腺癌的发生、发展中，一些关键基因发挥着很大的作用，寻找癌症发生、发展的内在分子机制有着至关重要的意义。

大规模癌症基因组项目已经揭示了广泛的肿瘤间和肿瘤内异质性^[2]。近期的单细胞组学研究，特别是通过单细胞RNA测序技术，极大地促进了对肿瘤细胞异质性、肿瘤免疫浸润细胞亚群在单细胞水平上的理解^[3]。然而，单细胞RNA测序技术仍有局限性，在组织解离成单细胞悬浮液后，会导致空间和形态学信息丢失。最近开发的空间转录组学(Spa-tial transcriptomics，ST)技术可以克服上述限制，ST提供了具有完整二维位置信息的高质量全基因组转录组数据。在最近的研究中，空间结构可以作为一种信息特征来改善细胞聚类，识别细胞亚型和重要基因^[4]。BayesSpace是一种完全贝叶斯统计方法，解决了在原始分辨率下无法检测到的组织结构，并识别了组织学分析无法达到的转录异质性^[5]。stLearn是一个综合分析三种数据类型的python库，可以在空间中重建组织内的细胞类型演变，并识别具有细胞间相互作用的区域^[6]。Zehra等^[7]将60例浸润性导管癌(Invasive ductal carcinoma，IDC)病例通过自动检测获得的结果与专家病理学家报告的结果进行比较发现，Ki-67染色的自动评分具有巨大潜力，在个性化医疗时代，病理学家可以在人工智能的支持下有效地给出精确的临床诊断。乳腺癌复发网络是一个深度学习框架，它会自动预测组织病理学载玻片的Oncotype DX复发风险，研究结果为以高度置信度自动对患者进行风险分层提供了强有力的证据^[8]。

本研究采用一种基于图卷积网络的细胞聚类方法，可以从空间基因表达数据中结合单个细胞的基因表达和复杂的全局空间信息。在DGI聚类之后，进行细胞簇之间基因表达的差异分析、基因本体(Gene ontology，GO)功能富集分析和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)通路富集分析，进一步揭示癌症发生、发展的分子机制^[9]。本文聚类方法在基于10×Genomics的乳腺癌数据集上进行了测试，实验结果表明，采用DGI模型对空间转录组数据进行聚类，与其他算法相比性能较好，能够发现具有更多区域连续性和更少噪声的区域。通过对乳腺癌数据集进行详细分析，发现肿瘤内部存在更多的异质性，从而为制定靶向治疗策略提供了理论基础^[10]。

1 材料与方法

1.1 数据集

使用两个来自10×Genomics平台(Datasets-10×Genomics)的ST数据集，第一个是人乳腺癌样本，包含3 798个细胞，36 601个基因。第二个是人宫颈癌样本，用来验证DGI算法的性能，包含2 781个细胞，17 943个基因。两个数据集中包括空间点的基因表达信息和对应位置坐标，其中乳腺癌组织样本注释文件可在https://zenodo.org/record/6925603#.ZEnf7s5BxsZ下载。

1.2 数据预处理

对于空间数据的细胞聚类，首先将数据转换成一个图，其中节点代表一个以基因表达谱为属性的细胞，边代表细胞之间的邻域关系。图由两个矩阵描述：表示图结构的邻接矩阵和表示节点属性的特征矩阵。关于细胞的空间信息可以用无向图表示，其中细胞由节点表示，边连接空间上彼此接近的一对细胞。如果图中节点i和j之间存在边，则A_{\mathrm{0}ij}=\mathrm{1}，否则A_{\mathrm{0}ij}=\mathrm{0}，A_{\mathrm{0}}为初始邻接矩阵。

为了平衡单个细胞的空间信息和基因表达之间的权重，引入超参数λ将初始邻接矩阵A_{\mathrm{0}}生成混合邻接矩阵A：

A=\lambda \times I+\left(\mathrm{1}-\lambda \right)A_{\mathrm{0}}(1)

其中：I\in R^{N\times N}是一个单位矩阵，在乳腺癌ST数据集上经过大量的计算和验证分析后得出，在0到1的范围内，λ的最佳取值为0.3。

对于ST数据集，使用空间基因表达数据来提取节点特征。首先删除具有低表达或低变异性的基因，主要包括删除数据集中表达少于三个位点的基因，每个细胞中平均表达水平低于0.02的基因以及方差低于0.05的基因，然后进行表达水平标准化：

\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}{\mathrm{n}}_{ij}=\frac{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}{\mathrm{t}}_{ij}}{{\sum }_j\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}{\mathrm{t}}_{ij}}\times \mathrm{10}\mathrm{000}(2)

其中：i指的是细胞，j指的是基因，通过主成分分析将基因表达维数降至200。

1.3 节点嵌入和聚类

DGI的输入是混合邻接矩阵A\in R^{N\times N}和一组节点特征X=\left\{x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_N\right\}，其中N是节点数，x_i\in R^F代表节点i的特征，F是节点特征的个数。

DGI的目标是学习一个编码器\epsilon，该编码器将节点特征和邻接矩阵映射到一个嵌入空间：

\epsilon \left(X,A\right)=H=\left\{{\vec{h}}_{\mathrm{1}},{\vec{h}}_{\mathrm{2}},\dots ,{\vec{h}}_N\right\}(3)

其中：H为节点在嵌入空间中的低维表示，每个节点i满足{\vec{h}}_i\in R^{F\mathrm{\text{'}}}，F\mathrm{\text{'}}为训练后的特征维数，所得到的节点特征的低维表示可以用于各种下游任务，例如节点聚类任务。

编码器通过不断聚合目标节点周边的邻居来完成特征学习，由四个图卷积层组成，用于聚集相邻节点上的特征，并以一个参数校正线性单元作为激活函数，参数校正线性单位函数为

(4)

其中：a是一个随机的超参数，可以用于反向传播的训练，与其他层同时优化，并且收敛速度快、错误率低。

全局特征S是通过读出函数从局部表示映射获得

S:\vec{s}=R\left(\epsilon \left(X,A\right)\right)(5)

结合局部和全局特征，引入判别器D:R^F\times R^F\to R，计算输入节点的特征和全局特征的互信息，使输入图节点的特征与全局特征的互信息接近1，其他图节点的特征与全局特征的互信息接近0。

对于单个图，负样本的生成需要一个随机变换函数C:R^{N\times F}\times R^{N\times N}\to R^{M\times F}\times R^{M\times M}，可以表述为C:\left(\tilde{X},\tilde{A}\right)=C\left(X,A\right)，最终目标函数为

\begin{array}{l}L=\frac{\mathrm{1}}{N+M}{\displaystyle \sum _{i=\mathrm{1}}^N}{E}_{X,A}\left[\mathrm{l}\mathrm{o}\mathrm{g}D\left({\vec{h}}_i,\vec{s}\right)\right]\\ +E_{\tilde{X},\tilde{A}}\left[\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{1}-D\left(\overrightarrow{\stackrel{̃}{h_j}},\vec{s}\right)\right)\right]\end{array}(6)

DGI通过最大化{\vec{h}}_i和\vec{s}之间互信息的近似表示，输出包含图结构信息的节点嵌入。DGI的核心思想在于通过最大化局部互信息来训练编码器。DGI输出每个节点训练后的特征，用主成分分析进行降维，采用K-means++聚类算法对top主成分进行聚类识别。

1.4 差异基因的鉴定与富集分析

为了分析每个簇不同的生物学功能，进而识别细胞亚型，使用Scanpy软件包进行差异基因筛选^[11]，Scanpy与Seurat相比能很好地解决内存不足的问题，在聚类时运算速度更快。根据差异显著性t，检验p值小于0.05以及差异倍数大于2的条件筛选出在每个簇中的差异显著基因^[12]。以相应数据集中的全基因表的前200个差异基因为背景，对指定簇进行GO项富集分析和KEGG通路分析，构建功能富集图谱^[13]。

2 结果

2.1 聚类分析

图1展示了DGI算法在乳腺癌和宫颈癌空间转录组学数据上的聚类结果。

本文使用的第一个ST数据集是人乳腺癌数据集，包括10×Genomics Visium空间转录组学和人工注释层，如图1(a)所示。该数据集包括21个区域和4种主要形态类型，即导管原位癌/小叶原位癌(Ductal carcinoma in situ/lobular carcinoma in situ，DCIS/LCIS)、健康组织、IDC和肿瘤周围低恶性区（肿瘤边缘）。为了衡量聚类标签与参考标签之间的一致性，采用轮廓系数(Silhouette coefficient，SC)和戴维森堡丁指数(Davies-bouldin index，DB)两种聚类内部评价指标，调整兰德系数(Adjusted rand index，ARI)、Fowlkes-Mallows指数(Fowlkes-mallows index，FMI)，标准化互信息(Normalized mutual information，NMI)和调整互信息(Adjusted mutual information，AMI) 4种聚类外部评价指标来比较聚类算法（越高越好）的性能。采用DGI模型的聚类SC值为0.32，DB值为1.39。ARI值为0.55，FMI值为0.58，NMI值为0.64，AMI值为0.63。DGI算法聚类结果如图1(b)所示，图中0~20表示聚类出的21个域，与图1(a)中的图例一一对应，可以看出聚类效果较好，在数值上接近注释层分割，并且可以找到边界明显平滑层的层次结构。使用人宫颈癌数据集对DGI算法进行验证，聚类结果如图1(c)所示，图中0~8表示聚类出的9个域，SC值为0.37，DB值为1.13，聚类效果较好。

为了验证DGI的性能，将DGI算法与其他算法进行比较发现，如表1所示，DGI算法的ARI、FMI值均略高于其他算法，NMI值略低SEDR，展现出良好的聚类效果。

2.2 标记基因分析

图2筛选了乳腺癌的10个标记基因，包括CDH1，CD24，EPCAM，ESR1，VIM，BRCA1，CD44，ERBB2，TP53和NF1。图中横轴表示基因名称，纵轴0~20表示聚类出的21个域，与图1(a)中的图例一一对应。BRCA1、CDH1和TP53是乳腺癌易感基因，其中BRCA1的致病突变与乳腺癌的高风险相关^[14]，CDH1致病性突变与雌激素受体阳性乳腺癌的风险增加相关，主要在簇8高表达^[15]。人类CD24基因位于6q21染色体，现已发现其与多种肿瘤的发生有关，其中CD24在正常乳腺组织中表达较少，在乳腺癌中具有较高的阳性表达率，促使乳腺癌的发生与转移。CD44，CD24在乳腺癌中具有比较高的表达率，其较高的表达率可能参与了乳腺癌的发生与发展过程^[16]。由此推测，将CD44，CD24作为乳腺癌患病风险的检测标志物颇具研究前景。PIK3CA突变是与ERBB2突变共同出现的最常见的基因突变^[17]，CDH1突变是第二常见。

2.3 差异表达基因筛选

本研究主要通过进行差异基因表达分析来比较簇4 (肿瘤边缘3)与簇8 (DCIS/LICS1)，检测显著的差异表达基因(|logFC|≥2，p值<0.05)^[18-19]，图3分析了聚类之后簇4和簇8的差异基因。火山图中红色点代表乳腺癌细胞中显著过表达的基因，蓝色点代表乳腺癌细胞中显著下调的基因。在簇4中有上调基因72个，下调基因8 883个。在簇8中有上调基因1 778个，下调基因1 956个。簇4中显著高表达基因有IGFBP5，MUC1和GATA3。IGFBP5在上皮生长中起抗凋亡等作用，与多种肿瘤的形成有关，并且其在不同肿瘤之间的作用有一定的差异性^[20]。GATA3在乳腺癌中的表达与乳腺癌的分化水平、转移能力相关，并且GATA3可作为乳腺癌预后的一个判断指标^[21]。簇8中显著高表达的基因有BACE2，HK2，APOE，CXCL14和APOC1。可以看出高表达基因在乳腺癌和癌旁组织中具有显著差异。

2.4 差异基因富集分析

通过R软件对筛选出的差异基因进行GO功能分析以及KEGG通路富集分析，如图4所示。结果显示，簇4上调基因GO功能主要富集在免疫球蛋白介导的免疫反应，B细胞介导的免疫反应，基于体细胞重组免疫受体的适应性免疫反应等各种免疫反应^[22]。下调基因GO功能主要与细胞质转运、氧化性磷酸化和腺嘌呤核苷三磷酸合成有关。KEGG通路分析显示这些基因主要富集在各种疾病上，如神经变性的途径-多种疾病、糖尿病心肌病等。簇8上调基因GO功能主要富集在细胞质翻译、氧化性磷酸化、有氧呼吸和生物合成与代谢过程。下调基因GO功能主要与B细胞、白细胞介导的免疫反应有关。KEGG通路分析显示这些基因主要富集在各种疾病与化学致癌作用。

3 结论

本文采用DGI模型分析乳腺癌空间转录组学数据，将空间位置作为一种信息特征来改善细胞聚类。与其他算法相比，DGI模型对乳腺癌细胞聚类的效果更好，可以识别出乳腺癌标记基因CD44和CD24，筛选出簇4 (肿瘤边缘3)与簇8 (DCIS/LICS1) 之间的差异表达基因，表明高表达基因在乳腺癌和癌旁组织具有显著差异，分析获得的基因可能为乳腺癌特异标志物的筛选以及精准治疗提供新的方向。

参考文献

[1] 龚珂, 肖钧方, 刘香婷, 等. 乳腺癌易感基因相关研究现状[J]. 现代肿瘤医学, 2020, 28(23): 4184-4190.

[2] WU R, GUO W B, QIU X Y, et al. Comprehensive analysis of spatial architecture in primary liver cancer[J]. Science Advances, 2021, 7(51): eabg3750.

[3] GOHIL S H, IORGULESCU J B, BRAUN D A, et al. Applying high-dimensional single-cell technologies to the analysis of cancer immunotherapy[J]. Nature Reviews Clinical Oncology, 2021, 18(4): 244-256.

[4] HEYDARI A A, SINDI S S. Deep learning in spatial transcriptomics: learning from the next next-generation sequencing[J]. Biophysics Reviews, 2023, 4(1): 011306.

[5] ZHAO E, STONE M R, REN X, et al. Spatial transcriptomics at subspot resolution with BayesSpace[J]. Nature Biotechnology, 2021, 39(11): 1375-1384.

[6] PHAM D, TAN X, XU J, et al. stLearn:integrating spatial location,tissue morphology and gene expression to find cell types,cell-cell interacti ons and spatial trajectories within undissociated tissues[DB/OL]. [2024-03-20].

[7] ZEHRA T, SHAMS M, AHMAD Z, et al. Ki-67 quantification in breast cancer by digital imaging AI software and its concordance with manual method[J]. Journal of the College of Physicians and Surgeons-Pakistan: JCPSP, 2023, 33(5): 544-547.

[8] SU Z Y, NIAZI M K K, TAVOLARA T E, et al. BCR-Net: a deep learning framework to predict breast cancer recurrence from histopathology images[J]. PloS One, 2023, 18(4): e0283562.

[9] VELIČKOVIĆ P, FEDUS W, HAMILTON W L, et al. Deep graph infomax[DB/OL]. [2024-03-20].

[10] ZHENG B W, FANG L. Spatially resolved transcriptomics provide a new method for cancer research[J]. Journal of Experimental & Clinical Cancer Research, 2022, 41(1): 179.

[11] WOLF F A, ANGERER P, THEIS F J. SCANPY: large-scale single-cell gene expression data analysis[J]. Genome Biology, 2018, 19(1): 15.

[12] 付兴前, 郭礼跃, 欧梁, 等. 基于数据挖掘和网络药理学的肝肾亏虚型膝骨关节炎的用药规律及机制分析[J]. 世界中医药, 2022, 17(5): 736-742.

[13] WU T Z, HU E Q, XU S B, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data[J]. Innovation, 2021, 2(3): 100141.

[14] CHANG H J, YANG U C, LAI M Y, et al. High BRCA1 gene expression increases the risk of early distant metastasis in ER⁺ breast cancers[J]. Scientific Reports, 2022, 12(1): 77.

[15] HAN T, JIANG S L, ZHENG H, et al. Interplay between c-Src and the APC/C co-activator Cdh1 regulates mammary tumorigenesis[J]. Nature Communications, 2019, 10(1): 3716.

[16] 张佳楠. 乳腺癌患者CD44, CD44v6及CD44⁺/CD24^-表达对生存预后影响的Meta分析[D]. 苏州: 苏州大学, 2016.

[17] 申越, 贺丽, 李华, 等. 2-脱氧葡萄糖对乳腺癌MCF7/ErbB2细胞增殖迁移和侵袭的影响[J]. 西部医学, 2018, 30(11): 1564-1569.

[18] XU C, JIN X Y, WEI S R, et al. DeepST: identifying spatial domains in spatial transcriptomics by deep learning[J]. Nucleic Acids Research, 2022, 50(22): e131.

[19] ZHANG Z, ZHENG Y L, BIAN X W, et al. Identification of key genes and pathways associated with oxidative stress in periodontitis[J]. Oxidative Medicine and Cellular Longevity, 2022, 2022: 9728172.

[20] DITTMER J. Biological effects and regulation of IGFBP5 in breast cancer[J]. Frontiers in Endocrinology, 2022, 13: 983793.

[21] 吴开良, 时倩倩, 杨壹羚, 等. GATA3在乳腺癌中的研究进展[J]. 中国现代医学杂志, 2021, 31(5): 34-38.

[22] 崔笑, 廖正凯. 非小细胞肺癌中免疫相关基因的预后作用[J]. 中国医药生物技术, 2021, 16(4): 322-331.

基金资助

内蒙古自治区自然科学基金项目(2019BS03025)