整合生物信息學(xué)分析篩選宮頸癌的關(guān)鍵預(yù)后生物標(biāo)志物

王小燕，李虎玲，林丹丹，張 晶, 王 凱

（新疆醫(yī)科大學(xué)1公共衛(wèi)生學(xué)院，2醫(yī)學(xué)工程技術(shù)學(xué)院，烏魯木齊 830017）

根據(jù)全球癌癥監(jiān)測機構(gòu)（globalcan）2020 年的數(shù)據(jù)，宮頸癌新發(fā)病例和死亡人數(shù)分別為60.4萬和34.2萬，發(fā)病率和死亡率排在第四位[1]。宮頸癌最常見的組織學(xué)亞型為鱗狀細(xì)胞癌和腺癌，分別約占所有宮頸癌的70%和25%[2]。一旦發(fā)展到轉(zhuǎn)移或復(fù)發(fā)階段，無法治愈，總生存時間（overall survival，OS）約為12個月[3]。因此，尋找新的預(yù)后生物標(biāo)志物和治療靶點可能有助于提高宮頸癌患者的生存率。

在已發(fā)表的報道中，關(guān)于宮頸癌發(fā)生發(fā)展中的分子標(biāo)志物的研究也取得了顯著性進展。例如，Wang 等[4]研究發(fā)現(xiàn)CDC7 基因的表達上調(diào)與宮頸癌的發(fā)生和發(fā)展相關(guān)，靶向這種生物標(biāo)志物可能會改善宮頸癌的早期診斷和治療；DudeaSimon 等[5]認(rèn)為血管內(nèi)皮生長因子（VEGF）參與ILK 信號傳導(dǎo)，能夠預(yù)測總體存活率，可能是影響宮頸癌預(yù)后的重要基因；此外，Zhang 等[6]研究表明MiR-378a-3p 下調(diào)與預(yù)后相關(guān)，可能是宮頸癌的潛在生物標(biāo)志物。然而，單一的標(biāo)志物在預(yù)測宮頸癌患者的預(yù)后時可能會存在一定的局限性。因此，整合預(yù)后生物標(biāo)志物在預(yù)測宮頸癌的不良預(yù)后方面具有重要意義。

隨著微陣列芯片技術(shù)和高通量測序技術(shù)的發(fā)展，新的預(yù)后生物標(biāo)志物和治療靶點不斷被研究發(fā)現(xiàn)。整合生物信息學(xué)方法[7]能夠高效的利用多方數(shù)據(jù)庫，因而在癌癥組學(xué)領(lǐng)域被廣泛應(yīng)用。在本研究中，利用生物信息技術(shù)和方法整合分析基因表達數(shù)據(jù)庫（Gene Expression Omnibus，GEO）和癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫中宮頸癌的轉(zhuǎn)錄組數(shù)據(jù)，以期識別出參與宮頸癌發(fā)生發(fā)展中的Hub 基因、生物學(xué)功能及信號通路；同時將多個Hub基因構(gòu)建一個標(biāo)識，利用標(biāo)識的風(fēng)險評分對患者的預(yù)后進行預(yù)測及評價，證實其可作為關(guān)鍵生物標(biāo)記物能夠更好地判斷宮頸癌的預(yù)后。

1 資料與方法

1.1 數(shù)據(jù)資料宮頸癌的mRNA 表達數(shù)據(jù)來源于GEO 和TCGA 數(shù)據(jù)庫。首先，從GEO 中選擇微陣列數(shù)據(jù)集GSE90738，包括10個宮頸腫瘤組織和10個匹配的宮頸癌患者的癌旁組織。根據(jù)補充文件GSE90738_ cervical_cancer_ mRNA_ processed.xlsx，提取20 個樣本對應(yīng)的mRNA 基因符號數(shù)據(jù)，根據(jù)基因平均表達量值最大去除重復(fù)基因，平均表達量小于1過濾基因的原則，共獲得16 367個基因用于下一步分析。其次，通過TCGA 獲取宮頸癌的mRNA 表達數(shù)據(jù)，共納入307 個宮頸癌樣本，包括304 個宮頸腫瘤組織和3 個正常宮頸組織，根據(jù)篩選原則，共獲得12 571 個基因的log2（FPKM+1）的基因表達數(shù)據(jù)用于后續(xù)分析。cBioPortal（http://www.cbioportal.org/）是一個開源資源平臺，可下載多種癌癥基因組數(shù)據(jù)集及臨床數(shù)據(jù)。從cBioPortal 下載了TCGA 數(shù)據(jù)庫中對應(yīng)的宮頸癌患者相應(yīng)的臨床信息。

