卵巢癌耐藥相關(guān)的miRNAs篩選及其預(yù)后意義研究

朱軼昊,楊奕辰,馬琳琳*

(1.上海健康醫(yī)學(xué)院醫(yī)學(xué)技術(shù)學(xué)院,中國上海 201318;2.上海杉達(dá)學(xué)院信息科學(xué)與技術(shù)學(xué)院,中國上海 201209)

卵巢癌是最致命的婦科惡性腫瘤,5年生存率僅為46%[1]。缺乏早期診斷的特異性生物標(biāo)志物和對(duì)化療藥物的耐藥是患者不良預(yù)后的重要原因[2]。雖然80%的卵巢癌患者最初對(duì)紫杉醇/卡鉑聯(lián)合化療敏感,但約有一半以上的患者在5年內(nèi)因獲得性耐藥而復(fù)發(fā)[3]。因此,識(shí)別和開發(fā)針對(duì)化療耐藥卵巢癌患者的特異性生物標(biāo)志物和潛在靶點(diǎn)至關(guān)重要。

微RNA(microRNA,miRNA)是一類長度約22個(gè)核苷酸的非編碼RNA分子,可以在轉(zhuǎn)錄后水平上調(diào)控基因的翻譯表達(dá)或mRNA的降解[4]。miRNA通過靶向mRNA的3′-UTR,實(shí)現(xiàn)對(duì)靶mRNA的降解或翻譯抑制[5～6],從而調(diào)控基因表達(dá)。越來越多的證據(jù)表明,miRNA能夠調(diào)節(jié)細(xì)胞的增殖、分化和應(yīng)激反應(yīng)。miRNA功能障礙與人類疾病,尤其是癌癥的形成和發(fā)展有關(guān)[7～8]。許多miRNAs在卵巢癌中異常表達(dá),對(duì)癌癥的發(fā)生發(fā)展及預(yù)后有重要影響。例如:miR-450a是卵巢癌中顯著下調(diào)的miRNAs之一,其過表達(dá)可抑制與上皮-間質(zhì)轉(zhuǎn)化相關(guān)的多個(gè)基因,減少腫瘤細(xì)胞的遷移和侵襲,增加腫瘤細(xì)胞的凋亡,在卵巢腫瘤異種移植模型中,miR-450a過表達(dá)可顯著抑制腫瘤生長[9];miR-27a被發(fā)現(xiàn)在卵巢癌組織和細(xì)胞系中顯著升高,抑制miR-27a的表達(dá)可抑制SK-OV-3和OVACAR-3卵巢癌細(xì)胞系的細(xì)胞周期進(jìn)展,進(jìn)而抑制其增殖[10]。此外,Zhou等[11]研究表明,低表達(dá)的miR-595與卵巢癌患者更短的生存期相關(guān)。本研究利用TCGA數(shù)據(jù)庫下載卵巢癌的轉(zhuǎn)錄組數(shù)據(jù)及對(duì)應(yīng)的臨床信息,通過差異分析篩選與卵巢癌耐藥相關(guān)的miRNAs,隨后通過單因素Cox比例風(fēng)險(xiǎn)回歸分析及Kaplan-Meier分析篩選出與預(yù)后相關(guān)的miRNAs,并通過構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型來預(yù)測(cè)卵巢癌患者的預(yù)后風(fēng)險(xiǎn),探討其與卵巢癌預(yù)后的相關(guān)性,以期為提高卵巢癌患者的個(gè)體治療提供參考。

1 材料與方法

1.1 數(shù)據(jù)的獲取

