膀胱癌m6A調節因子預后模型建立與分析

白洋洋，郭依琳，陳瑞廷，潘世杰，孫繼建

1. 河南省中醫院（河南中醫藥大學第二附屬醫院）泌尿外科（鄭州 450002） 2. 鄭州大學第二附屬醫院婦科腫瘤臨床醫學研究中心（鄭州 450014）

膀胱癌（bladder cancer, BC）是常見的男性泌尿系統惡性腫瘤之一。2020年全球BC 新發病例數超過55 萬，死亡病例數超過20 萬[1]。在中國，BC 發病率位居全部惡性腫瘤的第13 位，嚴重威脅男性的生命健康[2-3]。根據是否侵犯膀胱肌層，BC 可分為非肌層浸潤性和浸潤性。非肌層浸潤性BC 患者首選經尿道膀胱腫瘤切除術聯合膀胱內治療，術后生活質量較好，但約30%患者存在局部浸潤或遠處轉移可能[4]。對于浸潤性BC，約25%患者伴有淋巴結轉移，5%出現遠處轉移，預后較差[5]。近年來，隨著免疫治療的應用，BC患者預后有所改善，但由于腫瘤異質性，并非所有BC 患者均能從中獲益[6-7]。因此，探究BC 發生發展的分子機制對提高BC 患者治療反應率意義重大。

N6-甲基腺嘌呤（N6-methyladenosine, m6A）是在RNA 腺嘌呤殘基的第6 位含氮堿基位置上添加一個甲基基團，是真核生物中最豐富的表觀轉錄組學修飾[8]。研究發現，m6A 修飾主要由m6A調節因子，即編碼器-甲基轉移酶復合物、消碼器-去甲基化酶和讀碼器-m6A 讀取蛋白調控[9]。有研究指出，m6A 調節因子具有免疫調節功能，可影響腫瘤微環境（tumor microenvironment, TME）中免疫細胞浸潤，并對免疫治療療效產生影響[10]。目前，m6A 調節因子對BC 預后影響的相關研究較少。本研究利用預后相關的m6A 調節因子構建BC 預后預測模型，評估TME 中免疫細胞浸潤情況，進一步預測免疫治療和靶向治療的療效，為BC 的臨床治療提供參考。

1 資料與方法

1.1 RNA測序轉錄組數據和臨床數據獲取

利用癌癥基因組圖譜（The Cancer Genome Atlas, TCGA）數據庫（https://www.cancer.gov/ccg/research/genome-sequencing/tcga），獲取397 例BC組織的高通量轉錄組測序數據和對應的臨床病理特征數據。397 例BC 患者的臨床病理特征見表1。

表1 397例BC患者的臨床病理特征Table 1. Clinical and pathological characteristics of 397 patients of BC

1.2 預后相關的m6A調節因子篩選

檢索中國知網、Pubmed、Embase 數據庫，中文檢索詞包括m6A 調節因子、m6A 修飾，英文檢索詞包括m6A RNA methylation、m6A regulatory factors。共篩選到26 個m6A 調節因子，包括9 個編碼器（WTAP、ZC3H13、RBM15、METTL3、RBM15B、KIAA1429、METTL5、METTL14、METTL16）、14 個讀碼器（LRPPRC、HNRNPA2B1、IGF2BP1、IGF2BP2、IGF2BP3、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、ELAVL1、YTHDC1、YTHDC2、RBMX）和3 個消碼器（ALKBH5、FTO、ALKBH3）。通過R 軟件corrplot 包計算26 個m6A 調節因子之間表達的相關性，利用survival 包對26 個m6A調節因子進行單因素Cox 回歸分析，結果通過forestplot 包輸出，篩選預后相關的m6A 調節因子。

1.3 預后預測模型構建及驗證

