長鏈非編碼RNA預測頭頸部鱗狀細胞癌患者預后的研究

向 琳，徐曉晨，譚君武，杜 波*

(1.湖北民族大學附屬民大醫(yī)院,湖北 恩施445000；2.吉林大學第一醫(yī)院,吉林 長春130021)

頭頸部鱗狀細胞癌(head and neck squamous cell carcinoma ，HNSC)每年全球約新增75萬病例，死亡約40萬，嚴重威脅人類健康[1]，這類腫瘤包括起源于口腔、口咽、下咽、喉、鼻咽、腭舌和扁桃體的癌癥。這類患者的預后受多種因素的影響，預測其預后對臨床工作十分重要。長鏈非編碼RNA(long non-coding RNA，lncRNA)在HNSC的發(fā)生發(fā)展過程中發(fā)揮重要作用，影響腫瘤的生存、遷移和侵襲，可能是潛在的預后標志物[2]。既往已有l(wèi)ncRNA預測喉癌預后的相關研究[3]，本研究利用癌癥基因組圖譜(The Cancer Genome Atlas，TCGA)數(shù)據(jù)庫中HNSC表達數(shù)據(jù)預測患者5年生存率，篩選預測預后的關鍵lncRNA。

1 材料與方法

1.1 HNSC患者信息

