基于前列腺癌肿瘤相关成纤维细胞的单细胞转录组的生物信息学分析

Bioinformatics Analysis of Single-Cell Transcriptome Based on Prostate Cancer Cancer-Associated Fibroblasts

摘要: 本研究运用生物信息学分析方法对前列腺癌(PCa)相关的单细胞数据进行了深入研究。首先，我们从GEO数据库中筛选出了1035个与肿瘤相关成纤维细胞(CAFs)相关的差异表达基因。然后，我们选取了其中最显著的Top100个基因，通过蛋白质相互作用网络分析得到了10个核心差异表达基因。进一步利用GEPIA2工具进行生存预后分析，发现MYLK、MYH11、COL18A1和CALD1可能与前列腺癌的预后显著相关。这些结果为我们深入探究前列腺癌的发病机制提供了重要的信息，有助于制定个体化的治疗策略。

Abstract: This study utilized bioinformatics analysis methods to conduct in-depth research on single-cell data related to prostate cancer (PCa). Initially, we screened 1035 differentially expressed genes related to cancer-associated fibroblasts (CAFs) from the GEO database. Subsequently, we selected the most significant Top 100 genes and identified 10 core differentially expressed genes through protein-protein interaction network analysis. Further, using GEPIA2 tool, we performed survival prognosis analysis and found that MYLK, MYH11, COL18A1, and CALD1 may be significantly associated with the prognosis of prostate cancer. These results provide important information for further exploration of the pathogenesis of prostate cancer and contribute to the development of personalized treatment strategies.

1. 引言

前列腺癌(Prostate Cancer, PCa)是男性常见的恶性肿瘤之一 [1] [2] [3] [4] ，在全球范围内造成了广泛的关注和研究。尽管已经取得了一定的进展，但前列腺癌的发病机制和治疗方法仍然面临着挑战。近年来，单细胞RNA测序技术的发展为揭示肿瘤细胞异质性和细胞间相互作用提供了全新的视角 [5] [6] [7] 。

在前列腺癌组织中，除了肿瘤细胞外，还存在着多种其他细胞类型，如肿瘤相关成纤维细胞(CAFs) [8] [9] [10] [11] 。多种恶性肿瘤组织如前列腺癌、肺癌、乳腺癌、胃癌、结直肠癌和胰腺癌等均含有大量CAFs，尤其在乳腺癌和胰腺癌中CAFs含量相当高，甚至高达肿瘤质量的80%。CAFs为肿瘤微环境中的重要成分，对于肿瘤的发生、发展和转移具有重要影响。然而，迄今为止，关于前列腺癌中CAFs的分子特征及其与肿瘤细胞之间相互作用的了解仍然有限 [12] [13] [14] 。

因此，本研究旨在通过单细胞RNA测序，对前列腺癌组织中的CAFs进行全面的生物信息学分析，以深入探究CAFs在前列腺癌发展中的作用和机制。我们将重点关注CAFs的基因表达谱、信号通路激活情况、蛋白互作网络以及关键基因对生存的影响。通过这一研究，我们期望能够为前列腺癌的早期诊断、治疗靶点的发现以及个体化治疗策略的制定提供新的理论和实践基础。

2. 材料与方法

2.1. 单细胞转录组序列数据

为了获取前列腺癌组织中成纤维细胞的基因表达数据，我们利用了公共基因数据库Gene Expression Omnibus (GEO)中的多个相关数据集，包括GSE137829、GSE117403和GSE141445 [15] [16] [17] 。这些数据集提供了来自不同来源的前列腺癌和正常前列腺组织样本的基因表达信息，为我们进行生物信息学分析提供了重要的数据支持。具体来说，GSE137829数据集包含了6个前列腺癌样本，GSE117403包含了3个正常前列腺样本，而GSE141445则包含了13个前列腺癌样本。这些样本的丰富性使得我们能够进行对比分析，从而深入了解不同细胞类型的基因表达模式以及细胞群的特征。

2.2. PCA、tSNE和UMAP分析

