三陰性乳腺癌m6A相關lncRNA的免疫預后模型構建*

鮑淑梅，管浩欽，蘇瑩，劉沛，呂小毅?

（1.新疆大學軟件學院，新疆 烏魯木齊 830091；2.新疆大學信息科學與工程學院，新疆 烏魯木齊 830017）

0 引言

根據最新全球癌癥統計報告顯示乳腺癌（BRCA）已超過肺癌，成為世界上確診人數最多的癌癥，占全球報告病例的11.7%.2020年，全球僅女性乳腺癌就有230萬例，同時有超過68萬人死于乳腺癌[1].三陰性乳腺癌（TNBC）是一種特殊的乳腺癌分子亞型，約占乳腺癌的15%，其特征是雌激素受體（ER）、孕激素受體（PR）和人表皮生長因子受體2（HER2）的缺失[2].TNBC的復發和死亡主要發生在確診后的前5年，尤其是確診后的1～3年內.確診后的前5年，TNBC患者內臟器官轉移的發生率遠高于非TNBC患者且預后較差，5年生存率不到30%[3-4].由于TNBC的異質性和缺乏分子靶點，TNBC患者對內分泌治療或HER2靶向治療不敏感，導致化療仍然是TNBC患者的標準治療方法[5-6].因為腫瘤對化療的耐藥性，導致TNBC患者在接受化療后腫瘤迅速復發或轉移，所以TNBC仍然是一種預后較差和治療方法有限的疾病[7].N6-甲基腺苷（m6A）是真核細胞最豐富的核糖核酸（RNA）修飾，在各種生物過程和信使核糖核酸（mRNA）代謝中至關重要，例如RNA加工、轉運和穩定性[8-9].m6A修飾包括甲基化轉移酶、信號轉導和去甲基化酶，它們分別被稱為寫入者、讀取者和擦除者.此外，m6A修飾也是一個可逆的RNA表觀遺傳過程[10].RNA結構的變化可以影響各種細胞過程.因此，m6A調控的長鏈非編碼RNA（lncRNAs）的作用可能對癌細胞的增殖和遷移至關重要[11].

有研究表明，m6A修飾可調節腫瘤發生和腫瘤發展.例如，甲基轉移酶3（METTL3）的表達減少通過m6A甲基化介導的III型膠原蛋白α1鏈（COL3A1）上調促進三陰性乳腺癌的轉移[12].此外，YTH N6-甲基腺苷RNA結合蛋白F3（YTHDF3）通過以m6A依賴性方式穩定鋅指E-box綁定同源盒1（ZEB1）促進三陰性乳腺癌進展和轉移[13].最近，一些生物信息學研究表明，m6A調節器的失調與TNBC有關[14-15].m6A調節因子在lncRNA中的具體作用尚不明確；因此，了解m6A相關lncRNA在TNBC發展中的機制可能對預后目標有用.

本文首先從癌癥基因組圖譜（TCGA）數據集中提取了10 357個lncRNA和23個m6A基因的表達譜.接下來使用皮爾遜相關性分析確定了m6A相關的lncRNA.構建的模型是基于m6A的新型預后模型，該模型旨在預測TNBC患者的總生存期（OS）.使用公開的藥物敏感性數據庫，發現了針對這種m6A相關lncRNA特征的候選藥物，并探討了與免疫治療反應的關系.最后，建立了列線圖來預測TNBC患者的OS.

1 材料與方法

1.1 TNBC患者信息獲取

TCGA數據庫（https://cancergenomenihgov/）是目前最大的公開資源數據庫，也是國際上最常用的數據庫.該數據庫提供數以千計的與腫瘤樣本相關的基因表達、體細胞突變、基因甲基化等數據集.世界各地的癌癥研究人員可以在TCGA數據庫中免費獲取相關數據，識別可能在癌癥發展中發揮作用的重要基因.本文從TCGA數據庫中獲取了TNBC患者的RNA序列轉錄組數據、相關臨床信息和突變數據.篩選三陰性乳腺癌患者樣本的標準為：（1）組織病理學診斷為乳腺癌；（2）按TNBC診斷標準進行篩查免疫組織化學中ER、PR、HER2表達缺失患者.

1.2 m6A基因和m6A相關lncRNA的選擇

