楊 燕 賀 鑫 龔鵬舉 宋文靜 魏 蕾 張京偉

1武漢大學(xué)中南醫(yī)院甲狀腺乳腺外科 湖北 武漢 430071；2武漢大學(xué)基礎(chǔ)醫(yī)學(xué)院病理與病理生理學(xué)教研室 湖北 武漢 430071

乳腺癌是女性中最常見的惡性腫瘤，也是女性癌癥死亡的主要原因之一［1］。先前的研究表明，雌激素受體（estrogen receptor,ER）陽性乳腺癌約占乳腺癌的70%，雌激素受體α（ERα）在乳腺癌的發(fā)生和發(fā)展過程中起重要作用［2-4］。ERα是一種核受體，由幾個結(jié)構(gòu)域組成：募集轉(zhuǎn)錄調(diào)節(jié)因子的N端激活功能域1（AF1）和C端激活功能域2（AF2），與靶基因的雌激素反應(yīng)元件結(jié)合的DNA結(jié)合域（DBD）［5］。ERα的活性主要受雌激素調(diào)節(jié)，ERα的激活促進了許多調(diào)控致癌過程的靶基因的轉(zhuǎn)錄，包括細胞周期、細胞增殖和上皮/管腔細胞分化等過程［6，7］。內(nèi)分泌治療通過拮抗雌激素與ERα的結(jié)合，下調(diào)ERα或抑制雌激素的產(chǎn)生來抑制ERα信號通路的激活，結(jié)合手術(shù)和放射療法，顯著改善了ER陽性乳腺癌患者的臨床結(jié)局［8］。然而，仍有部分ER陽性的乳腺癌患者在內(nèi)分泌治療過程中產(chǎn)生獲得性耐藥，出現(xiàn)復(fù)發(fā)或轉(zhuǎn)移［9-11］。因此，迫切需要可以預(yù)測ER陽性乳腺癌患者預(yù)后的生物標(biāo)志物。

微小RNA（miRNA）是一類調(diào)節(jié)基因表達的短非編碼單鏈RNA，長度約為20～24個核苷酸，可通過與靶基因的3'-非翻譯區(qū)（3'-UTR）堿基配對抑制信使RNA（mRNA）的翻譯或降解mRNA從而調(diào)節(jié)靶基因的表達［12］。許多證據(jù)表明，miRNA在細胞增殖分化、血管生成和代謝等各種生物過程發(fā)揮重要作用［13-15］。多種惡性腫瘤發(fā)生發(fā)展也與miRNA的表達失調(diào)相關(guān)［16，17］，例如結(jié)腸癌［18，19］、肺癌［20］、胃癌［21］、鼻咽癌［22，23］等。有研究顯示，乳腺癌中也存在多種miRNA表達失調(diào)，且部分與ERα狀態(tài)有關(guān)［24］。例如，miRNA-519a的高表達與ER陽性乳腺癌患者較差的生存率相關(guān)，且miRNA-519a上調(diào)可導(dǎo)致ER陽性乳腺癌患者對他莫昔芬的耐藥。

因此，本研究旨在通過分析TCGA數(shù)據(jù)庫的miRNA測序數(shù)據(jù)構(gòu)建miRNA-臨床的預(yù)測ER陽性乳腺癌患者預(yù)后的模型，用于有效選擇高?；颊吆皖A(yù)測ER陽性乳腺癌患者的3年和5年生存率，從而有利于制定更有效的個體化治療決策。

1 資料與方法

