綜合分析和鑒定m6A RNA甲基化調(diào)節(jié)因子對前列腺癌進(jìn)展及預(yù)后的影響

張樂，寧靜華，張鑫，王振，劉虹汝，張鈺哲，2，3

（1.大理大學(xué)基礎(chǔ)醫(yī)學(xué)院病理學(xué)與病理生理學(xué)教研室，云南 大理 671000；2.云南抗病原藥用植物篩選重點(diǎn)實(shí)驗(yàn)室，云南 大理 671000；3.云南省昆蟲生物醫(yī)藥重點(diǎn)實(shí)驗(yàn)室，云南 大理 671000）

前列腺癌是中老年男性常見的惡性腫瘤，已成為威脅人類健康的公共問題［1］。目前，前列腺癌的治療方法包括手術(shù)、內(nèi)分泌治療、放化療、冷凍治療和免疫治療等［2］，但由于人們對前列腺癌的早期篩查不夠重視，通常確診時(shí)已為晚期或轉(zhuǎn)移。前列腺癌患者即使進(jìn)行積極的雄激素剝奪治療，但由于長時(shí)間的藥物治療，患者會(huì)產(chǎn)生耐藥性，從而演變?yōu)槿莸挚剐郧傲邢侔?］。因此，尋找更為有效的早期診斷標(biāo)志物和治療方案勢在必行。

N6-甲基腺苷 （N6-methyladenosine，m6A） 是最常見的RNA甲基化修飾。m6A修飾與腫瘤的發(fā)生、分化、增殖、侵襲和轉(zhuǎn)移有關(guān)［4-7］。本研究基于以往研究最廣泛報(bào)道的13種m6A RNA甲基化調(diào)節(jié)因子 （以下簡稱m6A調(diào)節(jié)因子），利用癌癥基因組圖譜（The Cancer Genome Atlas，TCGA） 數(shù)據(jù)庫中的RNA測序數(shù)據(jù)，系統(tǒng)分析了498例前列腺癌中這13種m6A調(diào)節(jié)因子的差異表達(dá)與臨床病理特征的關(guān)聯(lián)，評估m(xù)6A甲基化與前列腺癌進(jìn)展和生存的相關(guān)性，為從甲基化修飾角度理解前列腺癌的發(fā)生和發(fā)展提供參考。

1 材料與方法

1.1 數(shù)據(jù)采集

從TCGA 數(shù)據(jù)庫 （https：//cancergenome.nih.gov/）中獲取TCGA-PRAD數(shù)據(jù)，下載前列腺癌和癌旁樣本的基因表達(dá)數(shù)據(jù)。排除同一患者的重復(fù)樣本后，共有498例前列腺癌組織樣本和52例正常前列腺組織樣本。

1.2 m6A調(diào)節(jié)因子的選擇

根據(jù)之前的研究［8-10］，篩選出最廣泛報(bào)道的13種m6A調(diào)節(jié)因子：METTL3、METTL14、WTAP、KIAA1429、RBM15、ZC3H13、YTHDC1、YTHDC2、YTHDF1、YTHDF2、HNRNPC、FTO、ALKBH5。提取13種基因的表達(dá)矩陣和臨床樣本的信息。

1.3 生物信息學(xué)分析

1.3.1 差異基因分析：為了研究m6A調(diào)節(jié)因子在前列腺癌中的作用，應(yīng)用R語言 （4.2.1） limma包，分析13種m6A調(diào)節(jié)因子在前列腺癌和正常前列腺組織中的差異表達(dá)。

1.3.2 利用m6A調(diào)節(jié)因子在癌癥中的表達(dá)水平進(jìn)行無監(jiān)督聚類分析：為了確定m6A調(diào)節(jié)因子是否與前列腺癌的預(yù)后相關(guān)，將498例前列腺癌樣本進(jìn)行分組。使用 R語言Consensus Cluster Plus包執(zhí)行，根據(jù)集群和項(xiàng)目一致性原理，輸出圖形。集群數(shù)量通過累積分布函數(shù)曲線確定，為了驗(yàn)證根據(jù)13種m6A調(diào)節(jié)因子表達(dá)量進(jìn)行的聚類分組在所有基因中的準(zhǔn)確性，進(jìn)行主成分分析，驗(yàn)證聚類分組效果。

