N6-甲基腺苷相關lncRNAs是預測腎癌患者預后和免疫浸潤的潛在生物標志物①

肖勝英 閆志廣 曾福仁 邱 俊 盧義晨 朱小東

（湖南省人民醫院（湖南師范大學附屬第一醫院腫瘤科），長沙 410016）

據美國癌癥協會推測，2021 年將出現超過76 080 例腎癌（renal cell carcinoma，RCC）新病例，約13 780 例將死于該?。?］。隨著診斷技術、手術和分子治療的進步，腎癌患者的預后得到顯著改善。然而，晚期腎癌患者的預后仍然很差，五年生存率約為11%［2］。據報道，分子標志物的發現和應用可為癌癥患者提供預后價值［3］。

N6-甲基腺苷（N6-methyladenosine，m6A）修飾是在真核生物的mRNA、miRNA、長鏈非編碼RNA（long non-coding RNA，lncRNA）和其他RNA 分子上發現的最常見和最豐富的修飾之一。這些修飾是可逆的，可影響RNA 分子的轉錄、加工、翻譯和代謝［4］。m6A 修飾受m6A 調節器調控，包括“寫入器”（甲基轉移酶）、“讀取器”（識別蛋白）和“擦除器”（去甲基化酶）［5］。m6A 調節多種細胞的關鍵生物過程，包括干細胞更新、細胞分化和對DNA 損傷的反應［6］。m6A 調節因子的異常表達與腫瘤發生、免疫調節異常和惡性腫瘤進展有關［6-7］。有報道稱，甲基轉移酶樣因子14（METTL14）的表達水平與腎癌的病理和臨床分期呈負相關，但與總生存期（overall survival，OS）呈正相關［8］。因此，m6A 調節的lncRNAs可能對癌癥的增殖和轉移至關重要［4］。研究證明，m6A 調節因子的失調與腎癌的發生有關［9］。然而，m6A 調節因子對lncRNA 的作用機制尚不明確。因此，了解腎癌中m6A 相關lncRNAs 的具體機制有助于識別腎癌的治療靶點和與其預后相關的生物標志物。

本研究從癌癥基因組圖譜（TCGA）數據庫中提取了14 087 個lncRNAs 表達譜，通過生物信息學和統計學分析鑒定了m6A 相關的lncRNAs。結果表明，27個m6A相關lncRNAs在腎癌中具有預后價值。此外，本研究構建了一個新的預后模型，即9 個m6A相關的lncRNAs的預后簽名（m6A-LPS），該模型用于預測腎癌患者的OS。此外，本研究分析了不同亞組間腫瘤微環境（tumor microenvironment，TME）和免疫相關性的差異，以闡明免疫細胞浸潤及其與臨床預后的相關性。

1 資料與方法

1.1 資料 于2021 年4 月從TCGA 數據門戶（https：//portal.gdc.cancer.gov/）下載腎癌患者的RNA-seq 轉錄組數據和相應的臨床數據。共得到539 個腎癌樣本和72 個鄰近正常組織數據，納入標準：①經病理證實的腎癌；②同時獲得的關于mRNA表達譜和OS 的信息。最后納入530 例腎癌患者和相應臨床病理信息（年齡、生存時間、生存狀態、性別、分級和TNM 分期）進一步分析。使用caret 包將530 例腎癌患者隨機分配到訓練集（266 例患者）或驗證集（264 例患者）隊列。表1 為兩個隊列的臨床病理特征基線資料。根據已發表的文獻確定了23個m6A 相關基因，包括“寫入器”［METTL3、METTL14、METTL16、WTAP、VIRMA（KIA1499）、RBM15、RBM15B 和ZC3H13］、“擦除器”（FTO 和ALKBH5）和“讀取器”（YTHDC1、YTHDC2、FMR1、IGF2BP1、IGF2BP2、IGF2BP3、LRPPRC、YTHDF1、YTHDF2、YTHDF3、HNRNPC、HNRNPA2B1 和RBMX），并從TCGA數據集中提取相應表達矩陣。

