肺腺癌预后相关miRNA生物信息学筛选及其临床意义*

张皓旻， 杨 波， 陈红飞， 迟小华， 席义博， 陈熙勐， 郭 斌， 贺培凤， 卢学春△

(1. 山西医科大学管理学院, 太原 030001； 2. 中国人民解放军总医院南楼血液科 国家老年疾病临床医学研究中心, 北京 100853； 3. 中国人民解放军火箭军总医院药剂科, 北京 100800)

根据世界卫生组织(WHO)国际癌症研究机构(IARC)最新报道，肺癌是人类最常见的恶性肿瘤类型，其发病率(11.6%)和死亡率(18.4%)位居所有癌症类型之首[1]。在肺癌中，非小细胞性肺癌(non-small cell lung cancer, NSCLC)最常见，其中肺腺癌(lung adenocarcinoma, LA)又是主要的病理类型，超过40%的患者发现时已属晚期[2]。随着基因组学、表观基因组学的发展，现代医学已进入精准医学时代，通过对肺腺癌的生物信息学分析，已有大量关于肺腺癌发生发展、转移、耐药及预后相关分子标记物的研究报道[?]，且开发出分子靶向药物改变了肺腺癌的临床治疗模式，如表皮生长因子受体酪氨酸激酶抑制剂(epidermal growth factor receptor-tyrosine kinase inhibitor，EGFR-TKIs)、间变性淋巴瘤激酶(anaplastic lymphoma kinase, ALK)激酶抑制剂等[3,4]。由此可见，基于组学大数据深度挖掘是寻找肺腺癌诊断、治疗及预后相关靶点的重要手段。

目前已知，微小RNA(microRNA, miRNA)可参与肺腺癌的发生发展，它们通过调控靶基因及其相关信号通路，发挥与癌基因或抑癌基因类似的作用[5]。同时，有些由于某些基因会对患者预后产生影响，故而称之为预后基因，在miRNA 中同样如此。以往研究重点往往关注某个miRNA在肺腺癌病理机制中的作用，缺乏针对肺腺癌全miRNA组的系统研究[6,7]。本研究通过对癌症基因组图谱(The cancer genome atlas, TCGA)数据库中567例肺腺癌miRNA组学数据的生物信息学分析，探讨肺腺癌预后相关miRNA组学特征及其临床意义，以期为肺腺癌的诊断和治疗提供理论依据。

1 材料与方法

1.1 数据采集

TCGA[8]由美国国家肿瘤中心( National Cancer Institute，NCI) 和国家人类基因组研究所( National Human Genome Research Institute，NHGRI ) 共同发起，截至目前已收录50多种肿瘤类型。本研究从TCGA数据库收集并处理的公开肺腺癌miRNA表达谱数据567例，其中正常样本46例、肺腺癌患者样本521例，同时下载了504例肺腺癌患者的生存数据。

1.2 肺腺癌miRNA差异分析

用R语言DEseq程序包，对46份正常组织和521份肺腺癌组织样本进行差异miRNA分析并绘制热图，使用Origin软件绘制火山图[9]。差异miRNA筛选阈值为FDR<0.05，/LogFC/>1。

1.3 差异miRNA的生存分析

随后，运用R语言对得到的差异miRNA进行生存分析。首先，对差异miRNA进行单因素Cox风险回归分析，依据假设检验结果初步筛选对肺腺癌预后可能具有影响的miRNA；随后，对单因素Cox风险回归分析筛选出的miRNA进行多因素Cox比例风险回归分析，依据回归系数(β)及风险比(hazard ratios，HR)来识别保护(HR<1)和危险miRNA(HR>1)，同时确定肺腺癌中高风险miRNA(β>0)。所有统计检验均采用双侧检验[10]。

1.4 miRNA靶基因分析

