N6-甲基腺嘌呤調(diào)控因子對(duì)宮頸癌預(yù)后及免疫浸潤(rùn)的作用①

郭依琳 白洋洋 王 璐 徐 臻 王習(xí)亮 王武亮 （鄭州大學(xué)第二附屬醫(yī)院婦產(chǎn)科，鄭州 450014）

宮頸癌是女性生殖系統(tǒng)最常見的惡性腫瘤，其 中90%發(fā)生于發(fā)展中國家［1］。2020 年，全球?qū)m頸癌新發(fā)病例604 127 例，死亡病例341 831 例［2］。隨著宮頸癌篩查系統(tǒng)逐步開展及宮頸癌疫苗普及，早期宮頸癌預(yù)后較好，但仍有部分宮頸癌患者首診時(shí)已為晚期，失去了最佳治療時(shí)機(jī)。因此闡明宮頸癌發(fā)生發(fā)展的分子機(jī)制對(duì)改善宮頸癌患者預(yù)后具有重要臨床意義。

m6A 是真核生物最豐富的編碼和非編碼RNA甲基化修飾［3］。m6A 是位于腺嘌呤堿基第6 位氮原子的甲基化修飾，主要分布于終止密碼子附近和3′-非翻譯區(qū)，定位于保守的RRACH 序列（R=G 或A；H=A，C或U），在RNA加工成熟、剪接和轉(zhuǎn)錄調(diào)控中發(fā)揮重要作用［4-5］。m6A 受甲基轉(zhuǎn)移酶復(fù)合物（writers）、去甲基化酶（erasers）和m6A 讀取蛋白（readers）調(diào)控［6］。最近研究提出m6A 調(diào)控因子具有免疫調(diào)節(jié)功能，參與腫瘤免疫微環(huán)境（tumor immune microenvironment，TIME）的免疫調(diào)控，但其表達(dá)對(duì)宮頸癌免疫浸潤(rùn)和預(yù)后的研究未見報(bào)道。因此，本研究利用TCGA 和GEO 數(shù)據(jù)庫分析m6A 調(diào)控因子在宮頸癌中的表達(dá)，建立LASSO Cox回歸預(yù)后模型，利用ssGSEA 計(jì)算免疫細(xì)胞浸潤(rùn)程度，為宮頸癌患者臨床診療和預(yù)后改善提供依據(jù)。

1 資料與方法

1.1 TCGA 和GEO 數(shù) 據(jù) 獲 取 從TCGA 數(shù) 據(jù) 庫（http：//cancergenome.nih.gov/）下載306例宮頸癌組織和3 例正常組織RNA 測(cè)序轉(zhuǎn)錄組數(shù)據(jù)（FPKM 格式）及相應(yīng)臨床信息。從GEO 數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo/）下載55 例宮頸癌組織和17 例正常組織（GSE52903）芯片轉(zhuǎn)錄組數(shù)據(jù)及相應(yīng)臨床信息。采用R 軟件（4.1.2）“l(fā)imma”包將TCGA數(shù)據(jù)集中FPKM 格式的RNA 測(cè)序轉(zhuǎn)錄組數(shù)據(jù)轉(zhuǎn)換為TPM 格式，并與GSE52903 數(shù)據(jù)集合并，對(duì)合并數(shù)據(jù)進(jìn)行矯正。

1.2 m6A 調(diào)控因子選擇及相關(guān)性分析 選擇25個(gè)公認(rèn)的m6A 調(diào)控因子，包括8 個(gè)writers（METTL3、METTL5、METTL14、METTL16、WTAP、ZC3H13、RBM15、RBM15B）、14 個(gè)readers（YTHDC1、YTHDC2、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、LRPPRC、HNRNPA2B1、IGF2BP1、IGF2BP2、IGF2BP3、ELAVL1、RBMX）和3 個(gè)erasers（FTO、ALKBH5、ALKBH3）。采用R 軟件“l(fā)imma”包對(duì)25 個(gè)m6A 調(diào)控因子進(jìn)行Wilcox 秩和檢驗(yàn)比較宮頸癌組織和正常組織表達(dá)差異?！皃heatmap”包繪制m6A 調(diào)控因子表達(dá)和臨床特征熱圖。“corrplot”包比較m6A 調(diào)控因子間的相關(guān)性。