本文從TCGA數據庫中獲得了lncRNA和m6A基因的表達情況.根據前人研究，從TCGA中檢索到23個m6A基因的表達矩陣，包括甲基化轉移酶（METTL3、METTL14、METTL16、VIRMA、RBM15、RBM15B、ZC3H13和WTAP）、結合蛋白（IGFBP1、IGFBP2、IGFBP3、LRPPRC、FMR1、YTHDC1、YTHDC2、YTHDF1、YTHDF 2、YTHDF3、HNRNPA2B1、HNRNPC和RBMX）和去甲基化酶（ALKBH5和FTO）[16].通過皮爾遜相關性分析篩選了m6A相關的lncRNA，確定1 149個m6A相關的lncRNA.該過程使用了|Pearson R|的標準大于0.4且p＜0.001.

1.3 風險簽名的建立和驗證

整個TCGA數據集被隨機分為訓練集和測試集.利用訓練集構建m6A相關的lncRNA模型，并應用整個集和測試集來驗證所建立的模型.結合TCGA中的TNBC生存信息，本文采用單變量Cox回歸從TCGA數據集中的275個m6A相關lncRNAs中篩選出m6A預后相關的lncRNAs（p＜0.05）.使用R包glmnet進行LASSO回歸（使用通過10倍交叉驗證估計的懲罰參數）篩選候選基因，同時使用R包randomForestSRC構建RSF算法模型篩選候選基因.對LASSO與RSF算法篩選出的候選基因取交集，獲得5個m6A相關的lncRNAs與來自TCGA數據集的TNBC患者的OS明顯相關.應用多因素Cox回歸分析5個m6A相關的lncRNAs，最終建立了基于3個m6A相關的lncRNA風險模型.風險評分計算公式如下：Risk score=coef (lncRNA1)×expr (lncRNA1)+coef (lncRNA2)×expr (lncRNA2)+···+coef (lncRNAn)×expr (lncRNAn).其中：coef表示系數，coef (lncRNAn)是與生存相關的lncRNA的系數，expr (lncRNAn)是lncRNA的表達量.根據中位風險評分，建立了包括低風險組和高風險組在內的亞組[17].

1.4 差異基因富集分析

進行了基因本體論（GO）富集分析以確定高低風險組之間差異表達的基因所富集的相關生物學通路.此過程使用R包clusterProfiler.分析閾值由p值確定，p＜0.05表示相關生物功能顯著豐富[18].

1.5 免疫治療模式探索

使用R包maftools評估和匯總突變數據.根據腫瘤特異性突變基因測量腫瘤突變負荷（TMB）[19].使用腫瘤免疫功能障礙和排斥（TIDE）算法計算TIDE評分，通過TIDE評分預測免疫治療反應的效果[20].

1.6 PCA和KM生存分析

使用PCA根據3個m6A基因的表達模式，對全基因表達譜、23個m6A基因、275個m6A相關的lncRNA和風險模型相關的lncRNAs的高維數據進行有效降維、模型識別和分組可視化[21].使用KM生存分析評估高風險組和低風險組之間OS的多樣性.使用R包survMiner和survival完成實驗.

1.7 臨床治療中靶向m6A相關lncRNA模型的潛在化合物探索

為了在臨床上獲得用于TNBC治療的潛在化合物，本文在TNBC數據集中計算了從GDSC網站獲得的化合物的IC50.R包pRRhetic用于預測從GDSC網站獲得的化合物在TNBC患者中的IC50.

1.8 m6A相關lncRNA模型的獨立性

考慮到TNBC患者的其它臨床特征，如年齡、原發腫瘤的范圍和大?。═）分期、淋巴結播散情況（N）分期等，進行多變量和單變量Cox回歸分析以檢驗預后模型是否為獨立變量.

1.9 建立和證明預測列線圖

建立列線圖使用預測因子（年齡、風險評分、T分期、N分期）預測2年、3年和5年TNBC患者的OS.應用基于Hosmer-Lemeshow檢驗的校正曲線來說明實際結果和模型預測結果之間的一致性.

2 實驗結果

2.1 TNBC患者m6A相關lncRNA的鑒定

風險模型構建和分析的詳細工作流程如圖1所示.從TCGA數據庫中提取了23個m6A基因和10 357個lncRNA的矩陣表達.本文將m6A相關lncRNA定義為與大于或等于23個m6A基因之一顯著相關的lncRNA（|Pearson R|＞0.4和p＜0.001）.最后，使用圖2（a）中的?；鶊D可視化m6A-lncRNA共表達網絡，將275個m6A相關lncRNA識別為m6A相關的lncRNA.TCGA所有m6A基因與m6A相關lncRNA的相關性如圖2（b）所示.

圖1 工作流程

圖2 TNBC患者中m6A相關lncRNA的鑒定

