古嘉基 連泳欣 丘福滿

廣州醫(yī)科大學公共衛(wèi)生學院(廣州 511436)

卵巢癌(ovarian cancer)是女性在全球范圍內(nèi)第七位常見的惡性腫瘤，患者的五年生存率常低于45%[1]。卵巢腫瘤根據(jù)其組織學可分為：上皮性腫瘤，性索-間質(zhì)瘤和生殖細胞腫瘤。在卵巢的惡性腫瘤中，上皮性卵巢癌占90%，其中漿液性癌是上皮性卵巢癌的主要亞型，占70%以上[2]。腫瘤的分期和分級對卵巢癌的治療至關重要，加強對卵巢癌的檢測能夠有效提高患者的生存率。但卵巢癌常在晚期才會被診斷，一方面是因為長期無癥狀和體征，另一方面患者初步治療后臨床癥狀雖完全緩解，但復發(fā)率高[3]，且卵巢癌的風險評估和傳統(tǒng)臨床指征仍存在不足。因此，需開發(fā)更可靠的方法評估卵巢癌的預后并指導臨床的個體化治療。

N6-甲基腺苷(N6-methyladenosine, m6A)是高等真核生物中最常見的編碼和非編碼RNA甲基化修飾，影響mRNA代謝的各個階段,包括折疊、剪接、翻譯和降解等[4]。m6A主要和三類蛋白發(fā)生作用，分別是m6A甲基轉移酶、m6A去甲基酶、能與甲基化結合位點結合的特殊蛋白，其編碼基因分別為寫入基因(Writers)，擦除基因(Erasers)，讀取基因 (Readers)[5]。m6A與Writers相互作用，在RNA上加入甲基，Readers對修飾的RNA進行識別并產(chǎn)生不同功能，依靠Erasers的作用使得甲基化修飾過程動態(tài)可逆，并起到調(diào)節(jié)各基因表達的作用[6]。m6A甲基化基因參與多種細胞功能并與多種疾病相關聯(lián)，與癌癥關聯(lián)尤為密切，既可以通過促進癌基因表達和抑制抑癌基因的表達來促進癌癥發(fā)展，也可以通過抑制癌基因表達和促進抑癌基因表達來抑制癌癥發(fā)展[7]，具有雙重調(diào)節(jié)作用。迄今為止，研究證實異常的m6A甲基化將會影響肝癌、大腸癌、膀胱癌、腎透明細胞癌等多個癌癥的發(fā)生發(fā)展[8]，但目前針對卵巢癌的m6A甲基化研究還相對較少。

本研究采用生物信息學的方法分析m6A甲基化基因在卵巢癌組織與正常組織中的差異表達，使用相關統(tǒng)計學方法篩選候選的m6A甲基化基因并構建風險評估模型，旨在探討其潛在的預后預測價值，并分析相關生物學功能，為卵巢癌的靶向治療及預后評估提供科學依據(jù)。

1 資料與方法

1.1 數(shù)據(jù)集的獲取與處理