表1 TCGA訓練和驗證集中患者的臨床病理特征［例（%）］Tab.1 Clinicopathological features of patients in TCGA training and validation cohorts［n（%）］

1.2 方法

1.2.1 lncRNA 注釋 Genome Reference Consortium Human Build 38（GRCh38）的lncRNA 注釋文件來自GENCODE 網站（https：//www.gencodegenes.org/hu‐man/），用于注釋TCGA數據集中的lncRNAs。

1.2.2 生物信息學分析 首先，以Pearson相關分析挖掘數據集中的m6A 相關lncRNAs（|Pearsonr|>0.3和P<0.001）。與已發表的23 個m6A 相關基因中的一個或多個lncRNA 的表達值相關（|Pearsonr|>0.3和P<0.001）即為m6A 相關的lncRNA。進行單變量Cox 回歸分析以過濾數據集中的與m6A 預后相關lncRNA?！癵lmnet”R 包被用于最小絕對收縮和選擇算子（LASSO）Cox 回歸（懲罰參數通過10 倍交叉驗證估計）［10］，共發現9 個m6A 相關lncRNAs 可預測腎癌患者預后。以LASSO 回歸算法獲得回歸系數和風險評分方程：風險評分=系數×m6A相關lncRNA 表達水平總和。在該方程的基礎上，分別計算訓練集和驗證集中每個患者的風險評分。并以中位風險評分為分界點將患者分為高危組和低危組。使用“ConsensusClusterPlus”包闡明腎癌中m6A 相關lnc-RNAs 的生物學特性。然后，將腎癌患者分為不同亞型。使用JAVA 程序對MSigDB 的Hallmark 基因集“h.all.v6.2.symbols.gmt”進行基因集富集分析（GSEA），探索腎癌不同亞型間存活率差異的原因。兩種亞型間的富集途徑通過誤檢率<0.05 和歸一化富集評分（NES）確定。采用R 包“estimate package”計算每個患者的免疫評分［11］。通過CIBERSORT 網站（https：//cibersort.stanford.edu/）得到每個樣本中22 種免疫細胞類型比例，采用1 000 次排列算法，將CIBERSORTP<0.05 的樣本進行亞組間的免疫浸潤水平差異比較，亞組按照聚類亞型和風險評分分組。

1.3 統計學處理 R4.0.2 和SPSS24.0 軟件用于統計學檢驗。采用Mann-WhitneyU檢驗比較癌組織與正常組織中m6A 相關lncRNAs 的表達水平。兩個及兩個以上的亞組采用學生t檢驗和單因素方差分析比較。應用卡方檢驗比較訓練集和驗證集中的分類變量。兩組間差異比較采用log-rank 檢驗。Pearson 相關性檢驗用于分析亞型、風險評分、臨床病理特征、程序性死亡配體1（PD-L1）和免疫浸潤水平間的相關性。Cox 回歸模型用于單變量和多變量分析，以確定風險評分與其他臨床特征相結合的獨立預后價值。受試者工作特征曲線（ROC）用于評估m6A 相關lncRNAs 的預后簽名在三年和五年OS 的預測效率。P<0.05表示差異具有統計學意義。

2 結果

2.1 鑒定腎癌患者中m6A 相關的lncRNAs 使用“GENCODE”網站上下載的文件，基于ensemble ID，在TCGA 數據集中共識別14 087 個lncRNAs 用于后續分析。從TCGA數據集中提取了23個m6A相關基因的表達矩陣。通過Pearson相關分析尋找m6A相關的lncRNAs，得到239個與m6A顯著相關的lncRNAs，隨后將其與預后信息結合，進行單變量Cox 回歸分析（P<0.05）。最后，在TCGA 數據集中篩出27 個與腎癌患者OS顯著相關的m6A相關lncRNAs（圖1A），其表達水平在腎癌和正常組織中差異很大（圖1B）。AL008718.3、AC009948.2、AF117829.1、LINC01409、AL133243.3、AC084018.1、AC005253.1、AC090589.3、RUSC-1AS1、AC012615.6、RAD51-AS1、AL162586.1、AL928654.2、AC006435.2、AL136295.7、LINC00115、AC245052.4、AC114730.3、ARHGAP27P1-BPTFB1-KPNA2P3、AL135999.1、PTOV1-AS2、LINC00342 和AC00414在腎癌組織中的表達水平均高于正常鄰近組織（P<0.05）。AC018752.1、COL18A1-AS1、RPL34-AS1 和SNHG10 在腎癌組織中的表達水平明顯低于正常組織（P<0.05）。