1.3 LASSO 預(yù)后模型構(gòu)建及驗(yàn)證 采用R 軟件“Survival”包對(duì)25 個(gè)m6A 調(diào)控因子進(jìn)行單因素Cox回歸分析，“forestplot”包繪制森林圖。選取5個(gè)預(yù)后相關(guān)調(diào)控因子（P<0.05），使用LASSO Cox 回歸算法建立宮頸癌預(yù)后特征和風(fēng)險(xiǎn)評(píng)分，風(fēng)險(xiǎn)評(píng)分=Codfi×xi，N 為m6A 調(diào)控因子數(shù)量，Codfi 為系數(shù)，xi 為各基因Z-score 轉(zhuǎn)換的相對(duì)表達(dá)。該公式可計(jì)算每例宮頸癌患者風(fēng)險(xiǎn)評(píng)分，并將其分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，繪制相應(yīng)Kaplan-Meier 曲線。繪制ROC 曲線和計(jì)算AUC 檢驗(yàn)預(yù)后模型的預(yù)測(cè)效能。利用單變量和多變量Cox回歸分析評(píng)價(jià)預(yù)后模型和臨床特征對(duì)宮頸癌預(yù)后的價(jià)值。

1.4 m6A 調(diào)控因子在組織中的表達(dá)驗(yàn)證 人類蛋白圖譜數(shù)據(jù)庫（Human Protein Atlas，HPA，https：//www.proteinatlas.org/）分析預(yù)后相關(guān)m6A 調(diào)控因子（HNRNPC、LRPPRC、ELAVL1 和FMR1）在宮頸癌組織及配對(duì)的正常宮頸組織的表達(dá)。收集鄭州大學(xué)第二附屬醫(yī)院就診的15 例宮頸癌患者腫瘤組織和正常宮頸組織，qRT-PCR 檢測(cè)METTL3表達(dá)，本研究經(jīng)鄭州大學(xué)第二附屬醫(yī)院倫理委員會(huì)批準(zhǔn)（2021040），采用Trizol 法提取RNA，測(cè)定RNA 濃度和純度，逆轉(zhuǎn)錄為cDNA，得到的cDNA 以ACTIN基因?yàn)閮?nèi)參進(jìn)行PCR 反應(yīng)。PCR 反應(yīng)體系為SYBR?Premix 10 μl、PCR正向引物 （10 μmol/L） 0.5 μl、PCR反向引物（10 μmol/L）0.5 μl、DNA 模板1 μl、ddH2O 8 μl；PCR 反應(yīng)條件為95 ℃ 1 min；95 ℃ 15 s，57 ℃30 s，共40個(gè)循環(huán)。引物序列為：METTL3 F：5'-CCAGCACAGCTTCAGCAGTTCC-3'，R：5'-GCGTGGAGATGGCAAGACAGATG-3'。

1.5 m6A 調(diào)控因子聚類分型 采用R 軟件“ConsensusClusterPlus”包對(duì)25 個(gè)m6A 調(diào)控因子進(jìn)行樣本聚類分析。參數(shù)設(shè)定如下：iterations=50，resample rate=80%，分析方法為Pearson 相關(guān)性分析。選擇k=2 對(duì)樣本進(jìn)行分型，將樣本分為m6Acluster A組和m6Acluster B 組。主成分分析算法（principal component analysis，PCA）比較兩種聚類分型m6A 調(diào)控因子表達(dá)譜。

1.6 GSEA 分析 對(duì)聚類分型的m6Acluster A 和m6Acluster B 進(jìn)行GSEA 分析。從MSigDB 數(shù)據(jù)庫（http：//www.broadinstitute.org/gsea/msigdb/index.jsp）下載“c2.cp.kegg.v7.4.symbols.gmt”基因集。R 軟件“GSVA”包進(jìn)行GSEA 分析，調(diào)整后的P<0.05 為差異有統(tǒng)計(jì)學(xué)意義。采用“clusterProfiler”包對(duì)m6A調(diào)控因子進(jìn)行功能注釋，明確模型中潛在的生物學(xué)功能和參與的信號(hào)通路。