1.3.3 m6A調(diào)節(jié)因子預(yù)后風(fēng)險(xiǎn)模型的構(gòu)建：為了研究聚類分析與臨床病理特征的相關(guān)性，從TCGA中下載前列腺癌患者的臨床信息，使用limma包分析13種m6A調(diào)節(jié)因子在Cluster 1和Cluster 2中的表達(dá)情況，并研究這13種m6A調(diào)節(jié)因子與前列腺癌患者的臨床病理特征之間的相關(guān)性。隨后篩選出Cluster 1與 Cluster 2中的差異基因，對這些基因進(jìn)行后續(xù)分析。利用基因本體論 （gene ontology，GO） 以及京都基因和基因組數(shù)據(jù)庫 （Kyoto Encyclopedia of Genes and Genomes，KEGG） 富集分析，對差異基因進(jìn)行功能注釋。對差異基因進(jìn)行單因素Cox回歸分析，選擇與生存密切相關(guān)的基因，應(yīng)用LASSO Cox回歸算法構(gòu)建強(qiáng)大的預(yù)后特征，并計(jì)算每個(gè)基因的系數(shù)［11］。計(jì)算每例患者的風(fēng)險(xiǎn)評分，為每個(gè)基因評分的總和。根據(jù)前列腺癌患者風(fēng)險(xiǎn)評分的中位數(shù)，將前列腺癌分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，繪制受試者操作特征曲線，計(jì)算曲線下面積 （area under the curve，AUC），獲得預(yù)后特征的靈敏度和特異度。

1.4 風(fēng)險(xiǎn)因子的表達(dá)水平

進(jìn)行單因素Cox回歸分析，從297個(gè)差異表達(dá)基因中選出P< 0.01的基因，利用LASSO Cox算法選取誤差最小的點(diǎn)，獲得6個(gè)風(fēng)險(xiǎn)基因并構(gòu)建風(fēng)險(xiǎn)模型。

1.5 基因表達(dá)和免疫浸潤分析

使用癌癥基因組分析平臺(tái) （http：//bioinfo.life.hust.edu.cn/web/GSCALite/），對ROMO1在前列腺免疫微環(huán)境中的免疫浸潤情況進(jìn)行分析。為了探討ROMO1的表達(dá)與前列腺癌免疫細(xì)胞的關(guān)系，使用癌癥基因組分析平臺(tái)中的Immune工具，分析ROMO1表達(dá)與前列腺癌免疫細(xì)胞浸潤水平的關(guān)系。

2 結(jié)果

2.1 前列腺癌中m6A調(diào)節(jié)因子的關(guān)聯(lián)性分析

與正常前列腺組織相比，前列腺癌患者中METTL3、METTL14、WTAP、KIAA1429、RBM15、ZC3H13、YTHDC1、YTHDC2、YTHDF1、YTHDF2、HNRNPC、FTO、ALKBH5的表達(dá)水平較高 （圖1A、1B），超過半數(shù) （11/13） 的m6A調(diào)節(jié)因子在前列腺癌中高表達(dá)。m6A調(diào)節(jié)因子與前列腺癌的TNM分期、年齡以及生存狀態(tài)相關(guān)，尤其與T3、M0、N0和年齡≤65歲患者的臨床特征顯著相關(guān) （圖1C）。

圖1 前列腺癌中m6A調(diào)節(jié)因子的分布Fig.1 Distribution of m6A regulators in prostate cancer

2.2 m6A調(diào)節(jié)因子的共識(shí)聚類鑒定2種不同的前列腺癌亞型

為了提供更為準(zhǔn)確的分型，通過498例前列腺癌組織樣本中13個(gè)m6A調(diào)節(jié)因子的mRNA 表達(dá)共識(shí)，將前列腺癌分為幾個(gè)集群。聚類指數(shù) k=2 時(shí)，是聚類間最大差異的最佳點(diǎn)。從m6A調(diào)節(jié)因子的表達(dá)相似性來看，選擇一致累積分布函數(shù) （cumulative distribution function，CDF） 坡度較小的曲線對應(yīng)的k值為最佳。結(jié)果發(fā)現(xiàn)，k=3時(shí)，似乎具有較小的CDF值（圖 2A、2B），說明將前列腺癌組織樣本分為3個(gè)聚類最佳。但將其分為3組后，組間相關(guān)性較高，樣本數(shù)量較少，因此分為2組 （Cluster 1和Cluster 2，圖2C）。為了驗(yàn)證聚類的分組效果以及說明利用m6A調(diào)節(jié)因子進(jìn)行的聚類分組在所有基因中的適用性，進(jìn)行主成分分析。結(jié)果顯示，可以將Cluster 1與 Cluster 2區(qū)分開來 （圖2D），說明利用m6A調(diào)節(jié)因子的分類方法合理且準(zhǔn)確。

圖2 m6A調(diào)節(jié)因子一致聚類的識(shí)別Fig.2 Identification of consistent clusters of m6A regulators

2.3 通過共識(shí)聚類確定的類別與臨床病理特征密切相關(guān)

