基于GEO數(shù)據(jù)庫建立上皮性卵巢癌分子預(yù)后風(fēng)險評分模型

王興國，馬 剛，董 健，徐智陽，劉淑娟(空軍軍醫(yī)大學(xué)西京醫(yī)院婦產(chǎn)科，陜西 西安 700；解放軍第90部隊醫(yī)院，吉林 延邊 000；空軍軍醫(yī)大學(xué)西京醫(yī)院消化內(nèi)科，陜西 西安 700)

卵巢癌是婦科常見的惡性腫瘤之一，在全球范圍內(nèi)發(fā)病率居婦科惡性腫瘤第3位，死亡率居第2位[1]。上皮性卵巢癌是卵巢癌最常見的組織學(xué)類型，手術(shù)和化療是上皮性卵巢癌的主要治療方法[2]。由于卵巢癌早期缺乏典型癥狀，且無有效的早期篩查手段，以致75%的患者在確診時已處于晚期，5年生存率僅為39%[3]。隨著手術(shù)、化療、靶向及免疫治療等手段的聯(lián)合應(yīng)用，近年來卵巢癌患者的生存獲益改善明顯，但其預(yù)后仍不盡人意[4]。通常，臨床分期、腫瘤細(xì)胞減滅術(shù)的實施情況和對化療藥物的敏感性是預(yù)測卵巢癌患者預(yù)后的基本因素[5]。然而，上皮性卵巢癌是一種具有復(fù)雜分子特征和遺傳物質(zhì)改變的高度異質(zhì)性疾病，以致臨床特征相似病例的生存率和治療反應(yīng)差異很大，故探索基于分子水平的預(yù)后評估顯得尤為重要[6]。

隨著生物信息學(xué)的發(fā)展，挖掘在線數(shù)據(jù)庫已成為當(dāng)下熱門。然而，在基因表達(dá)譜數(shù)據(jù)中，預(yù)測變量遠(yuǎn)遠(yuǎn)大于樣本含量且各變量之間常具有強(qiáng)相關(guān)性，呈現(xiàn)高維度和共線性。最小化絕對收縮和選擇算子(least absolute shrinkage and selection operator，LASSO)回歸算法可在系數(shù)的絕對值之和上增加一個約束條件來對高維資料進(jìn)行降維，實現(xiàn)變量選擇的目的[7]。逐步回歸的基本思想是逐個引入新變量，每引入一個新變量時考慮是否剔除已選變量，直至不再引入新變量。這種方法所得的回歸方程的變量較少，并保留了影響最顯著的重要變量[8]。LASSO回歸聯(lián)合逐步回歸的優(yōu)勢在于可篩選出相關(guān)度更高的預(yù)后變量。

本研究基于高通量基因表達(dá)(Gene Expression Omnibus，GEO)數(shù)據(jù)庫中的GSE26712數(shù)據(jù)集，利用LASSO回歸聯(lián)合逐步回歸向后法對上皮性卵巢癌預(yù)后相關(guān)基因進(jìn)行篩選與模型構(gòu)建，為評估上皮性卵巢癌患者的預(yù)后及個體化治療方案的制定提供了一個較為準(zhǔn)確的工具。

1 對象與方法

1.1 對象