從TCGA(https://portal.gdc.cancer.gov/)數(shù)據(jù)庫中下載卵巢癌的基因轉錄組測序(RNA-seq)數(shù)據(jù)并獲得相應臨床患者的資料，從GTEx(Genotype-Tissue Expression，https://commonfund.nih.gov/gtex)數(shù)據(jù)庫中下載癌旁組織的RNA-seq，共下載373例具有完整生存信息的卵巢癌患者的癌組織和180例正常組織的測序信息。

1.2 m6A基因的選擇及差異表達分析

使用R(version 3.5.2)軟件進行統(tǒng)計分析和繪圖。根據(jù)相關m6A甲基化與腫瘤的研究[9]，選擇23個m6A甲基化基因(METTL3， METTL14，WTAP， RBM15， ZC3H13， METTL16， VIRMA， YTHDC1， YTHDC2， YTHDF1， YTHDF2， HNRNPC， YTHDF3， HNRNPA2B1， RBMX， IGF2BP1， IGF2BP2， IGF2BP3， EIF3A， FMR1， FTO， ALKBH5，ALKBH3)并使用Wilcoxon秩和檢驗(P<0.05)進行組間差異表達分析。

1.3 m6A基因風險模型的構建與檢驗

使用Least Absolute Shrinkage and Selection Operator(LASSO)分析篩選候選m6A甲基化基因，該過程主要使用“glmnet”和“survival” R包。進一步使用逐步Cox回歸分析構建卵巢癌風險預測模型并繪制森林圖。根據(jù)回歸系數(shù)加權建立卵巢癌的預后特征和風險評分公式，公式如下：

公式中N指基因數(shù)量，Expi是每個基因的表達量，Coei是多因素Cox回歸系數(shù)。計算每位患者的風險評分，并參照風險評分中位數(shù)將患者分為高風險組和低風險組，繪制相應的Kaplan-Meier曲線?；贚og-rank檢驗對高風險組和低風險組進行分析，分別繪制1、3、5年的ROC曲線并計算曲線下面積(AUC)，以檢驗風險評分模型的預測效能。進一步利用風險評分結合患者的年齡和腫瘤分級與分期，以觀察其是否能夠提高預測的準確度與價值，并使用多因素Cox回歸分析風險評分因素對卵巢癌預后的效應。

1.4 調(diào)控網(wǎng)絡的構建

應用Pearson相關分析篩選與風險評分模型中m6A基因相關的潛在靶基因(abs|r|>0.7，P<0.05)，應用Cytoscape (version 3.7.2)軟件進行調(diào)控網(wǎng)絡的可視化。

1.5 GO功能注釋和KEGG通路分析

應用“clusterProfile” R包對風險評估模型中m6A基因相關的潛在靶基因(abs|r|>0.7，P<0.05)進行GO功能富集和KEGG通路分析，以明確模型中基因潛在的生物學功能和參與的信號通路。

2 結 果

2.1 m6A基因在卵巢癌組織的差異表達分析

本研究373例卵巢癌患者的臨床信息顯示：患者年齡中位數(shù)為59歲，腫瘤惡性程度較高(見表1)。373例癌組織及180例正常組織樣本m6A甲基化基因的差異表達分析結果顯示：20個基因的

表1TCGA患者的臨床病理資料 n=373

差異表達具有統(tǒng)計學意義(P<0.05)，其中FMR1，HNRNPC， IGF2BP1， IGF2BP2， IGF2BP3， RBM15， VIRMA， YTHDF1， YTHDF2， YTHDF3在癌組織中上調(diào)，EIF3A， FTO， HNRNPA2B1， METTL14， METTL16， METTL3， RBMX， WTAP， YTHDC1， YTHDC2在癌組織中下調(diào)。(見圖1)。

2.2 確定3個m6A甲基化基因并應用于構建卵巢癌的風險特征

使用LASSO過濾篩選得到10個m6A甲基化基因(HNRNPA2B1，ZC3H13，WTAP，VIRMA，YTHDC2，HNRNPC，F(xiàn)MR1，IGF2BP1，RBMX，ALKBH5)(見圖2)。對篩選得到的10個基因進一步使用逐步Cox回歸分析，進而保留3個m6A甲基化基因(HNRNPA2B1，ZC3H13，WTAP)，其中HNRNPA2B1為保護基因，其HR為0.59(95%CI為0.45～0.77)(P<0.01)，而ZC3H13和WTAP為危險基因，其HR分別為1.29(95%CI為1.05～1.57)(P=0.01)和1.43(95%CI為1.15～1.78)(P<0.01)(見圖3)。

圖1 癌組織與正常組織m6A甲基化基因差異表達

圖2 LASSO模型篩選結果注：A為LASSO模型中最合適的log(Lambda)值；B為LASSO模型中m6A基因名及其回歸系數(shù)值

圖3 逐步Cox回歸分析結果

2.3 風險預測模型的構建與評價

基于風險評分模型計算患者風險評分：風險評分=(- 0.41×HNRNPA2B1表達水平)+(0.29×ZC3H13表達水平)+(0.43×WTAP表達水平)。根據(jù)風險評分中位數(shù)將患者分為高風險組和低風險組(見圖4)。Kaplan-Meier曲線及Log-rank檢驗顯示高風險組預后較低風險組差(P=0.001 9)(見圖5)。多因素Cox回歸分析顯示校正年齡、臨床分級和分期后，風險評分對卵巢癌生存預后的危害效應仍具有統(tǒng)計學意義(HR=2.643,P<0.01)(見圖6A)，并依據(jù)風險評分、年齡、臨床分級和分期繪制相應預后生存列線圖(見圖6B)。ROC曲線的AUC可用于評估風險評分模型預測效能，該風險評分模型1年、3年、5年的AUC分別為0.62、0.58、0.62，預測效能一般。當風險評分模型結合患者年齡、臨床分級和分期后，1年、3年、5年的AUC提高至0.74、0.64、0.64，提示該模型結合患者臨床資料將提高卵巢癌預后預測的準確性和價值(見圖6C)。

圖4 風險評分模型注：A為納入患者風險評分；B為隨訪終點高、低風險組生存狀態(tài)；C為風險預測模型中的m6A基因表達與風險等級

圖5 高、低風險組患者生存曲線

圖6 風險預測模型的評價注：A為多因素Cox回歸分析結果；B為列線圖預測患者生存；C為風險評分模型結合臨床病理資料預測AUC值(t代表生存月數(shù))

2.4 m6A甲基化基因與mRNA轉錄調(diào)控網(wǎng)絡的構建及功能分析

使用Pearson相關篩選與風險預測模型中3個m6A甲基化基因相關的mRNA(abs|r|>0.7，P<0.05)，最終構建其與919個相關mRNA的調(diào)控網(wǎng)絡(見圖7)。為研究預測模型的功能作用，將|r|>0.7的相關基因進行GO功能富集和KEGG通路分析(見圖8)，發(fā)現(xiàn)RNA剪接、RNA代謝的調(diào)控、RNA定位、RNA及核酸轉運等生物學功能明顯富集，而KEGG通路分析中發(fā)現(xiàn)靶基因參與RNA轉運、剪接體、mRNA監(jiān)視、蛋白水解、核糖體合成、細胞周期等信號途徑。上述結果提示，模型中的m6A甲基化基因在卵巢癌的發(fā)生發(fā)展中具有重要作用，可能通過多種癌癥通路調(diào)控卵巢癌細胞的惡性生物學進而影響患者預后。

圖7 m6A甲基化基因與相關mRNA調(diào)控網(wǎng)絡的構建

3 討 論

m6A甲基化作為最常見的RNA修飾，在多種癌癥中都發(fā)揮著不同的作用，包括腫瘤的發(fā)生、侵襲、增殖、轉移等[7]，m6A甲基化對腫瘤的發(fā)生發(fā)展是一把雙刃劍，因此需迫切了解m6A甲基化對各種癌癥的影響。本研究從TCGA數(shù)據(jù)庫及GTEx數(shù)據(jù)庫中下載了373例具有完整生存信息的卵巢癌患者的癌組織樣本及180例正常組織樣本測序信息，組間差異表達結果表明m6A甲基化基因在卵巢癌中起到重要作用。在LASSO回歸篩選出10個顯著的m6A甲基化基因后，進一步使用逐步Cox回歸分析保留了3個基因(HNRNPA2B1，ZC3H13，WTAP)并構建風險評分模型作為預測卵巢癌預后的指標，Kaplan-Meier曲線提示高風險組預后較低風險組差，校正患者臨床資料后風險評分模型仍具有統(tǒng)計學意義。ROC曲線的AUC值提示風險評估模型具備一定預測效能，模型結合患者臨床資料后其預測效能將進一步提升。

圖8 風險評分模型功能分析注 ：A為關聯(lián)基因的GO分析結果；B為關聯(lián)基因的KEGG通路分析結果

本研究模型中的HNRNPA2B1作為保護基因，能夠改善卵巢癌患者預后。HNRNPA2B1主要介導初級microRNA的加工及選擇性剪接[10],并在microRNA分選及包裝過程中也起到重要作用[11]。HNRNPA2B1對腫瘤的影響主要體現(xiàn)在其長鏈非編碼RNA(lncRNA)的調(diào)控[12],例如作為lncRNA的miR503HG可以與HNRNPA2B1特異性結合，抑制肝癌細胞的NF-κB信號通路，從而抑制腫瘤轉移[13]。另外，本研究的模型中的ZC3H13和WTAP作為危險基因，將提高卵巢癌患者的死亡風險，其中ZC3H13主要負責調(diào)節(jié)細胞核中RNA的m6A甲基化及細胞的自我更新[14]。有趣的是，ZC3H13可通過滅活Ras-ERK信號通路從而抑制結直腸癌的侵襲和增殖[15]，但其在肝細胞癌高風險患者中上調(diào)[16]，這提示在不同癌癥中相同的m6A甲基化基因可能存在不同的作用，在本研究中ZC3H13將導致預后不良。WTAP則通常與METTL3及METTL14相互作用，負責調(diào)節(jié)轉錄、RNA加工相關基因的表達和選擇性剪接[17]，并主要通過參與RNA加工和細胞增殖發(fā)揮致癌作用[18]。在肝癌[19]、膀胱癌[20]、膽管癌[21]、胃癌[22]中WTAP均顯示明顯的過量表達并與預后不良顯著相關，這表明大多數(shù)癌癥中WTAP的高表達將會促進癌癥發(fā)展并導致不良預后，在本研究中WTAP的高表達同樣可能導致卵巢癌預后不佳。

綜上，基于TCGA數(shù)據(jù)庫和GTEx數(shù)據(jù)庫，我們分析了卵巢癌組織與正常組織m6A甲基化基因的差異表達情況，隨后建立了一個三基因聯(lián)合風險評分預測模型并構建相關調(diào)控網(wǎng)絡及進行GO功能注釋和KEGG通路分析。我們發(fā)現(xiàn)該風險評估模型具備一定的預后預測效能，該模型可能是卵巢癌預后潛在的預測因子及治療靶標，然而其具體的分子機制和相關信號通路仍有待進一步研究與探索。