李元幸，牛曉辰，常晶晶，雒海瑕，王偉

宮頸癌是第二大導(dǎo)致女性死亡的癌癥[1]，宮頸鱗狀細(xì)胞癌（cervical squamous carcinoma，CESC）是最常見的臨床病理類型。盡管近年來宮頸癌的預(yù)防和診治都有了長足發(fā)展，但短期內(nèi)預(yù)防措施覆蓋面有限，診治水平也受地域局限[2]，宮頸癌的預(yù)后仍然不容樂觀，其5 年總生存率反呈下降趨勢[3]。因此，探索宮頸癌的分子機(jī)制，尋找與宮頸癌預(yù)后相關(guān)的潛在生物標(biāo)志物和治療靶標(biāo)仍然是臨床的重要工作目標(biāo)。

高通量測序的快速發(fā)展使各類基因生物標(biāo)志物應(yīng)用于臨床成為可能。微小RNA（microRNA，miRNA）是一種小的非編碼RNA，可以調(diào)節(jié)多種信使RNA（messenger RNA，mRNA）的表達(dá)，從而影響各種細(xì)胞生物學(xué)行為，包括轉(zhuǎn)錄、修飾、染色體重塑和信號(hào)傳導(dǎo)等。越來越多的證據(jù)表明，miRNA 和mRNA 的異常表達(dá)與多種癌癥的發(fā)生、發(fā)展有關(guān)，并且可能成為癌癥的預(yù)后生物標(biāo)志物和治療靶標(biāo)[4-5]。然而，目前miRNA 與mRNA 在宮頸鱗癌中的研究還較為有限。隨著生物信息學(xué)研究技術(shù)的發(fā)展，現(xiàn)可通過網(wǎng)絡(luò)公開的癌癥數(shù)據(jù)庫信息來輔助篩選并驗(yàn)證各類調(diào)控分子在癌癥中的遺傳信息[6]。癌癥和腫瘤基因圖譜（The Cancer Genome Atlas，TCGA）計(jì)劃是一種涵蓋多中心腫瘤基因數(shù)據(jù)的大型在線數(shù)據(jù)庫。本研究通過TCGA 數(shù)據(jù)庫中CESC 的相關(guān)樣本數(shù)據(jù)，探討miRNA 與mRNA 的相互作用關(guān)系，即miRNAmRNA 關(guān)系對，為將來CESC 分子機(jī)制、預(yù)后標(biāo)志物和治療靶點(diǎn)的研究奠定基礎(chǔ)。

1 資料與方法

1.1 數(shù)據(jù)獲取與整理通過TCGA 數(shù)據(jù)庫（https://portal.gdc.cancer.gov/）下載CESC 的3 級mRNA 與miRNA 測序信息和臨床數(shù)據(jù)。利用perl 語言（perl 5.30.2，https://www.perl.org/）完成數(shù)據(jù)的提取與合并，整理為R 語言可識(shí)別的矩陣文件。

1.2 差異分析使用R project 3.6.3 的edgeR 包對mRNA 與miRNA 表達(dá)數(shù)據(jù)分別進(jìn)行差異分析。利用pheatmap 包分別繪制表達(dá)上調(diào)或下調(diào)的前20 個(gè)基因的熱圖。

1.3 miRNA 預(yù)后模型構(gòu)建與評價(jià)使用survival 包進(jìn)行miRNA 表達(dá)和CESC 預(yù)后的單因素與多因素Cox回歸分析，篩選有統(tǒng)計(jì)學(xué)差異的miRNA 并構(gòu)建miRNA 預(yù)后風(fēng)險(xiǎn)評分模型，計(jì)算公式為：風(fēng)險(xiǎn)值（risk score）=風(fēng)險(xiǎn)基因表達(dá)量1×coef1＋風(fēng)險(xiǎn)基因表達(dá)量2×coef2＋……＋風(fēng)險(xiǎn)基因表達(dá)量n×coefn（coef 為風(fēng)險(xiǎn)系數(shù)）。

