基于免疫相關(guān)IncRNA建立晚期頭頸部鱗癌預(yù)后模型

吳嘉雯，唐加山

(南京郵電大學(xué)理學(xué)院，江蘇 南京 210023)

頭頸部癌癥起源于上呼吸道和消化道的各種解 剖部位，居世界惡性腫瘤第六位。2020年全球范圍 內(nèi)頭頸部癌癥新增病例約93萬(wàn)人，死亡病例約47 萬(wàn)人[1]；其中，約90%的病例被歸類為HNSCC。盡管近年來(lái)放療、手術(shù)和輔助化療等治療策略取得了 進(jìn)展，但是大約有三分之二的HNSCC患者在診斷 時(shí)已經(jīng)處于局部晚期(III期和IV期)，預(yù)后較差[2]。 因此，建立可以準(zhǔn)確預(yù)測(cè)晚期HNSCC患者預(yù)后的模型，對(duì)指導(dǎo)臨床診斷和治療有重要意義。

IncRNA被定義為大于等于200個(gè)核昔酸的RNA,目前還沒(méi)有證據(jù)表明它們可以被翻譯成肽。隨著研究的深入,IncRNA的異常已經(jīng)被證明具有抑制或促進(jìn)腫瘤的作用，在腫瘤發(fā)展中發(fā)揮著不可或缺的作用[3-4]。此外,最近的研究表明，IncRNA在抗原呈遞、免疫激活和免疫細(xì)胞浸潤(rùn)等癌癥免疫的不同階段具有重要作用[5]，因此免疫相關(guān)IncRNA引起了相當(dāng)大的關(guān)注。免疫相關(guān)IncRNA作為新興的癌癥生物標(biāo)志物已被用于多種癌癥的診斷和生存預(yù)測(cè),例如肝癌[6]、肺癌[7]、胃癌[8]等。同樣，當(dāng)前已有相關(guān)研究表明了免疫相關(guān)IncRNA在HNSCC中具有重要的預(yù)后價(jià)值[9-10]。然而，目前缺少關(guān)于晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA的研究。本文旨在識(shí)別與晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA,并基于這些IncRNA建立預(yù)后模型以改善晚期HNSCC的預(yù)后預(yù)測(cè)。

1 材料與方法