miRWalk(http：//mirwalk.uni-hd.de/)是一个用于预测miRNA靶基因的开源在线分析工具[11]。我们首先应用miRWalk数据库对上述能够影响肺腺癌患者生存时间的miRNA进行靶基因预测。随后，依据靶基因在DAVID数据库中的聚类结果，进一步分析miRNA功能。最后，利用Cytoscape软件绘制miRNA-mRNA调控网络[12]。

2 结果

2.1 筛选肺腺癌差异miRNA

利用来自TCGA数据库的46份正常组织和521份肺腺癌组织样本，筛选出46个差异表达miRNA，其中上调19个、下调27个，差异基因的层次聚类图如图1A所示。随后利用Origin绘制火山图，结果如图1B所示。差异miRNA的层次聚类图显示，肿瘤和正常组织形成了明显的类团，各组组内数据关联度强，可以很好地区分正常组织与癌组织这两组样本(图1)。

Fig.1Volcano map and hierarchical cluster heat map of differential miRNA

A： A hierarchical clustering heat map; B： Volcano map of differential miRNA

2.2 差异miRNA生存分析

为分析差异miRNA对肺腺癌的预后影响，选取504个带有随访的肺腺癌病例进行生存分析。首先，利用单因素Cox风险回归模型进行生存分析，初步筛选出7个可能对肺腺癌预后产生影响miRNA，且具有统计学及意义(P<0.05，表1)。随后，利用多因素Cox风险回归模型对初步筛选出的7个miRNA进行生存分析，进一步确定高风险miRNA，结果显示(表2)：hsa-mir-21、hsa-mir-378e与预后不良有关(β>0，HR>1)，且为高风险miRNA，生存曲线具有统计学意义(图2A)。绘制受试者工作特征曲线，即ROC(receiver operating characteristic curve，ROC)曲线，显示AUC=0.618，结果较为准确(图2B)。hsa-mir-101、hsa-let-7c、hsa-mir-142、hsa-mir-200a为保护性miRNA(β<0，HR<1)，提示其表达与肺腺癌预后良好有关。所有预后相关miRNA热图见图2C。

Tab.1Analysis of single factor Cox regression analysis for the effect of microRNA expression on survival time of patients with lung adenocarcinoma

miRNAHRZ-scorePExpressionhsa-mir-148a0.785301-3.216940.001296uphsa-mir-1010.735497-2.985850.002828downhsa-mir-378e12.15492.8405660.004503downhsa-mir-211.2800672.3078110.02101uphsa-mir-200a0.855878-2.2620.023697uphsa-mir-1420.881663-2.166020.03031uphsa-let-7c0.858666-2.096930.036down

miRNA： microRNA; HR： Hazard ratio

Tab.2Analysis of multiple Cox regression analysis for the effect of microRNA expression on survival time of patients with lung adenocarcinoma

miRNAβHRse(β)zPExpressionhsa-mir-210.17111.18660.11391.50.01331uphsa-mir-101-0.22530.79830.1182-1.910.0466downhsa-let-7c-0.12520.88230.087-1.440.01504downhsa-mir-142-0.18610.83020.0654-2.840.0044uphsa-mir-200a-0.20170.81730.0702-2.870.0041uphsa-mir-378e2.21259.1390.88722.490.0126down

miRNA： microRNA; β： Regression coefficient; HR： Hazard ratio

2.3 miRNA靶基因分析

首先，利用miRWalk数据库进行靶基因预测，选择miRWalk数据库内包含的全部数据库，设置各数据库连词符为“AND”即输出所有数据库共有靶基因，设置阈值P<0.01筛选靶基因，分析结果详见表3，挑选各miRNA排名前20靶基因进行图形化展示，详见图3A。随后，利用DAVID数据库进行聚类分析，从中挑选肺腺癌相关通路，分析各miRNA与肺腺癌的潜在联系。最后，利用Cytoscape软件绘制miRNA—基因—通路互作网络，详见图3B。

Fig.2Survival curves, ROC curves, and prognosis related miRNA heat maps of prognostic-associated miRNAs

A： Risk survival curve; B： The ROC curve; C： The prognosis related miRNA heat map

Tab.3Target genes analysis of miRNA which affecting the survival time of lung adenocarcinoma

miRNAExpressionTarget genehsa-mir-21↑265hsa-mir-142↑368hsa-mir-200a↑317hsa-mir-101↓173hsa-let-7c↓326hsa-mir-378e↓144

Fig.3miRNA target gene analysis

A： miRNA and its top 20 target genes; B： miRNA-target gene-pathway network

3 讨论

微小RNA(microRNA，miRNA)是长度约为～22 nt的小的非编码RNA。它们通过与靶mRNA的互补碱基配对在转录后水平调节基因表达，使mRNA降解，阻断翻译[13]。近年来，miRNA与肺腺癌相关性的报道日渐增多。各国研究人员为方便研究开发了大量用于预测miRNA-mRNA相互作用的工具，所有这些工具都使用不同的算法。因此，合理选择靶基因预测工具是miRNA靶基因预测的关键[14]。miRWalk(http：//mirwalk.uni-hd.de/)是一个用于预测miRNA靶基因的开源在线分析工具[11]。涵盖了人、小鼠和大鼠等物种的所有目前已知miRNA结合位点信息。 此外，它还整合了许多新功能。例如，同时对来自不同数据库的miRNA靶基因信息的比较分析、miRNA及其靶基因在通路中的相互作用网络、以及分析miRNA及其靶基因在疾病和遗传中的相互作用，同时还包含实验信息(例如，细胞系，疾病，miRNA加工蛋白)。此类工具的出现为肺腺癌miRNA研究的开展提供了便利。有研究发现，miRNA能通过调控上皮细胞或维持间质细胞特征相关的靶基因以及相关信号通路来调控EMT，从而促进肺腺癌的转移[15]。又如，Misono等人发现[16]，miR-145可能与肺腺癌患者的预后不良有关。此外，还有学者发现miR-205可能与EGFR患者的酪氨酸激酶抑制剂(tyrosine kinase inhibitors，TKIs)耐药有关[17]。上述研究结果表明，miRNA可能参与调控肺腺癌发生、发展的全过程，同时其表达紊乱还与肺腺癌转移、耐药、预后不良有关[18,19]。这些都提示miRNA可以成为肺腺癌早期诊断及预后判断的分子标志物，针对肺腺癌的全miRNA组研究有望帮助临床医生发现新的肺腺癌治疗策略。

本研究通过对来自TCGA数据库的567例肺腺癌患者的miRNA组数据进行差异分析，得到46个差异miRNA。随后，利用单因素Cox风险回归模型和多因素Cox风险回归模型进行生存分析筛选出6个与肺腺癌预后不良相关的风险miRNA，可作为肺腺癌预后独立危险因素，高表达与患者生存期缩短显著相关(P<0.05)。通过对6个风险miRNA进行功能富集分析发现，他们均与肺腺癌的发生发展有关，主要体现在对细胞周期、代谢、细胞间相互作用、转录调控产生的影响。

已有大量文献证实，hsa-miR-21与肺腺癌预后不良有关[20]。Kunita等人[21]利用体外细胞实验发现hsa-miR-21上调可诱导成纤维细胞转分化为癌症相关成纤维细胞，从而支持肺腺癌进展，与本研究结果一致。此外，本研究中发现hsa-miR-21所调控的靶基因大量的富集在免疫相关通路中，例如Toll样受体信号通路、NF-Κb信号通路、趋化因子信号通路以及细胞因子与细胞因子受体相互作用等，提示hsa-miR-21可能通过刺激肿瘤微环境中促炎因子的释放，促进肺腺癌的发展，其分子机制有待进一步验证。此外，有国内学者发现在乳腺癌中hsa-mir-21可以通过影响DNA甲基化水平影响乳腺癌患者预后[22]。Kaduthanam等人[23]发现，hsa-miR-142表达显著增加肺腺癌患者术后复发的风险，可作为肺腺癌术后风险评估的标志物，与本研究结果一致。本研究发现hsa-miR-142能够影响代谢相关通路，猜测其可能与肺腺癌代谢环境改变有关，同时发现hsa-miR-142能够通过调控KRAS基因，而KRAS基因的异常表达与KRAS驱动肺腺癌关系密切。Lin等人[24]发现hsa-miR-200a可能与肺腺癌预后不良有关，与本研究结果一致。Yang等人[25]发现DNMT1/miR-200a/GOLM1信号通路对肺腺癌细胞的增殖具有调控作用。本研究发现，hsa-miR-200a在肺腺癌样本中表达量增加，且hsa-miR-200a调控靶基因大量富集于免疫、代谢相关通路，尤其是AMPK信号通路，这提示hsa-miR-200a对肺腺癌细胞增殖的调控，不仅与其表达变化及靶基因结合有关，还与肺腺癌细胞条件有关，有待进一步验证。Li等人[26]发现hsa-miR-101下调与肺腺癌患者预后不良有关，与本研究结果一致。本研究发现hsa-miR-101与药物代谢相关通路有关，推测hsa-miR-101可能在肺腺癌治疗耐药上发挥作用，具体分子机制有待进一步验证。Cui等人[27]发现，hsa-let-7c与肺腺癌化疗耐药有关。本研究发现，hsa-let-7c可以调控ERBB2基因，而ERBB2基因是表皮生长因子受体(EGFR)的编码基因，因而推测其与EGFR驱动肺腺癌关系密切，具体机制有待进一步验证。目前暂无文献报道hsa-miR-378e与肺腺癌的关系。本研究中尽管在对其调控靶基因的通路富集分析中未发现具有统计学意义的信号通路，但其调控靶基因大都与肺腺癌密切相关，例如FGF12[28]、RBX1[29]、SLC2A1[30]等，因此hsa-miR-378e有望成为新的肺腺癌预后标志物，其在肺腺癌中的具体作用有待进一步研究。

目前已知，miRNA不仅可以通过调节靶mRNA的表达来发挥作用，还可以通过影响甲基化进而影响基因表达。同时，也可受到其他包括lncRNA在内的非编码RNA调控。因此，开展肺腺癌全miRNA组学研究，有助于明晰肺腺癌的转录调控机制，为后续开展肺腺癌表观遗传学研究提供理论依据。也可为肺腺癌早期诊断标志物的发现；转移、复发、耐药及预后不良相关的分子靶点的发现提供方向和指导。最终为肺腺癌现有临床方案的优化提供理论依据，为难治肺癌患者尽早选择有效方案提供指导。

综上所述，通过应用来自TCGA的带有患者完整生存信息的肺腺癌miRNA组数据，发现6种miRNA与肺腺癌预后有关，其中包括4种保护性miRNA(hsa-mir-142、hsa-mir-200a、hsa-mir-101、hsa-let-7c)，2种风险miRNA(hsa-mir-21、、hsa-mir-378e)，未来有望作为肺腺癌预后的生物标志物继续研究。其中，hsa-mir-378e与肺腺癌预后不良的详细研究尚未见报道，有望成为新的肺腺癌预后标志物，本研究为临床肺腺癌诊疗提供新的方向和理论依据。