对参与肺腺癌发生的核心基因的综合生信分析

Integrative Bioinformatics Analysis Identified Hub Genes in Association with Development of Lung Adenocarcinoma

摘要: 目的：侵袭性的肺腺癌(LUAD)是引起肺癌死亡的主要原因之一。因此，鉴定重要的LUAD相关基因及进一步分析其预后意义对于LUAD患者的生存率至关重要。方法：采用加权基因共表达网络分析(WGCNA)和差异基因表达分析方法，从TCGA-LUAD数据库和GEO的GSE32863筛选出有差异的共表达基因，对其进行功能富集分析和蛋白质相互作用网络(PPI)分析。此外，通过应用Cytoscape的CytoHubba插件来识别12个核心基因进行生存分析和肿瘤分期相关性的分析。结果：从TCGA和GEO数据库中共提取了358个差异共表达基因。这些基因在GO分析中主要富集于细胞外结构组织，细胞–细胞连接和DNA结合转录酶激活活性。在KEGG分析中，主要富集于药物代谢–细胞色素P450。此外，在PPI网络中鉴定了12个核心基因。在LUAD患者中，ADCY4、VIPR1和TGFBR2的表达水平与临床分期和整体生存率(OS)相对应。结论：ADCY4、VIPR1和TGFBR2可能在肿瘤发生中起重要作用，因此它们可作为LUAD的预后生物标志物和治疗靶点。

Abstract: Objective: Lung cancer-related death is mainly caused by lung adenocarcinoma (LUAD), an aggressive malignant tumor. Therefore, the identification of important LUAD-related genes and the further analysis of its prognostic significance are critical for the survival of LUAD patients. Method: Weighted Gene Co-expression Network Analysis (WGCNA) and differential gene expression analysis methods were adopted to screen out TCGA-LUAD database and the gene expression profiles of GSE32863 from GEO. Functional annotation analysis and protein-protein interaction (PPI) network were conducted on differential co-expression genes. Furthermore, survival analysis was carried out on twelve hub genes that were identified by applying the CytoHubba plugin of Cytoscape. Results: A total of 358 differential co-expression genes were extracted from the database of TCGA and GEO. These genes were mainly enriched in extracellular structure organization, cell-cell junction and DNA-binding transcription activator activity. In the KEGG analysis, the main pathways were Drug metabolism-cytochrome P450. Moreover, in a PPI network, the 12 hub genes were identified. The expression level of ADCY4, VIPR1, and TGFBR2 was corresponded with clinical stages and overall survival (OS) in LUAD patients. Conclusion: ADCY4, VIPR1 and TGFBR2 may play an important role in the mechanism of the tumorigenesis, so they will serve as prognostic biomarkers and therapeutic targets of LUAD in the future.

1. 引言

肺癌作为全球发病率最高的恶性肿瘤，具有较高的死亡率 [1]。肺癌可分为非小细胞肺癌(non-small cell lung carcinoma, NSCLC)和小细胞肺癌类型。非小细胞肺癌的主要亚型是肺腺癌(lung adenocarcinoma, LUAD)、肺鳞状细胞癌和大细胞癌 [2]。最近的研究表明，全球LUAD的发病率高于肺鳞癌 [3]。分子靶向治疗和免疫治疗使LUAD患者的生存率逐渐提高。例如，酪氨酸激酶抑制剂被用作晚期LUAD患者表皮生长因子受体基因突变的一线临床治疗方案 [4]。此外，肺癌免疫治疗逆转了免疫检查点，主要包括程序性死亡蛋白-1和程序性死亡配体-1，在特定肺癌患者中观察到明显的治疗效果 [5]。然而，LUAD患者的预后仍然较差，5年生存率仍处于较低水平 [6] [7]。由于缺乏对LUAD具体分子机制的研究，更多靶向治疗路径并未被发现。因此，准确的预测LUAD的靶向指标被认为具有显著改善LUAD预后和治疗的能力。

