肺腺癌动态5节点ceRNA调控网络预测及机理分析

Predictive and Mechanistic Analysis of Dynamic 5-Node ceRNA Regulatory Network in Lung Adenocarcinoma

摘要: 目的：构建肺腺癌发生发展过程中的动态五节点ceRNA调控网络，挖掘核心基因，为肺腺癌诊断及预后提供新思路。方法：从TCGA及GEO数据库获得肺腺癌mRNA、lncRNA、miRNA、circRNA、TF表达数据，将患者样本根据临床分期分为癌旁样本、早期样本(stage I期)、晚期样本(stage II、III、IV期)，并将癌旁与早期、早期与晚期分别进行差异分析，将两组差异结果取交集，基于ChipBase、HOCOMOCO V11、AnimalTFDB、GTRD、TransmiR、TRRUST、CircBank、Starbase、miR2Disease、miRecords、miRTarBase和TarBase、LncBase、LncLocator数据库获得调控关系对，构建五节点ceRNA调控网络，对网络中的靶基因进行GO富集以及构建PPI网络挖掘核心基因。结果：构建了随分期动态变化的LUAD 5节点ceRNA调控网络，网络中的靶基因主要富集在脂肪酸代谢和突触成熟等生物过程中，最后获得与肺腺癌发生发展有关的8个核心基因NEFL、RBP4、FGA、SLC2A1、ALB、AFP、SLC7A5、DKK1。结论：调控网络中靶基因富集的相关通路以及8个核心基因NEFL、RBP4、FGA、SLC2A1、ALB、AFP、SLC7A5、DKK1为肺腺癌发生发展过程的机制分析、诊断及预后提供新思路。

Abstract: Objective: To construct a dynamic 5-node ceRNA regulatory network during the development of lung adenocarcinoma, to mine the core genes, and to provide new ideas for lung adenocarcinoma diagnosis and prognosis. Methods: Lung adenocarcinoma mRNA, lncRNA, miRNA, circRNA and TF expression data were obtained from TCGA and GEO databases, and the patient samples were divided into paraneoplastic samples, early samples (stage I), and advanced samples (stage II, III, IV) according to the clinical staging and the differences between the paraneoplastic and the early, early and the late stage were analyzed separately, and the differences between the two groups were analyzed. The results were taken as intersection, based on ChipBase, HOCOMOCO V11, AnimalTFDB, GTRD, TransmiR, TRRUST, CircBank, Starbase, miR2Disease, miRecords, miRTarBase, and TarBase, LncBase, LncLocator database to obtain regulatory pairs, construct 5-node ceRNA regulatory network, GO enrichment of target genes in the network as well as construction of PPI network to mine core genes. Results: A LUAD 5-node ceRNA regulatory network that changes dynamically with staging was constructed, and the target genes in the network were mainly enriched in biological processes such as fatty acid metabolism and synaptic maturation, and finally eight core genes related to lung adenocarcinoma development were obtained, NEFL, RBP4, FGA, SLC2A1, ALB, AFP, SLC7A5, and DKK1. Conclusion: The pathways involved in the enrichment of target genes in the regulatory network, as well as the eight core genes NEFL, RBP4, FGA, SLC2A1, ALB, AFP, SLC7A5, and DKK1, provide new insights for the mechanism analysis, diagnosis, and prognosis of lung adenocarcinoma.

1. 引言

肺癌是世界范围内致死率居于首位的癌症 [1] ，大约85%的癌症病例被诊断为非小细胞肺癌，其中肺腺癌(lung adenocarcinoma, LUAD)是最常见的类型 [2] 。尽管针对LUAD的早期诊断和治疗有了很大改善，据美国肺协会(American Lung Association)的2023年《肺癌状况》报告揭示，国际范围内整体肺癌五年生存率从2015年到2019年上增加了22% [3] 。但由于肿瘤的快速进展转移，患者的长期预后仍然较差，LUAD患者的5年总生存率仅有5%~18% [4] 。然而，LUAD复杂的分子生物学机制尚未完全阐明，潜在的治疗靶点有限，需要进一步探索潜在的分子机制并确定新的LUAD治疗靶点或生物标志物。