頭頸部鱗狀細胞癌(HNSC)患者RNA-seq數(shù)據(jù)和臨床資料由TCGA網(wǎng)站(https://tcga-data.nci.nih.gov/tcga/)下載。臨床資料包括總體生存率(overall survival，OS)、年齡、性別、腫瘤分級和AJCC分期。

1.2 HNSC樣本表達譜

HNSC患者HTSeq-counts數(shù)據(jù)與GENCODE數(shù)據(jù)庫資料(https://www.gencodegenes.org/human/，gencode.v22)匹配添加注釋信息。使用RPKM對lncRNAs和mRNAs表達水平進行標準化。lncRNAs篩選標準：Ⅰ.轉錄位置不在蛋白編碼區(qū)域；Ⅱ.在Gencode中存在注釋信息；Ⅲ.至少在一半HNSC樣本中表達；Ⅳ.平均RPKM>0.1。

1.3 差異表達分析

使用edgeR包以log2|fold change|>1 和調整后P<0.001作為閾值計算差異表達lncRNAs。

1.4 關鍵lncRNA篩選

在訓練集中使用單因素COX分析計算差異表達lncRNAs與總體生存率(OS)之間的關系，然后使用多因素COX分析隨機生存森林法篩選關鍵lncRNAs，風險評分(risk scores，RS)等于每個lncRNA的Coeffcient系數(shù)乘以其表達量之和，大于中位值定義為高風險，反之為低風險，以P-value <0.05表示差異有統(tǒng)計學意義。

1.5 生存分析

使用Kaplan-Meier曲線log-rank檢驗計算兩組的生存差異，多因素COX分析和分層分析評估關鍵lncRNAs和臨床特征的關系。使用受試者工作特征(receiver operating characteristic，ROC)曲線評估五年生存預測的敏感性和特異性。

1.6 功能富集

使用Spearman相關性分析篩選lncRNA-mRNA共表達的蛋白編碼基因，相關系數(shù)>0.40，P<0.01表示顯著相關。功能富集使用GO分析，通路富集使用KEGG分析,以P-value <0.001,Q-value<0.01作為閾值。所有分析均使用R(version 3.6.2)軟件。

2 結果

2.1 關鍵lncRNAs

通過與臨床信息匹配共篩選出475例樣本，隨機分為訓練集(n=238)和測試集(n=237)，共篩選出493個滿足閾值的差異表達lncRNAs。在訓練集共篩選出5個關鍵lncRNAs與樣本預后明顯相關(表1)，其中1個(RP11-865I6.2)為正系數(shù)，表示其高表達與短生存期相關，4個(RP11-417L19.2、RP11-567M16.1、RP11-44K6.2、FALEC)為負系數(shù)，表示其高表達與長生存期相關。

表1 HNSC訓練集中5個lncRNAs與總體生存率的關系

2.2 訓練集中關鍵lncRNAs與總體生存率

訓練集中Kaplan-Meier曲線顯示高風險患者(n=119)預后顯著差于低風險患者(n=119)(P<0.0001，圖1A)。高風險患者3年生存率38.61%、5年生存率23.59%、8年生存率0%，低風險患者3年生存率79.80%、5年生存率63.78%、8年生存率39.25%。使用ROC曲線評估5個lncRNAs預測HNSC患者預后的作用，其ROC曲線下(area under curv，AUC)面積為0.774(圖1B)。單因素COX分析顯示5個lncRNAs風險評分顯著和患者預后相關(表2)。高風險患者死亡率顯著高于低風險患者(圖1C)，熱圖顯示RP11-865I6.2在高風險組上調，RP11-417L19.2、RP11-567M16.1、RP11-44K6.2和FALEC在高風險組下調(圖1C)。

2.3 測試集、整集中關鍵lncRNAs與總體生存率

測試集中高風險患者(n=132)OS明顯差于低風險患者(n=105)(P=0.0042，圖2A)。整集中高風險患者(n=251)生存期明顯短于低風險患者(n=224)(P<0.0001，圖2B)。測試集高風險患者3年生存率51.02%、5年生存率37.93%、8年生存率15.17%；低風險患者3年生存率61.19%、5年生存率52.35%、8年生存率52.35%。整集中高風險患者3年生存率45.11%、5年生存率31.61%、8年生存率11.08%；低風險患者3年生存率72.09%、5年生存率58.66%、8年生存率47.16%。測試集和整集AUC分別為0.651(圖2C)和0.712(圖2D)。

表2 不同數(shù)據(jù)集COX分析結果

圖1 HNSC訓練集中5個lncRNAs風險評分模型預測總體生存率

2.4 關鍵lncRNAs的獨立性分析

COX分析顯示5個lncRNAs危險評分、年齡與預后相關(表2)。按中位年齡61歲將樣本分為兩組，結果顯示在不同年齡分層中，高風險患者OS顯著短于低風險患者(圖3A-3C)。低風險中不同年齡組患者預后無差異(P=0.091，圖3D)。

圖2 HNSC測試集與整集中5個lncRNAs風險評分模型預測總體生存率

圖3 5個lncRNAs在不同年齡患者的分層分析

2.5 關鍵lncRNAs的功能分析

共篩選出728個蛋白編碼基因，GO分析顯示5個關鍵lncRNAs相關蛋白編碼基因在242個GO項中顯著富集，其中BP 219項，CC 13項，MF 10項，KEGG分析顯示在33個通路中顯著富集。功能富集主要集中在白細胞粘附、淋巴細胞分化、淋巴和T細胞激活(圖4A)，通路富集主要集中在抗原處理和呈遞、T細胞分化、細胞粘附分子(圖4B)。

圖4 5個lncRNAs的功能富集分析

3 討論

長鏈非編碼RNA(lncRNA)調節(jié)腫瘤發(fā)生、遷移和侵襲。已有眾多研究發(fā)現(xiàn)lncRNA通過下游基因靶向調節(jié)HNSC細胞增殖與入侵[4-6]，許多研究表明lncRNA可以作為癌癥預后的獨立預測因子[7,8]。Xiong D 等[9]的研究發(fā)現(xiàn) LINC00958 和 HOXC13-AS可作為HNSC患者的診斷標志物。lncRNA是癌和其他頭頸部腫瘤的重要生物標志物[10]。但是，lncRNA在預測HNSC患者預后中的作用尚不明確。

本研究中，我們利用TCGA數(shù)據(jù)庫，在訓練集中通過單因素和多因素COX分析構建了一個基于5-lncRNAs的風險模型，利用中位風險值將患者分為高風險組和低風險組。在測試集和全集中驗證了5-lncRNAs模型，與低危患者相比，高?；颊呖傮w生存率低。將5-lncRNAs模型與臨床信息相結合做分層分析，證實5-lncRNAs模型獨立于臨床因素，這增加了預后預測的可靠性。RP11-865I6.2上調，與患者預后呈負相關，RP11-417L19.2、RP11-567M16.1、RP11-44K6.2和FALEC下調，與患者預后呈正相關。我們的研究中，共篩選出5個lncRNAs作為預測HNSC患者預后的關鍵因子，據(jù)我們所知，之前沒有關于這5個lncRNAs的報告，表明它們是在本研究中新發(fā)現(xiàn)的。Li J等[11]研究發(fā)現(xiàn)了10個lncRNAs，Xing L等[12]發(fā)現(xiàn)4個lncRNAs，Zhang Z等[13]發(fā)現(xiàn)3個lncRNAs，Yang B等[14]發(fā)現(xiàn)8個lncRNAs，這些lncRNAs均與HNSC患者預后明顯相關。這可能是分析過程中數(shù)據(jù)處理方法的差異，樣本量的不同，樣本來源數(shù)據(jù)庫的不同以及不同lncRNA之間的差異導致患者處于不同的風險水平所導致。本研究利用一個穩(wěn)健的隨機生存森林分析來篩選與預后相關的關鍵lncRNA，確定了一個5-lncRNAs風險模型，而且沒有和以上研究相重復的lncRNA。我們希望更多的類似研究揭示預測HNSC患者預后的關鍵因子，Pan Y等[15]研究發(fā)現(xiàn)RP11-865I6.2、RP11-366H4.1、HOTTIP、RP11-275N1.1可作為HNSC患者的預后預測基因，這些重復較多的基因可能為進一步研究HNSC的分子機制和生物標志物提供新的靶點和理論基礎，而且，采用生物信息學工具對HNSC中l(wèi)ncRNA進行表達分析時必須使用不同的方法對結果進行驗證。我們的研究提示這個5個lncRNAs獨立于臨床因素。事實上，包含基因信息、臨床病理分期的變量相較于單純的TNM分期變量能更加準確的估計喉癌總體生存率[3]。既往有研究發(fā)現(xiàn)在中國人中l(wèi)ncRNA WWTR1-AS1過表達與不良預后相關[16]，我們的結果中不包括上述基因，可能是由于人種的差異導致其潛在的分子機制不同。富集分析顯示與5個lncRNAs相關蛋白編碼基因主要集中在細胞免疫及分子粘附方面。這與以前腫瘤免疫以及細胞粘附分子的相關研究結果相似[17,18]。

本研究存在一定的局限和不足。首先，本研究主要集中在數(shù)據(jù)挖掘和分析兩個方面，這些都是基于統(tǒng)計學方法，研究結果沒有通過進一步的實驗驗證，沒有試驗分析其潛在機制，此外，我們只分析驗證了TCGA數(shù)據(jù)集中5個lncRNAs的預測能力，沒有其他數(shù)據(jù)庫lncRNA表達數(shù)據(jù)用于進一步驗證。其次，由于不同的lncRNA檢測方法可能導致不同的結果，因此必須對lncRNA的檢測、量化和轉錄活性的測定過程進行標準化。最后，lncRNA通過復雜的調控網(wǎng)絡調節(jié)腫瘤過程，涉及到不同種類的順式和反式調控元件，在廣泛的生物過程中發(fā)揮著重要的調控作用，需要進一步的綜合分析5個lncRNA在HNSC中的作用。因此，我們建議未來的研究應增加多變量預測模型以提高HNSC患者預后預測的準確性。

本研究中，我們證明了lncRNAs在HNSC患者中的預測價值，提示5個lncRNAs有助于預測臨床結果，并且是獨立預測HNSC患者生存率的有效預后生物標志物。