1.7 免疫浸潤(rùn)分析 從Pornpimol Charoentong 獲得TIME 浸潤(rùn)免疫細(xì)胞基因集［7］。該基因集包括多種人類免疫細(xì)胞亞型，如活化的CD8+T細(xì)胞、活化的樹突狀細(xì)胞、巨噬細(xì)胞、自然殺傷T 細(xì)胞、調(diào)節(jié)性T細(xì)胞等。ssGSEA計(jì)算聚類分型的兩組樣本TMIE中免疫細(xì)胞浸潤(rùn)相對(duì)豐度。

2 結(jié)果

2.1 宮頸癌中m6A 調(diào)控因子表達(dá)和相關(guān)性 剔除重復(fù)及不完整數(shù)據(jù)后，共收集宮頸癌樣本356例，正常樣本20 例。在376 例樣本中對(duì)25 個(gè)m6A 調(diào)控因子進(jìn)行分析，結(jié)果顯示，12 個(gè)m6A 調(diào)控因子表達(dá)在宮頸癌樣本和正常樣本中存在差異，其中RBM15、IGF2BP2、IGF2BP3、FMR1、YTHDF1、METTL3 和YTHDF2 表達(dá)上調(diào)，ZC3H13、METTL14、YTHDC1、FTO、ALKBH3 表達(dá)下調(diào)（P<0.05，圖1A、B）。計(jì)算25個(gè)m6A調(diào)控因子表達(dá)的相關(guān)性，發(fā)現(xiàn)m6A調(diào)控因子表達(dá)均呈正相關(guān)，其中METTL14和YTHDC1正相關(guān)性最為顯著（r=0.76，圖1C、D）。

2.2 m6A 調(diào)控因子對(duì)宮頸癌預(yù)后的評(píng)估價(jià)值 為研究m6A 調(diào)控因子對(duì)宮頸癌預(yù)后的影響，對(duì)25 個(gè)m6A 調(diào)控因子進(jìn)行單變量Cox 回歸，結(jié)果顯示FMR1（P=0.017，HR：0.680，95%CI：0.496～0.933）、HNRNPC（P=0.016，HR：2.065，95%CI：1.145～3.725）、LRPPRC（P=0.011，HR：1.711，95%CI：1.133～2.584）、ELAVL1（P=0.028，HR：0.569，95%CI：0.343～0.942）和METTL3（P=0.038，HR：1.551，95%CI：1.026～2.346）是獨(dú)立的預(yù)后因素（圖2A）。對(duì)5 個(gè)m6A 調(diào)控因子進(jìn)行Kaplan-Meier 生存分析，結(jié)果顯示HNRNPC、LRPPRC 和METTL3 均為危險(xiǎn)調(diào)控因子，其高表達(dá)與較差的總生存期（overall survival，OS）顯著相關(guān)；而ELAVL1 和FMR1 為保護(hù)調(diào)節(jié)因子，其高表達(dá)與較好的OS顯著相關(guān)（圖2B～F）。

2.3 宮頸癌風(fēng)險(xiǎn)預(yù)后模型構(gòu)建及組織樣本驗(yàn)證

選 擇HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1 建立風(fēng)險(xiǎn)預(yù)后模型，通過LASSO Cox 回歸模型計(jì)算宮頸癌患者風(fēng)險(xiǎn)評(píng)分（圖3A、B）。HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1 的系數(shù)分別為0.686 4、0.383 6、0.428 5、-0.880 4 和-0.434 4。根據(jù)中位風(fēng)險(xiǎn)評(píng)分將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。高風(fēng)險(xiǎn)組OS 與低風(fēng)險(xiǎn)組OS 相比明顯降低（P=2.554e-08，圖3C）。ROC 曲線評(píng)估預(yù)后特征的特異度和敏感度，AUC=0.734，表明預(yù)測(cè)性能良好（圖3D）。為驗(yàn)證篩選的5個(gè)m6A 調(diào)控因子，從HPA數(shù)據(jù)庫分析HNRNPC、LRPPRC、ELAVL1 和FMR1在正常組織和宮頸癌組織中的免疫組化表達(dá)，發(fā)現(xiàn)LRPPRC和HNRNPC在宮頸癌組織中的表達(dá)高于正常組織，ELAVL1和FMR1在正常組織中的表達(dá)高于宮頸癌組織（圖4A～H）。qRT-PCR 檢測(cè)收集的15對(duì)宮頸癌和正常組織中METTL3 表達(dá)，發(fā)現(xiàn)宮頸癌組織中METTL3表達(dá)更高（P<0.000 1，圖4I）。