2.2 根據TNBC患者的m6A相關lncRNA構建和驗證風險模型

本文使用單變量Cox回歸分析從TCGA訓練集中的275個m6A相關lncRNA中篩選出m6A相關預后lncRNA.TCGA數據集中的11個m6A相關lncRNA與OS顯著相關（圖3(a)）.

圖3 基于m6A相關lncRNA的TNBC患者風險模型

LASSO分析是多元回歸分析的常用方法.該方法的應用不僅提高了統計模型的預測精度和可解釋性，而且可以同時進行變量選擇和正則化.該方法廣泛應用于相關性較差、預測值突出的高維數據中特征的最優選擇，避免過擬合.因此，LASSO可以有效地辨別最可用的預測標記并產生預后指標來預測臨床結果.RSF是一種基于決策樹集合的生存分析機器學習方法，已經證明其對嘈雜的變量數據具有較好的魯棒性.此外生存樹既能檢測非線性關系的變化（變量和預測目標值之間），又能檢測變量之間相互作用.所以適用于從包含生存狀態的基因表達數據中篩選出關鍵基因[22].之后，對LASSO和RSF選出的候選基因取交集，獲得5個m6A相關的lncRNA用于后續的多變量Cox分析（圖3(b～d)）.接下來使用多變量Cox回歸分析來篩選基因構建模型.最終3個m6A相關lncRNA用于構建風險模型以評估TNBC患者的預后風險.

根據預后風險等級的中值，將TNBC樣本分為低風險組和高風險組.二者之間的風險等級分布如圖4（a）所示，兩個不同風險組患者的生存狀態和生存時間如圖4（b）所示.每個患者的3個m6A相關lnc-RNA的相對表達標準如圖4（c）所示.生存分析表明，低風險組的OS長于高風險組（p＜0.01）（圖4(d)）.

圖4 TCGA訓練集中3個m6A相關lncRNA風險模型的預后價值

為了測試這個已建立模型的預后能力，本文使用統一公式計算了測試組和整個組中每個患者的風險評分.圖5描述了風險等級的分布、生存狀態和生存時間的模式，以及測試集（圖5(a～c)）和整個集（圖5(e～g)）中m6A相關lncRNA的表達.對測試集和整個集進行的KM生存分析顯示與TCGA訓練集的結果沒有差異：具有較高風險評分的TNBC患者的OS比具有較低風險評分的患者更差（p＜0.05）（圖5(d)和圖5(h)），表明m6A相關lncRNA模型水平影響TNBC患者的預后.

圖5 TCGA測試和整個組中3個m6A相關lncRNA風險模型的預后價值

2.3 PCA進一步驗證了m6A相關lncRNA模型的分組能力

基于全基因表達譜、23個m6A基因、23個m6A相關lncRNA以及按3個m6A相關lncRNA表達譜分類的風險模型，進行PCA檢驗低危組和高危組的差異（圖6）.圖6（a～c）顯示高風險和低風險人群的分布相對分散.然而，根據本文模型獲得的結果表明，低風險和高風險群體具有不同的分布（圖6(d)）.這些結果表明預后特征可以區分低風險和高風險人群.

圖6 高風險組和低風險組之間的主成分分析

2.4 使用m6A相關lncRNA模型估計腫瘤免疫微環境和癌癥免疫治療反應

基于來自108個TNBC樣本的m6A相關lncRNA模型，進一步分析了TNBC中高低風險組之間幾種免疫細胞、通路或功能的活性.低危組和高危組在免疫指標的表達上表現出顯著差異（圖7(a～b)）.為探索基于m6A模型的潛在分子機制，本文進行了GO富集分析，揭示了許多免疫相關生物過程的參與（圖7(c)）.接下來研究了m6A相關lncRNA模型與免疫治療生物標志物之間的相關性.不出所料，發現高風險組比低風險組更有可能對免疫療法產生反應，表明這種基于m6A的分類器指數可以作為預測TIDE的指標（圖7(d)）.使用R包maftools分析和總結突變數據.根據變異效應預測因子對突變進行分層.圖7（e）和圖7（f）顯示了高風險和低風險亞組之間突變頻率最高的前20個驅動基因.

圖7 在TCGA整個數據集中使用m6A相關lncRNA模型估計腫瘤免疫微環境和癌癥免疫治療反應

2.5 鑒定靶向m6A相關lncRNA模型的新型候選化合物