1.2 篩選差異基因采用R軟件“l(fā)imma”包標(biāo)準(zhǔn)化矩陣數(shù)據(jù)并分別鑒定GSE90738 和TCGA 數(shù)據(jù)中宮頸腫瘤組織與宮頸正常組織間的差異基因。以校正后的P<0.05和|logFC|>1為差異基因篩選條件，盡可能消除假陽性結(jié)果。使用“ggplot2”和“ pheatmap”包分別繪制差異基因的火山圖和熱圖。為了消除不同數(shù)據(jù)平臺上不同表達量類型造成的背景誤差，本研究使用Venny 2.1.0（https://bioinfogp.cnb.csic.es/tools/ven?ny）映射篩選出共享的差異基因，并以GSE90738數(shù)據(jù)為參考繪制共享差異基因熱圖。

1.3 功能和通路富集分析及蛋白質(zhì)互作網(wǎng)絡(luò)構(gòu)建

利用R 軟件clusterProfiler 包中的enrichGO 和en?richKEGG 函數(shù)進行共享差異基因的GO 富集分析和KEGG通路分析。以P<0.05為閾值鑒定顯著GO的生物學(xué)過程和KEGG 通路。STRING（search tool for theretrieval of interacting genes）是一個免費在線生物分析數(shù)據(jù)庫，提供已知和預(yù)測的蛋白質(zhì)-蛋白質(zhì)相互作 用 網(wǎng) 絡(luò)（protein-protein interaction,PPI）。應(yīng) 用STRING 數(shù)據(jù)庫，構(gòu)建共享差異基因的PPI 網(wǎng)絡(luò)。然后，通過Cytoscape 軟件對PPI 網(wǎng)絡(luò)進行可視化，并應(yīng)用其MCODE 插件對網(wǎng)絡(luò)進行聚類，從網(wǎng)絡(luò)圖中找到Hub基因組塊和基因進入下一步分析。

1.4 Hub 基因的篩選及驗證使用R 軟件“survival”包進行單因素Cox 比例風(fēng)險回歸分析選擇與宮頸癌患者總體生存期相關(guān)的預(yù)后基因。使用R 語言“glm?net”包對單因素Cox 顯著性分析結(jié)果P<0.05 的預(yù)后基因進行多元逐步Cox 回歸分析，將篩選出的預(yù)后基因作為Hub 基因。GEPIA2（http://gepia2.cancer-pku.cn/# analysis）是一個在線分析數(shù)據(jù)庫，能夠?qū)CGA和GTEx 項目共9 736 個腫瘤樣本、8 587 個正常樣本的RNA-seq 表達數(shù)據(jù)進行分析。利用GEPIA2 數(shù)據(jù)庫中的Boxplots 工具進一步驗證Hub 基因在宮頸腫瘤組織與正常組織之間的表達水平。