在TCGA數(shù)據(jù)庫下載前列腺癌患者的臨床信息，對患者的分型和臨床病理特征進(jìn)行分析，發(fā)現(xiàn)大部分m6A調(diào)節(jié)因子在Cluster 1中高表達(dá) （圖3）。與Cluster 2相比，Cluster 1與診斷時(shí)的年齡 （≤65歲）、TNM分期 （T3、M0、N0）、生存狀態(tài) （存活） 顯著相關(guān)，聚類結(jié)果與前列腺癌的惡性程度密切相關(guān)。

圖3 Cluster 1、Cluster 2 與前列腺癌臨床病理特征的關(guān)系Fig.3 Relationship between Clusters 1 and 2 and the clinicopathological features of prostate cancer

2.4 聚類間差異基因分析

從Cluster 1和Cluster 2中篩選出297個(gè)差異基因。GO分析結(jié)果表明，差異基因主要富集在惡性腫瘤的細(xì)胞基質(zhì)黏附、黏附連接組織、線粒體呼吸鏈復(fù)合體組裝、磷脂酰絲氨酸代謝過程、上皮細(xì)胞遷移的正向調(diào)節(jié)、細(xì)胞連接組織、氧化磷酸化作用等經(jīng)典途徑 （圖4A、4B）。KEGG結(jié)果顯示，差異基因富集于腫瘤的蛋白多糖、細(xì)胞外基質(zhì)受體相互作用、轉(zhuǎn)化生長因子-β信號通路、細(xì)胞骨架作用的調(diào)節(jié)、黏著斑、黏著連接 （圖4C、4D）。

圖4 利用生物過程GO分析和KEGG通路分析對Cluster 1亞群中表達(dá)較高的基因進(jìn)行功能注釋Fig.4 Functional annotation of highly expressed genes in the Cluster 1 subgroup using bioprocess GO and KEGG pathway analyses

2.5 m6A調(diào)節(jié)因子的風(fēng)險(xiǎn)特征與預(yù)后價(jià)值

從297個(gè)顯著差異表達(dá)基因中，選出在單因素Cox回歸分析中P< 0.01的基因，繪制森林圖 （圖5A）。風(fēng)險(xiǎn)比 （hazard ratio，HR） >1時(shí)，說明該基因是高風(fēng)險(xiǎn)基因；HR<1時(shí)，說明該基因是低風(fēng)險(xiǎn)基因。結(jié)果表明，SNORD60、SNHG25、SNORD104、C4orf48、OAS3、NME3、ROMO1、SNHG9、TMEM238、MRPL41、SNHG19這些基因均HR>1，說明以上基因表達(dá)水平較高時(shí)，前列腺癌患者的生存期較差。而SMIM22和MYOF均HR<1，說明這2種基因表達(dá)水平較高時(shí)，前列腺癌患者有更好的生存期。后續(xù)研究中，由這些基因代替m6A調(diào)節(jié)因子。

圖5 6個(gè)m6A調(diào)節(jié)因子的風(fēng)險(xiǎn)特征Fig.5 Risk profiles of six m6A regulators

進(jìn)行多因素Cox回歸分析，利用LASSO Cox算法得到交叉驗(yàn)證的誤差，選取誤差最小的點(diǎn)并獲得6個(gè)風(fēng)險(xiǎn)基因 （圖5B），構(gòu)建風(fēng)險(xiǎn)模型來預(yù)測樣本的風(fēng)險(xiǎn)得分。為了研究風(fēng)險(xiǎn)模型預(yù)后作用，根據(jù)LASSO Cox算法得到的風(fēng)險(xiǎn)評分，利用風(fēng)險(xiǎn)評分的中位數(shù)，將TCGA數(shù)據(jù)集中的前列腺癌患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，進(jìn)行生存分析，結(jié)果顯示，6個(gè)基因 （OAS3、MYOF、SMIM22、SNORD60、ROMO1、MRPL41） 構(gòu)建的風(fēng)險(xiǎn)模型對前列腺癌的臨床預(yù)后有較好參考價(jià)值 （P= 1E-07）。提示高風(fēng)險(xiǎn)組前列腺癌患者的生存期較差 （圖5C）。

2.6 風(fēng)險(xiǎn)評分與前列腺癌的臨床病理特征密切相關(guān)