為了確定針對本文用于治療TNBC患者的lnc-RNA模型的潛在藥物，使用pRRhetic算法根據每個樣本的GDSC數據庫中可用的半數最大抑制濃度（IC50）來估計治療反應.篩選出了3種化合物在高低風險組之間估計的IC50存在顯著差異，高風險組對這些化合物都更敏感.圖8顯示了可用于TNBC患者進一步分析的3種化合物.

圖8 三種新型候選化合物在高低風險組之間的IC50

2.6 m6A相關lncRNA預后風險模型評估

本文進行了單變量和多變量Cox回歸分析，以此評估3個m6A相關lncRNA的風險模型是否具有TNBC的獨立預后特征.在單變量Cox回歸分析中，風險評分的風險比（HR）和95%置信區間（CI）分別為1.050和1.011～1.090（p＜0.011）.在多變量Cox回歸分析中，HR為1.092，95%CI為1.002～1.189（p＜0.044）（圖9），表明3種m6A相關lncRNA的風險模型與臨床病理參數無關，如種族、年齡和腫瘤/淋巴結/轉移（TNM）分期等.

圖9 TCGA整個組TNBC中m6A相關lncRNA和臨床特征的預后風險模型評估

2.7 預后列線圖的構建和評估

構建了包含風險等級和臨床風險特征的列線圖來預測2年、3年和5年的OS發生率（圖10(a)）.相關實驗結果顯示，2年、3年和5年OS率與預測率表現出理想的一致性（圖10(b～d)）.對于2年、3年和5年OS率，受試者工作特征曲線（ROC）的AUC值分別為0.984、0.805和0.853（圖10(e)）.

圖10 預后列線圖的構建和評估

3 分析與討論

與其它乳腺癌相比，TNBC是一種更具侵襲性的亞型，與其它亞型相比存活率較低.TNBC患者中只有30%～45%達到了與其它乳腺癌亞型相似的完全病理反應和生存率.在預后方面，由于缺乏有效的治療靶點，TNBC患者的預后效果一般不佳.因此，越來越多的研究專注于識別非編碼RNA的特征，以預測TNBC患者的生存和免疫治療反應.

本文從TCGA數據集中鑒定了275個m6A相關的lncRNA，以探索m6A相關的lncRNA的預后功能.TCGA數據集證實了11個m6A相關lncRNA的預后價值，其中3個用于構建m6A相關lncRNA模型以預測TNBC患者的OS.其中RP11-126K1.6與卵巢癌中基于鉑的化學耐藥性相關[23].此外，其它lncRNA首次被揭示與癌癥相關.隨后，根據中間風險評分將TNBC患者分為高風險組和低風險組，高風險組的臨床結果明顯較差.多變量Cox回歸分析表明，m6A相關lncRNA模型是OS的獨立風險因素.本文還建立了列線圖，顯示2年、3年和5年OS的觀察率和預測率之間的完美一致性.ROC分析表明，該模型在TNBC的2年、3年和5年OS預測率顯示出較好的敏感性和特異性.最后，觀察到的2年、3年和5年OS預測率顯示出極好的一致性.基于與OS獨立相關的3個m6A相關lncRNA的風險模型相當準確，該預測模型可以為后續研究識別新的生物標志物.

TMB是體細胞編碼突變的總數，與觸發抗腫瘤免疫的新抗原的出現有關[24].最近的研究表明，TMB是預測程序性細胞死亡配體1（PD-L1）治療反應的有效生物標志物[25].本文發現TMB在高低風險組之間具有差異.

此外，越來越多的研究使用了TIDE預測評分，這是一種為免疫治療預測開發的計算框架，其預測功能已成功驗證[26].本文中，TIDE算法的預測表明，高風險亞型患者對免疫療法有更好的反應.根據以上結果，推斷該預測模型可能為腫瘤治療提供可靠的免疫生物標志物.此外，本文為TNBC中m6A相關lncRNA的分子生物學機制提供了新的見解.

在臨床上，病理分期是影響TNBC預后的決定性因素[27].然而，同一分期的TNBC患者往往有不同的臨床結局，提示目前的分期系統在提供可靠的預測和反映TNBC的異質性方面是不準確的[28].因此，應探索潛在的預測和治療生物標志物.本文建立的m6A相關lncRNA模型為TNBC患者的預后預測提供了一種新方法.該結果也為未來研究lncRNAs m6A修飾的過程和機制提供了思路.

4 結論

本文為TNBC患者的預后預測提供了線索，并可能有助于闡明m6A調節的lncRNA的過程和機制.此外，該預后模型在2年、3年和5年的OS預測中表現出良好的性能.這些發現將為今后TNBC的預后和免疫治療提供一定的借鑒.