1.1 下載與處理原始數(shù)據(jù)2019年11月15日從TCGA官方網(wǎng)站獲取乳腺癌患者miRNA表達譜和臨床數(shù)據(jù)，miRNA表達譜數(shù)據(jù)中總共包括1 096例乳腺癌組織和104例正常乳腺組織，臨床數(shù)據(jù)中包括1 174例乳腺癌患者臨床病理信息。納入標(biāo)準(zhǔn)：（1）病理學(xué)確診為乳腺癌。（2）既有miRNA的樣品測序數(shù)據(jù)又有完整的臨床預(yù)后資料（生存狀態(tài)及生存時間）。排除標(biāo)準(zhǔn)：（1）雌激素受體表達情況不明確。（2）生存時間未知或生存時間≤0 d。最后共納入730例乳腺癌患者，其中555例雌激素受體陽性患者，175例雌激素受體陰性患者。從GEO中下載數(shù)據(jù)集GSE37405（GPL 13703），該數(shù)據(jù)集包含60例ER陽性乳腺癌患者的miRNA測序數(shù)據(jù)，患者均接受他莫昔芬治療，其中乳腺癌復(fù)發(fā)患者和未復(fù)發(fā)患者各30例。

1.2 篩選ER陽性乳腺癌特異性差異miRNA以|log2FC|≥2，F(xiàn)DR＜0.05為篩選條件，使用R3.6.0軟件中的edgeR包分別篩選出ER陽性乳腺癌差異表達的miRNA以及ER陰性乳腺癌中差異表達的miRNA，并使用ggplot2和pheatmap包繪制火山圖和熱圖。使用VennDiagram包繪制韋恩圖篩選出ER陽性乳腺癌特異性上調(diào)和下調(diào)的差異miRNA。

1.3 構(gòu)建風(fēng)險評分公式和miRNA-臨床預(yù)后模型篩選得到的ER陽性乳腺癌特異性差異miRNA進行Lasso回歸分析，篩選出與生存相關(guān)miRNA，根據(jù)每個miRNA的回歸系數(shù)，構(gòu)建基于miRNA表達的風(fēng)險評分公式。繪制受試者工作曲線（ROC）評估m(xù)iRNA風(fēng)險評分的預(yù)測性能，根據(jù)ROC曲線選擇最佳臨界風(fēng)險評分將ER陽性乳腺癌患者分為低風(fēng)險組和高風(fēng)險組，比較兩組患者之間生存的差異及miRNA的表達情況。進一步通過單因素和多因素Cox回歸分析結(jié)合臨床病理因素和miRNA風(fēng)險評分構(gòu)建預(yù)測ER陽性乳腺癌患者預(yù)后的預(yù)后模型，并使用rms包繪制列線圖。

1.4 評估風(fēng)險評分公式和miRNA-臨床預(yù)后模型的預(yù)測性能繪制風(fēng)險評分公式和miRNA-臨床預(yù)后模型的ROC曲線并計算曲線下面積（AUC）評估模型的預(yù)測能力。繪制校準(zhǔn)曲線評估模型的準(zhǔn)確性。

1.5 靶基因預(yù)測和功能富集分析通過Starbase預(yù)測miRNA的候選靶基因。選擇至少3個數(shù)據(jù)庫重疊的基因作為miRNA的靶基因。使用Cytoscape可視化miRNA與靶基因之間的相互作用網(wǎng)絡(luò)。最后，使用DAVID 6.8進行靶基因的GO功能富集和KEGG通路富集分析。

1.6 統(tǒng)計學(xué)分析正態(tài)分布的數(shù)值變量采用均數(shù)±標(biāo)準(zhǔn)差描述，t檢驗用于組間比較。偏態(tài)分布的數(shù)值變量采用中位數(shù)（四分位數(shù)間距）描述，Mann-WhitneyU檢驗用于組間比較。分類變量采用頻數(shù)（構(gòu)成比）描述，χ2檢驗或Fisher精確檢驗用于組間比較。Lasso回歸分析用于篩選與預(yù)后相關(guān)miRNA并構(gòu)建風(fēng)險評分公式。使用單因素Cox回歸和多因素Cox回歸分析構(gòu)建miRNA-臨床預(yù)測模型，并使用R中rms包繪制列線圖。survival ROC包繪制ROC曲線評估了列線圖預(yù)測能力，校準(zhǔn)曲線用于評估模型的準(zhǔn)確性。使用Kaplan-Meier方法繪制生存曲線，并使用Log-rank檢驗比較生存率。P＜0.05被認為具有統(tǒng)計學(xué)意義。應(yīng)用R3.6.0進行數(shù)據(jù)分析。