圖1 m6A 相關lncRNAs 的表達及其在腎癌患者預后中的作用Fig.1 Expressions of m6A-related lncRNAs and their role in prognosis of RCC patients

2.2 TCGA 數據集中m6A-LPS 的構建和驗證 進一步分析530 例腎癌患者和相應的臨床病理信息。為構建m6A-LPS 以預測腎癌患者的OS，基于TCGA訓練集中的27 個預后相關的m6A 相關lncRNAs進行LASSO Cox 分析并生成了包含9 個m6A 相關lncRNAs 的m6A-LPS，即AC018752.1、COL18A1-AS1、RPL34-AS1、AL008718.3、AC009948.2、AL133243.3、AL162586.1、AL136295.7和LINC00342。使用LASSO算法得到的系數計算TCGA 訓練集和驗證集的風險評分，公式如下：風險評分=1.659 2×AL008718.3 表達水平?0.460 4×AC018752.1 表達水平?1.044 3×COL18A1-AS1表達水平+0.637 1×AC009948.2表達水平?1.8806×RPL34-AS1 表達水平+0.461 1×AL133243.3 表達水平+0.193 2×AL162586.1 表達水平+0.118×LINC00342 表達水平?0.085 8×AL136295.7 表達水平。根據中位風險評分將患者分為高風險組和低風險組。圖2A、B 顯示了TCGA訓練集和驗證集中9 個m6A-LPS 的風險評分、熱圖、OS、OS 狀態和表達譜分布。在高危人群中，與m6A相關的高風險基因AC009948.2、LINC00342、AL008718.3、AL162586.1、AL133243.3和AL136295.7高表達，而保護性m6A相關基因如AC018752.1、COL18A1-AS1 和RPL34-AS1 在低風險人群中表達上調。Kaplan-Meier生存曲線顯示，在TCGA 訓練集和驗證集中，具有高風險評分的腎癌患者臨床結果較差（OS率低，OS時間短，P<0.000 1，圖2C、D）。結果表明，m6A-LPS 對TCGA 訓練集中腎癌患者的OS有較好的預測能力（三年ROC AUC=0.737，五年AUC=0.779，圖3A、B）。在TCGA驗證集中，9 個風險簽名的三年和五年AUC 值分別為0.703 和0.724（圖3C、D）。以上結果表明m6A-LPS 對腎癌患者具有穩健的OS預測能力。

圖2 m6A 相關lncRNAs 風險評分及其對腎癌患者預后的影響Fig.2 m6A-related lncRNAs risk score and their influence on prognosis of RCC patients

圖3 ROC曲線Fig.3 ROC curve

2.3 m6A-LPS的分層分析 為更好地評估m6A-LPS的預后能力，進行分層分析以確認m6A-LPS 是否保留其預測每個亞組OS的能力。與低風險患者相比，高風險腎癌患者在G1~2 和G3~4、Ⅰ~Ⅱ期和Ⅲ~Ⅳ期亞組的OS較差（圖4A~D）。還證明m6A-LPS具有預測年齡≤65 歲或>65 歲患者（圖4E、F）及不同性別患者（圖4G、H）OS 的能力。這些數據表明m6A-LPS可作為腎癌患者的潛在預測因子。