利用最小絕對值選擇與收縮算子（least absolute selection and shrinkage operator, LASSO）Cox 回歸，篩選出5 個預后相關m6A 調控因子（P＜0.05），構建預后預測模型。每例BC 患者的風險評分等于5 個m6A 調控因子中每個基因表達量乘以系數之和（Coefi為風險系數，Xi 為基因表達量）。計算公式如下：風險評分=分別計算397例BC患者的風險評分，根據中位風險評分分為高風險組和低風險組，比較兩組間5年生存差異。繪制受試者工作特征（receiver operating characteristic curve, ROC）曲線并計算ROC 曲線下面積（area under ROC curve, AUC），評估所構建預后預測模型的預測價值。利用單因素Cox 回歸和多因素Cox 回歸評估預后預測模型和臨床病理特征對BC 預后的影響。

1.4 基因集富集分析和免疫浸潤分析

在上述構建的BC 預后預測模型中，通過MSigDB 數據庫（http://www.broadinstitute.org/gsea/msigdb/index.jsp）獲取京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes,KEGG）相關通路基因集，利用R 軟件GSVA 包對高風險組和低風險組進行富集分析，明確兩組KEGG 相關信號通路。利用聚類分析和單樣本基因集富集分析（single sample gene set enrichment analysis, ssGSEA）[12]評估高風險組和低風險組中23 種免疫細胞類型的頻率差異，并計算免疫細胞的相對浸潤豐度。利用估計惡性腫瘤組織中基質和免疫細胞（estimation of stromal and immune cells in malignant tumors using expression data,ESTIMATE）方法評估高風險組和低風險組中每個樣本中免疫、基質成分及兩者之和在TME 中的比例，即Immune 分數、Stromal 分數和Estinate 分數。

1.5 基因相關性分析

利用R 軟件limma 包，比較高風險組和低風險組在免疫檢查點相關基因PD-L1和CTLA- 4[13]、靶向治療相關基因HER-2[14]、BRAF[15]、KRAS[16]和EGFR[17]的表達差異，繪制箱式圖，以P＜0.05為差異具有統計學意義。

2 結果

2.1 預后相關m6A調節因子篩選

本研究納入397 例BC 組織樣本。利用R 軟件分析26 個m6A 調節因子在BC 中表達水平之間的相關性，結果顯示除IGF2BP2 與YTHDC1、IGF2BP2 與METTL3，以及ALKBH3 與FMR1 外，大部分m6A 調節因子之間的表達存在顯著正相關性，其中HNRNPC 與HNRNPA2B1 的正相關性最強（r=0.66），見圖1-A、圖1-B。為進一步篩選與BC 預后相關的m6A 調節因子，參考Feng 等[18]的研究，通過單變量Cox 回歸發現YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 是BC 獨立的預后因素，其中IGF2BP3、LRPPRC、FTO 和ALKBH3 為危險調節因子，其高表達與較差的預后有關，而YTHDC1 為保護調節因子，其高表達與較好的預后有關，見圖2。

圖1 BC中m6A調節因子表達相關性Figure 1. Correlation of m6A regulator expression in BC注：A. 26個m6A調節因子表達的相關性熱圖；B. 26個m6A調節因子表達的相關性圈圖。

圖2 BC中m6A調節因子單因素Cox回歸分析Figure 2. Univariate Cox regression analysis of m6A regulators in BC

2.2 構建BC預后預測模型

利用上述篩選的YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3，通過LASSO Cox 回歸方法建立BC 預后預測模型，見圖3-A。在模型中，YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 的系數分別為-0.08、0.12、0.03、0.19和0.03。根據模型計算每例BC 患者的風險評分，將患者分為高風險組和低風險組。低風險組的OS顯著高于高風險組（P＜0.001），見圖3-B。利用ROC 曲線分析預后模型預測BC 患者5年生存率的準確性，結果顯示AUC 值為0.665。

圖3 基于m6A調節因子構建BC風險預后模型Figure 3. Prognostic risk model of BC based on m6A regulators注：A. LASSO Cox回歸模型構建預后預測模型；B. 高風險組和低風險組生存分析。

2.3 預后預測模型的臨床價值

繪制YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 在高風險組和低風險組中的表達量及與臨床特征的熱圖，結果顯示其表達量在臨床分期中的差異有統計學意義，見圖4。單因素Cox 回歸和多因素Cox 回歸分析結果顯示，年齡、臨床分期和風險評分與OS 相關（P＜0.001），均為BC 獨立的預后影響因素，見表2。