用survminer 與survival ROC 包進(jìn)行Kaplan-Meier 生存分析并繪制生存曲線，計(jì)算受試者工作特征（ROC）曲線下面積（AUC）進(jìn)行模型評價(jià)。根據(jù)預(yù)后模型計(jì)算所有樣本的風(fēng)險(xiǎn)值，由低到高排序后根據(jù)中位數(shù)將患者分為高、低風(fēng)險(xiǎn)2 組，輸出2 組中模型的生存曲線，驗(yàn)證模型的準(zhǔn)確性。利用survival 包對所有樣本的生存時(shí)間、生存狀態(tài)、年齡、世界衛(wèi)生組織（WHO）分級、TNM 分期及本研究構(gòu)建的預(yù)測模型風(fēng)險(xiǎn)評分進(jìn)行單因素與多因素的獨(dú)立預(yù)后分析，評估構(gòu)建預(yù)測模型的風(fēng)險(xiǎn)值是否能作為患者的獨(dú)立預(yù)后因子（評價(jià)標(biāo)準(zhǔn)為P.adj＜0.01）。

1.4 miRNA-mRNA調(diào)控網(wǎng)絡(luò)的構(gòu)建利用3 種生物信息學(xué)算法（miRDB、miRTarBase 與TargetScan 數(shù)據(jù)庫）對納入模型的miRNA 進(jìn)行靶基因預(yù)測，篩選同時(shí)被2 個(gè)及以上數(shù)據(jù)庫收錄的靶基因，并與1.2 中差異表達(dá)的mRNA 取交集，明確其靶關(guān)系以及上下調(diào)表達(dá)關(guān)系，并通過Cytoscape 3.7.2 將miRNA-mRNA調(diào)控網(wǎng)絡(luò)可視化。

1.5 miRNA-mRNA 關(guān)系對的篩選用同樣方法對mRNA數(shù)據(jù)進(jìn)行生存分析，篩選與患者生存預(yù)后密切相關(guān)的mRNA。為了進(jìn)一步篩選miRNA-mRNA調(diào)控關(guān)系對，根據(jù)miRNA 對mRNA 的功能發(fā)揮抑制作用這一生物學(xué)基礎(chǔ)[7]制定篩選標(biāo)準(zhǔn)：“高表達(dá)、生存率低”的miRNA 對應(yīng)的mRNA“低表達(dá)、生存率低”；“低表達(dá)、生存率低”的miRNA 對應(yīng)的mRNA“高表達(dá)、生存率低”。

1.6 功能富集分析為了進(jìn)一步明確miRNA-mRNA調(diào)控網(wǎng)絡(luò)所參與的主要功能，對其中的mRNA 進(jìn)行基因本體功能（gene ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析。以P＜0.05 且偽發(fā)現(xiàn)率＜0.05 作為富集的判斷標(biāo)準(zhǔn)。

2 結(jié)果

2.1 樣本特征和RNA 差異分析從TCGA 數(shù)據(jù)庫進(jìn)行檢索分析后，共納入255 例CESC 患者的腫瘤組織測序信息和2 例正常宮頸組織測序信息。與正常宮頸組織相比，CESC 腫瘤組織共有3 167 個(gè)差異表達(dá)的mRNA（上調(diào)1 030 個(gè)、下調(diào)2 137 個(gè)），115 個(gè)差異表達(dá)的miRNA（上調(diào)39 個(gè)、下調(diào)76 個(gè)）。上調(diào)與下調(diào)最顯著的前20 個(gè)mRNA 與miRNA 見圖1。

圖1 CESC 中差異表達(dá)的RNA 聚類熱圖