2.4 m6A-LPS 是腎癌患者的獨立預后因素 應用單變量和多變量Cox 分析評估m6A-LPS 是否是腎癌患者的獨立預后因素。根據TCGA 數據集中的腎癌患者數據，單變量Cox 分析顯示m6A-LPS 與OS 顯著相關［風險比（HR）：2.261，95%置信區間（CI）：1.884~2.713，P<0.001，圖4I］，而多變量Cox 分析進一步表明m6A-LPS 是OS 的獨立預測因子（HR：2.058，95%CI：1.655~2.560，P<0.001，圖4J）。因此，m6A-LPS 是驗證集中腎癌患者OS 的獨立預測因子（單變量分析HR：1.255，95%CI：1.149~1.371，P<0.001；多變量分析HR：1.206，95%CI：1.082~1.344，P<0.001，圖4K、L）。以上結果表明m6A-LPS是一個獨立的預后指標，可能有助于腎癌的臨床預后評估。

圖4 模型驗證的生存曲線和多變量和單變量獨立預后分析Fig.4 Survival curve for model validation and multivariate and univariate analysies of independent prog?nostic analysis

2.5 m6A 相關預后lncRNAs 的共識聚類與腎癌患者的特征和生存期顯著相關 根據m6A 相關預后lncRNAs 的表達水平與模糊聚類檢驗比例的相似性，k=2 在k=2~9 間具有最好的聚類穩定性。根據m6A相關預后lncRNAs的表達水平，將530例腎癌患者分為兩個亞型，即聚類1（n=313）和聚類2（n=217），見圖5A。27個m6A相關預后lncRNAs在聚類1中的表達水平大多低于聚類2（圖5B）。聚類1的OS比聚類2 長（P<0.001，圖5C）。分析27 個m6A 相關預后lncRNAs 與PD-L1 間的相關性，結果表明，AL008718.3、AL133243.3、SNHG10、LINC00115、AF117829.1、AC084018.1、AC005253.1 和PD-L1 的表達水平呈正相關（P<0.05，圖5D）。

圖5 TCGA 隊列中Cluster1/2 亞型腎癌的不同臨床病理特征和生存率以及m6A 預后lncRNAs 與靶基因的表達和關系Fig.5 Differential clinicopathological features and survival of RCC in Cluster1/2 subtypes in TCGA cohort and expression and relationship of m6A prognostic lncRNAs and target genes

2.6 m6A 相關預后lncRNAs 的共識聚類及其22 種免疫細胞浸潤水平 評估聚類1 與下調的m6A 相關預后lncRNAs 表達和聚類2 與上調的m6A 相關預后lncRNA 表達間的22 種免疫細胞浸潤水平差異（圖6A）。聚類1 顯示幼稚B 細胞和巨噬細胞M0 浸潤水平更高（圖6A、B），而聚類2 中CD8 T 細胞和濾泡輔助T 細胞浸潤水平更高（圖6A、C、D）。隨后進行GSEA 以闡明導致兩個亞組間免疫細胞浸潤差異的潛在調節機制。結果表明，腫瘤的惡性特征，包括哺乳動物雷帕霉素靶標（mTOR）信號（NES=1.827，標準化P=0.007）、血管內皮生長因子（VEGF）信號（NES=1.672，標準化P=0.019）和Notch 信號（NES=1.739，標準化P=0.037）在聚類2中活躍（圖6E）。因此，兩個亞組間免疫細胞浸潤差異可能與mTOR、VEGF 和Notch 信號通路的活躍有關。

圖6 免疫細胞浸潤和TME分析Fig.6 Analysis of immune cell infiltration and TME

2.7 與腎癌分級、分期和免疫評分相關的預后風險評分 評估臨床特征和風險評分間的關系。圖7A為TCGA數據集中高風險和低風險組中9個m6A-LPS的表達水平。其中AC018752.1、COL18A1-AS1 和RPL34-AS1在高危組中的表達水平通常低于低危組。AC009948.2、LINC00342、AL008718.3、AL162586.1、AL133243.3和AL136295.7在低危組中表達水平較低。M 分期（P<0.001）、T 分期（P<0.001）、TNM 分期（P<0.001）、分級（P<0.001）、聚類亞型（P<0.001）和免疫評分（P<0.001）在高風險組和低風險組間差異有統計學意義。進一步探索風險評分、亞型、分類和免疫評分間的關系。風險評分隨著組織學分級和分期的增加而增加（P<0.05，圖7B、C）。聚類1的風險評分明顯低于聚類2（P<0.001，圖7D）。與低免疫評分組相比，高免疫評分組的風險評分較高（P<0.001，圖7E）。提示腎癌患者的風險評分與分級、分期、聚類亞型和免疫評分顯著相關。