圖4 基于5個m6A調節因子構建的預后模型與臨床特征的關系熱圖Figure 4. Heatmap of the relationship between prognostic model based on five m6A regulators and clinical characteristics注：***P＜0.001，**P＜0.01。

表2 BC預后預測模型與臨床特征的關系Table 2. Relationship between BC prognostic prediction model and clinical characteristics

2.4 基因集富集分析

基因集富集分析結果顯示，高風險組在趨化因子、NOD 樣受體、嘌呤代謝、丙酮酸代謝等信號通路富集，而低風險組在藥物代謝細胞色素P450、亞油酸代謝、α-亞麻酸代謝等信號通路富集，見圖5。

2.5 免疫浸潤分析

采用ssGSEA 評估高風險組和低風險組中23種免疫細胞類型的頻率差異，結果發現高風險組免疫細胞浸潤程度顯著高于低風險組，特別是激活的CD4+T 細胞、激活的CD8+T 細胞、激活的B細胞、嗜酸性粒細胞、激活的樹突細胞、Th1 細胞、Th17 細胞、中性粒細胞等，提示高風險組和低風險組具有不同的免疫細胞浸潤特征，見圖6-A。ESTIMATE 計算高風險組和低風險組中每個樣本中免疫、基質成分及兩者之和在TME 中的比例，發現高風險組Immune 分數、Stromal 分數和Estinate 分數均高于低風險組，見圖6-B。上述結果表明高風險組的TME 為免疫炎癥型，具有豐富的免疫細胞浸潤。

圖6 BC預后預測模型的免疫浸潤分析Figure 6. Immune infiltration analysis of prognostic risk model in BC注：A. ssGSEA分析；B. ESTIMATE分析；***P＜0.001。

2.6 基因相關性分析

進一步分析免疫檢查點相關基因PD-L1和CTLA-4、靶向治療相關基因HER-2、BRAF、KRAS和EGFR在高風險組和低風險組之間的表達差異，結果發現PD-L1、CTLA-4、EGFR和KRAS在高風險組中表達更高，HER-2和BRAF在低風險組表達更高，差異均有統計學意義（P＜0.05），見圖7。以上結果提示高風險組可能對PD-L1 抑制劑、CTLA-4 抑制劑、EGFR 抑制劑和KRAS 抑制劑更敏感，而低風險組可能對HER2 抗體和BRAF 抑制劑更敏感。

圖7 高低風險組間免疫檢查點相關基因和靶向治療相關基因表達差異Figure 7. Differences in the expression of immune checkpoint related genes and targeted therapy related genes between high and low risk group注：A. PD-L1基因；B. CTLA-4基因；C. KRAS基因；D. EGFR基因；E. HER-2基因；F. BRAF基因。

3 討論

BC 是一種來源于膀胱黏膜的惡性腫瘤。手術和放化療是其主要的治療方式，然而即使接受標準治療，仍有50%～70%患者出現復發或轉移，晚期和轉移性BC 的治療效果不理想[5]。目前BC仍缺乏特異性預后標志物，因此，探究新型預后預測標志物并建立預后預測模型對治療BC 具有重要的臨床價值。

m6A 調節因子的異常表達可能影響腫瘤的發生發展[9]。隨著m6A 修飾研究領域的不斷突破，越來越多的新m6A 調節因子被發現。本研究分析了26 個m6A 調控因子表達的相關性，發現BC中m6A 調節因子之間的表達存在顯著相關性。進一步篩選與BC 預后相關的m6A 調節因子，單變量Cox 回歸結果顯示YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 是BC 獨立的預后因素。IGF2BP3 和IGF2BP1 同屬IGF2BP 家族，是m6A閱讀蛋白的一種。有研究證實IGF2BP3 在BC 組織內的高表達，可通過調節HMGB1 信號通路促進BC 的發展，進而影響BC 患者預后[19]。周子健等通過臨床標本驗證了BC 組織內IGF2BP1 的高表達可提示BC 患者預后不良，轉染IGF2BP1 siRNA 后BC 細胞的增殖和生長被顯著抑制[20]。脂肪量和肥胖相關蛋白FTO 在多種癌癥組織中表達，刺激腫瘤細胞代謝，誘導腫瘤細胞浸潤和進展。蘇甲林通過檢測24 例BC 組織及癌旁正常組織FTO 的表達，發現與對應癌旁組織相比，BC組織中FTO 相對低表達，且與病理分期分級呈負相關[21]。

