張樂(lè)，寧?kù)o華，張?chǎng)?，王振，劉虹汝，張鈺哲?，3

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

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

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

1 材料與方法

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

從TCGA 數(shù)據(jù)庫(kù) （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語(yǔ)言 （4.2.1） limma包，分析13種m6A調(diào)節(jié)因子在前列腺癌和正常前列腺組織中的差異表達(dá)。

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

1.3.3 m6A調(diào)節(jié)因子預(yù)后風(fēng)險(xiǎn)模型的構(gòu)建：為了研究聚類(lèi)分析與臨床病理特征的相關(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中的差異基因，對(duì)這些基因進(jìn)行后續(xù)分析。利用基因本體論 （gene ontology，GO） 以及京都基因和基因組數(shù)據(jù)庫(kù) （Kyoto Encyclopedia of Genes and Genomes，KEGG） 富集分析，對(duì)差異基因進(jìn)行功能注釋。對(duì)差異基因進(jìn)行單因素Cox回歸分析，選擇與生存密切相關(guān)的基因，應(yīng)用LASSO Cox回歸算法構(gòu)建強(qiáng)大的預(yù)后特征，并計(jì)算每個(gè)基因的系數(shù)［11］。計(jì)算每例患者的風(fēng)險(xiǎn)評(píng)分，為每個(gè)基因評(píng)分的總和。根據(jù)前列腺癌患者風(fēng)險(xiǎn)評(píng)分的中位數(shù)，將前列腺癌分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，繪制受試者操作特征曲線(xiàn)，計(jì)算曲線(xiàn)下面積 （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á)和免疫浸潤(rùn)分析

使用癌癥基因組分析平臺(tái) （http：//bioinfo.life.hust.edu.cn/web/GSCALite/），對(duì)ROMO1在前列腺免疫微環(huán)境中的免疫浸潤(rùn)情況進(jìn)行分析。為了探討ROMO1的表達(dá)與前列腺癌免疫細(xì)胞的關(guān)系，使用癌癥基因組分析平臺(tái)中的Immune工具，分析ROMO1表達(dá)與前列腺癌免疫細(xì)胞浸潤(rùn)水平的關(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），超過(guò)半數(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í)聚類(lèi)鑒定2種不同的前列腺癌亞型

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

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

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

在TCGA數(shù)據(jù)庫(kù)下載前列腺癌患者的臨床信息，對(duì)患者的分型和臨床病理特征進(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)，聚類(lèi)結(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 聚類(lèi)間差異基因分析

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

圖4 利用生物過(guò)程GO分析和KEGG通路分析對(duì)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í)，說(shuō)明該基因是高風(fēng)險(xiǎn)基因；HR<1時(shí)，說(shuō)明該基因是低風(fēng)險(xiǎn)基因。結(jié)果表明，SNORD60、SNHG25、SNORD104、C4orf48、OAS3、NME3、ROMO1、SNHG9、TMEM238、MRPL41、SNHG19這些基因均HR>1，說(shuō)明以上基因表達(dá)水平較高時(shí)，前列腺癌患者的生存期較差。而SMIM22和MYOF均HR<1，說(shuō)明這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)模型來(lái)預(yù)測(cè)樣本的風(fēng)險(xiǎn)得分。為了研究風(fēng)險(xiǎn)模型預(yù)后作用，根據(jù)LASSO Cox算法得到的風(fēng)險(xiǎn)評(píng)分，利用風(fēng)險(xiǎn)評(píng)分的中位數(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)模型對(duì)前列腺癌的臨床預(yù)后有較好參考價(jià)值 （P= 1E-07）。提示高風(fēng)險(xiǎn)組前列腺癌患者的生存期較差 （圖5C）。

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

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

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

2.7 ROMO1在前列腺癌中的免疫浸潤(rùn)情況

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

圖7 ROMO1在前列腺癌免疫微環(huán)境中的免疫浸潤(rù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ì)胞的免疫浸潤(rùn)相關(guān)，其有望獨(dú)立作為前列腺癌臨床預(yù)后的指標(biāo)。

ROMO1可以通過(guò)干擾活性氧的產(chǎn)生，影響細(xì)胞的狀態(tài)［12］。通過(guò)釋放細(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ì)胞浸潤(rùn)水平和免疫途徑相關(guān)，與本研究結(jié)果一致，說(shuō)明了本研究結(jié)果的準(zhǔn)確性。本研究發(fā)現(xiàn)，ROMO1與前列腺癌的CD8+T細(xì)胞間存在免疫浸潤(rùn)關(guān)系，這也許為前列腺癌的免疫療法提供了新的治療方案。在機(jī)體正常情況下，體內(nèi)的T細(xì)胞可以通過(guò)識(shí)別異常的癌細(xì)胞，啟動(dòng)抗腫瘤的免疫反應(yīng)，T細(xì)胞越多，人體發(fā)動(dòng)免疫反應(yīng)抵抗癌細(xì)胞的效果就越好。CD8+T細(xì)胞對(duì)腫瘤細(xì)胞具有較強(qiáng)的殺傷功能。因此，從免疫治療的角度考慮，可能通過(guò)調(diào)節(jié)ROMO1的表達(dá)導(dǎo)致前列腺癌免疫細(xì)胞浸潤(rùn)的改變，從而影響前列腺癌微環(huán)境的改變，進(jìn)而導(dǎo)致前列腺癌異質(zhì)性，從而達(dá)到精準(zhǔn)化個(gè)體治療的目的。

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