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

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

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

1 材料與方法

1.1 HNSC患者信息

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

1.2 HNSC樣本表達(dá)譜

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

1.3 差異表達(dá)分析

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

1.4 關(guān)鍵lncRNA篩選

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

1.5 生存分析

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

1.6 功能富集

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

2 結(jié)果

2.1 關(guān)鍵lncRNAs

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

表1 HNSC訓(xùn)練集中5個lncRNAs與總體生存率的關(guān)系

2.2 訓(xùn)練集中關(guān)鍵lncRNAs與總體生存率

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

2.3 測試集、整集中關(guān)鍵lncRNAs與總體生存率

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

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

圖1 HNSC訓(xùn)練集中5個lncRNAs風(fēng)險評分模型預(yù)測總體生存率

2.4 關(guān)鍵lncRNAs的獨立性分析

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

圖2 HNSC測試集與整集中5個lncRNAs風(fēng)險評分模型預(yù)測總體生存率

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

2.5 關(guān)鍵lncRNAs的功能分析

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

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

3 討論

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

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

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

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