圖7 TCGA訓練集中預后風險評分與臨床病理特征和免疫評分相關Fig.7 Prognostic risk scores correlated with clinicopatho?logical features and immunoscore in TCGA training cohort

2.8 22 種免疫細胞的浸潤豐度與風險評分的關系 分析22 種免疫細胞的浸潤水平與風險評分間的關系，以評估9 個m6A-LPS 對腎癌免疫微環境的影響。風險評分與漿細胞、CD4 記憶T 細胞、CD8 T細胞、調節性T 細胞和濾泡輔助性T 細胞的浸潤水平呈正相關（P<0.001，圖8A~E）。風險評分和單核細胞浸潤水平（P<0.001）、靜息肥大細胞（P<0.001）、CD4 記憶靜息T 細胞和巨噬細胞M2 呈負相關（P<0.001，圖8F~I）。以上結果表明m6A-LPS 與腎癌TME有關。

圖8 22種免疫細胞類型的浸潤豐度與風險評分間的關系Fig.8 Relationship between risk score and infiltration abundances of 22 immune cells

3 討論

眾所周知，病理分期是臨床上預測腎癌預后的關鍵依據［12］。然而，相同分期的患者預后不同，表明目前的分期系統在反映腎癌的異質性方面不甚準確［13］。因此，應探索潛在的預測和治療生物標志物。本研究探討了m6A 相關lncRNAs 對來自TCGA數據集的530 例患者的預后意義。在TCGA 數據集中共有27 個m6A 相關lncRNAs 表現出預后價值，其中9個用于建立m6A-LPS預測腎癌患者的OS。根據中位風險評分，腎癌患者被分為低風險組和高風險組。高風險組表現出較差的臨床結果。高風險組和低風險組還與不同的聚類亞型、免疫評分、等級和TNM 分期相關。多變量Cox 回歸分析證實m6ALPS 為OS 的獨立危險因素。此外，本研究中的模型可獨立于影響患者預后的其他臨床預后因素。因此，該模型可應用于不同的臨床亞組。

m6A 修飾mRNA 和lncRNAs 使其具有多種功能，調控mRNA 的剪接、穩定、核輸出和翻譯，是真核生物中最豐富的修飾。此外，與m6A 相關的lncRNAs 研究也引起眾多癌癥領域專家的關注［14］。ALKBH5 通過降低lncRNA 核旁斑點組裝轉錄本1的甲基化加速胃癌的轉移和侵襲［15］。DKK1 的lncRNA 激活調節因子通過與熱休克蛋白家族A（Hsp70）成員1A 和Y-box 結合蛋白1 形成三元復合物穩定m6A 甲基化，從而促進頭頸部鱗狀細胞癌進展［16］。lncRNAs 的m6A 修飾可介導基因轉錄抑制、改變結構、調節lncRNAs 的穩定性影響腫瘤的發生發展［17-19］。lncRNAs 和m6A 是腎癌發生的重要調節因子［9］。目前已初步了解lncRNA 和m6A 在腎癌中的功能，然而，lncRNA 和m6A 在腎癌進展中的病理作用機制尚不明確，其預后生物標志物的研究甚少。本研究基于m6A 相關的lncRNA 構建了腎癌的獨立預后模型。