生物系统内部癌症的恶化机制十分复杂，并且癌症发生发展是一个动态变化的过程。越来越多的研究表明非编码RNA (non-coding RNA, ncRNA)的异常表达在癌症的发生发展中也扮演了重要角色 [5] 。研究非编码RNA在癌症中作用的重大发现是关于一个内源性竞争RNA (competing endogenous RNA, ceRNA)假说的提出，Pandolfi等于2011年在cell上发表了研究同时提出了“ceRNA假说”。该假说的理论核心为：微小RNA (microRNA, miRNA)是一种内源性的小分子非编码RNA，通常由20~25个核苷酸组成，大部分在细胞质中。在编码蛋白的基因转录本上的信使RNA (messenger RNA, mRNA)的3’ UTR非编码区存在多种miRNA的应答元件(miRNA response element, MRE)，miRNA可通过MRE与mRNA结合，导致mRNA降解或者抑制其翻译，从而导致靶基因表达下调。而在细胞中除mRNA之外，还存在另外一些RNA分子，例如长非编码RNA (long noncoding RNA, lncRNA)、环状RNA (circular RNA, circRNA)也存在MRE，当lncRNA、circRNA与mRNA存在相同的MRE时，他们之间构成了竞争相同种类miRNA的关系 [6] 。也就是说，lncRNA、circRNA通过MRE这个桥梁，参与调控了mRNA的表达水平，从而调控细胞功能。同时，转录因子(TF)作为基因调控网络和细胞信号传导的中枢，参与正常组织和肿瘤组织生长发育过程中的转录调控 [7] 。在转录水平，mRNA、lncRNA、circRNA和miRNA均受到转录因子TF的调控 [8] [9] [10] [11] 。TF和mRNA之间的相互作用是通过基因转录调控机制实现的，转录因子可以调节基因的转录过程，从而影响mRNA的合成 [12] [13] [14] [15] 。人们普遍认为，转录失调在大量癌症的癌变、转移、预后和药物依赖中同样起着关键作用 [16] 。

随着ceRNA假说的提出，越来越多的研究人员利用基因组学数据分析ceRNA网络探究癌症的潜在分子机制，或者构建以TF为主导的转录调控网络寻找与疾病有关的失调网络寻找潜在生物标志物 [17] [18] [19] 。然而，当前大部分基于ceRNA的研究忽略了转录阶段TF与ceRNA靶向结合关系，miRNA、lncRNA、circRNA、TF和mRNA在LUAD中结合癌症分期的特异性全转录组调控机制有待进一步阐明。本研究整合TF、mRNA、miRNA、lncRNA、circRNA构建表达量随分期动态变化的LUAD ceRNA五节点调控网络，揭示LUAD发生发展的分子机制，挖掘有助于LUAD诊断及预后的关键基因，为LUAD早期的诊断和治疗提供新思路。

2. 材料与方法

2.1. 数据获取