我们使用了R软件中的Seurat V4工具包 [18] [19] [20] 。我们进行了PCA分析，该分析通过计算主成分之间的方差来确定数据中的主要变化方向。PCA将高维度的数据转换为低维度空间，同时保留了数据中最重要的变化信息。这使得我们可以更好地区分不同细胞类型，并在更低维度的空间中表示数据。我们运用tSNE和UMAP算法对经过PCA降维后的数据进行进一步的降维和聚类分析。UMAP是一种非线性降维算法，能够有效地捕捉数据中的非线性结构，并将其映射到二维或三维空间中。这样，我们可以更直观地观察不同细胞类型之间的关系和分布模式。UMAP的优势在于可以在保留数据结构的同时，减少数据的维度，使得数据更易于可视化和解释 [21] [22] [23] [24] 。

2.3. 细胞类型注释

我们利用marker基因对聚类结果进行细胞类型注释。通过对已知的细胞类型特异性marker基因进行分析，我们能够将每个聚类识别为特定的细胞类型 [25] [26] [27] 。marker基因是与特定细胞类型或状态相关联的基因，在特定细胞类型中表达水平显著高于其他细胞类型。通过识别和分析这些marker基因，我们可以确定每个聚类所代表的细胞类型，从而对单细胞RNA测序数据进行更精确的注释和解读。

2.4. 差异基因

在完成UMAP分析后，我们使用了|avg_log2FC| > 0.25这个标准，其中avg_log2FC代表每个基因在不同聚类之间的平均对数折叠变化 [28] [29] [30] [31] 。这个标准帮助我们筛选出在不同细胞类型之间表达差异显著的基因，因为大于0.25的对数折叠变化意味着基因的表达在不同细胞类型之间有明显的差异。我们使用了P < 0.05的显著性水平作为筛选标准，这意味着我们只考虑在统计上显著的差异基因。P值是用来衡量差异是否显著的统计指标，较小的P值表示差异更显著。因此，我们选取了P < 0.05的显著性水平，以确保我们识别的差异基因具有统计学意义。

2.5. GO和KEGG通路富集分析

使用R软件包首先，我们利用Clusterprofiler进行基因本体(gene ontology，GO)分析 [5] [32] [33] [34] 。基因本体分析是一种用于揭示差异基因生物学功能和调控机制的方法。通过GO分析，我们可以将差异基因分为三个主要的功能分类：分子功能(Molecular Function)、细胞组成(Cellular Component)和生物过程(Biological Process)。这样的分类有助于我们理解差异基因在细胞中的功能作用和相互关系。我们利用Clusterprofiler进行京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG)通路分析。KEGG通路分析可以帮助我们确定差异基因富集的信号通路，从而揭示这些基因在细胞信号传导、代谢途径和疾病发生中的重要作用。这种分析能够帮助我们深入理解差异基因的功能调控网络，为进一步的机制研究和临床应用提供重要的指导。

2.6. 蛋白互作网络分析