6個(gè)風(fēng)險(xiǎn)基因與前列腺癌患者的年齡、分級、分期、T分期、M分期、N分期以及生存狀態(tài)存在相關(guān)性 （圖6A）。通過計(jì)算AUC，表明風(fēng)險(xiǎn)特征 （6個(gè)基因） 預(yù)測的準(zhǔn)確率，AUC值越高，模型預(yù)測生存率越高。結(jié)果表明，風(fēng)險(xiǎn)評分可以預(yù)測前列腺癌的生存率 （AUC=0.727） （圖6B）。單因素Cox回歸分析結(jié)果顯示，T分期、N分期、風(fēng)險(xiǎn)評分均可以作為高風(fēng)險(xiǎn)因素與風(fēng)險(xiǎn)評分相關(guān) （圖6C）。多因素Cox回歸分析結(jié)果顯示，風(fēng)險(xiǎn)評分可作為獨(dú)立預(yù)后因子 （圖6D）。以上證據(jù)表明，風(fēng)險(xiǎn)評分與前列腺癌臨床病理特征密切相關(guān)，并且SNORD60、ROMO1、MRPL41與前列腺癌的惡性程度相關(guān)，特別是在高風(fēng)險(xiǎn)的前列腺癌中ROMO1表達(dá)水平更高。

圖6 高低風(fēng)險(xiǎn)組中風(fēng)險(xiǎn)評分和臨床病理特征的關(guān)系Fig.6 Relationship between risk scores and clinicopathological characteristics in high- and low-risk groups

2.7 ROMO1在前列腺癌中的免疫浸潤情況

結(jié)果發(fā)現(xiàn)，ROMO1與前列腺癌的CD8+T細(xì)胞的免疫浸潤相關(guān)性強(qiáng) （圖7），這也揭示了ROMO1與前列腺癌免疫浸潤有一定關(guān)系。

圖7 ROMO1在前列腺癌免疫微環(huán)境中的免疫浸潤情況Fig.7 Immune infiltration by ROMO1 in the prostate cancer immune microenvironment

3 討論

本研究結(jié)果顯示，m6A調(diào)節(jié)因子的表達(dá)與前列腺癌的臨床病理特征高度相關(guān)；使用6種基因構(gòu)建的臨床風(fēng)險(xiǎn)模型可以作為預(yù)后指標(biāo)；ROMO1與高風(fēng)險(xiǎn)組前列腺癌的惡性程度相關(guān)性較高；ROMO1與前列腺癌CD8+T細(xì)胞的免疫浸潤相關(guān)，其有望獨(dú)立作為前列腺癌臨床預(yù)后的指標(biāo)。

ROMO1可以通過干擾活性氧的產(chǎn)生，影響細(xì)胞的狀態(tài)［12］。通過釋放細(xì)胞色素c激活細(xì)胞凋亡途徑，最終導(dǎo)致細(xì)胞死亡，從而引發(fā)癌癥。從正常細(xì)胞轉(zhuǎn)化為癌細(xì)胞，ROMO1作為一種活性氧調(diào)節(jié)蛋白，已在多種腫瘤細(xì)胞中被發(fā)現(xiàn)，并在腫瘤的發(fā)生、發(fā)展中起作用［13］。有研究［14］發(fā)現(xiàn)，ROMO1在前列腺癌組織中高表達(dá)，且與前列腺癌的免疫細(xì)胞浸潤水平和免疫途徑相關(guān)，與本研究結(jié)果一致，說明了本研究結(jié)果的準(zhǔn)確性。本研究發(fā)現(xiàn)，ROMO1與前列腺癌的CD8+T細(xì)胞間存在免疫浸潤關(guān)系，這也許為前列腺癌的免疫療法提供了新的治療方案。在機(jī)體正常情況下，體內(nèi)的T細(xì)胞可以通過識(shí)別異常的癌細(xì)胞，啟動(dòng)抗腫瘤的免疫反應(yīng)，T細(xì)胞越多，人體發(fā)動(dòng)免疫反應(yīng)抵抗癌細(xì)胞的效果就越好。CD8+T細(xì)胞對腫瘤細(xì)胞具有較強(qiáng)的殺傷功能。因此，從免疫治療的角度考慮，可能通過調(diào)節(jié)ROMO1的表達(dá)導(dǎo)致前列腺癌免疫細(xì)胞浸潤的改變，從而影響前列腺癌微環(huán)境的改變，進(jìn)而導(dǎo)致前列腺癌異質(zhì)性，從而達(dá)到精準(zhǔn)化個(gè)體治療的目的。

綜上所述，本研究結(jié)果系統(tǒng)展示了m6A調(diào)節(jié)因子在前列腺癌中的表達(dá)、潛在功能和預(yù)后價(jià)值。通過分析篩選出的目標(biāo)基因ROMO1與前列腺癌惡性程度高度相關(guān)，有望作為前列腺癌臨床預(yù)后的獨(dú)立指標(biāo)，為去勢抵抗性前列腺癌的治療提供理論基礎(chǔ)。