从TCGA数据库(https://portal.gdc.cancer.gov/)下载LUAD mRNA、lncRNA和成熟体miRNA表达谱。从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)下载LUAD circRNA数据集GSE101586、GSE101684及人类circRNA ID映射文件。

2.2. 数据预处理

将TCGA数据库miRNA、lncRNA、mRNA表达谱临床样本中不包括分期信息的样本剔除，剩余样本分为癌旁样本(normal)、早期样本(Stage I)、晚期样本(Stage II & III & IV)。

在GEO数据集中，对circRNA表达量进行标准化处理，用分位值法检查样本的数据分布，如果99%分位数大于100或最大值与最小值的差大于50且25%分位数大于0则对该样本数据进行log2对数转换。其次根据soft文件对样本进行注释并将Arraystar Human circRNA microarray V2 ID映射人类circRNA标准化7位ID。

2.3. 差异分析

使用R“Deseq2”包对miRNA、lncRNA、mRNA表达谱分别进行早期与癌旁、晚期与早期差异分析，阈值为|log2FC| > 0.58，校正后P < 0.05 (FC：fold change，差异倍数)，将差异分析结果取交集，获得差异表达mRNA (differentially expressed mRNA, DEmRNA)、差异表达miRNA (differentially expressed miRNA, DEmiRNA)和差异表达lncRNA (differentially expressed lncRNA, DElncRNA)。

使用R“limma”包对GSE101586、GSE101684数据集circRNA癌症样本和癌旁样本进行差异分析，阈值为|logFC| > 0.58，校正后P < 0.05，将两个数据集差异结果取交集，获得差异表达circRNA (differentially expressed circRNA, DEcircRNA)。

2.4. 五节点动态ceRNA调控网络构建

基于已经发表的经实验鉴定的数据库HOCOMOCO V11、AnimalTFDB手动整理人类转录因子数据集，并从DEmRNA中选出充当TF的基因。基于ChipBase实验数据库中筛选转录起始位点附近−30 KB到10 kb区域内的DElncRNA-TF；基于GTRD预测数据库及TransmiR实验数据库中取交集获得DEmiRNA-TF；在GTRD预测数据库和TRRUST实验数据库中取交集获得DEmRNA-TF。

基于CircBank数据库和Starbase数据库筛选DEmiRNA-DEcircRNA相互作用。其中Circbank数据库采用miRanda和Targetscan算法对140790条human circRNA以及1917条human miRNA的结合位点进行预测，筛选其中预测分数 > 200的相互作用对。使用miR2Disease、miRecords、miRTarBase、TarBase、Starbase 5个实验数据库筛选DEmiRNA-DEmRNA。其中miR2Disease、miRecords、miRTarBase和TarBase 4个数据库的数据主要来自高通量和低通量等生物实验，其中选择miR2Disease和miRecords中的所有靶标对，选择miRTarBase中强实验验证的靶标对，选择TarBase数据库中直接验证的靶标对。使用LncBase实验数据库和Starbase预测数据库筛选DEmiRNA-DElncRNA相互作用。

由于lncRNA具有组织特异性，根据定位不同产生不同的作用机制。对lncRNA进行亚细胞定位首先在UCSC数据库中下载筛选出的lncRNA的序列，将序列输入lncLocator数据库中预测lncRNA亚细胞定位。

利用cytoscape整合以上相互作用，构建转录和转录后水平表达量随分期动态变化的TF、mRNA、miRNA、lncRNA、circRNA五节点ceRNA调控网络。

2.5. 基因本体论(Gene Ontology, GO)富集分析

利用cytoscape中的cluego插件对靶基因进行GO富集分析，筛选其中P < 0.05的通路。

2.6. 蛋白互作(Protein-Protein Interaction Network, PPI)网络构建

利用String构建五节点调控网络中靶基因的蛋白互作PPI网络，筛选MCC排名前8的基因作为核心基因。

3. 结果

3.1. 数据获取及预处理

下载TCGA包含16,869个lncRNA、19,937个mRNA (539个LUAD样本，59个癌旁样本)以及2213个成熟体miRNA (521个LUAD样本，46个癌旁样本)，GSE101586包含5490个circRNA (LUAD样本5个，癌旁样本5个)，预处理后，TCGALUAD数据集包含16,414个lncRNA、19,494个mRNA (296个LUAD早期样本，235个LUAD晚期样本，58个癌旁样本)以及2212个成熟体miRNA (283个LUAD早期样本，231个LUAD晚期样本，46个癌旁样本)。

下载GSE101586包含5490个circRNA (LUAD样本5个，癌旁样本5个)，GSE101684包含9114个circRNA (LUAD样本4个，癌旁样本4个)。

3.2. 差异基因筛选

癌旁与早期差异miRNA共541个，其中300个上调241个下调，差异lncRNA共6699个，其中5150个上调1549个下调。DEmRNA共8140个，其中5303个上调、2837个下调。

早期与晚期差异miRNA共41个，其中8个上调23个下调。差异lncRNA共2115个，其中371个上调1744个下调。mRNA共1149个，其中300个上调、849个下调。

GSE101684中差异circRNA共2125个，其中1124个上调，1001个下调。GSE101586中差异circRNA共182个，其中154个上调，28个下调。

取交集后，在早晚期两组差异分析结果中均差异表达23个DEmiRNA、1441个DElncRNA和880个DEmRNA，45个DEcircRNA。在差异基因中筛选出27个TF (图1)。

Figure 1. Differential analysis volcano plot; (A) paraneoplastic vs. early differential expression of miRNA; (B) early vs. late differential expression of miRNA; (C) paraneoplastic vs. early differential expression of lncRNA; (D) early vs. late differential expression of lncRNA; (E) paraneoplastic vs. early differential expression of mRNAF; (F) early vs. late differential expression of mRNA; (G) GSE101684 differentially expresses circRNA; (H) GSE101586 differentially expresses circRNA, the blue color in the plot represents down-regulation, the red color represents up-regulation

图1. 差异分析火山图；(A) 癌旁与早期差异表达miRNA；(B) 早期与晚期差异表达miRNA；(C) 癌旁与早期差异表达lncRNA；(D) 早期与晚期差异表达lncRNA；(E) 癌旁与早期差异表达mRNA；(F) 癌旁与早期差异表达mRNA；(G) GSE101684差异表达circRNA；(H) GSE101586差异表达circRNA，图中蓝色代表下调，红色代表上调

3.3. 构建五节点动态ceRNA调控网络

Figure 2. 5-node dynamic ceRNA regulatory network; (A) paraneoplastic vs. early stage differential regulatory network; (B) early vs. late stage differential regulatory network, where triangles represent TFs, squares represent mRNAs, V-shape represents miRNAs, ovals represent circRNAs, and diamonds represent lncRNAs

图2. 五节点动态ceRNA调控网络；(A) 癌旁与早期差异调控网络；(B) 早期与晚期差异调控网络，其中三角形代表TF，方形代表mRNA，V型代表miRNA，椭圆形代表circRNA，菱形代表lncRNA

基于差异表达结果，在数据库中得到转录关系对：4对DEmiRNA-TF，其中1个lncRNA、4个TF；29对DEmiRNA-TF，其中7个miRNA、9个TF；87对DEmRNA-TF，其中32个mRNA、9个TF。转录后关系对：31对DEcircRNA-miRNA，其中15个miRNA18个circRNA；DEmRNA-miRNA共58对，其中miRNA 13个、mRNA 45个；DElncRNA-miRNA共70对相互作用，其中43个lncRNA、12个miRNA。亚细胞定位结果显示，4个lncRNA定位于细胞质中。

整合的转录水平和转录后水平的ceRNA调控网络中共50个节点174对调控关系，其中miRNA 9个、circRNA 2个、lncRNA 4个、mRNA 30个、TF 11个。最终得到早晚期表达量动态变化的五节点ceRNA网络(图2)，图中蓝色节点代表下调，红色节点代表上调，颜色深浅代表差异表达logFC的大小，节点大小代表度的大小。

从以上分期网络中可见随疾病分期显著动态表达的关键调控因子和基因。如ALB在早期表达量相比于癌旁显著上调(logFC = 5.8, adj.P < 0.01)，而在晚期的表达量显著下调(logFC = −3.2, adj.P < 0.01)。HNF4A在早期表达量相比于癌旁显著上调(logFC = 2.7, adj.P < 0.01)，晚期的表达量持续上调(logFC = 0.7, adj.P < 0.01)。RBP4在早期表达量相比于癌旁显著下调(logFC = −2.6, adj.P < 0.01)，而在晚期的表达量上调(logFC = 0.7, adj.P < 0.01)。ERBB4在早期表达量相比于癌旁显著下调(logFC = −1.4, adj.P < 0.01)，晚期的表达量持续下调(logFC = 0.8, adj.P < 0.01)。

3.4. 基因本体论(Gene Ontology, GO)富集分析

靶基因主要富集在亮氨酸运输、脂肪酸代谢、突触成熟，以及乳腺肺泡发育等生物过程中(图3)。

图中不同颜色代表富集不同的生物学过程，其中紫色代表脂肪酸代谢过程，深蓝色代表色氨酸运输，黄色代表乳腺肺泡发育，粉色代表突触成熟，圆圈大小代表富集生物种类中包含的通路数量，圆圈越大代表富集到的通路数量越多。

Figure 3. GO enrichment analysis

图3. GO富集分析

3.5. 蛋白相互作用(Protein-Protein Interaction Network, PPI)网络的构建

通过String构建核心基因蛋白互作网络(图4)，最终得到MCC排名前8的关键靶基因，分别为NEFL、RBP4、FGA、SLC2A1、ALB、AFP、SLC7A5、DKK1。

颜色越深代表MCC值越大。

Figure 4. PPI network interoperability diagram

图4. PPI网络互作图

4. 讨论

通过构建的表达量随分期动态变化的LUAD 5节点ceRNA调控网络，得到的靶基因主要富集在营养物质代谢和运输、突触成熟、以及乳腺肺泡发育等生物学过程，这些过程都和LUAD发生发展有紧密联系。代谢重编程是恶性肿瘤的标志之一，主要表现为糖酵解增强、谷氨酰胺代谢活跃及脂质代谢异常 [20] 。肺癌中，肿瘤细胞糖酵解产生的乳酸诱导微环境中调节性T细胞上的程序性死亡蛋白-1 (programmed death-1, PD-1)表达上调，导致PD-1免疫阻断疗法的失效 [21] 。而抑制葡萄糖转运蛋白阻碍LUAD糖酵解，可减缓LUAD发展 [22] 。哺乳动物雷帕霉素靶点(mammalian target of rapamycin, mTOR)是细胞生命活动和代谢调控中的关键因素 [23] 。它包含两种复合体，mTORC1和mTORC2，其中mTORC1在氨基酸充足时激活，促进细胞生长通过增加合成代谢和抑制分解代谢。研究发现SAR1B通过控制亮氨酸依赖的mTORC1信号来影响肿瘤发生，特别是在人LUAD中，SAR1B常缺失，导致亮氨酸敏感性消失和mTORC1的结构性激活。通过实验模型，研究证明SAR1A和SAR1B的缺失不仅促进皮下和肺内肿瘤生长，而且促进肿瘤的代谢重编程和增殖 [24] 。此外，神经元还可以和肿瘤细胞形成肿瘤–神经突触。这些信号传导机制通常会激活促进肿瘤生长的典型致癌信号通路 [25] 。下丘脑–垂体轴及其产生的激素α-MSH可以介导肿瘤诱导的髓系造血(MDSCs)和免疫抑制，研究人员通过肿瘤模型发现，非小细胞肺癌(NSCLC)患者血清中α-MSH浓度显著升高并与外周血中的MDSCs比例呈正相关 [26] 。另外，乳腺癌与肺癌的发生存在某些共同机制，如基因、晚期糖基化终产物受体、雌激素及烟草等。乳腺癌放疗则通过一些直接和间接作用在一定程度上促进第二原发性肺癌的发生 [27] 。

通过构建的网络的靶基因构建的PPI网络，获得了NEFL、RBP4、FGA、SLC2A1、ALB、AFP、SLC7A5、DKK1这8个核心基因，它们都与肿瘤发生发展有关。

神经丝轻链多肽(NEFL)是细胞骨架的关键成分。NEFL基因被证明为一种抑癌基因，位于人类染色体8p21，而这一区域富含抑癌基因 [28] 。研究人员通过对108例肺癌患者的NEFL的表达进行免疫组化检测和随访观察发现，NEFL的表达与肺癌淋巴结转移有关系，并且其阳性表达与患者预后呈正相关 [29] 。

视黄醇结合蛋白4 (RBP4)是一种在肝脏合成的RBP家族分泌分子 [30] 。流行病学研究表明，RBP4与多种恶性肿瘤的发生和发展有关，在肝癌、胰腺癌、急性白血病等肿瘤中高表达 [31] [32] 。胡等人将256例确诊的NSCLC病例和256例年龄和性别匹配的健康对照者病例进行对照研究，发现NSCLC患者血清RBP4水平明显高于健康对照组(36.05 ± 8.28 vs 29.54 ± 7.71 μg/mL, P < 0.05)。血清RBP4水平升高与NSCLC风险增加相关(P-value = 0.001) [33] 。

纤维蛋白原α链(FGA)，是纤维蛋白原的α成分 [34] 。研究人员使用CRISPR/Cas9基因组编辑在两个LUAD细胞系A549和H1299中进行实验，发现敲除FGA可促进LUAD细胞的生长、迁移和侵袭，但服用FGA可抑制LUAD细胞在肺部的表达。在体外LUAD细胞和体内异种移植肿瘤模型中进行功能分析显示，FGA通过凋亡和EMT介导的肿瘤生长和转移调控还参与了整合素-AKT信号通路。研究结果表明，FGA对LUAD细胞的生长和转移起着抑制作用 [35] 。

溶质运载家族2促进葡萄糖转运体成员1 (SLC2A1)是细胞能量代谢途径中的一种重要蛋白质。在正常组织中，SLC2A1仅限于在红细胞和血脑屏障的内皮细胞上表达 [36] 。最近，SLC2A1被证明是恶性肿瘤细胞中葡萄糖转运的关键限速因子，并在几种不同类型的人类癌症中过度表达，包括肝癌、胰腺癌、子宫内膜癌、乳腺癌和肺癌 [37] [38] [39] 。一项研究表明，SLC2A1表达的升高与LUAD患者的预后密切相关，同时，敲除SLC2A1可抑制LUAD细胞的增殖、集落形成、侵袭和葡萄糖利用 [40] 。

溶质载体家族7成员5 (SLC7A5)是色氨酸转运蛋白 [41] 。一项基于血清代谢变化和体外实验的研究发现，在A549肺癌细胞中敲除SLC7A5后，色氨酸的细胞水平显著降低，糖酵解活性的减弱，表明SLC7A5在促进肺癌的糖酵解过程中起重要作用 [42] 。

Dickkopf相关蛋白1 (DKK1)是Wnt/B-catenin信号通路的典型抑制剂，属于DKK家族 [43] 。DKK1作为一种分泌型糖蛋白，主要通过与Wnt蛋白竞争性结合LRP5/6受体，拮抗Wnt/β-catenin信号通路活性，调节细胞增殖、分化及癌变，影响细胞凋亡、瘤细胞侵袭和转移 [44] 。Wnt拮抗剂启动子高甲基化导致的基因沉默发生在LUAD的早期阶段，并随着恶性肿瘤的发展而不断累积 [45] 。

血清白蛋白(ALB)水平是血液学指标中最常用的一项评估指标 [46] 。血浆ALB水平降低会引起体内环境紊乱，导致机体代谢异常，多项研究表明血浆ALB降低与肿瘤不良预后相关 [47] 。血浆白蛋白可以反映人体营养状态及凝血功能，机体的营养状态及凝血功能可以用于判断恶性肿瘤患者的预后 [48] [49] [50] ，近年来ALB被广泛应用于判断恶性肿瘤患者的预后 [51] [52] [53] 。

甲胎蛋白(AFP)是胚胎期血浆蛋白的主要成分，其功能与结构都与白蛋白有很多相似之处(属于类白蛋白家族) [54] 。AFP在免疫应答中发挥作用，认为其具有免疫抑制的作用，是T细胞和B细胞的免疫调节剂，主要表现为抑制母体对胚胎发育的免疫应答以及肿瘤患者对肿瘤的免疫应答 [55] 。AFP是当前诊断肝癌常用的方法和最重要的肿瘤指标，约60%~70%肝癌患者AFP升高，其中约18%肝癌患者AFP低浓度升高 [56] 。

5. 结论

综上所述，本研究基于生物信息学构建了LUAD动态五节点ceRNA调控网络，通过数据处理、差异分析、数据库预测与筛选得到了关键调控因子与关键基因，并通过富集分析及PPI蛋白互作网络构建筛选出8个核心基因，初步探索LUAD发生发展的机制，此结果有助于为LUAD诊断与治疗提供新思路。然而，由于目前能收集到的circRNA数据有限，后续的研究应该结合包含circRNA临床分期样本的表达数据构建调控网络，进一步挖掘分子机制。同时，后续的研究应该进一步结合相关临床数据和实验对上述研究结果进行验证。

NOTES

^*通讯作者。