随着基因组技术的发展，利用生物信息学的基因表达谱分析来探索癌症的分子机制和发现准确的生物标志物变成了一种主流方式 [8]。可以采用加权基因共表达网络分析(Weighted Gene Co-expression Network Analysis, WGCNA)的系统生物学算法将基因与临床性状进行关联分析，并识别这些基因的潜在生物学功能 [9]。WGCNA的重要作用是构建高度相关基因的共表达模块，并评估相关模块与临床性状之间的关联 [10]。此外，我们还利用差异基因表达分析来探讨基因组调控的分子机制，并发现肿瘤与正常组织之间表达水平的变化 [11]。差异基因表达分析有助于发现潜在的癌症相关生物标志物。因此，将WGCNA和差异基因表达分析得到的结果结合起来可用于筛选出肿瘤高度相关基因。

本研究采用WGCNA和差异基因表达分析方法，对TCGA和GEO数据库中的LUAD相关mRNA表达谱进行了分析，以区分不同的共表达基因。通过蛋白质–蛋白质相互作用(protein-protein interaction, PPI)分析、功能富集分析和生存分析，进一步探讨了LUAD的发展。我们的研究提供了一种新的方法来了解LUAD的潜在分子机制与差异共表达基因的分析。

2. 材料和方法

2.1. 数据采集