2 結(jié)果

2.1 患者的基線特征共555名ER陽性乳腺癌患者納入研究，用隨機抽樣法按1∶1比例分為訓(xùn)練集278例和驗證集277例，訓(xùn)練集用于構(gòu)建預(yù)測模型，驗證集用于驗證預(yù)測模型效能。兩組基本情況比較差異均無統(tǒng)計學(xué)意義（P＞0.05），兩組具有可比性。所有患者的5年總生存率為91.7%。

2.2 篩選ER陽性乳腺癌與生存相關(guān)的特異性差異miRNA從TCGA數(shù)據(jù)庫中獲得了乳腺癌的miRNA測序數(shù)據(jù)，包括555例ER陽性乳腺癌，175例ER陰性乳腺癌和104正常組織樣本。使用R3.6.0軟件中的edgeR包分別篩選出ER陽性乳腺癌差異表達的miRNA和ER陰性乳腺癌差異表達的miRNA（|log2FC|≥2，錯誤發(fā)現(xiàn)率FDR＜0.05）。根據(jù)兩組上調(diào)的miRNA和下調(diào)的miRNA分別繪制韋恩圖，選擇ER陽性乳腺癌非交集部分miRNA為ER陽性乳腺癌特異性miRNA，其中5個miRNA在ER陽性乳腺癌特異性上調(diào)，17個miRNA在ER陽性乳腺癌特異性下調(diào)。為了進一步篩選與生存相關(guān)的miRNA，在訓(xùn)練集中對22個miRNA進行了Lasso回歸，使用交叉驗證以建立風(fēng)險模型，最終篩選出4個與生存相關(guān)的miRNA，分別為miR-331、miR-615、miR-653、miR-887（見圖1）。

圖1 ER陽性與ER陰性乳腺癌下調(diào)的miRNA的韋恩圖（A）和上調(diào)的miRNA的韋恩圖（B）以及Lasso回歸分析篩選miRNA(C)和交叉驗證結(jié)果(D)

2.3 風(fēng)險評分公式的建立和miRNA-臨床預(yù)后模型的構(gòu)建根據(jù)Lasso回歸得到相應(yīng)的系數(shù)，建立立風(fēng)險評分公式：風(fēng)險評分=0.098×miR-331exp+0.122×miR-615exp+0.102×miR-653exp+0.113×miR-887exp。根據(jù)此方程，計算每個患者的風(fēng)險評分，并根據(jù)ROC曲線的最佳臨界風(fēng)險評分（風(fēng)險得分為1.757）將訓(xùn)練集患者分為低風(fēng)險組（n=189）的和高風(fēng)險組（n=89），驗證集也分為低風(fēng)險組（n=73）和高風(fēng)險組（n=204）。Kaplan-Meier生存分析表明，在訓(xùn)練集（P=0.000 26）和驗證集（P＜0.000 1）中，高風(fēng)險組的預(yù)后比低風(fēng)險組差（圖2）。隨后，將風(fēng)險評分作為變量，結(jié)合其他臨床病理特征進行單變量和多變量Cox回歸分析來識別與預(yù)后相關(guān)的危險因素，結(jié)果顯示TNM分期（HR=5.068，95%CI：1.799～14.278，P=0.002）、風(fēng)險評分（HR=1.191，95%CI：1.030～1.377，P=0.018）和術(shù)后是否放療（HR=0.411，95%CI：0.176～0.965，P=0.041）是預(yù)后的獨立危險因素（表1）。將上述危險因素用于構(gòu)建miRNA-臨床預(yù)后模型并繪制列線圖，在列線圖中，每個變量的得分相加計算總分，通過總分評估患者的3年和5年生存率（圖3）。

圖2 訓(xùn)練集（A）和驗證集（B）中高低風(fēng)險評分組患者的生存曲線

表1 影響ER陽性乳腺癌預(yù)后的臨床病理因素的單因素和多因素Cox回歸分析