圖4 組織中5個(gè)m6A調(diào)控因子表達(dá)驗(yàn)證Fig.4 Validation of five m6A regulators expressions in tissues

2.4 風(fēng)險(xiǎn)預(yù)后模型在宮頸癌中的臨床應(yīng)用價(jià)值

繪制5 個(gè)m6A 調(diào)控因子在高、低風(fēng)險(xiǎn)組中的表達(dá)及臨床特征熱圖，結(jié)果發(fā)現(xiàn)其表達(dá)在國際婦產(chǎn)科聯(lián)盟（Federation International of Gynecology and Obstetrics，F(xiàn)IGO）分期（P<0.05）和生存狀態(tài)（P<0.001）中差異有統(tǒng)計(jì)學(xué)意義（圖5A）。對(duì)風(fēng)險(xiǎn)預(yù)后模型分別進(jìn)行單變量和多變量Cox 回歸模型分析，發(fā)現(xiàn) FIGO 分期、風(fēng)險(xiǎn)預(yù)后評(píng)分與OS 相關(guān)（P<0.001），可獨(dú)立預(yù)測(cè)宮頸癌患者預(yù)后（圖5B、C）。

圖5 風(fēng)險(xiǎn)預(yù)后模型在宮頸癌中的臨床應(yīng)用價(jià)值Fig.5 Clinical value of prognostic risk signature in cervical cancer

2.5 宮頸癌中存在的兩種m6A 修飾模式 一致性聚類分析顯示，Delta 面積在k>3 時(shí)增大不明顯（圖6A）。分型矩陣顯示k=2 時(shí)無樣品特別小的分組，分組內(nèi)相關(guān)性高，分組間相關(guān)性低，聚類穩(wěn)定性最佳（圖6B）。因此，基于m6A調(diào)控因子的表達(dá)確定了兩種m6A 修飾模式，命名為m6Acluster A（n=224）和m6Acluster B（n=152，圖6C）。生存分析表明，兩種m6A 修飾模式在OS 上差異無統(tǒng)計(jì)學(xué)意義（P=0.299，圖6D）。熱圖顯示25 個(gè)m6A 調(diào)控因子在m6Acluster A中的表達(dá)高于m6Acluster B（圖6E）。

圖6 25個(gè)m6A調(diào)節(jié)劑介導(dǎo)的m6A甲基化修飾模式Fig.6 m6A methylation modification patterns mediated by 25 m6A regulators

2.6 兩種m6A 修飾模式生物學(xué)功能差異和免疫浸潤(rùn)分析 GSVA 富集分析提示兩種m6A 修飾模式存在功能通路差異。m6Acluster A 與RNA 降解、細(xì)胞周期、核苷酸修復(fù)（切除修復(fù)和錯(cuò)配修復(fù)）、剪接體等KEGG 通路相關(guān)，而m6Acluster B 則與PPAR 通路、細(xì)胞色素P450代謝、花生四烯酸代謝、鈣信號(hào)通路等KEGG 通路相關(guān)（圖7A）。為研究?jī)煞Nm6A 修飾模式對(duì)宮頸癌TIME 的影響，評(píng)估了m6A 調(diào)控因子表達(dá)上調(diào)的m6Acluster A 和表達(dá)下調(diào)的m6Acluster B 免疫細(xì)胞的浸潤(rùn)水平，結(jié)果發(fā)現(xiàn)，m6Acluster B在免疫細(xì)胞浸潤(rùn)中顯著富集，包括激活的B細(xì)胞、激活的CD8+T 細(xì)胞、嗜酸性粒細(xì)胞、未成熟的B 細(xì)胞等，兩種修飾模式顯示在TIME 中存在不同免疫細(xì)胞浸潤(rùn)特征（圖7B）。