近年來，有研究利用m6A 調節因子建立頭頸鱗癌[22]、肝癌[23]、結直腸癌[24]的預后預測模型，并驗證發現其具有較好的預測價值，可以指導后續臨床治療。Wu 等分析了23 個m6A 調節因子與臨床病理特征和預后的相關性，提示m6A 調節因子可以作為BC 預后評估的潛在分子靶點[25]。Chen 等分析了13 個m6A 調節因子在BC 中的表達情況，并選擇FTO、YTHDC1 和WTAP 構建預后預測模型[26]。Zheng 等利用FTO 和YTHDF2評估BC 預后預測模型的應用價值[27]。本研究納入了更多m6A 調控因子，更全面地探討了m6A修飾對BC 預后預測的價值。通過Cox 回歸模型選擇5 個m6A 調節因子（YTHDC1、IGF2BP3、FTO、ALKBH3和LRPPRC）建立BC預后預測模型，ROC 曲線檢測風險評分預測BC 患者5年生存率的準確性，結果顯示AUC 值為0.665。本研究與Zheng 等[27]建立的模型（5年AUC=0.614）相比，5年AUC 值更高，進一步分析發現Zheng 等[27]利用單因素Cox 篩選預后相關m6A 調節因子時設置篩選標準為P＜0.15，最后納入預后模型構建的YTHDF2（P=0.112）對BC 患者預后的影響較小，本研究將篩選標準設置為P＜0.05，以確保參與模型構建的每個m6A 調節因子都與BC 預后相關。另外，按照風險評分將BC 患者分為高風險組和低風險組，經過單因素Cox 回歸及多因素Cox 回歸發現年齡、臨床分期和風險評分與OS 相關，是BC 獨立的預后影響因素。

m6A 調節因子具有免疫調節功能，可影響TME 中免疫細胞浸潤，并對免疫治療療效產生影響[10]。利用ssGSEA 和ESTIMATE 法比較高風險組和低風險組中免疫細胞浸潤水平的差異，結果發現，兩組具有不同的免疫細胞浸潤特征，高風險組具有豐富的免疫細胞浸潤，其TME 為免疫炎癥型，即熱腫瘤對免疫治療療效更好[28]。

靶向治療和免疫治療的快速發展極大程度提高了BC 患者的OS，并改善了其預后。但由于腫瘤異質性，并非所有的BC 患者均能從靶向治療和免疫治療中獲益。為了探討建立的預后預測模型對靶向治療和免疫治療的指導作用，本研究選擇了免疫檢查點相關基因PD-L1和CTLA- 4、靶向治療相關基因HER-2、BRAF、KRAS和EGFR，比較高風險組和低風險組中這些基因表達的差異，結果發現PD-L1、CTLA-4、EGFR和KRAS在高風險組中表達更高，HER-2和BRAF在低風險組中表達更高，提示高風險組可能對PD-L1 抑制劑、CTLA-4 抑制劑、EGFR 抑制劑和KRAS 抑制劑更敏感，而低風險組可能對HER2抗體和BRAF 抑制劑更敏感，間接預測了高低風險組患者對免疫治療和靶向治療的療效差異。本研究存在一定局限性，僅分析了TCGA 數據庫，未進行臨床試驗驗證。

綜上所述，本研究利用5 個m6A 調節因子在BC 中構建預后預測模型并驗證其臨床應用價值，評估預后預測模型中免疫細胞浸潤程度，該模型可能成為評估BC 患者是否適合免疫治療和靶向治療的參考指標，為臨床醫師的治療選擇提供依據。