卵巢癌相關(guān)的miRNAs表達(dá)數(shù)據(jù)及相應(yīng)的臨床數(shù)據(jù)均來源于TCGA數(shù)據(jù)庫(https://portal.gdc.cancer.gov/)。下載的基因表達(dá)數(shù)據(jù)包含499例卵巢癌患者,臨床數(shù)據(jù)包含年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲、耐藥與否及腫瘤殘留大小等信息[12]。

1.2 差異表達(dá)miRNAs篩選

利用R軟件對(duì)下載的卵巢癌數(shù)據(jù)進(jìn)行預(yù)處理,將counts數(shù)轉(zhuǎn)化為表達(dá)值,并將樣本ID轉(zhuǎn)換為miRNAs基因名。根據(jù)卵巢癌耐藥與否,我們分別將卵巢癌樣本分為耐藥組與非耐藥組,通過R軟件的“edgeR”包篩選出與卵巢癌耐藥相關(guān)的mi-RNAs,篩選標(biāo)準(zhǔn)為:|log2FC|≥1(FC:fold change),P≤0.05[13]。

1.3 與卵巢癌總體生存相關(guān)miRNAs基因集的篩選

為了篩選與卵巢癌總體生存相關(guān)的miRNAs基因集,利用 R 軟件的“survival”和“survminer”對(duì)卵巢癌耐藥相關(guān)的miRNAs分別進(jìn)行單因素Cox比例風(fēng)險(xiǎn)回歸分析及Kaplan-Meier分析,P≤0.05被認(rèn)為與卵巢癌的總體生存相關(guān),最終取兩種分析方法的交集進(jìn)行下一步分析。

1.4 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的構(gòu)建

利用R軟件中“survival”包的多因素Cox回歸分析構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型。采用前向逐步回歸法篩選多因素模型納入變量,以縮小用于構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型的miRNAs基因集數(shù)。根據(jù)風(fēng)險(xiǎn)評(píng)分的中位數(shù)將卵巢癌患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。采用Kaplan-Meier分析比較高風(fēng)險(xiǎn)組及低風(fēng)險(xiǎn)組患者的總體生存期,并采用log-rank檢驗(yàn)進(jìn)行評(píng)價(jià)。采用時(shí)間依賴性ROC曲線分析復(fù)發(fā)預(yù)測(cè)模型的敏感性和特異性。多因素分析用于評(píng)估風(fēng)險(xiǎn)評(píng)分是否獨(dú)立于年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲及腫瘤殘留大小等臨床病理參數(shù),成為卵巢癌患者預(yù)后的獨(dú)立風(fēng)險(xiǎn)因子[14]。

1.5 列線圖模型的構(gòu)建

利用R軟件的“rms”軟件包繪制列線圖模型。用于構(gòu)建列線圖的變量包括年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲、腫瘤殘留大小和風(fēng)險(xiǎn)評(píng)分。校準(zhǔn)曲線可以展示列線圖模型在全部預(yù)測(cè)概率范圍內(nèi)的校準(zhǔn)能力[15]。為評(píng)估列線圖中實(shí)際生存和預(yù)測(cè)生存的一致性,用校準(zhǔn)曲線評(píng)價(jià)列線圖模型的預(yù)測(cè)能力。

1.6 統(tǒng)計(jì)學(xué)方法

所有分析均采用SPSS 23.0和R 3.5.3進(jìn)行。所有統(tǒng)計(jì)檢驗(yàn)均為雙側(cè)檢驗(yàn),P值小于0.05被認(rèn)為具有統(tǒng)計(jì)學(xué)意義。符合正態(tài)分布的連續(xù)變量采用獨(dú)立t檢驗(yàn)進(jìn)行組間差異比較,而偏態(tài)分布的連續(xù)變量則采用Mann-Whitney U檢驗(yàn)進(jìn)行比較。

2 結(jié)果

2.1 卵巢癌耐藥相關(guān)的miRNAs

根據(jù)卵巢癌耐藥與否將患者分為耐藥組及非耐藥組,利用R軟件的“edgeR”包進(jìn)行差異分析,篩選出與卵巢癌耐藥相關(guān)的差異miRNAs共84個(gè)(P＜0.05)。

2.2 卵巢癌預(yù)后相關(guān)的miRNAs

為了探討卵巢癌耐藥相關(guān)的miRNAs與患者預(yù)后的關(guān)系,我們首先利用單變量Cox比例風(fēng)險(xiǎn)回歸分析法進(jìn)行篩選,結(jié)果顯示:16個(gè)miRNAs與卵巢癌患者的總體生存相關(guān)(表1,P＜0.05)。隨后對(duì)這16個(gè)基因進(jìn)行Kaplan-Meier分析,結(jié)果表明,共有12個(gè)基因與卵巢癌患者的預(yù)后相關(guān)(圖1,P＜0.05)。我們將這12個(gè)基因用于下一步的分析。

表1 單因素Cox比例風(fēng)險(xiǎn)回歸分析篩選出的與卵巢癌預(yù)后相關(guān)的miRNAsTable 1 miRNAs related to the prognosis of ovarian cancer selected by univariate Cox proportional risk regression analysis

2.3 卵巢癌預(yù)后風(fēng)險(xiǎn)評(píng)分模型

采用前向逐步回歸法在單變量Cox比例風(fēng)險(xiǎn)回歸分析和Kaplan-Meier分析的基礎(chǔ)上繼續(xù)縮小用于構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型的miRNAs基因集數(shù)。結(jié)果表明,經(jīng)過前向逐步回歸法,我們最終篩選出8個(gè)miRNAs(表2)。將篩選的基因置入多變量Cox比例風(fēng)險(xiǎn)回歸分析,根據(jù)這些miRNAs的回歸系數(shù)構(gòu)建以下公式:Risk10=-653.47×hsa-miR-6786-5p+282.93×hsa-miR-7155-3p-0.06×hsamiR-30e-3p+4 972.34×hsa-miR-4502-8 158.13×hsa-miR-3927-3p+427.09×hsa-miR-3662-0.09×hsa-miR-30d-5p+3 121.64×hsa-miR-548ai,構(gòu)建預(yù)后風(fēng)險(xiǎn)回歸模型。根據(jù)模型每個(gè)患者將得到一個(gè)復(fù)發(fā)風(fēng)險(xiǎn)評(píng)分。根據(jù)風(fēng)險(xiǎn)評(píng)分的中位值將患者分為高風(fēng)險(xiǎn)組及低風(fēng)險(xiǎn)組,Kaplan-Meier生存曲線提示高風(fēng)險(xiǎn)組患者具有較差的預(yù)后結(jié)果(P=0.029,圖2A),ROC曲線表明該風(fēng)險(xiǎn)模型在預(yù)測(cè)復(fù)發(fā)風(fēng)險(xiǎn)上具有較好的敏感性及特異性(圖2B)。卵巢癌患者的風(fēng)險(xiǎn)評(píng)分、生存狀態(tài)及基因表達(dá)的聚類熱圖分別見圖2C、圖2D和圖2E。

表2 用于構(gòu)建風(fēng)險(xiǎn)評(píng)分模型的miRNAs及其回歸系數(shù)Table 2 miRNAs used to construct prognostic risk score model and their regression coefficients

圖2 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的構(gòu)建(A)Kaplan-Meier曲線;(B)ROC曲線;(C)卵巢癌患者風(fēng)險(xiǎn)評(píng)分;(D)卵巢癌患者生存狀態(tài)分析;(E)卵巢癌患者基因表達(dá)聚類熱圖。Fig.2 Construction of prognostic risk score models(A)Kaplan-Meier survival curve;(B)ROC curve;(C)The risk scores of ovarian cancer patients;(D)Analysis of survival status of ovarian cancer patients;(E)The gene expression heat map of ovarian cancer patients.

2.4 預(yù)后風(fēng)險(xiǎn)評(píng)分是預(yù)測(cè)卵巢癌患者預(yù)后的獨(dú)立風(fēng)險(xiǎn)因子

為了探索構(gòu)建的預(yù)后風(fēng)險(xiǎn)評(píng)分模型對(duì)卵巢癌患者預(yù)后風(fēng)險(xiǎn)的預(yù)測(cè)是否獨(dú)立于患者的年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲及腫瘤殘留大小等臨床病理參數(shù),我們針對(duì)風(fēng)險(xiǎn)評(píng)分及臨床病理參數(shù)進(jìn)行了單因素及多因素分析。結(jié)果表明,風(fēng)險(xiǎn)評(píng)分(HR=1.41,95%CI=1.30～1.53;P=5.57E-17)及FIGO分期(HR=0.55,95%CI=0.39～0.78;P=0.000 7)均與卵巢癌患者的總體預(yù)后相關(guān)。其中,風(fēng)險(xiǎn)評(píng)分(HR=1.37,95%CI=1.09～1.73;P=0.008)為預(yù)測(cè)卵巢癌患者預(yù)后的獨(dú)立風(fēng)險(xiǎn)因子(表3)。

表3 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的單因素和多因素Cox回歸分析Table 3 Univariate and multivariate Cox regression analysis of prognostic risk score model

2.5 預(yù)后風(fēng)險(xiǎn)評(píng)分與卵巢癌患者臨床病理參數(shù)的相關(guān)性

進(jìn)一步探索預(yù)后風(fēng)險(xiǎn)評(píng)分與卵巢癌患者年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲及腫瘤殘留大小等臨床病理參數(shù)之間的關(guān)系。結(jié)果表明,預(yù)后風(fēng)險(xiǎn)評(píng)分與卵巢癌患者的FIGO分期密切相關(guān)(P＜0.000 1,表4),而與年齡、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲及腫瘤殘留大小等病理參數(shù)無關(guān)(P＞0.05,表4)。

表4 預(yù)后風(fēng)險(xiǎn)評(píng)分與卵巢癌患者臨床病理參數(shù)之間的相關(guān)性Table 4 The correlation between the prognostic risk score and clinicopathologic parameters of ovarian cancer patients

2.6 列線圖模型

最后,為了聯(lián)合風(fēng)險(xiǎn)評(píng)分及年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲、腫瘤殘留大小等臨床病理參數(shù)對(duì)患者的預(yù)后進(jìn)行有效預(yù)測(cè),我們構(gòu)建了列線圖模型。結(jié)果顯示,該模型可以充分利用各個(gè)預(yù)后變量對(duì)卵巢癌患者的死亡風(fēng)險(xiǎn)進(jìn)行精準(zhǔn)預(yù)測(cè)(圖3)。校準(zhǔn)曲線表示,列線圖模型在預(yù)測(cè)復(fù)發(fā)風(fēng)險(xiǎn)上具有較好的校準(zhǔn)能力,即模型預(yù)測(cè)風(fēng)險(xiǎn)與患者實(shí)際風(fēng)險(xiǎn)具有較好的一致性(圖4)。

圖3 列線圖模型Fig.3 Nomogram model

圖4 列線圖的校準(zhǔn)圖(A)卵巢癌患者1年生存率預(yù)測(cè);(B)卵巢癌患者2年生存率預(yù)測(cè);(C)卵巢癌患者3年生存率預(yù)測(cè)。Fig.4 Calibration diagram of nomogram(A)Prediction of 1-year survival rate for ovarian cancer patients;(B)Prediction of 2-year survival rate for ovarian cancer patients;(C)Prediction of 3-year survival rate for ovarian cancer patients.

3 總結(jié)與討論

卵巢癌的主要化療方法是以鉑類為基礎(chǔ)的聯(lián)合化療,其已被證明是最有效的化療方案。然而,在鉑類為基礎(chǔ)的聯(lián)合化療中,疾病常因耐藥而復(fù)發(fā)[16]。因此,尋找耐藥相關(guān)的腫瘤標(biāo)記物,并探索其對(duì)卵巢癌預(yù)后的影響,對(duì)卵巢癌患者的精準(zhǔn)治療至關(guān)重要。miRNA是一類由內(nèi)源性基因編碼的長度約為22個(gè)核苷酸的非編碼RNA分子,參與調(diào)控動(dòng)植物轉(zhuǎn)錄后的基因表達(dá)[17]。Shin等[18]研究表明,miRNA可以參與腫瘤的發(fā)生發(fā)展,調(diào)節(jié)腫瘤細(xì)胞對(duì)放療和化療的敏感性。因此,本研究篩選出與卵巢癌耐藥相關(guān)的miRNAs,并從中篩選出與預(yù)后相關(guān)的miRNAs構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型,以預(yù)測(cè)卵巢癌患者的復(fù)發(fā)風(fēng)險(xiǎn)。

首先,我們從TCGA數(shù)據(jù)庫中下載卵巢癌轉(zhuǎn)錄組數(shù)據(jù)及臨床數(shù)據(jù),根據(jù)卵巢癌耐藥與否將患者分為耐藥組及非耐藥組,并進(jìn)行差異分析,篩選出與卵巢癌耐藥相關(guān)的基因。隨后,對(duì)篩選出來的基因進(jìn)行單因素Cox比例風(fēng)險(xiǎn)回歸分析及Kaplan-Meier分析,進(jìn)一步篩選出與預(yù)后相關(guān)的基因。最后,利用這些基因構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型。Kaplan-Meier生存曲線分析表示,與低風(fēng)險(xiǎn)組患者相比,高風(fēng)險(xiǎn)組患者具有較差的預(yù)后。ROC曲線表示,該風(fēng)險(xiǎn)模型在預(yù)測(cè)復(fù)發(fā)風(fēng)險(xiǎn)上具有較好的敏感性及特異性。進(jìn)一步的Cox回歸分析提示,該預(yù)后風(fēng)險(xiǎn)評(píng)分是預(yù)測(cè)卵巢癌患者復(fù)發(fā)風(fēng)險(xiǎn)的獨(dú)立影響因子。此外,我們還聯(lián)合風(fēng)險(xiǎn)評(píng)分及年齡、FIGO分期、病理分級(jí)、淋巴結(jié)轉(zhuǎn)移、脈管侵襲、腫瘤殘留大小等臨床病理參數(shù)構(gòu)建了列線圖模型,結(jié)果顯示,該模型可以對(duì)卵巢癌患者的死亡風(fēng)險(xiǎn)進(jìn)行精準(zhǔn)預(yù)測(cè)。

在用于構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型的8個(gè)miRNAs中,目前僅3個(gè)是有相關(guān)文獻(xiàn)報(bào)道過的,即hsa-miR-30e-3p、hsa-miR-3662 和 hsa-miR-30d-5p。Koutsaki等[19]研究表明 hsa-miR-30e-3p的表達(dá)與卵巢癌的上皮-間質(zhì)轉(zhuǎn)化相關(guān)。hsa-miR-3662在實(shí)體瘤肺腺癌中過表達(dá),且在可手術(shù)腺癌的診斷中發(fā)揮著重要的作用[20]。然而,hsa-miR-3662被報(bào)道在肝癌組織和細(xì)胞系低表達(dá),且低表達(dá)的hsa-miR-3662與腫瘤大小、腫瘤多樣性、Edmondson分級(jí)和腫瘤-淋巴結(jié)轉(zhuǎn)移分期相關(guān)[21]。Shi等[22]研究表明hsa-miR-30d-5p的高表達(dá)與卵巢癌良好的預(yù)后水平相關(guān)。

綜上所述,我們通過一系列生物信息學(xué)方法篩選出8個(gè)與卵巢癌預(yù)后相關(guān)的miRNA標(biāo)志物。多因素Cox比例風(fēng)險(xiǎn)回歸分析表明,8個(gè)miRNAs是卵巢癌預(yù)后的獨(dú)立風(fēng)險(xiǎn)因素。結(jié)合預(yù)后風(fēng)險(xiǎn)評(píng)分及卵巢癌的臨床病理參數(shù)構(gòu)建的列線圖模型可有效預(yù)測(cè)卵巢癌的生存率。該預(yù)測(cè)方法在一定程度上可為未來卵巢癌預(yù)后相關(guān)研究提供方向。