乳腺癌预后的相关基因生物信息学分析

左 　然 任晓虎 吴德生李 　萍 吴 　雯 谢 　妮袁建辉* 让蔚清*

（1．南华大学公共卫生学院，湖南　衡阳　421001；2．深圳市疾病预防控制中心毒理研究所，广东 　深 圳 　518055；3．深圳市第二人民医院转化医学研究院，广东 　深 圳 　518000）

乳腺癌是全世界女性最常见的侵袭性肿瘤之一，其发病率高居女性肿瘤榜首，对女性健康造成极大威胁[1-2]。乳腺癌的发病机制目前尚不明确，遗传、雌激素水平、生长因子、细胞因子、激酶和非编码RNA等均在肿瘤细胞生长中起重要作用[3-4]。microRNA是一类重要的内源性非编码RNA，主要通过转录后水平调控基因表达[5]。microRNA和表观遗传调控过程之间存在复杂的反馈网络并形成表观遗传学microRNA调控回路[6]。随着对microRNA功能研究的不断涌现，越来越多的研究表明microRNA与肿瘤之间存在密切关系。有报道发现，microRNA-214-3p通过抑制靶基因MELK表达从而抑制肝细胞癌进展[7]；而miR-137下调和肝癌淋巴结转移被发现与不良预后密切相关[8]；另外microRNA-200可影响管腔祖细胞增殖进而促进乳腺癌生长和转移[9]。尽管microRNA及其潜在靶向调控的基因已被广泛应用于肿瘤研究，但类似的分子事件在乳腺癌预后中鲜有涉及。本研究通过挖掘肿瘤基因组图谱计划(the cancer genome atlas，TCGA)中的乳腺癌基因表达谱和microRNA测序数据，寻找乳腺癌预后相关的关键性分子事件。

本研究中，我们利用van Iterson等[10]提出的全局算法整合基因表达谱和microRNA测序数据，并预测给定临床背景下的microRNA靶基因。该算法用到的统计学检验可以评估Pearson相关系数出现假阳性的情况，同时基于全局检验通过Lasso回归预测数据库中所有可能的目标mRNA[11-12]。为了进一步分析microRNA及相关靶基因在乳腺癌中起到的作用，我们利用通用基因集富集算法(general applicable gene set enrichment，GAGE)[13]分析microRNA及其靶基因涉及的异常变化通路。进一步证实所筛选的microRNA在乳腺癌组织和细胞中能负向调节其靶基因，并使用Cox回归风险模型比例风险模型评估其作为预后标志物的潜力。本研究分别对乳腺癌相关microRNAs和基因的表达数据进行分析。

1　材料与方法

1.1　数据准备