LUAD的基因表达数据来源于TCGA数据库和GEO数据库。使用TCGA biolinks包 [12] 下载的TCGA数据中包含535个LUAD样本和44个正常组织的信息。在研究中，所有基因保存在cpm (count per million) ≥ 1。我们在edgeR包 [13] 中使用函数rpkm筛选总共15，127个具有RPKM值的基因进行下一次分析。从Illumina HumanWG-6 v3.0 expression beadchip (GPL6884)下载GSE32863的基因表达谱。GSE32863包括来自LUAD患者的58个肺腺癌和58个邻近的非肿瘤肺组织。使用Perl语言的merge脚本(https://www.perl.org/)将RNA-seq数据文件和临床数据文件分别合并到一个矩阵文件中。

2.2. 共表达网络中功能模块的识别

利用R中的WGCNA包构建TCGA-LUAD和GSE32863数据的共表达网络 [9]。对于与外部样本性状相关的模块，WGCNA对样本间高度相关基因的模块进行了分析。为了建立一个大规模的网络，根据pickSoftThreshold函数，我们选择了soft powers β = 3 and20。然后，使用以下公式：a_ij = |S_ij| β (a_ij：基因i与基因j之间的邻接矩阵，S_ij：Pearson相关性分析计算的相似性矩阵，β：soft power值)创建邻接矩阵。对拓扑重叠矩阵(TOM)和相应的异种(1-TOM)进行了变换。然后，通过构建基于1-TOM矩阵的基因层次聚类，将相似的表达基因分为不同的基因共表达模块。根据先前的研究，我们计算了模块–特性关联来识别功能模块 [14]。因此，高相关系数模块与临床性状相关，并被筛选出用于后续分析。在前面的研究中对这些方法有过具体的描述 [14]。

2.3. 差异表达基因(DEGS)筛选和与功能模块的相互作用

利用limma包 [15] 对RNA测序和微阵列数据进行差异表达分析。将limma包分别应用于TCGA-LUAD和GSE32863数据集中，筛选出LUAD与正常组织之间差异表达的基因(DEGS)。为了控制假发现率，采用Benjamini-Hochberg法对p值进行了调整。识别DEGS的阈值设置在|logFC| ≥ 1.0和adj. p < 0.05。最后，为了识别潜在的预后基因，我们使用VennDiagra包构建DEGS与共表达基因之间的重叠基因。

2.4. 潜在预后基因的富集分析

GO分析包括三个领域：分子功能(MFs)、细胞成分(CCs)和生物过程(BPS)。利用KEGG分析信号通路中异常表达的交叉基因的潜在功能。ClusterProfiler包 [16] 用于探索潜在预后基因的功能(筛选标准为adjusted p < 0.05)。

2.5. PPI的构建及核心基因的筛选

利用STRING在线工具(https://string-db.org)构建了这些候选中枢基因的蛋白质–蛋白质网络(PPI) [17]。在本研究中，我们使用Cytoscape (v3.7.2)将综合分数 ≥ 0.4作为截止值来建立可视化的网络模型 [18]。由CytoHubba (Cytoscape中的一个插件)计算Degree algorithm后，在共表达网络中查找集线器节点 [19]。随后，我们选择Degree algorithm位于前12 (连接度为 ≥ 6)的基因作为核心基因。

2.6. 核心基因的临床相关性和预后价值

利用TCGA数据集的临床信息，我们分析了在癌症分期中核心基因的相对表达。在TCGA数据库中运用survival包中的Kaplan-Meier单变量生存分析对核心基因进行生存分析。

3. 结果

3.1. 共表达基因模块的鉴定

我们使用R上的WGCNA包从TCGA-LUAD和GSE32863数据集构建基因共表达网络。总共确定了TCGA-LUAD中的15个模块和GSE32863中的7个模块，并排除了一个灰色模块，因为它没有被分配到集群中。然后绘制模块–性状关系热图，评价两个临床性状与每个模块之间的相关性。结果表明，TCGA-LUAD的蓝绿色模块和GSE32863的棕色模块与正常肺组织的相关性最高(蓝绿色模块：r = 0.79，p = 5e−127，棕色模块：r = 0.97，p = 9e−71)。此外，与正常组织相比，蓝绿色模块和棕色模块中的基因在LUAD组织中的表达有所下调。

3.2. DEG和共表达网络间重叠基因的提取

我们通过limma包(cut-off criteria of |logFC| ≥ 1.0, adj. p < 0.05)发现TCGA-LUAD数据集中的3583个DEGS和GSE32863数据集中的957个DEGS在肿瘤组织中表达失调。在TCGA-LUAD数据集中的蓝绿色模块和GSE32863中的棕色模块中，分别有5313个和1699个共表达基因。最终，我们提取了两个模块中重叠的358个基因来验证共表达网络。

3.3. 功能富集分析

为了进一步探讨358个重叠基因的潜在功能，采用cluster profiler包进行GO和KEGG富集分析。生物过程(BP)的结果表明，这些基因主要富集于细胞外结构组织、细胞–底物粘附和对酸性化学物质的反应中。这些基因在细胞成分(CC)中主要参与细胞–细胞连接和细胞外基质胶原等部分。最重要的是，分子功能(MF)分析证明，这些基因与DNA结合转录激活酶活性，糖胺聚糖结合和生长因子结合相关。在KEGG分析中，这些基因主要作用在药物代谢–细胞色素P450、醚脂代谢和花生四烯酸代谢等途径。

3.4. PPI网络的构建及轮毂基因的鉴定

我们应用STRING数据库在重叠基因之间建立PPI网络，并从PPI网络中提取的具有连接度排序的集线器基因的信息。最后，选择GNG11、ADRB2、ADCY4、TGFBR2、IL6、GPC3、VIPR1、GRK5、CAV1、RAMP3、RAMP2和CALCRL等12个最高分基因作为核心基因。

3.5. 核心基因表达在LUAD患者中的临床相关性

在筛选出12个枢纽基因后，基于TCGA中的数据集，我们通过在R中应用cliCor包验证了ADCY4的表达与肿瘤分期有很强的正相关关系(p < 0.05；图1(a))。在VIPR1的表达中I期与IV期有显着性差异(p = 0.033，图1(e))。最后，I期和III期TGFBR2的表达有显着性差异(p = 0.024，图1(i))。

3.6. 12个核心基因的预后值和蛋白质表达

我们应用R生存包的Kaplan-Meier分析法对12个中枢基因进行生存分析(图2)。Kaplan-Meier分析表明，ADCY4、VIPR1和TGFBR2的表达水平较低与LUAD患者总体生存率(overallsurvival, OS)较差有关(ADCY4: p = 0.003, VIPR1: p = 0.019, TGFBR2: p = 0.002) (图2(a)、图2(e)和图2(i))。

Figure 1. Twelve central genes expressed in different tumor stages based on TCGA-LUAD database

图1. TCGA-LUAD数据库中十二个中枢基因在不同肿瘤分期的表达

Figure 2. Analysis of the overall survival (OS) for the expression of 12 hub genes based on TCGA database

图2. TCGA数据库对12个中枢基因表达的总体生存(OS)分析

4. 讨论

虽然分子靶向治疗和免疫治疗已逐渐应用于LUAD患者，但LUAD的预后仍然较差。不同的癌症亚型呈现不同的分子靶点。因此，探索更好的生物标志物对LUAD的具体预后和进展至关重要。在本研究中，通过对与LUAD相关的DEGS和共表达基因的综合生物信息学分析，我们最终在TCGA和GSE32863数据库中获得了358个具有相同表达趋势的基因。GO和KEGG通路富集分析表明，差异共表达基因主要富集于与细胞增殖、侵袭和迁移有关的通路中。此外，根据Cytoscape中CytoHubba插件的连接度排名，我们筛选出了前12个LUAD相关基因，包括GNG11、ADRB2、ADCY4、TGFBR2、IL6、GPC3、VIPR1、GRK5、CAV1、RAMP3、RAMP2和CALCRL。其中ADCY4、VIPR1和TGFBR2的表达与LUAD的临床分期有关。同时，ADCY4、VIPR1和TGFBR2的下调表达与LUAD中的不良OS密切相关。最后，我们对ADCY4、VIPR1和TGFBR2进行了生存和免疫组织化学分析。

特异性腺苷酸环化酶(adenylyl cyclases，ACs)可激活cAMP信号通路。由于特定AC亚型的组织分布不同，cAMP可以选择性地调节效应蛋白。因此，ACs往往能整合细胞信号，并被视为这些信号通路中的关键酶 [20]。例如，ACs中的ADCY高表达可导致白血病细胞的生长、刺激细胞凋亡和c-Myc表达增加 [21]。此外，一项研究表明，ADCY3作为胃癌的癌基因，可以被DNA甲基化调控 [22] [23]。同时发现在乳腺癌中，ADCY4与G蛋白偶联受体和cAMP信号通路密切相关。且ADCY4的表达与乳腺癌患者的OS和肿瘤分期有关 [24]。最近的研究报道，VIPR1可以抑制多种癌症的发展，包括髓母细胞瘤 [25]，前列腺癌 [26] 和淋巴母细胞瘤 [27]，提示VIPR1在抑制癌细胞生长发育中起着重要作用。在肝细胞癌中，VIPR1的DNA甲基化显著抑制基因转录，VIPR1的低表达与较差的预后和组织学分化有关 [28]。据报道，与正常组织相比，VIPR1在肺腺癌组织中表达下调 [29]。肺癌细胞H1299中VIPR1的表达低于正常肺上皮细胞BEAS-2B。此外，VIPR1对H1299细胞的生长、迁移和侵袭也有明显的抑制作用 [29]。TGFBR2可抑制各种癌症的恶性进展，如胶质母细胞瘤 [30]，宫颈癌 [31]，膀胱癌 [32] 等。在TGF-β信号通路中，TGF-β配体与TGF-BR2的结合促进了TGFBR1的磷酸化和激活 [33]。随后，活化的TGFBR1激活Smad蛋白，可抑制癌症中的EMT [33] [34] [35]。Li的研究表明，TGFBR2的低表达与肺癌细胞(A549)增殖、迁移和侵袭能力的上调有关 [36]。

在我们的研究中也存在一些局限。虽然LUAD与正常组织之间的潜在诊断基因是通过综合生物信息学分析来确定的，但对于每个LUAD患者来说，很难做到非常准确。此外，ADCY4、VIPR1和TGFBR2对LUAD患者预后产生影响的具体分子机制应通过基础实验加以验证。

5. 结论

本研究旨在挖掘与LUAD相关的生存相关基因，从而进一步开发一种新的靶向基因治疗策略。通过对与LUAD相关的DEGS和共表达基因的综合生物信息学分析，我们最终获得了358个具有相同表达趋势的基因。进一步的分析表明，包括ADCY4、VIPR1和TGFBR2在内的三个核心基因与肿瘤的临床分期和整体生存率显著相关。因此，ADCY4、VIPR1和TGFBR2可能在肿瘤发生中起重要作用，并作为LUAD的预后生物标志物和治疗靶点。