圖3 基于4個miRNA表達的預(yù)測ER陽性乳腺癌患者3年和5年生存率的列線圖

2.4 評估風(fēng)險評分和miRNA-臨床預(yù)后模型的預(yù)測效能繪制風(fēng)險評分和基于miRNA-臨床預(yù)后模型的ROC曲線并計算AUC值以評估模型的預(yù)測能力。在訓(xùn)練集和驗證集中，miRNA-臨床預(yù)后模型預(yù)測3年生存率的AUC分別為0.768和0.909，預(yù)測5年生存率的AUC分別為0.849和0.860（圖4）。

圖4 訓(xùn)練集（A）和驗證集（B）評估m(xù)iRNA-臨床預(yù)后模型3年和5年預(yù)測能力的ROC曲線

2.5 miRNA在他莫昔芬治療患者中的表達從基因公共表達數(shù)據(jù)庫（GEO）中下載數(shù)據(jù)集GSE37405（GPL 13703），該數(shù)據(jù)集包含30例他莫昔芬治療后未復(fù)發(fā)患者及30例他莫昔芬治療后復(fù)發(fā)患者的miRNA測序數(shù)據(jù)，分析4個miRNA在兩組間的表達水平，結(jié)果顯示，與未復(fù)發(fā)組相比，復(fù)發(fā)患者miR-331的表達水平顯著升高，而兩組之間的miR-615、miR-653和miR-887的水平?jīng)]有明顯差異（圖5）。miR-331可能和ER陽性乳腺癌他莫昔芬耐藥相關(guān)。

圖5 他莫昔芬治療后未復(fù)發(fā)組和復(fù)發(fā)組中4個miRNA的表達水平

2.6 mi-331的靶基因的GO功能富集分析和KEGG通路富集分析GO功能富集分析表明，miR-331的靶基因細胞組成（圖6A）主要富集在細胞核和核質(zhì)，生物過程（圖6B）主要參與在RNA聚合酶Ⅱ轉(zhuǎn)錄的終止、核轉(zhuǎn)錄mRNA poly（A）尾巴的縮短和內(nèi)皮細胞增殖的調(diào)控，分子功能（圖6C）主要參與蛋白質(zhì)結(jié)合。KEGG通路分析（圖6D）結(jié)果表明，miR-331的靶基因主要參與mRNA監(jiān)視途徑，Rap1信號通路和MAPK信號通路相關(guān)。

圖6 miR-331的GO功能富集和KEGG通路富集

3 討論

高通量測序結(jié)合生物信息學(xué)對基因組學(xué)數(shù)據(jù)進行分析，可以從分子角度探索與惡性腫瘤診斷，治療與預(yù)后密切相關(guān)的標(biāo)志物。通過分子標(biāo)志物對患者生存率的預(yù)測有助于制定個體化的臨床決策。許多研究［25-29］表明，miRNA與ER陽性乳腺癌的發(fā)生和進展有關(guān)。Mulrane等［30］研究發(fā)現(xiàn)miR-187在乳腺癌中的高表達可能與更具侵略性的表型相關(guān)，且miR-187與乳腺癌預(yù)后不良相關(guān)，是預(yù)測患者生存率的獨立預(yù)測因子。Okuda等［31］研究表明，miR-7可以通過下調(diào)KLF4基因的表達進而抑制乳腺癌細胞的遠處轉(zhuǎn)移。有研究者［32］發(fā)現(xiàn)ER陽性患者的腫瘤中miRNA-375水平顯著高于ER陰性乳腺癌患者，且miRNA-375在他莫昔芬耐藥細胞中下調(diào)，其重新表達可以抑制MTDH基因表達從而恢復(fù)腫瘤細胞對他莫昔芬敏感性，同時也抑制腫瘤細胞侵襲和逆轉(zhuǎn)EMT樣特性。此外，Gong等［33］建立了一個基于miRNA的模型來預(yù)測激素受體陽性HER2陰性乳腺癌患者的復(fù)發(fā)風(fēng)險。因此，在多種癌癥中基于miRNA的預(yù)測模型對預(yù)測患者的生存或復(fù)發(fā)具有重要意義。傳統(tǒng)的生物學(xué)方法驗證miRNA和患者的預(yù)后之間的關(guān)聯(lián)既耗時又昂貴，隨著生物學(xué)數(shù)據(jù)的積累，我們可以利用大量數(shù)據(jù)分析miRNA和腫瘤的相關(guān)性并開發(fā)預(yù)測模型來預(yù)測腫瘤患者的預(yù)后及治療反應(yīng)，通過另一數(shù)據(jù)集的驗證檢驗?zāi)Ｐ偷目煽啃?，建立基于miRNA的預(yù)測模型既對個體化治療方案的選擇具有重要的參考價值。也可為miRNA與腫瘤的相關(guān)性進一步的實驗驗證提供理論依據(jù)。