从TCGA数据库中(http://cancergenome.nih.gov/)下载经Lowess标准化和对数转化的乳腺癌3级基因表达数据和microRNA测序数据。将所有的基因符号转换为Entrez Gene ID后，整理成矩阵并导入结构调查语言(structure query language，SQL)数据库中。

1.2　乳腺癌中的基因与microRNA表达的定量分析

将来自120对癌-癌旁组织的基因表达谱数据以及204对癌-癌旁组织的microRNA测序数据整理、优化后，利用R语言Limma软件包(Limma package R 3.1.1)通过整合线性模型、经验贝叶斯理论与传统t检验对基因表达谱数据进行定量及统计分析。利用Holm-Bonferroni方法控制多次统计学检验可能产生的假阳性，将假阳性率(false discovery rate，FDR)控制在0.01。差异表达基因的筛选标准为：FDR值(即矫正后的P值)小于0.01，差异表达绝对值变化大于4。对于异常改变的microRNA，筛选标准为FDR值(即矫正后的P值)小于0.01，差异表达绝对值变化大于1.5[14-15]。

1.3　乳腺癌中microRNAs潜在调控靶基因的预测

收集388例乳腺癌和388例配对癌旁组织的基因表达谱和microRNA表达数据。从miRTarBase(http://mirtarbase.mbc.nctu.edu.tw/)，MSigDBC3 motif Gene set(http://www.broadinstitute.org/)，PITA(http://genie.weizmann.ac.il/pubs/mir07/mir07_data.html)，Microcosm(http://www.ebi.ac.uk/enright-srv/microcosm/htdocs/targets/v5/)和TargetScan(http://www.targetscan.org/) 数据库获得microRNA靶基因。将来自388名患者基因谱表达数据和microRNA测序数据整理为矩阵格式。利用van Iterson等提出的基于线性(或广义线性)回归模型的全局无效假设的microRNA靶基因预测算法分析microRNA的潜在靶基因，该方法旨在提供一种有效的方法来评估多种疾病之间的定量比较，并预测每个潜在的microRNA-mRNA的关联性[10]。该算法在R软件microRNA-mRNA包中实现(http://www.humgen.nl/Microarray Analysis Group.html)。选择同时在2个或以上数据库中预测到的microRNA靶基因作为基础数据库。将样本随机分为训练集(所有样品的2/3)和测试集(所有样品的1/3)，把microRNA预测靶基因对引入线性回归模型的全局无效假设方法应用于训练集，找出关联并在测试集中进行验证。

1.4　基因集富集分析

传统定量分析关注单个基因显著变化却忽视了多个基因之间的协同作用，因此通用基因集富集算法(GAGE)是对特定基因集合(如特定信号通路)的变化进行分析，而不考虑单个基因改变。将基因表达谱数据整理成矩阵(基因表达为行，样本分组为列)导入R软件，通过GAGE包进行分析。基因集来自MSigDB C2 Curated Gene set(http://www.broadinstitute.org/)。不同于传统的基因集富集，GAGE算法允许基因同时向不同方向(上调或下调)变化，因此更适合于通路分析。Benjamini和Hochberg方法来控制多重测试的假阳性率，选择校正后的P值(即Q值)＜0.05基因集进一步分析。基因集的大小(包括在基因集中的基因的数目)设定在50～500个之间，选择差异变化大于1倍标准差的基因集。

1.5　Cox回归风险模型筛选乳腺癌潜在评估预后的分子靶标基因

利用R survival包(R package survival)分析microRNA及潜在靶基因在乳腺癌预后中的作用。以microRNA测序数据及基因表达谱数据作为主要自变量，生存时间为因变量，相关临床资料作为次要自变量，使用Cox回归风险模型进行迭代分析，Kaplan-Meier进行统计学检验。

2　结果

2.1　乳腺癌中异常表达的基因和microRNA

对基因表达谱和microRNA测序数据的分析结果显示，17 533个基因中只有344个基因表达呈显著改变(矫正后p＜0.01，差异表达绝对值变化大于4)，见图1A和表1。而microRNA中，135个在乳腺癌组织中表达发生异常改变化(矫正后p＜0.01，差异表达绝对值变化大于1.5)，见图1B和表2。

图1　mRNA和microRNA表达谱的火山图

表1　microRNA及其潜在靶基因实验验证

2.2　乳腺癌中基于表达数据的microRNAs潜在靶基因预测

对338对乳腺癌-癌旁组织的microRNA及基因表达谱数据进行关联性研究发现，31个microRNA与6 429个预测靶基因显著相关。其中有部分microRNA对基因表达呈正相关，尽管研究报道发现microRNA可诱导上调靶mRNA对细胞周期阻滞[16-17]，亦可促进结合核糖体mRNA基因的5´UTR和增强翻译[18]，但目前尚无证据表明在人群中microRNA上调能引起mRNA表达增加，因此正相关可能由于microRNA对基因的间接调控作用。然而为保证结果的可靠性，我们仅选择与靶基因表达呈负相关的microRNA进行验证。利用另外60对癌-癌旁的microRNA和基因表达谱数据，对14个microRNA和67个负相关的潜在靶基因进行验证。结果显示，其中6个microRNA表达上调，其潜在靶基因表达下调，而8个microRNA表达下调，其潜在调控靶基因上调。

2.3　microRNA及其潜在调控靶基因的通路富集

为了进一步探讨microRNA及相关靶基因在乳腺癌中可能的作用机制，基因集富集算法对microRNA测序及基因表达谱数据进行富集分析，发现了8个乳腺癌相关microRNA及其调控的31个关键潜在调控靶基因参与139条乳腺癌富集通路(见图2)。

表2　乳腺癌相关临床资料的生存分析

2.4　乳腺癌潜在预后标志物

生存曲线显示，乳腺癌患者从初始诊断起的6年，生存率明显下降(2 000 d，图3A)。为避免样本变异和未记录样本的选择偏倚，选择有完整记录的样本(有临床结局和治疗方案)进行分析。由于某些类别的肿瘤分期(IV，Tis和X)样本量非常少，因此未被纳入生存分析。Cox回归风险模型显示，肿瘤边缘状态、更年期状态和淋巴结数量等相关临床资料对患者预后无明显影响。而放射治疗显示肿瘤边缘状态呈保护作用(HR=0.432，P=0.058，图3B)。进一步深入分析发现，分泌卷曲相关蛋白1(SFRP1)高表达呈保护作用(HR=0.9，P=0.015，图3C)，has-miR-342-5p既没有保护作用，也没有危害效应(HR=0.99，P=0.144，图3D)。然而，交互作用研究显示，has-miR-342-5p抑制靶基因SFRP1时提示不良预后(见表3)。另外，hasmiR-342-5p对生长特异型抑制基因2(GASP2)的抑制作用则代表不良预后(Cox多因素分析，P=0.016，HR=18.88)。

图3　乳腺癌患者的生存曲线

表3　乳腺癌潜在预后标志物

3　讨论

在本研究中，我们将microRNA与其靶基因的关联性研究与通路富集分析整合后筛选出有明显改变的microRNAs及其潜在靶基因。此外，我们对所有预测的microRNA-靶基因中，发现9对microRNA-靶基因已在其他研究中被实验验证(表1)。其中7个microRNA及其30个潜在靶基因被发现在139个乳腺癌异常通路中发挥关键作用。Cox回归风险模型分析揭示了来自这些microRNA及其靶基因参与的分子事件可作为乳腺癌重要的预后特征。使用Cox比例风险分析，SFRP1、CLDN19及KRT5在患者生存中被发现有保护作用。已知SFRP1是与乳腺癌细胞生长负相关的肿瘤抑制因子[19]。据报道，乳腺癌中SFRP1的失活与不良预后有关[20]。此外，SFRP1仅在2 400～3 500 d内表示出显著的危害效应，而在3 500 d后危害效应降低。然而，由于miR-342-5p的调节，SFRP1的保护作用受到了影响。

目前已有一些研究表明microRNA与乳腺癌的相关性，在乳腺癌的发生过程中起到促癌或抑癌的作用。miR-342-5p是一种多功能的microRNA，参与细胞增殖、分化、信号转导等许多重要的生物学过程。据报道miR-342-5p在由MCF-7细胞构成的三维组织模型中通过靶向DNA结合抑制剂4(ID4)和DNA甲基转移酶1(DNMT1)来调节CEACAM1诱导的管腔形成[21]。此外，miR-342-5p还通过调节Notch信号传导参与神经元干细胞分化[22]，并可作为阿尔茨海默症的潜在血浆生物标志物[23]。其他研究表明，miR-342-5p可能通过下调阿尔茨海默症转基因模型小鼠中的锚蛋白G来促成轴突病变[24]。

SFRP1是富含半胱氨酸结构域的SFRP家族成员，可能被Wnt结合。同时SFRP1还介导视网膜中感光细胞的极性[20]。CLDN19是一个重要的膜蛋白，在细胞间紧密连接中起关键作用[25]，CLDN19的异常表达可促进癌前病变的细胞向恶性方向发展[26]。也有报道CLDN19 在乳腺癌中表达异常降低[27]，与我们的研究结果一致。

在乳腺癌中microRNA靶向调控特定基因作为预后关键分子事件的研究还相对较少。因此我们的研究不仅揭示了miR-342-5p、SFRP1、CLDN19以及KRT5在乳腺癌发生发展中具有的重要生物学功能，同时还表明miR-342-5p可能通过抑制乳腺癌中的SFRP1、KRT5以及CLDN19来提示不良预后，这为miR-342-5p相关的分子事件作为乳腺癌临床监测的潜在指标提供了科学依据。

[1] SIEGEL R L，MILLER K D，JEMAL A. Cancer statistics，2017[J]. CA Cancer J Clin，2017，67(1)：7-30.

[2] CHEN W，ZHENG R，ZHANG S，et al. Cancer incidence and mortality in China，2013[J]. Cancer Lett，2017，401：63-71.

[3] JORDAN C T，GUZMAN M L. Mechanisms controlling pathogenesis and survival of leukemic stem cells[J]. Oncogene，2004，23(43)：7178-7187.

[4] TSATSANIS C，SPANDIDOS D A. The role of oncogenic kinases in human cancer (Review)[J]. Int J Mol Med，2000，5(6)：583-590.

[5] CHEN K，RAJEWSKY N. The evolution of gene regulation by transcription factors and microRNAs[J]. Nat Rev Genet，2007，8(2)：93-103.

[6] MOUTINHO C，ESTELLER M. MicroRNAs and epigenetics[J].Adv Cancer Res，2017，135：189-220.

[7] LI Y，LI Y，CHEN Y，et al. MicroRNA-214-3p inhibits proliferation and cell cycle progression by targeting MELK in hepatocellular carcinoma and correlates cancer prognosis[J].Cancer Cell Int，2017，17：102.

[8] CUI S， SUN Y， LIU Y， et al. MicroRNA137 has a suppressive role in liver cancer via targeting EZH2[J]. Mol Med Rep，2017，16(6)：9494-9502.

[9] SANCHEZ-CID L， PONS M， LOZANO J J， et al.MicroRNA-200， associated with metastatic breast cancer，promotes traits of mammary luminal progenitor cells[J].Oncotarget，2017，8(48)：83384-83406.

[10] VAN ITERSON M，BERVOETS S，DE MEIJER E J，et al.Integrated analysis of microRNA and mRNA expression：adding biological significance to microRNA target predictions[J].Nucleic Acids Res，2013，41(15)：e146.

[11] GOEMAN J J，MANSMANN U. Multiple testing on the directed acyclic graph of gene ontology[J]. Bioinformatics，2008，24(4)：537-544.

[12] LI X，GILL R，COOPER N G，et al. Modeling microRNA-mRNA interactions using PLS regression in human colon cancer[J]. BMC Med Genomics，2011，4：44.

[13] LUO W，FRIEDMAN M S，SHEDDEN K，et al. GAGE：generally applicable gene set enrichment for pathway analysis[J].BMC Bioinformatics，2009，10：161.

[14] TOMCZAK K，CZERWINSKA P，WIZNEROWICZ M. The Cancer Genome Atlas (TCGA)：an immeasurable source of knowledge[J]. Contemp Oncol (Pozn)，2015，19(1A)：68-77.

[15] LEE H，PALM J，GRIMES S M，et al. The cancer genome atlas clinical explorer：a web and mobile interface for identifying clinical-genomic driver associations[J]. Genome Med，2015，7：112.

[16] VASUDEVAN S，TONG Y，STEITZ J A. Switching from repression t o a ctivation：microRNAs c an u p-regulate t ranslation[J].Science，2007，318(5858)：1931-1934.

[17] ZHANG R，SU B. Small but influential：the role of microRNAs on gene regulatory network and 3'UTR evolution[J]. J Genet Genomics，2009，36(1)：1-6.

[18] OROM U A，NIELSEN F C，LUND A H. MicroRNA-10a binds the 5'UTR of ribosomal protein mRNAs and enhances their translation[J]. Mol Cell，2008，30(4)：460-471.

[19] KLOPOCKI E，KRISTIANSEN G，WILD P J，et al. Loss of SFRP1 is associated with breast cancer progression and poor prognosis in early stage tumors[J]. Int J Oncol，2004，25(3)：641-649.

[20] VEECK J， NIEDERACHER D， AN H， et al. Aberrant methylation of the Wnt antagonist SFRP1 in breast cancer is associated with unfavourable prognosis[J]. Oncogene，2006，25(24)：3479-3488.

[21] WENG C， NGUYEN T， SHIVELY J E. miRNA-342 Regulates CEACAM1-induced lumen formation in a threedimensional model of mammary gland morphogenesis[J]. J Biol Chem，2016，291(32)：16777-16786.

[22] GAO F，ZHANG Y F，ZHANG Z P，et al. miR-342-5p regulates neural stem cell proliferation and differentiation downstream to Notch signaling in mice[J]. Stem Cell Reports，2017，8(4)：1032-1045.

[23] LUGLI G，COHEN A M，BENNETT D A，et al. Plasma exosomal miRNAs in persons with and without Alzheimer disease：Altered expression and prospects for biomarkers[J].PLoS One，2015，10(10)：e0139233.

[24] SUN X，WU Y，GU M，et al. miR-342-5p decreases ankyrin G levels in Alzheimer's disease transgenic mouse models[J]. Cell Rep，2014，6(2)：264-270.

[25] TSUKITA S， FURUSE M. Pores in the wall： claudins constitute tight junction strands containing aqueous pores[J]. J Cell Biol，2000，149(1)：13-16.

[26] SINGH A B，SHARMA A，DHAWAN P. Claudin family of proteins and cancer：an overview[J]. J Oncol，2010，2010：541957.

[27] TOKES A M，SZASZ A M，JUHASZ E，et al. Expression of tight junction molecules in breast carcinomas analysed by array PCR and immunohistochemistry[J]. Pathol Oncol Res，2012，18(3)：593-606.