1.1 材料HNSCC患者的RNA測(cè)序(RNA Sequencing,RNA-seq)數(shù)據(jù)來(lái)自癌癥基因組圖譜(The Cancer Genome Adas,TCGA；https://portal,gdc.canc-er.gov/)中的HNSCC項(xiàng)目組，并從中提取了IncRNA表達(dá)數(shù)據(jù);HNSCC患者的臨床病理信息來(lái)自UCSC Xena(https://xenabrowser.net/)。將HNSCC患者的IncRNA表達(dá)數(shù)據(jù)和臨床數(shù)據(jù)進(jìn)行整理和合并，使得整理后的樣本均有對(duì)應(yīng)的表達(dá)數(shù)據(jù)和臨床數(shù)據(jù)，并刪除生存信息缺失和生存時(shí)間小于30天的樣本。最后根據(jù)臨床信息提取出臨床III期和IV期的HNSCC樣本。

1.2 方法

1.2.1 免疫相關(guān)IncRNA的提取 本文選取edgeR方法來(lái)對(duì)IncRNA進(jìn)行差異表達(dá)分析,R軟件中的edgeR是基于負(fù)二項(xiàng)分布的統(tǒng)計(jì)方法[11]，根據(jù)樣本的IncRNA表達(dá)量，選擇FDR<0.05和|log2FC|≥1作為閾值篩選出在腫瘤樣本和非腫瘤樣本之間有明顯差異表達(dá)的IncRNA。另外從免疫學(xué)數(shù)據(jù)庫(kù)(Immunology Database and Analysis Portal,ImmPort)(https：//www.immport.org/home)網(wǎng)站獲得免疫相關(guān)基因列表，進(jìn)行免疫基因與IncRNA的Pearson相關(guān)性分析，以相關(guān)系數(shù)|R|>0.4且P<0.001為篩選條件得到HNSCC免疫相關(guān)IncRNA。

1.2.2 關(guān)鍵預(yù)后免疫相關(guān)IncRNA的篩選 首先在訓(xùn)練集中對(duì)上述確定的免疫相關(guān)IncRNA進(jìn)行單因素Cox回歸分析，以P值小于0.05為標(biāo)準(zhǔn)確定與晚期HNSCC患者總生存期(overall survival,OS)相關(guān)的免疫IncRNA。然后基于這些預(yù)后相關(guān)的免疫IncRNA建立Lasso-Cox回歸模型,并進(jìn)行10折交叉驗(yàn)證確定最優(yōu)模型,進(jìn)一步篩選出與晚期HNSCC患者預(yù)后密切相關(guān)的免疫IncRNA。同樣，在訓(xùn)練集中對(duì)確定的免疫相關(guān)IncRNA進(jìn)行Coxboost分析來(lái)選擇與晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA。基于R語(yǔ)言中的“Coxboost”包構(gòu)建Coxboost模型，并使用cv.CoxBoost函數(shù)進(jìn)行5折交叉驗(yàn)證來(lái)選擇最優(yōu)提升步數(shù)。最后比較篩選得到的兩組與晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA,確定用于建立預(yù)后模型的關(guān)鍵免疫IncRNA。

1.2.3 預(yù)后模型的構(gòu)建 基于篩選出來(lái)的關(guān)鍵預(yù)后免疫相關(guān)IncRNA,利用隨機(jī)生存森林算法建立晚期HNSCC患者的預(yù)后模型,并計(jì)算每個(gè)患者的預(yù)后風(fēng)險(xiǎn)值。繪制模型的5年OS的時(shí)間依賴性受試者工作特征曲線(time depesndent receiver operating characteristic curve,timeROC)曲線，計(jì)算使Youden指數(shù)最大時(shí)的風(fēng)險(xiǎn)值作為閾值，并根據(jù)風(fēng)險(xiǎn)值的閾值將訓(xùn)練集和測(cè)試集中的患者分為高風(fēng)險(xiǎn)組(風(fēng)險(xiǎn)值大于閾值)和低風(fēng)險(xiǎn)組(風(fēng)險(xiǎn)值小于等于閾值)。

1.3 統(tǒng)計(jì)學(xué)分析本文所有的統(tǒng)計(jì)分析和可視化均基于R軟件(4.1.2版)進(jìn)行,P值小于0.05被認(rèn)為具有統(tǒng)計(jì)學(xué)意義,P值小于0.01認(rèn)為具有顯著差異。

2 結(jié)果

2.1 數(shù)據(jù)處理從UCSC Xena網(wǎng)站下載了612名HNSCC患者的臨床信息，并將臨床信息和RNA-seq數(shù)據(jù)進(jìn)行匹配，得到539個(gè)匹配樣本;其中包括495個(gè)腫瘤樣本和44個(gè)正常組織樣本。接著去除生存信息缺失和生存時(shí)間小于30 d的樣本，并從中提取III期和IV期樣本，最終得到362個(gè)晚期HNSCC樣本。

根據(jù)基因注釋文件，從HNSCC的RNA-seq數(shù)據(jù)中得到15 878個(gè)IncRNA的表達(dá)數(shù)據(jù)，通過(guò)差異表達(dá)分析(FDR<0.05且|log2FC|≥1)得到1 729個(gè)差異表達(dá)的IncRNA。根據(jù)ImmPort數(shù)據(jù)庫(kù)下載的免疫相關(guān)基因信息，從基因表達(dá)數(shù)據(jù)中提取了1 279個(gè)免疫相關(guān)基因的表達(dá)數(shù)據(jù)。采用Pearson相關(guān)分析篩選與免疫相關(guān)基因相關(guān)的IncRNA,結(jié)果確定了988個(gè)免疫相關(guān)IncRNA(|R|>0.4,P<0.001)。

2.2 關(guān)鍵預(yù)后相關(guān)免疫IncRNA的篩選及預(yù)后模型的構(gòu)建按照7∶3的比例將樣本隨機(jī)分為訓(xùn)練集和測(cè)試集。在訓(xùn)練集中，我們對(duì)包含988個(gè)免疫相關(guān)IncRNA的生存數(shù)據(jù)進(jìn)行了單因素Cox回歸分析，以確定與預(yù)后相關(guān)的免疫IncRNA；根據(jù)Wald檢驗(yàn)的P值小于0.05為標(biāo)準(zhǔn)，篩選與預(yù)后相關(guān)的免疫IncRNA,結(jié)果得到159個(gè)免疫相關(guān)IncRNA。再基于這159個(gè)免疫相關(guān)IncRNA建立Lasso-Cox回歸模型，通過(guò)10折交叉驗(yàn)證的Lasso-Cox回歸分析結(jié)果如圖1所示。圖1A中的橫坐標(biāo)為A的對(duì)數(shù)值，縱坐標(biāo)為方差值，最上方為L(zhǎng)asso-Cox回歸分析后的IncRNA的數(shù)量，選擇使模型達(dá)到最小損失時(shí)的懲罰參數(shù)λ(λ=0.118)，得到8個(gè)系數(shù)非零的預(yù)后相關(guān)免疫IncRNA,分別是LINC01305、RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1。

圖1 Lasso-Cox回歸分析篩選關(guān)鍵預(yù)后相關(guān)免疫IncRNA

同樣，在訓(xùn)練集上建立Coxboost模型，進(jìn)行10次提升之后，得到7個(gè)系數(shù)非零的預(yù)后相關(guān)免疫lncRNA,為L(zhǎng)INC01305、RP11-30P6.6、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1。

對(duì)篩選得到的兩組預(yù)后相關(guān)免疫IncRNA進(jìn)行比較,可以看到基于Cox回歸和Lasso-Cox回歸篩選出來(lái)的與預(yù)后相關(guān)的免疫IncRNA集合中僅比基于Coxboost分析篩選得到的預(yù)后相關(guān)的免疫lncRNA集合多了IncRNA RP11-890B15.2。表明了UNC01305、RP11-30P6.6、RP11-65M17.3、RP11- 497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1這7個(gè)免疫相關(guān)IncRNA與預(yù)后具有密切的聯(lián)系。為了確定RP11-890B15.2對(duì)于預(yù)后的重要性，我們基于RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3、RP11-445F12.1、LINC01305和RP11-30K9.5這8個(gè)IncRNA,利用隨機(jī)生存森林構(gòu)建預(yù)后模型，并根據(jù)VIMP值對(duì)于這8個(gè)預(yù)后相關(guān)免疫IncRNA進(jìn)行重要性排序。通過(guò)網(wǎng)格搜索計(jì)算并選擇選擇使隨機(jī)生存森林的袋外錯(cuò)誤率最小的mtry和node size組合作為該模型的最優(yōu)參數(shù)。根據(jù)計(jì)算結(jié)果可以得到，在mtry=45,node size=1時(shí)構(gòu)建的隨機(jī)生存森林模型的袋外錯(cuò)誤率達(dá)到最低，此時(shí)的袋外錯(cuò)誤率為33.3%,并且模型的袋外錯(cuò)誤率在1 000棵生存樹(shù)時(shí)已經(jīng)趨于穩(wěn)定(見(jiàn)圖2A)。計(jì)算變量VIMP值并根據(jù)VIMP值對(duì)IncRNA的重要性進(jìn)行排序，從圖2B和表1中可以看出這8個(gè)免疫相關(guān)IncRNA對(duì)于預(yù)后模型的重要性均是正值，提示了這8個(gè)免疫相關(guān)IncRNA對(duì)于患者預(yù)后都具有一定的影響，因此將這8個(gè)免疫相關(guān)IncRNA全部納入預(yù)后模型。

圖2 隨機(jī)生存森林構(gòu)建預(yù)后模型

表1 8個(gè)免疫相關(guān)IncRNA的重要性值

2.3 預(yù)后模型的驗(yàn)證在訓(xùn)練集上,根據(jù)上述基于隨機(jī)生存森林建立的預(yù)后模型計(jì)算每個(gè)患者的風(fēng)險(xiǎn)值，并根據(jù)風(fēng)險(xiǎn)值的閾值(cutoff=33.07)將訓(xùn)練集中風(fēng)險(xiǎn)值大于閾值的患者歸為高風(fēng)險(xiǎn)組，反之則歸為低風(fēng)險(xiǎn)組;訓(xùn)練集和測(cè)試集中患者風(fēng)險(xiǎn)值的分布、相應(yīng)的生存狀態(tài)和8個(gè)免疫相關(guān)IncRNA的表達(dá)水平的熱圖分別如圖3和圖4所示，從圖中可以看出高風(fēng)險(xiǎn)組發(fā)生死亡的人數(shù)多于低風(fēng)險(xiǎn)組，且高風(fēng)險(xiǎn)組患者的生存時(shí)間少于低風(fēng)險(xiǎn)組。熱圖顯示,UNC01305和RP11-30K9.5在低風(fēng)險(xiǎn)組中高表達(dá),而RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3和RP11-445F12.1在高風(fēng)險(xiǎn)組中高表達(dá)。Kaplan-Meier生存分析結(jié)果顯示，低風(fēng)險(xiǎn)組患者總生存率顯著高于高風(fēng)險(xiǎn)組患者(P<0.001)(圖5A)。在測(cè)試集中進(jìn)行了相同的分析，并且得到了類似的結(jié)果(結(jié)果見(jiàn)圖5B)。另外,在總數(shù)據(jù)集中，低風(fēng)險(xiǎn)組患者的中位生存期為2 900 d而高風(fēng)險(xiǎn)組患者的中位生存期882 d，可見(jiàn)高風(fēng)險(xiǎn)組相比低風(fēng)險(xiǎn)組的總體生存較差(見(jiàn)圖5C)。

圖3 訓(xùn)練集風(fēng)險(xiǎn)值分析

圖4 測(cè)試集風(fēng)險(xiǎn)值分析

圖5 高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組患者的Kaplan Meier生存曲線

另外，我們繪制了timeROC曲線并根據(jù)AUC值來(lái)評(píng)估模型的預(yù)后價(jià)值，計(jì)算得到訓(xùn)練集上1年、2年和3年OS的ROC曲線的AUC值分別為0.75、0.76和0.77(圖6A)；測(cè)試集上1年、2年和3年ROC曲線的AUC值分別為0.61、0.66和0.64(圖6B)；總數(shù)據(jù)集上1年、2年和3年ROC曲線的AUC值分別為0.71、0.73和0.73(圖6C)。

圖6 timeROC曲線

t-分布隨機(jī)鄰域嵌入(t-SNE)算法是一種降維技術(shù),經(jīng)常用于將高維數(shù)據(jù)映射到二維空間從而使其可視化。我們使用t-SNE探索了訓(xùn)練集(圖7A)和測(cè)試集(圖7B)中高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組患者的分布,從圖中可以直觀地看到，具有不同風(fēng)險(xiǎn)的患者被很好地分為兩個(gè)集群。

圖7 訓(xùn)練集和測(cè)試集的t分布隨機(jī)鄰域嵌入分析

2.4 預(yù)后模型在不同亞組中的驗(yàn)證為了進(jìn)一步評(píng)估所構(gòu)建的預(yù)后模型的可靠性，我們?cè)谡麄€(gè)數(shù)據(jù) 集中進(jìn)行了分層分析。根據(jù)年齡小于70歲和大于等于70歲劃分為不同的年齡組;分別在不同的年齡組內(nèi)對(duì)局低風(fēng)險(xiǎn)組患者進(jìn)行Kaplan-Meier生存分析并利用Log-rank檢驗(yàn)進(jìn)行高低風(fēng)險(xiǎn)組之間生存率的比較。根據(jù)Kaplan-Meier曲線可知,在小于70歲的分組中高風(fēng)險(xiǎn)組患者生存率顯著低于低風(fēng)險(xiǎn)組患者的生存率(P<0.0001,結(jié)果見(jiàn)圖8A)；根據(jù)Kap-lan-Meier方法估算得到的高風(fēng)險(xiǎn)組患者的中位生存期為2.5年,低風(fēng)險(xiǎn)組患者中位生存期為7.9年;可知低風(fēng)險(xiǎn)組患者的預(yù)后較高風(fēng)險(xiǎn)組患者更好。在年齡大于等于70歲的分組中,高風(fēng)險(xiǎn)組患者的中位生存期為2.3年，低風(fēng)險(xiǎn)組中位生存期為5.7年，高風(fēng)險(xiǎn)組患者生存率低于低風(fēng)險(xiǎn)組患者生存率，且具有統(tǒng)計(jì)學(xué)意義(P=0.011,見(jiàn)圖8B)。同樣,在不同性別組內(nèi)進(jìn)行Kaplan-Meier生存分析，根據(jù)分析結(jié)果可知，男性患者組(P<0.000 1,見(jiàn)圖8C)和女性患者組(P=0.034,見(jiàn)圖8D)的高、低風(fēng)險(xiǎn)組之間的生存率存在差異且具有統(tǒng)計(jì)學(xué)意義;其中在男性患者組中,高風(fēng)險(xiǎn)組患者的中位中位生存期為2.6年，低風(fēng)險(xiǎn)組患者中位生存期為7.9年;女性患者組中的高風(fēng)險(xiǎn)組患者的中位生存期為2.3年。綜上可知,該預(yù)后模型在不同的年齡亞組和性別亞組中均具有良好的預(yù)測(cè)能力，表明了該模型是可靠的。

3 討論

當(dāng)前的數(shù)據(jù)表明，近年來(lái)HNSCC的發(fā)病率和死亡率在不斷地上升，盡管診斷和治療手段的進(jìn)步使得早期患者得到了改善，但是大部分患者在確診時(shí)處于晚期[12]，耐藥性難以克服，使其預(yù)后不良。近年來(lái),lncRNA成為腫瘤研究的熱點(diǎn)，人們發(fā)現(xiàn)lncRNA在腫瘤的發(fā)生和發(fā)展中起著重要作用，lncRNA的異常表達(dá)被認(rèn)為是影響腫瘤活性的關(guān)鍵因素;其中與免疫相關(guān)的IncRNA被證明與多種惡性腫瘤的預(yù)后有關(guān)[13-14]。然而，目前缺少關(guān)于晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA的研究。因此，本文的研究目的是通過(guò)對(duì)TCGA數(shù)據(jù)庫(kù)中HNSCC晚期患者的相關(guān)數(shù)據(jù)進(jìn)行分析，確定與晚期HNSCC預(yù)后相關(guān)的免疫IncRNA,用于晚期HNSCC患者的預(yù)后風(fēng)險(xiǎn)評(píng)估。

本研究中，使用單因素Cox回歸、Lasso-Cox回歸和Coxboost分析,確定了8個(gè)與晚期HNSCC預(yù)后相關(guān)的免疫IncRNA(LINC01305、RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1),基于8個(gè)預(yù)后相關(guān)免疫IncRNA利用隨機(jī)生存森林建立了晚期HNSCC患者的生存預(yù)后風(fēng)險(xiǎn)模型,并根據(jù)風(fēng)險(xiǎn)值的閾值將晚期HNSCC患者劃分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。根據(jù)8個(gè)IncRNA的表達(dá)量熱圖可知，IncRNA LINC01305和RP11-30K9.5在低風(fēng)險(xiǎn)組中高表達(dá),提示IncRNA LINC01305和RP11-30K9.5的高表達(dá)量可能與患者較好的預(yù)后有關(guān);而RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3和RP11-445F12.1在高風(fēng)險(xiǎn)組中高表達(dá)，提示它們的高表達(dá)量可能與患者較差的預(yù)后有關(guān)。根據(jù)Kaplan-Meier曲線可以看到，高風(fēng)險(xiǎn)組患者的生存率顯著低于低風(fēng)險(xiǎn)組患者的生存率(P<0.01),其中根據(jù)基于整個(gè)數(shù)據(jù)集的Kaplan-Meier分析結(jié)果可知，高風(fēng)險(xiǎn)組患者的3年生存率約為41%,中位生存期約為2.4年；而低風(fēng)險(xiǎn)組患者的3年生存率約為85%,中位生存期約為7.9年。高風(fēng)險(xiǎn)組患者的3年生存率比低風(fēng)險(xiǎn)組患者的3年生存率高了44%,可見(jiàn)高風(fēng)險(xiǎn)組的患者預(yù)后較低風(fēng)險(xiǎn)組更差。同時(shí)，在訓(xùn)練集、測(cè)試集和總數(shù)據(jù)集中，模型3年OS的timeROC曲線的AUC值分別為0.77、0.64和0.73,表明該模型具有較好的預(yù)測(cè)能力。為了進(jìn)一步驗(yàn)證預(yù)后風(fēng)險(xiǎn)模型在不同臨床特征分層中的適用性，分別在不同性別組和不同年齡組中對(duì)高低風(fēng)險(xiǎn)組患者的生存率進(jìn)行了分析，分析結(jié)果表明,基于8個(gè)免疫相關(guān)IncRNA的預(yù)后風(fēng)險(xiǎn)模型不受性別和年齡因素的影響，在不同的分層中的高低風(fēng)險(xiǎn)組患者生存率的差異均具有顯著性。

另外通過(guò)查閱文獻(xiàn)發(fā)現(xiàn),IncRNA LINC01305可以促進(jìn)宮頸癌的進(jìn)展，且LINC01305的高表達(dá)與宮頸癌患者的低生存率相關(guān)[15]；且有研究證明了UNC01305可以促進(jìn)食管鱗狀細(xì)胞癌的轉(zhuǎn)移和增殖。迄今為止還沒(méi)有相關(guān)報(bào)告證明這8個(gè)免疫相關(guān)IncRNA與HNSCC的關(guān)系，它們?cè)贖NSCC中的作用機(jī)制還需要進(jìn)行進(jìn)一步的探索。

這項(xiàng)研究也存在一定的局限性。首先,本文僅基于TCGA-HNSCC數(shù)據(jù)進(jìn)行分析，樣本量較少可能存在分析偏差;另外,這是一項(xiàng)純粹的統(tǒng)計(jì)分析和數(shù)據(jù)挖掘研究，構(gòu)建的預(yù)后風(fēng)險(xiǎn)評(píng)估模型在臨床預(yù)測(cè)中是否可行還有待進(jìn)一步驗(yàn)證。

綜上所述，本文篩選出8個(gè)與晚期HNSCC患者預(yù)后相關(guān)的免疫IncRNA，并基于這8個(gè)免疫相關(guān)IncRNA成功構(gòu)建了預(yù)后模型,經(jīng)驗(yàn)證表明該預(yù)后模型能較好的預(yù)測(cè)晚期HNSCC患者的預(yù)后情況?？赡転槲磥?lái)的研究和臨床實(shí)踐提供新見(jiàn)解。