1.5 構(gòu)建預(yù)后Hub 基因風(fēng)險標(biāo)識為進一步評估Hub 基因的預(yù)后價值，對Hub 基因進行多因素Cox 分析，并構(gòu)建Hub基因的風(fēng)險評分模型。風(fēng)險評分模型中的每個Hub 基因系數(shù)來自于多因素Cox 分析對應(yīng)的變量系數(shù)。Hub 基因風(fēng)險標(biāo)識（Hub Genes Risk Signature，HGRS）計算如下：HGRS = （βHubgenei* EX?PHubgenei）。以HGRS 的中位數(shù)為cutoff 值，將患者分為高、低風(fēng)險組，采用Kaplan-Meier曲線進行生存分析，log-rank 檢驗P<0.05 認(rèn)為差異有統(tǒng)計學(xué)意義。SangerBox 是一個免費的在線數(shù)據(jù)分析平臺（http://www.sangerbox.com/tool），使 用SangerBox 工 具 對HGRS 進行可視化分析，并繪制ROC 曲線評價HGRS的預(yù)測能力。

2 結(jié)果

2.1 差異表達分析在GSE90738 數(shù)據(jù)集中，共鑒定出差異基因1 282個，包括上調(diào)基因777個，下調(diào)基因505 個（圖1A）。在TCGA 數(shù)據(jù)庫中下載的宮頸癌數(shù)據(jù)集中，共鑒定出差異基因2 203 個，包括上調(diào)基因1 003 個，下調(diào)基因1 200 個（圖1B）。GSE90738 數(shù)據(jù)集和TCGA 數(shù)據(jù)集中的差異基因取交集后，獲得486個共享差異基因，包括上調(diào)基因319個，下調(diào)基因167個（圖2）。以GSE90738 數(shù)據(jù)為參考繪制486 個共享差異基因熱圖（圖1C）。

圖1 差異基因的鑒定

圖2 共享差異基因的韋恩圖

2.2 GO 富集分析和KEGG 信號通路分析通過“clusterProfiler”包對486 個共享差異基因進行GO 富集分析和KEGG 信號通路分析。在生物過程（BP）、細(xì)胞組分（CC）、分子功能（MF）這三個生物學(xué)方面的GO 富集分析中，最顯著的前10 項進行分析展示。結(jié)果顯示，在BP 中，共享差異基因主要富集在“染色體分離”，“核分裂”，“有絲分裂核分裂”，“細(xì)胞器裂變”，“核染色體分離”，“DNA 復(fù)制”，“姊妹染色單體分裂”，“細(xì)胞周期G1/S 期轉(zhuǎn)變”，“有絲分裂姐妹染色單體分離”，“有絲分裂細(xì)胞周期的G1/S 轉(zhuǎn)變”等生物過程。在宮頸癌中，共享差異基因主要富集在“染色體區(qū)域”，“染色體著絲粒區(qū)域”，“濃縮染色體”，“染色體濃縮，著絲粒區(qū)域”，“紡錘體”等細(xì)胞組分。在MF中，共享差異基因主要富集在“催化活性，作用于DNA”，“DNA 復(fù)制起始結(jié)合”，“DNA 解旋酶的活動”，“微管結(jié)合”，“微管蛋白結(jié)合”，“解旋酶的活動”等分子功能（圖3A）。根據(jù)KEGG信號通路分析，共享差異基因主要參與“細(xì)胞周期”，“癌癥中的微小RNA”，“PI3K-Akt信號通路”，“細(xì)胞衰老”，“人乳頭瘤病毒感染”，“Epstein-Barr病毒感染”等路徑過程（圖3B）。

圖3 共享差異基因的GO和KEGG富集分析

2.3 PPI 網(wǎng)絡(luò)構(gòu)建及Hub 基因模塊分析使用STRING 在線數(shù)據(jù)庫構(gòu)建由445 個節(jié)點9 605 個連接組成的共享差異基因的PPI 網(wǎng)絡(luò)圖（圖4A）。應(yīng)用MCODE 插件的默認(rèn)參數(shù)設(shè)置，共得到19 個模塊。本研究選擇MCODE1 作為Hub 基因模塊，因評分最高達到100.393 分，由118 個節(jié)點和5 873 個連接組成（圖4B），且均為上調(diào)共享差異基因。

圖4 共享差異基因的PPI網(wǎng)絡(luò)和Hub基因模塊分析