2.2 miRNA 預(yù)后模型的構(gòu)建對網(wǎng)絡(luò)中的miRNA分別進(jìn)行單因素與多因素Cox回歸分析，最終確定4 個(gè)miRNA 可納入預(yù)后風(fēng)險(xiǎn)評分模型，分別是hsamiR-505-5p、hsa-miR-142-3p、hsa-miR-3613-5p和hsa-miR-532-5p。預(yù)后風(fēng)險(xiǎn)評分模型計(jì)算公式為：風(fēng)險(xiǎn)值（risk score）=hsa-miR-505-5p 表達(dá)量×（-0.031 21）+hsa-miR-142-3p 表達(dá)量×（-0.000 06）+hsa-miR-3613-5p 表達(dá)量×0.00424+hsa-miR-532-5p 表達(dá)量×（-0.000 28）。

2.3 miRNA 預(yù)后模型評價(jià)hsa-miR-505-5p、hsa-miR-142-3p 與hsa-miR-532-5p 的生存曲線顯示，低表達(dá)組生存率低（見圖2A、B、C）；hsa-miR-3613-5p 生存曲線顯示，高表達(dá)組生存率低（見圖2D）。根據(jù)預(yù)后模型計(jì)算所有樣本的風(fēng)險(xiǎn)值，根據(jù)中位數(shù)分為高風(fēng)險(xiǎn)組與低風(fēng)險(xiǎn)組，高風(fēng)險(xiǎn)組生存率低于低風(fēng)險(xiǎn)組，證實(shí)模型準(zhǔn)確（見圖2E）。進(jìn)一步對所有樣本的生存時(shí)間、生存狀態(tài)、年齡、WHO 分級、TNM 分期及本研究構(gòu)建的預(yù)測模型風(fēng)險(xiǎn)評分進(jìn)行單因素與多因素預(yù)后分析并繪制ROC 曲線，風(fēng)險(xiǎn)值是患者生存預(yù)后的獨(dú)立預(yù)測因子，且隨著風(fēng)險(xiǎn)值增大，患者生存率逐漸降低，其AUC 為0.872（見圖2F），進(jìn)一步證明了此miRNA 預(yù)后模型的可靠性。

圖2 納入預(yù)后模型的miRNA 不同表達(dá)者的生存曲線及預(yù)后模型評價(jià)

2.4 miRNA-mRNA調(diào)控網(wǎng)絡(luò)的構(gòu)建利用miRDB、miRTarBase 與TargetScan 數(shù)據(jù)庫對納入模型的4 個(gè)miRNA 進(jìn)行靶基因預(yù)測，共預(yù)測到23 469個(gè)靶基因，其中同時(shí)被2 個(gè)及以上數(shù)據(jù)庫收錄的靶基因有1 674 個(gè)，與2.1 中差異表達(dá)的mRNA 取交集，共得到111 個(gè)靶基因，并由此構(gòu)建了miRNAmRNA調(diào)控網(wǎng)絡(luò)，見圖3。

圖3 CESC 的miRNA-mRNA調(diào)控網(wǎng)絡(luò)

2.5 miRNA-mRNA 關(guān)系對的篩選生存分析篩選出9 個(gè)與CESC 預(yù)后相關(guān)的差異表達(dá)mRNA，見圖4。分別是Ⅰ型血小板結(jié)合蛋白基序的解聚蛋白樣金屬蛋白酶3（ADAMTS3）、載脂蛋白B mRNA 編輯酶催化多肽3B（APOBEC3B）、染色質(zhì)結(jié)構(gòu)域蛋白8（CBX8）、整合素ɑ8（ITGA8）、微管剪切蛋白60 亞基Ⅰ型類似物（KATNAL1）、痩素受體疊加轉(zhuǎn)錄蛋白（LEPROT）、蛋白酪氨酸磷酸酶受體B（PTPRB）、G蛋白信號(hào)傳導(dǎo)調(diào)節(jié)因子5（RGS5）和SEC23 同源物A（SEC23A）。根據(jù)miRNA-mRNA 關(guān)系對篩選標(biāo)準(zhǔn)，最終有4 個(gè)關(guān)系對具有重要生物學(xué)意義，通過miRNA 的表達(dá)變化調(diào)控下游mRNA 水平的變化，進(jìn)而影響患者的生存預(yù)后情況，見表1。