在本研究中，我們從TCGA數(shù)據(jù)庫下載miRNA表達數(shù)據(jù)，通過差異表達分析和Lasso回歸發(fā)現(xiàn)了4種ER陽性乳腺癌特異性表達上調(diào)的miRNA（miR-331、miR-615、miR-653、miR-887）并建立了miRNA表達的風(fēng)險評分公式。根據(jù)ROC曲線的最佳臨界風(fēng)險評分將患者分為低風(fēng)險組和高風(fēng)險組，風(fēng)險評分與ER陽性乳腺癌患者的總生存率顯著相關(guān)。通過Cox回歸分析風(fēng)險評分和臨床病理因素與預(yù)后的關(guān)系，構(gòu)建了miRNA-臨床預(yù)測模型用于預(yù)測ER陽性乳腺癌患者的預(yù)后。該模型預(yù)測3年和5年生存率的AUC值分別為0.768和0.849，表明其具有良好的預(yù)測性能。模型中包含的4個miRNA中，部分miRNA已有試驗證實與乳腺癌相關(guān)。Jiang等［34］發(fā)現(xiàn)miR-331在乳腺癌細胞中上調(diào)，從而促進乳腺癌細胞的增殖、遷移和侵襲，miR-331的過表達與淋巴結(jié)轉(zhuǎn)移、TNM分期和預(yù)后不良有關(guān)。Lei等［35］發(fā)現(xiàn)miR-615-3p在乳腺癌細胞和組織中特別是在轉(zhuǎn)移的乳腺癌細胞和組織中顯著升高，miR-615-3p通過抑制PICK 1基因表達促進了腫瘤細胞的EMT過程。因此，我們構(gòu)建的miRNA-臨床的列線圖可能是ER陽性乳腺癌患者生存預(yù)測的重要工具，有助于制定個性化治療策略。隨后我們利用GEO數(shù)據(jù)庫中GSE37405數(shù)據(jù)集分析了4個miRNA在他莫昔芬治療后復(fù)發(fā)和未復(fù)發(fā)兩組之間的表達，發(fā)現(xiàn)miR-331在他莫昔芬治療后復(fù)發(fā)組中表達明顯高于未復(fù)發(fā)組，提示miR-331可能與ER陽性乳腺癌患者發(fā)生他莫昔芬耐藥相關(guān)。進一步進行了miR-331靶基因的GO和KEGG途徑的分析，結(jié)果顯示miR-331可能在蛋白質(zhì)結(jié)合、轉(zhuǎn)錄、mRNA監(jiān)視通路、MAPK信號通路中起關(guān)鍵作用。

綜上所述，我們構(gòu)建的預(yù)測ER陽性乳腺癌患者生存的miRNA-臨床風(fēng)險模型，該模型結(jié)合了4個miRNA的表達情況和臨床危險因素，可預(yù)測ER陽性乳腺癌患者的3年和5年生存率，且建立的風(fēng)險評分公式可有效地識別ER陽性乳腺癌中的高?；颊摺Ｒ虼?，該預(yù)測模型對制定個性化治療決策有重要意義。但在該模型廣泛應(yīng)用于臨床實踐之前，仍需要進行多中心、大規(guī)模、前瞻性研究以驗證該預(yù)測工具。