2.4 Hub 基因的篩選和驗證通過單因素Cox 比例風(fēng)險回歸模型分析Hub 基因模塊中的118 個基因，發(fā)現(xiàn)16 個基因（P<0.05）與患者總生存期顯著相關(guān)。將以上16 個基因進行多元逐步Cox 回歸分析，最終得到4 個與預(yù)后相關(guān)的Hub 基因，為CENPM、ANLN、CHAF1A 和HELLS。利用GEPIA2 工具驗證4 個Hub基因的表達水平。結(jié)果表明，4 個Hub 基因均在腫瘤組織中高表達（圖5），這與4 個Hub 基因均為上調(diào)基因的結(jié)果一致。

圖5 4個Hub基因在宮頸腫瘤組織和正常組織中的表達比較

2.5 預(yù)后Hub 基因風(fēng)險標(biāo)識的構(gòu)建多元逐步Cox回歸分析結(jié)果如表1 所示，CENPM、CHAF1A、HELLS為保護性因素，ANLN 為危險因素，其中CENPM、CHAF1A 和ANLN 為影響宮頸癌患者預(yù)后的獨立影響因素。根據(jù)4 個Hub 基因的回歸系數(shù)β 和基因表達量構(gòu)建HGRS，HGRS=（-0.458）* CENPM+（0.561）*ANLN+（-0.558）*CHAF1A +（-0.504）*HELLS。根據(jù)HGRS 值的中位數(shù)-3.083，患者被分為了136 人的高風(fēng)險組和137人的低風(fēng)險組。圖6A、6B、6C分別展示了在患者中風(fēng)險值的分布、生存時間和生存結(jié)局的分布及4 個Hub 基因Z-score 值熱圖。ROC 曲線顯示，HGRS 的1、3、5 年曲線下面積（Area under the curve，AUC）分別為0.67（95%CI：0.53～0.81）、0.72（95%CI：0.64～0.81）、0.76（95%CI：0.66～0.85）（圖6D）。圖6E 表明低風(fēng)險組的總體生存時間明顯高于高風(fēng)險組（P<0.001）。圖7 展示了HGRS 預(yù)測患者的生存狀況具有一定的穩(wěn)健性。

圖6 預(yù)后Hub基因風(fēng)險標(biāo)識HGRS的生存預(yù)測及效果評價

圖7 HGRS在5年內(nèi)各時間點上的AUC值及95%置信區(qū)間

表1 4個Hub基因的總體貢獻程度

3 討論

人乳頭瘤病毒（HPV）感染是宮頸癌的主要危險因素，其高危亞型幾乎導(dǎo)致所有宮頸癌[8]。微陣列芯片技術(shù)和高通量測序技術(shù)已廣泛用于研究癌癥的基因改變和確定疾病特異性預(yù)后生物標(biāo)志物和治療靶點。因此，本研究進行了基于微陣列和高通量測序的轉(zhuǎn)錄組分析，以確定宮頸癌的異常調(diào)節(jié)基因。

本研究中，結(jié)合GEO 數(shù)據(jù)庫中宮頸癌GSE90738數(shù)據(jù)集和TCGA 數(shù)據(jù)集的顯著差異基因，共得到486個共享差異基因，包括319 個上調(diào)基因和167 個下調(diào)基因。GO 功能富集分析表明，共享差異基因主要富集在“染色體分離”、“染色體區(qū)域”、“催化活性，作用于DNA”；KEGG 通路富集分析表明，共享差異基因主要參與“細(xì)胞周期”，“癌癥中的微小RNA”，“PI3KAkt 信號通路”，“細(xì)胞衰老”，“人乳頭瘤病毒感染”，“Epstein-Barr 病毒感染”等信號路徑過程。研究表明，細(xì)胞周期調(diào)控缺陷，是癌癥發(fā)病機制的基本特征[9]。此外，越來越多的證據(jù)表明，大量MicroRNAs在宮頸癌組織中異常表達，在腫瘤發(fā)生、進展和轉(zhuǎn)移中發(fā)揮著不可替代的作用[10]。研究表明PI3K-Akt 信號通路通過多條途徑介導(dǎo)化療耐藥過程，包括凋亡相關(guān)蛋白表達、ABC轉(zhuǎn)運、NF-κB、mTOR信號等[11]。