本研究展示了腎癌中m6A 相關lncRNAs 的表達情況、預后價值和差異的免疫浸潤豐度。篩選了9 個m6A 相關lncRNAs 構建模型以預測腎癌患者的OS。其中RPL34-AS1 是結直腸癌和胃癌的抑癌基因［20］。RPL34-AS1 在食管癌中通過下調RPL34 表達、在乳頭狀甲狀腺癌中通過競爭性結合miR-3663-3p/RGS4、在宮頸癌中通過MDM2-p53 途徑抑制細胞的增殖、遷移和侵襲［21-23］。本研究結果也提示RPL34-AS1 是腎癌的抑癌基因。COL18A1-AS1在腎透明細胞癌中是抑癌基因，并可能作為癌癥診斷和治療的預后生物標志物［24-25］。膽管癌患者中COL18A1-AS1 高表達可提高存活率［26］。這些報道為本研究結果提供支持，即COL18A1-AS1可能為腎癌的抑癌基因，有望成為其診斷和治療的預后生物標志物。AL162586.1是一種關鍵的免疫相關lncRNA，其表達與N 期呈正相關，可能通過剪接途徑影響前列腺癌患者預后［27］。本研究結果也認為AL162586.1可能為腎癌患者的致癌基因。此外，其他lncRNAs是首次被發現。這些lncRNAs可能是預測腎癌患者預后的潛在生物標志物。因此，本研究結果可能有助于識別與m6A 相關的lncRNAs，為深入了解其在腎癌中發生、發展的潛在作用提供見解。

進一步分析樣本中不同免疫細胞的浸潤水平以確定免疫細胞浸潤和TME 在腎癌中的作用。與聚類2 相比，聚類1 中的患者具有更長的OS。B 細胞、幼稚細胞和巨噬細胞M0的浸潤水平在聚類1中更高，而聚類2中CD8 T細胞和濾泡輔助T細胞的浸潤水平更高。GSEA 結果表明，腫瘤的惡性功能特征包括mTOR、VEGF 和Notch 信號通路，在聚類2 中顯著富集。ZHOU 等［28］觀察到METTL14 表達下調可上調AKT 和mTOR 表達水平，與腎癌患者的不良預后相關。越來越多的證據表明，腫瘤中VEGF 水平升高可能導致適應性和先天免疫反應的抑制，血清或腫瘤的升高與腎癌患者預后不良有關［29］。mTOR 和VEGF 在致癌過程中發揮關鍵作用，包括細胞存活、增殖和預后［30］。Notch 信號通路影響TME 并調節癌癥生物學過程，包括癌癥的轉移和進展［31］。Notch3 調節細胞周期進程和缺氧誘導因子2α，后者參與腎癌的增殖［32］。這些研究支持本文結果，即具有更短OS 的聚類2 顯著富集mTOR、VEGF和Notch 信號通路，表明mTOR、VEGF 和Notch 信號通路的mRNA可作為m6A甲基化修飾的靶標。因此推測m6A 和mTOR、VEGF 和Notch 信號通路的甲基化共同參與了腎癌兩個亞組間TME 的差異調節。隨后進行風險評分與免疫細胞的相關分析，以評估其間關系。CD4 記憶T 細胞、CD8 T 細胞、Treg 和濾泡輔助T細胞與風險評分呈正相關。TME中浸潤的活化CD4 記憶T 細胞、CD8 T 細胞、Treg 和濾泡輔助性T細胞的減少可能不利于腎癌患者的預后。本推論與何天基等［33］的結論大致相同，即大量CD4+T 細胞和瘤內CD8+T 細胞與短生存期和高腫瘤分級相關。TME 中Tregs 的誘導可促進腫瘤生長，且TM 中Tregs 的分化與癌癥亞型患者總體生存率低有關。Tregs可抑制抗原呈遞細胞的成熟、細胞因子的分泌以及細胞毒性顆粒酶和穿孔素的產生，還可促進血管生成、腫瘤生長、增殖和腫瘤向轉移性疾病的轉變［34-35］。以上研究結論同樣支持本研究結果，因此，通過敲除或阻斷Tregs、活化的CD4 記憶T 細胞和CD8 T細胞的相關信號通路或基因可能改善TME的抗腫瘤活性和免疫抑制。

總之，本研究建立了一個m6A 相關的lncRNA模型，為預測腎癌患者預后提供了生物標志物和治療靶點，并為未來對lncRNA 的m6A修飾機制的研究提供見解。然而，本研究也存在一些缺陷和局限性，m6A 相關lncRNA 的生物學機制尚未完全闡明。因此，課題組將在后續工作中收集臨床數據并擴大樣本量，通過外部試驗驗證模型的準確性。