圖4 與CESC 預(yù)后相關(guān)的9 個(gè)mRNA 不同表達(dá)者的生存曲線

表1 miRNA-mRNA調(diào)控關(guān)系對

2.6 功能富集對miRNA-mRNA調(diào)控網(wǎng)絡(luò)中的111 個(gè)mRNA 所屬基因進(jìn)行GO 富集分析，共富集到415 個(gè)GO 功能詞條，包括374 個(gè)生物學(xué)過程（biological process，BP）、23 個(gè) 細(xì) 胞 成 分（cellular component，CC）和 18 個(gè) 分 子 功 能（molecular function，MF），主要富集在BP 中，如蛋白定位調(diào)控細(xì)胞外周環(huán)境、原胚層形成、中胚層形態(tài)形成與分化。KEGG 通路富集分析顯示，差異表達(dá)的mRNA 主要富集在磷脂酰肌醇3 激酶-蛋白激酶B（PI3KAkt）信號(hào)通路與人乳頭瘤病毒（HPV）感染。按照P值由小到大分別篩選了各組前10 個(gè)GO 功能詞條和通路并進(jìn)行可視化。見圖5。

圖5 mRNA 功能富集氣泡圖

3 討論

CESC 現(xiàn)有的預(yù)后評價(jià)指標(biāo)主要依據(jù)腫瘤的大小、分期與淋巴結(jié)轉(zhuǎn)移情況，醫(yī)師的主觀判斷與診治水平對其預(yù)后評估影響較大[8]，缺乏可靠的特異性生物標(biāo)志物。尋找客觀的CESC 預(yù)后因子，對監(jiān)測腫瘤復(fù)發(fā)、提高患者生存質(zhì)量、CESC 的機(jī)制研究和靶向藥物研發(fā)都有著重要意義。miRNA-mRNA調(diào)控關(guān)系對是非編碼RNA 領(lǐng)域研究的經(jīng)典線性調(diào)控模式，在癌癥預(yù)后預(yù)測中表現(xiàn)出巨大潛力。因此，本研究基于TCGA 數(shù)據(jù)庫，利用基因表達(dá)差異分析、生存分析和預(yù)后模型等方法，最終得到與CESC 預(yù)后緊密相關(guān)的4 個(gè)miRNA-mRNA調(diào)控關(guān)系對，以期為臨床診療和研究提供依據(jù)。