本研究構(gòu)建了由4 個Hub 基因構(gòu)成的風(fēng)險標(biāo)識，能夠區(qū)分高風(fēng)險組患者和低風(fēng)險組患者，且兩組患者的總生存期存在明顯差異。4 個Hub 基因均已被證實與多種癌癥的預(yù)后密切相關(guān)。著絲粒蛋白M（centromere protein M，CENPM）是近年發(fā)現(xiàn)的促癌分子，它編碼一種動力蛋白，在細(xì)胞分裂過程中與紡錘體微管結(jié)合，調(diào)節(jié)染色體的分離[12]。Xiao 等[13]的研究證實，CENPM 與肝癌進展密切相關(guān)，CENPM 的上調(diào)通過多種機制促進肝癌的發(fā)生，可作為肝癌的新的可能生物標(biāo)志物和治療靶點。染色體組裝因子1 單位A（chromatin assembly factor 1，subunit A，CHAF1A）是一種高度保守的組蛋白伴侶分子，可調(diào)控細(xì)胞生長、胚胎發(fā)育以及DNA 修復(fù)[14]。Chen 等[15]通過GEO中的胃癌數(shù)據(jù)進行外部驗證，證實了一個由CHAF1A和RMI1 構(gòu)成的預(yù)后標(biāo)識可以有效預(yù)測胃癌患者的總體生存率。Han 等[16]的研究表明CHAF1A 可作為宮頸癌的潛在診斷和預(yù)后生物標(biāo)志物之一，與本研究結(jié)果類似。HELLS 是一種染色質(zhì)重塑因子，研究報道其在肝癌、胰腺癌、肺癌中高表達，通過介導(dǎo)多個抑癌基因的沉默，增強癌細(xì)胞的增殖和遷移，從而導(dǎo)致更差的患者預(yù)后[17-19]。Liu 等[20]的研究結(jié)果表明ANLN 由SP2調(diào)控，通過PI3K/AKT 和MAPK 信號通路促進結(jié)直腸癌細(xì)胞增殖。Xia 等[21]的研究與本研究ANLN 的結(jié)果吻合，再次說明過表達ANLN 組的患者預(yù)后與低表達ANLN 組的患者預(yù)后存在顯著差異，提示ANLN 可能是一種潛在的腫瘤致癌基因，可以作為預(yù)測宮頸癌患者預(yù)后的生物標(biāo)志物。本研究中，采用Kaplan-Meier 曲線和ROC 曲線分析證明由4 個Hub 基因構(gòu)成的風(fēng)險標(biāo)識能準(zhǔn)確的區(qū)分高風(fēng)險組和低風(fēng)險組，并能較好地預(yù)測患者的生存情況。由此可見，4 個Hub 基因作為宮頸癌的生物標(biāo)志物具有較高的診斷價值。

本研究存在的局限性如下：這是一項回顧性研究，后期還需進行前瞻性實驗驗證。其次，本研究限于273 例宮頸癌患者，臨床資料不全面，可能導(dǎo)致分析結(jié)果出現(xiàn)偏差。再者，由這4 個Hub 基因組成的基因標(biāo)志物的預(yù)測能力需要更大的樣本量進一步研究和驗證標(biāo)志物的有效性。

綜上，本研究篩選出4個可能在宮頸癌發(fā)生發(fā)展中起到重要作用的Hub 基因。這些基因可能作為診斷宮頸癌的潛在分子生物標(biāo)志物。此外，由4 個Hub基因組成的新的標(biāo)志物可以進一步更好地預(yù)測患者的生存結(jié)局，為宮頸癌患者的臨床治療決策提供有效的建議。