圖7 兩種m6A修飾模式生物學(xué)功能差異和免疫浸潤(rùn)分析Fig.7 Biological functional differences and immune infiltration of two m6A modification patterns

3 討論

宮頸癌已成為女性第四大常見癌癥死亡原因，嚴(yán)重威脅女性健康［8］。既往研究證實(shí)人乳頭瘤病毒（human papillomavirus，HPV）感染與宮頸癌密切相關(guān)［9］。HPV 篩查和疫苗接種是預(yù)防宮頸癌的有效策略［10］。手術(shù)及放化療是早期無轉(zhuǎn)移宮頸癌的主要治療手段。但傳統(tǒng)治療方法在晚期或復(fù)發(fā)性宮頸癌的治療中未取得令人滿意的效果，存在預(yù)后差、生活質(zhì)量低等問題［11］。m6A 甲基化是最常見的RNA 修飾，在腫瘤增殖、侵襲和轉(zhuǎn)移中起重要作用，但作用機(jī)制尚不明確［12］。因此本研究通過TCGA 和GEO 公共數(shù)據(jù)庫下載宮頸癌轉(zhuǎn)錄組表達(dá)數(shù)據(jù)和相應(yīng)臨床數(shù)據(jù)，通過生物信息學(xué)方法和組織驗(yàn)證實(shí)驗(yàn)探究m6A 甲基化修飾與宮頸癌發(fā)生發(fā)展及TIME 免疫細(xì)胞浸潤(rùn)的關(guān)系，為宮頸癌患者臨床診療和預(yù)后改善提供依據(jù)。

本研究首先分析了25 個(gè)m6A 甲基化調(diào)控因子在宮頸癌組織和正常宮頸組織中的表達(dá)和相關(guān)性，結(jié)果發(fā)現(xiàn)在宮頸癌組織中有7個(gè)顯著上調(diào)的m6A調(diào)控因子（RBM15、IGF2BP2、IGF2BP3、FMR1、YTHDF1、METTL3和YTHDF2）和5個(gè)顯著下調(diào)的m6A調(diào)控因子（ZC3H13、METTL14、YTHDC1、FTO、ALKBH3）。研究發(fā)現(xiàn)METTL14 可作為癌基因或抑癌基因在腫瘤中發(fā)揮不同作用，為m6A 修飾甲基轉(zhuǎn)移酶相關(guān)腫瘤藥物研發(fā)提供了思路，但目前未見其在宮頸癌中作用的研究［13］。IGF2BP 包括IGF2BP1、IGF2BP2 和IGF2BP3，是m6A閱讀蛋白的一種。研究發(fā)現(xiàn)IGF2BP在正常和應(yīng)激條件下均可增強(qiáng)mRNA 穩(wěn)定性和翻譯，導(dǎo)致Myc等致癌產(chǎn)物積累，在宮頸癌和肝癌細(xì)胞中敲除IGF2BP 可抑制腫瘤細(xì)胞增殖和侵襲［14］。FTO 是第一個(gè)被發(fā)現(xiàn)的m6A 去甲基化酶，通過調(diào)控m6A 修飾E2F1 和Myc 轉(zhuǎn)錄本，在宮頸癌細(xì)胞增殖和遷移中發(fā)揮重要致癌作用［15］。YTH 家族蛋白包括YTHDF1、YTHDF2、YTHDC1與YTHDC2。研究發(fā)現(xiàn)YTHDF1和YTHDF2在宮頸癌中表達(dá)增加，可促進(jìn)宮頸癌細(xì)胞增殖和糖酵解［16-17］。

近年基于m6A 調(diào)控因子的預(yù)后模型被用于評(píng)估和預(yù)測(cè)各種腫瘤患者預(yù)后，指導(dǎo)后續(xù)個(gè)體化治療方 案 制 定［18-20］。本 研 究 選 擇5 個(gè)m6A 調(diào) 控 因 子HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1，通過LASSO Cox 回歸建立了預(yù)后風(fēng)險(xiǎn)模型，利用ROC分析預(yù)測(cè)模型的特異度和敏感度，結(jié)果顯示AUC=0.734，表明預(yù)測(cè)性能良好。Kaplan-Meier 生存曲線顯示高風(fēng)險(xiǎn)組OS 明顯低于低風(fēng)險(xiǎn)組。單變量及多變量Cox回歸分析證明該模型可作為獨(dú)立的預(yù)后因子。利用HPA 數(shù)據(jù)庫和qRT-PCR 檢測(cè)5 個(gè)m6A 調(diào)控因子在臨床組織中的表達(dá)，與生物信息學(xué)預(yù)測(cè)結(jié)果一致，表明該模型可作為預(yù)測(cè)宮頸癌患者生存結(jié)局的工具。