使用在线STRING数据库(https://string-db.org/)对差异基因进行蛋白互作(Protein-protein interaction, PPI)网络分析是一种常用的方法 [32] [34] [35] [36] ，可以帮助我们揭示差异基因在细胞中的相互作用关系，进而理解其在生物学过程中的功能和调控机制。在进行PPI网络分析时，我们设置了相互作用的最高置信度阈值为大于0.7。这个阈值可以确保我们筛选出高可信度的蛋白互作关系，从而得到更可靠的网络信息。我们利用Cytoscape软件(版本3.10.1)对蛋白互作网络进行可视化分析。Cytoscape是一种用于网络分

析和可视化的强大工具，能够帮助我们更直观地理解PPI网络的结构和特征。通过在Cytoscape中加载STRING数据库中预测的PPI网络数据，我们可以进行节点布局、网络分析和功能注释等操作，从而更全面地理解差异基因之间的关联关系和功能调控网络。我们利用Cytoscape中的CytoHubba插件进行最大克里尔中心度和度的分析，以识别在PPI网络中具有关键作用的基因。最大克里尔中心度和度是衡量节点在网络中重要性的指标，通过识别这些关键基因，我们可以进一步深入探究其在细胞信号传导、代谢途径和疾病发生中的重要作用，为后续的功能研究和临床转化提供重要线索和指导。

2.7. 生存分析

根据TCGA-PRAD的表达信息和临床信息，我们使用了基因表达谱动态分析(Gene Expression Profiling Interactive Analysis, GEPIA)数据库(http://gepia.cancer-pku.cn/)来收集RNA表达谱数据 [37] [38] [39] [40] 。我们从GEPIA数据库中获取了TCGA-PRAD患者的RNA表达数据，并针对感兴趣的基因进行表达水平的分析。为了将患者样本分为高表达组和低表达组，我们通常会以中位数表达作为表达阈值。这种方法可以帮助我们将患者样本分成表达水平高和低的两组，以便后续的生存分析。我们利用Kaplan-Meier生存曲线来评估总无病生存期(Disease-Free Survival, DFS)的差异。生存分析会计算高表达组和低表达组的生存曲线，并计算风险比(Hazard Ratio, HR)、95%置信区间(Confidence Interval, CI)以及log-rank检验的p值。这些统计指标能够帮助我们评估基因表达水平与患者生存情况之间的关联性，确定基因在前列腺癌患者中的生物学和临床意义。

3. 结果

3.1. PCA、tSNE和UMAP分析

对数据进行基本的过滤时，我们设置了最小基因数为500，最多基因数为7000，以确保数据的质量和可信度。接下来，我们采用了harmony方法将所有前列腺癌(PCa)单细胞水平数据整合，通过PCA分析得到一系列主成分(Principal Component, PC)，如图1(a)所示。图1(a)，图1(b)展示了20个基因的PCA和热图分析结果，呈现了基因表达的变化趋势。使用FindClusters功能对细胞进行聚类时，我们设置了参数resolution = 0.5，以获得合适的聚类结果。接着，利用RuntSNE功能进行tSNE聚类，参数为dims = 1:20，最终得到了31个细胞簇，如图1(c)所示。这些细胞簇反映了PCa样本中不同细胞类型或亚群的分布和特征。

3.2. 细胞类型注释和差异基因分析

为了进一步理解这些细胞簇的性质和功能，我们通过marker基因对多细胞类型进行了注释，结果如图1(d)所示。marker基因的注释帮助我们确定了每个细胞簇所代表的细胞类型或状态，从而深入探索前列腺癌(PCa)组织中细胞类型的分布和相互关系。根据以往文献报道中的细胞标记基因对这31个细胞簇进行细胞类型的鉴定，见图1(d)。区分各细胞类型的典型标记基因在31个细胞簇中的表达情况见图2。经过鉴定，该数据集得到的细胞类型包括肿瘤相关成纤维细胞(Cancer-Associated Fibroblasts, CAFs)、内皮细胞(Endothelial Cells)、上皮细胞(Epithelial Cells)、髓细胞(Myeloid Cells)、增殖细胞(Proliferating Cells)、平滑肌细胞(Smooth Muscle Cells)、T细胞(T Cells)。我们单独提取了每种细胞类型的差异基因，并使用FindAllMarkers功能对PCa进行了差异分析。FindAllMarkers函数的参数为默认参数，得到了差异基因的列表，并对P值进行筛选。其中，CAFs中包含了1035个差异基因，这些差异基因可能在肿瘤相关成纤维细胞中发挥重要的生物学功能，值得进一步研究和探索。

Figure 1. (a): the first 20 genes of PC1, PC2, PC3 and PC4; (b): Heat map of the expression of the first 20 genes of PC1, PC2, PC3 and PC4; (c): tSNE of 31 clustering groups; (d): Cell type identification results

图1. (a): PC1、PC2、PC3和PC4的前20个基因；(b)：PC1、PC2、PC3和PC4的前20个基因的表达情况热图；(c)：tSNE聚类得到31个细胞簇；(d)：细胞类型注释结果

Figure 2. Violin plot of canonical marker genes used for distinguishing cell types

图2. 区分各细胞类型的典型标记基因的小提琴图

3.3. GO和KEGG通路富集分析

将单细胞转录组注释中CAFs中的Top100个显著差异表达基因。将这些基因进行GO和KEGG通路富集分析。GO富集柱状图展示了BP、CC、MF排名前6的GO术语，富集结果涉及到生物学过程(biological process, BP)、细胞组分(cellular component, CC)和分子功能(molecular function, MF)，如图3(a)所示。差异基因主要的生物学过程主要包括Ameboidal-型细胞迁移(ameboidal-type cell migration)、细胞外基质组织(extracellular matrix organization)。差异基因主要的细胞组分主要包括胶原蛋白-含细胞外基质(collagen-containing extracellular matrix)、内质网膜泡(endoplasmic reticulum lumen)。差异基因主要的分子功能主要包括细胞外基质结构组成(extracellular matrix structural constituent)等。对差异基因进行KEGG通路分析，如图3(b)所示，结果表明它们主要富集在肌细胞中的细胞骨架(Cytoskeleton in muscle cells)、血管平滑肌收缩(Vascular smooth muscle contraction)等。

Figure 3. (a): Results of GO enrichment process of differentially expressed genes; (b): Differentially expressed gene pathway annotation results

图3. (a)：差异表达基因GO富集过程结果；(b)：差异表达基因通路注释结果

3.4. 差异表达基因蛋白相互作用网络的构建与分析

我们将100个差异表达基因导入在线分析数据库STRING，并设定置信度0.7作为判断相互作用是否有意义的标准，同时去除网络中无连接的节点，构建了PPI网络。接着，我们将STRING构建的PPI网络导入Cytoscape软件进行了可视化分析，该网络包含了100个节点和87条边(见图4(a))。使用Cytoscape中的CytoHubba插件对PPI网络进行分析，我们发现了最显著的相互作用模块。这个模块网络由MCC得分最高的前10个基因组成，并用从黄色到红色渐变的颜色表示，颜色越红代表MCC分数越高，更倾向于是关键蛋白(见图4(b))。这10个基因分别是ACTA2、MYL9、CALD1、TPM1、TPM2、MYH11、TAGLN、MYLK、COL18A1、COL6A2。这些基因在PPI网络中的关联性和重要性被MCC分数所体现，MCC分数高的基因更可能在网络中扮演关键的角色，可能具有重要的生物学功能和调控作用。

Figure 4. Protein-protein interaction (PPI) networks in 100 genes and the selection of 10 hub genes

图4. 100个基因的PPI图和10个关键基因图

3.5. TCGA数据库中验证枢纽基因的临床意义

在TCGA数据库中选择了PRAD病例，并从中筛选出了关键基因，用于无病存活率(disease-free survival, DFS)生存分析。接着，利用在线工具GEPIA2对显著模块中的10个差异表达基因进行了预后价值分析。结果显示，其中4个基因MYLK、MYH11、COL18A1和CALD1的表达情况与前列腺癌患者的DFS显著相关(见图5)。这些分析结果为我们提供了有关前列腺癌的生存预后和潜在的关键基因的重要信息。

Figure 5. Survival analysis results

图5. 生存分析结果

4. 讨论

本研究通过生物信息学方法对前列腺癌(PCa)相关的单细胞数据进行了深入研究，并取得了一系列有意义的结果。首先，在数据预处理阶段，我们通过设置最小基因数和最多基因数的过滤条件，确保了数据的质量和可信度。接着，我们利用PCA、tSNE和UMAP等方法对数据进行了降维和聚类分析，成功地识别出了31个细胞簇，反映了PCa样本中不同细胞类型的分布和特征。

在细胞类型注释和差异基因分析方面，我们通过marker基因对细胞簇进行了注释，确定了每个细胞簇所代表的细胞类型。经过鉴定，我们发现数据集包括肿瘤相关成纤维细胞(CAFs)、内皮细胞、上皮细胞、髓细胞、增殖细胞、平滑肌细胞、T细胞等多种细胞类型。这些结果为我们深入探索PCa组织中细胞类型的分布和相互关系提供了重要的基础。

在进一步的分析中，我们单独提取了CAFs中的差异基因，并进行了GO和KEGG通路富集分析。结果显示，这些差异基因主要涉及到细胞迁移、细胞外基质组织和细胞骨架等生物学过程和通路，揭示了CAFs在PCa发展中的潜在作用和调控机制。

此外，我们构建了差异表达基因蛋白相互作用网络，并通过网络分析发现了最显著的相互作用模块。这些模块中的关键基因可能在PCa中扮演关键的调控角色，具有重要的生物学功能。进一步地，我们验证了其中一些关键基因在TCGA数据库中的临床意义，发现MYLK、MYH11、COL18A1和CALD1等基因表达情况与PCa患者的预后显著相关。

综上所述，本研究通过细胞分型、差异基因分析、通路富集和蛋白相互作用网络等多个层面，深入探究了PCa组织中细胞类型的特征和关键调控因子，并为该疾病的治疗和预后评估提供了重要的理论依据和研究方向。

基金项目

吉林省科技厅青年成长科技项目(批准号：20230508165RC)和吉林省自然科学基金资助(批准号：20240101200JC)。

NOTES

^*通讯作者。