上述miRNA-mRNA調(diào)控關(guān)系對還未在CESC中展開研究。miRNA 是重要的靶基因調(diào)控因子，通過切斷或阻止翻譯而下調(diào)mRNA 表達(dá)，其作用近乎貫穿腫瘤發(fā)生、發(fā)展全過程。但目前hsa-miR-505-5p 和hsa-miR-142-3p 與CESC 的相關(guān)研究較少。研究報(bào)道m(xù)iR-505-5p 與CESC 轉(zhuǎn)移相關(guān)的臨床病理特征密切相關(guān)，還發(fā)現(xiàn)miR-505-5p 過表達(dá)可以抑制宮頸癌細(xì)胞的轉(zhuǎn)移和上皮-間質(zhì)轉(zhuǎn)化過程，從而抑制宮頸癌的進(jìn)展和轉(zhuǎn)移[9]。另外，也有研究證實(shí)hsa-miR-142-3p 在鼻咽癌、食管鱗狀細(xì)胞癌組織中差異表達(dá)，可作為潛在的預(yù)后標(biāo)志物[10-11]。mRNA 在基因編碼蛋白質(zhì)的轉(zhuǎn)錄過程中至關(guān)重要，在腫瘤發(fā)生、發(fā)展中起主要作用。研究顯示，CBX8可以作為肝癌、結(jié)腸癌的預(yù)后因子[12-13]。Zhang 等[14]通過全基因組測序發(fā)現(xiàn)，CDX8 在食管鱗狀細(xì)胞癌組織中局灶性擴(kuò)增現(xiàn)象，可以作為食管鱗狀細(xì)胞癌早期診斷標(biāo)志物和治療靶標(biāo)。ADAMTS3 在細(xì)胞外基質(zhì)組裝與降解中有重要作用，已經(jīng)被證明廣泛參與血管和淋巴管生成、細(xì)胞遷移等生物學(xué)過程[15-16]，因其與腫瘤復(fù)發(fā)、轉(zhuǎn)移密切相關(guān)，有作為腫瘤預(yù)后因子的潛力。有研究報(bào)道PTPRB 通過與受體酪氨酸激酶的結(jié)合和去磷酸化，從而抑制腫瘤細(xì)胞增殖，可作為大腸癌、非小細(xì)胞肺癌的預(yù)后因子[17-18]。SEC23A 是SEC23 亞家族的成員，可通過影響外殼蛋白復(fù)合物Ⅱ（COPⅡ）的裝配來調(diào)節(jié)腫瘤微環(huán)境從而影響預(yù)后，是皮膚黑色素瘤的預(yù)后標(biāo)志物[19]。由上可見，CBX8、ADAMTS3、PTPRB、SEC23A 均與癌癥的發(fā)生、發(fā)展和預(yù)后密切相關(guān)，但這些mRNA 在CESC 中的研究尚屬空白，有待進(jìn)一步深入探討。而且，將差異表達(dá)的mRNA 與miRNA 數(shù)據(jù)相結(jié)合，有利于提供更多信息，進(jìn)一步提高對CESC預(yù)后的預(yù)測準(zhǔn)確性。

綜合生物學(xué)進(jìn)程的富集分析顯示，miRNAmRNA調(diào)控關(guān)系對主要參與以下生物學(xué)功能：蛋白定位調(diào)控細(xì)胞外周環(huán)境、原胚層形成、中胚層形態(tài)形成與分化。下游關(guān)鍵通路主要富集于PI3K-Akt 通路與HPV 感染。PI3K-Akt 通路是腫瘤領(lǐng)域經(jīng)典的致癌信號(hào)通路。在信號(hào)傳導(dǎo)復(fù)合物的作用下，PI3K 活化并發(fā)生級聯(lián)反應(yīng)，導(dǎo)致原癌基因Akt 的異?；罨瑥亩偈拱┌Y發(fā)生[20-21]。HPV 也已被證明是宮頸癌發(fā)生、發(fā)展的主要因素[22]。上述富集通路再次從功能層面佐證了本研究篩選出的與CESC 預(yù)后相關(guān)的miRNA-mRNA調(diào)控關(guān)系對的可靠性。后續(xù)研究也可借鑒本研究篩選出的富集通路，作為miRNA-mRNA調(diào)控關(guān)系對的下游調(diào)控通路進(jìn)行深入探討。

綜上，本研究基于miRNA-mRNA 關(guān)系對和Cox回歸模型，利用TCGA 數(shù)據(jù)庫綜合分析了與CESC預(yù)后緊密相關(guān)的生物標(biāo)志物。然而上述生物標(biāo)志物在CESC 中的研究還涉及較少，有待進(jìn)一步驗(yàn)證其表達(dá)及功能。本研究初步篩選了CESC 的預(yù)后標(biāo)志物，對CESC 的機(jī)制研究、靶向藥物研發(fā)以及監(jiān)測CESC 患者預(yù)后都有著重要意義。隨著更多功能性RNA 的發(fā)現(xiàn)，宮頸癌預(yù)后相關(guān)的基因研究范圍有待擴(kuò)展，期待宮頸癌遺傳基因圖譜的進(jìn)一步揭密。