本研究中，HNRNPC、LRPPRC 和METTL3 均為m6A 的危險(xiǎn)調(diào)控因子，能增加宮頸癌患者死亡風(fēng)險(xiǎn)。多項(xiàng)研究已證實(shí)METTL3 在宮頸癌中高表達(dá)，可促進(jìn)PDK4和HK2基因m6A 修飾，進(jìn)而增加其mRNA 穩(wěn)定性，促進(jìn)宮頸癌發(fā)生發(fā)展［16，21］。LRPPRC是一種RNA 結(jié)合蛋白，在肺腺癌、乳腺癌、胃癌等腫瘤中表達(dá)上調(diào)，抑制腫瘤細(xì)胞凋亡，促進(jìn)其侵襲和轉(zhuǎn)移［22］。HNRNPC 能夠促進(jìn)乳腺癌細(xì)胞和前列腺癌細(xì)胞增殖和腫瘤生長(zhǎng)［23-24］。ELAVL1 和FMR1 為m6A 調(diào)控因子的保護(hù)因子改善宮頸癌患者預(yù)后。ELAVL1 又稱HuR基因，與卵巢癌、尿路上皮癌、食管癌、結(jié)直腸癌等多種腫瘤侵襲和不良預(yù)后相關(guān)［25］。FMR1 通常與YTHDF2 共同參與RNA 剪接、核輸出和翻譯等，同時(shí)被認(rèn)為與脆性X綜合征、卵巢早衰和自然流產(chǎn)等有關(guān)［26-27］。

為進(jìn)一步研究宮頸癌中m6A 甲基化修飾情況，對(duì)25個(gè)m6A甲基化調(diào)節(jié)因子進(jìn)行一致性聚類，將其分為m6Acluster A 和m6Acluster B。生存分析發(fā)現(xiàn)兩組OS 差異無統(tǒng)計(jì)學(xué)意義，但25 個(gè)m6A 調(diào)控因子在m6Acluster A 中的表達(dá)明顯高于m6Acluster B。GSEA 富集分析探討m6Acluster A 和m6Acluster B兩組間生物學(xué)功能和信號(hào)通路差異，m6Acluster A與RNA 降解、細(xì)胞周期、核苷酸修復(fù)（切除修復(fù)和錯(cuò)配修復(fù)）、剪接體等KEGG通路相關(guān)，而m6Acluster B則與PPAR 通路、細(xì)胞色素P450 代謝、花生四烯酸代謝、鈣信號(hào)通路等KEGG 通路相關(guān)。ssGSEA 分析兩組間腫瘤微環(huán)境中免疫細(xì)胞浸潤(rùn)水平，發(fā)現(xiàn)m6Acluster B 在免疫細(xì)胞浸潤(rùn)中顯著富集，明顯強(qiáng)于m6Acluster A，提示m6Acluster B 存在免疫激活和免疫細(xì)胞浸潤(rùn)，可能對(duì)應(yīng)免疫炎癥型的“熱腫瘤”，對(duì)免疫治療更敏感［28］。

綜上，本研究系統(tǒng)證明了m6A 調(diào)控因子在宮頸癌中的表達(dá)改變和預(yù)后價(jià)值，構(gòu)建了宮頸癌風(fēng)險(xiǎn)預(yù)后模型并驗(yàn)證其應(yīng)用價(jià)值。同時(shí)，基于25個(gè)m6A調(diào)控因子確定了兩種不同的m6A 修飾模式，并與TIME免疫細(xì)胞浸潤(rùn)特征相關(guān)，為指導(dǎo)免疫治療策略提供了依據(jù)。