在GEO公共數(shù)據(jù)庫(https：//www.ncbi.nlm.nih.gov/gds)檢索并獲取GSE26712的基因表達(dá)譜數(shù)據(jù)。GSE26712數(shù)據(jù)集基于GPL96平臺(HG-U133A)Affymetrix Human Genome U133A Array，其中包含185例原發(fā)性上皮性卵巢癌樣本和10例正常卵巢表面上皮樣本的基因表達(dá)信息。U133A Array是一種覆蓋33 000個充分表征人類基因的單陣列，可用于探索疾病的發(fā)生發(fā)展。185例未經(jīng)治療的晚期(Ⅲ/Ⅳ期)高級別(G2/G3級)原發(fā)性上皮性卵巢癌患者均接受腫瘤細(xì)胞減滅術(shù)聯(lián)合術(shù)后含鉑類化療，其中滿意減瘤術(shù)患者92例，不滿意減瘤術(shù)患者93例。

1.2 方法

1.2.1 差異表達(dá)基因(differentially expressed genes，DEGs)的篩選 通過利用R軟件中“l(fā)imma”包分析芯片表達(dá)矩陣數(shù)據(jù)從而篩選出DEGs。篩選標(biāo)準(zhǔn)：P<0.05，|log2(FC)|>0.585。將log2(FC)>0.585的基因作為上調(diào)差異基因(up-regulated DEGs，UDEGs)，log2(FC)<-0.585的基因作為下調(diào)差異基因(down-regulated DEGs，DDEGs)。

1.2.2 預(yù)后特征基因的篩選及評估 基于R軟件，按照7∶3將患者隨機(jī)分為訓(xùn)練集(129例)和驗證集(56例)。利用“rms”“survival”包對訓(xùn)練集進(jìn)行單因素Cox回歸分析，篩選出卵巢癌患者預(yù)后相關(guān)的基因(P<0.05)。然后將預(yù)后相關(guān)基因與DEGs取交集進(jìn)而獲得重疊候選基因。通過“glmnet”包對重疊候選基因進(jìn)行LASSO回歸分析，來縮減預(yù)后基因的數(shù)量。采用逐步回歸向后法進(jìn)一步篩選與預(yù)后顯著相關(guān)的基因，根據(jù)赤池信息準(zhǔn)則(Akaike information criterion，AIC)獲得最優(yōu)特征基因集?；跀?shù)據(jù)集，分析預(yù)后特征基因在卵巢癌組織與正常組織中的表達(dá)差異。最后，利用cBioPortal(http：//cbioportal.org)數(shù)據(jù)庫在線分析預(yù)后特征基因改變的情況。

1.2.3 預(yù)后模型的構(gòu)建與評估 通過多因素Cox比例風(fēng)險回歸確定預(yù)后特征基因的回歸系數(shù)，并構(gòu)建風(fēng)險評分模型，最后，利用“nomogram”函數(shù)實現(xiàn)風(fēng)險模型的可視化。受試者工作特征(receiver operating characteristic，ROC)曲線和Kaplan-Meier(KM)生存分析用于評估模型預(yù)測的準(zhǔn)確性。公式為：

1.2.4 預(yù)后模型的驗證 為了更好地驗證模型的效能，使用“caret”包的十折交叉驗證[9]，即將全部樣本隨機(jī)劃分為10個樣本量相同的子集，依次將其中1個作為驗證集，其余9個為訓(xùn)練集。使用“riskRegression”包繪制ROC曲線。

1.2.5 統(tǒng)計學(xué)分析 利用R(4.1.1)完成數(shù)據(jù)統(tǒng)計分析。采用KM法進(jìn)行生存分析，兩組生存曲線比較采用Log-rank檢驗。Cox比例風(fēng)險回歸模型用于單變量和多變量分析。P<0.05表示差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 篩選DEGs

經(jīng)篩選后得到3 397個DEGs，其中包括1 490個UDEGs，1 907個DDEGs。

2.2 卵巢癌預(yù)后特征基因的篩選與評估

對訓(xùn)練集的13 237個基因進(jìn)行單因素Cox回歸分析，共計155個基因與預(yù)后顯著相關(guān)(P<0.05)，然后將預(yù)后相關(guān)基因與DEGs取交集，最后獲得16個重疊候選基因，其中風(fēng)險比(hazard ratio，HR)>1的基因12個，HR<1的基因4個(圖1)。

UDEGs：上調(diào)差異基因；DDEGs：下調(diào)差異基因。圖1 16個重疊候選基因韋恩圖

通過LASSO回歸分析進(jìn)一步篩選預(yù)后相關(guān)基因，選擇α=1，運(yùn)用十折交叉驗證法實現(xiàn)模型內(nèi)部驗證。在λmin為0.051 148 63時，進(jìn)而篩選出10個預(yù)后基因。最后，采用逐步回歸法，依據(jù)AIC最小值(718.25)判斷確定了5個與卵巢癌預(yù)后相關(guān)的特征基因，即IGFBP4、IGF2、TLR2、DIAPH2、AADAC。

然后，對比腫瘤組織與正常組織之間的5個基因表達(dá)譜發(fā)現(xiàn)，IGF2和IGFBP4在卵巢癌中顯著上調(diào)，而TLR2、DIAPH2和AADAC顯著下調(diào)(P<0.005，圖2)。

圖2 5個預(yù)后特征基因在卵巢癌組織與正常卵巢組織中的表達(dá)

最后，在cBioPortal數(shù)據(jù)庫中搜索了5個基因的遺傳改變類型，擴(kuò)增(Amplification)在AADAC和DIAPH2中常見，深度缺失(Depth Deletion)在IGFBP4和TLR2中常見(圖3)。

圖3 5個預(yù)后特征基因改變視覺圖

2.3 預(yù)后模型的構(gòu)建與評估

對5個預(yù)后基因進(jìn)行多因素Cox回歸分析(圖4)，并根據(jù)5個基因的Cox系數(shù)構(gòu)建風(fēng)險評分模型：risk score=(0.205 44)×IGFBP4+(0.096 32)×IGF2+(-0.796 39)×TLR2+(-0.562 80)×DIAPH2+(-0.218 78)×AADAC。

根據(jù)風(fēng)險評分公式，以風(fēng)險評分的中位數(shù)(risk score=-4.78)為截斷值將訓(xùn)練集劃分為高風(fēng)險組與低風(fēng)險組(圖5A)，高風(fēng)險組中卵巢癌患者死亡例數(shù)較低風(fēng)險組多，低風(fēng)險組中卵巢癌患者生存時間超過5年的患者例數(shù)較高風(fēng)險組多(圖5B)。對訓(xùn)練集中的高、低風(fēng)險組患者進(jìn)行KM生存分析，高風(fēng)險組患者預(yù)后顯著低于低風(fēng)險組患者，差異有統(tǒng)計學(xué)意義(圖5C)。進(jìn)一步通過1、3、5年ROC曲線評估預(yù)后模型的靈敏度和特異度，在訓(xùn)練集中1、3和5年總生存ROC曲線下面積(area under the curve，AUC)分別為0.813、0.876、0.895(圖5D)。

A：風(fēng)險評分的分布和中值；B：高低風(fēng)險評分組患者的生存點圖；C：高低風(fēng)險評分組患者的生存曲線；D：預(yù)后模型的ROC曲線。圖5 訓(xùn)練集中5個預(yù)后特征基因模型的預(yù)后分析

為了達(dá)到預(yù)后模型可視化的效果，基于訓(xùn)練集數(shù)據(jù)構(gòu)建列線圖預(yù)后模型來預(yù)測患者1、3和5年預(yù)后(圖6)。

圖6 預(yù)測卵巢癌患者總生存率的列線圖

2.4 預(yù)后模型的交叉驗證

通過十折交叉驗證分析，在驗證集中，預(yù)測1、3、5年總生存率的AUC分別為0.906、0.920、0.739，提示該預(yù)測模型對評估上皮性卵巢癌患者的預(yù)后具有較好的預(yù)測能力(圖7)。

圖7 十折交叉驗證下的ROC曲線

3 討論

上皮性卵巢癌復(fù)發(fā)率高、預(yù)后差的現(xiàn)狀是目前臨床的難點，因此，構(gòu)建準(zhǔn)確簡單的預(yù)后評估方法十分重要。本研究利用GSE26712篩選出16個與預(yù)后相關(guān)的基因，然后通過LASSO回歸與逐步回歸共確定了5個預(yù)后特征基因，即IGFBP4、IGF2、TLR2、DIAPH2和AADAC。

胰島素樣生長因子(insulin-like growth factor，IGF)家族包括兩種配體(IGF1、IGF2)、兩類受體(IGF1R、IGF2R)和6種胰島素樣生長因子結(jié)合蛋白(insulin-like growth factor binding proteins，IGFBPs)。IGF2通過與其受體和IGFBPs結(jié)合，在胚胎發(fā)育、代謝紊亂和腫瘤發(fā)生發(fā)展中發(fā)揮重要作用[10]。IGFBP4作為IGFBPs中的一員，其功能是抑制IGF1、IGF2與其相應(yīng)的受體結(jié)合，進(jìn)而促進(jìn)或抑制IGF[11]。研究表明，IGF2與卵巢癌的進(jìn)展、預(yù)后和化療耐藥反應(yīng)有關(guān)。與正常組織相比，卵巢癌組織中IGF2表達(dá)顯著增加，且高表達(dá)的IGF2是上皮性卵巢癌患者的獨立預(yù)后因素[12]。IGF2在紫杉醇耐藥卵巢癌細(xì)胞系中的表達(dá)顯著高于敏感細(xì)胞系，進(jìn)一步研究發(fā)現(xiàn)，抑制IGF通路可恢復(fù)對紫杉醇的敏感性[13]。研究發(fā)現(xiàn)，IGFBP4在多種腫瘤中發(fā)揮作用。肺癌血清IGFBP4的表達(dá)水平顯著高于健康對照組，且可作為潛在的腫瘤標(biāo)志物[14]；IGFBP4在卵巢癌組織中表達(dá)上調(diào)，且與總生存期、無進(jìn)展生存期相關(guān)[15]。TLR2作為TLR家族的成員，通過識別病原相關(guān)分子模式，參與天然免疫防御和腫瘤發(fā)生[16]。在慢性淋巴細(xì)胞白血病中，TLR2低表達(dá)與預(yù)后顯著相關(guān)，提示TLR2可能成為慢性淋巴細(xì)胞白血病患者臨床預(yù)后的潛在生物標(biāo)志物[17]。DIAPH2基因編碼的蛋白在生理環(huán)境中參與細(xì)胞運(yùn)動與黏附。研究表明，失調(diào)的DIAPH2可增強(qiáng)腫瘤細(xì)胞的遷移能力并促進(jìn)其轉(zhuǎn)移[18]。DIAPH2在喉鱗狀細(xì)胞癌患者中低表達(dá)，其低表達(dá)預(yù)示著患者預(yù)后不良，其可能參與喉鱗狀細(xì)胞癌的轉(zhuǎn)移，可作為喉鱗狀細(xì)胞癌轉(zhuǎn)移的預(yù)測指標(biāo)[19]。AADAC是一種絲氨酸酯酶，主要存在于肝微粒體，在水解各種藥物和異源物質(zhì)中起著至關(guān)重要的作用[20]。AADAC的過表達(dá)抑制卵巢癌細(xì)胞的惡性進(jìn)展[21]。

隨后，基于5個預(yù)后特征基因的回歸系數(shù)，構(gòu)建風(fēng)險評分模型并研究其在上皮性卵巢癌患者中的預(yù)后價值，發(fā)現(xiàn)高風(fēng)險組患者的預(yù)后顯著差于低風(fēng)險組患者。最后，采用“riskRegression”函數(shù)對預(yù)測模型進(jìn)行診斷性分析，結(jié)果表明，基于IGFBP4、IGF2、TLR2、DIAPH2和AADAC構(gòu)建的預(yù)后模型具有較高的診斷價值。在訓(xùn)練集中，預(yù)測1、3、5年總生存率的AUC分別為0.813、0.876、0.895；在交叉驗證集中，預(yù)測1、3、5年總生存率的AUC分別為0.906、0.920、0.739。因此，本研究構(gòu)建的風(fēng)險評分模型可幫助臨床醫(yī)生預(yù)測上皮性卵巢癌患者的預(yù)后，并用于輔助臨床決策。

本研究還存在著一些局限性。首先，本文采用的GEO數(shù)據(jù)庫缺少每例患者具體臨床預(yù)后指標(biāo)，如年齡、組織學(xué)類型、術(shù)后殘余病灶大小、重要器官轉(zhuǎn)移、腫瘤標(biāo)志物(CA125、HE4)檢測結(jié)果。因此，模型預(yù)測的準(zhǔn)確性還需進(jìn)一步提高。其次，本研究以在線數(shù)據(jù)為基礎(chǔ)進(jìn)行分析，但結(jié)論尚需實際臨床數(shù)據(jù)進(jìn)一步驗證。

綜上所述，本研究通過GEO數(shù)據(jù)庫篩選出了5個上皮性卵巢癌患者預(yù)后特征基因，基于這5個基因構(gòu)成的風(fēng)險模型能夠很好地對上皮性卵巢癌患者預(yù)后進(jìn)行高、低分險分組，為卵巢癌患者的預(yù)后預(yù)測及個性化治療方案制定提供了一定的依據(jù)。IGFBP4、IGF2、TLR2、DIAPH2和AADAC有可能成為上皮性卵巢癌的預(yù)后標(biāo)志物，有助于上皮性卵巢癌患者個性化治療及臨床預(yù)后判斷。