頭頸部鱗癌免疫相關lncRNA的生物信息學分析

吳嘉雯，唐加山

(南京郵電大學 理學院，南京 210023)

0 前 言

頭頸部癌癥位居世界惡性腫瘤第六位，2020年全球頭頸部癌癥新增病例約93萬人，死亡病例約47萬人[1]，其中約90%的病例被歸類為頭頸部鱗癌(HNSCC)。[2]盡管當前在治療方式上取得了重大進展，但HNSCC患者的5年生存率并沒有得到明顯的提高。[3]近年來隨著腫瘤免疫微環(huán)境領域的發(fā)展，有研究表明HNSCC患者的預后與腫瘤組織免疫微環(huán)境密切相關。[4，5]這提示我們可以識別免疫相關的預后指標，以改善HNSCC患者的預后并指導他們的治療。

長鏈非編碼RNA(ln cRNA)被定義為大于等于200個核苷酸的RNA，目前還沒有證據(jù)表明它們可以被翻譯成肽。隨著研究的深入，ln cRNA已經(jīng)被證明在腫瘤的發(fā)生和發(fā)展中起著不可或缺的作用。[6,7]研究表明，其中與免疫相關的ln cRNA在癌癥的預后中具有獨特的價值，免疫相關ln cRNA作為新興的癌癥生物標志物已被用于多種癌癥的診斷和生存預測。例如，Shen等人[8]在乳腺癌中識別了11種免疫相關ln cRNA并構建預后標志物，并且證明了這個標志物和免疫細胞亞型的浸潤有關；Li等人[9]在肺腺癌中識別了7個與免疫相關的ln cRNA并構建穩(wěn)健的預測模型，以提高肺腺癌的預后預測。同樣在HNSCC中，Wu等人[10]整合了與臨床數(shù)據(jù)相關的ln cRNA、microRNA和mRNA，并建立了一種新的免疫相關RNA特征來預測HNSCC患者的存活率；Yin等人[11]基于13個免疫相關ln cRNA對來建立HNSCC的預后模型，但在這項研究中僅基于訓練集數(shù)據(jù)進行了模型的構建和驗證。因此，本文旨在提出一種新的基于免疫相關ln cRNA的生存預后模型來預測HNSCC的預后，并在訓練集和測試集上對構建的模型進行驗證。

本文獲取了HNSCC患者的ln cRNA表達數(shù)據(jù)和相應的臨床信息，將樣本隨機分為訓練集和測試集?；谟柧毤?，篩選得到HNSCC中具有預后價值的免疫相關ln cRNA，建立預后風險評分模型；并在訓練集和測試集上對模型的性能進行了驗證。結果表明，基于7個免疫相關ln cRNA構建的預后模型可作為一個獨立的預后因素用于HNSCC的生存預后預測，并且這7個ln cRNA可能成為HNSCC的潛在治療靶點。

1 材料與方法

1.1 數(shù)據(jù)來源及處理

HNSCC患者的轉錄組數(shù)據(jù)來自癌癥基因圖譜數(shù)據(jù)庫(The Cancer Genome Atlas，TCGA)(https://portal.gdc.cancer.gov/)中的HNSCC項目組，并從中提取了lncRNA表達數(shù)據(jù)；HNSCC患者的臨床病理信息來自UCSC Xena(https://xenabrowser.net/)。將HNSCC患者的lncRNA表達數(shù)據(jù)和臨床數(shù)據(jù)進行整理和合并，使得整理后的樣本均有對應的表達數(shù)據(jù)和臨床數(shù)據(jù)，并刪除生存信息缺失和生存時間小于30天的樣本。

1.2 免疫相關lncRNA的提取

R軟件中的edgeR是基于負二項分布的統(tǒng)計方法，根據(jù)樣本的ln cRNA表達量，選擇FDR<0.05和|log2FC|≥1作為閾值篩選出在腫瘤樣本和非腫瘤樣本之間有明顯差異表達的ln cRNA。另外從免疫學數(shù)據(jù)庫(Immunology Database and Analysis Portal，ImmPort)(https://www.immport.org/home)網(wǎng)站獲得免疫相關基因列表，采用Pearson相關系數(shù)計算并確定與免疫基因相關的lncRNA。

1.3 預后模型的構建

首先，對免疫相關lncRNA進行單因素Cox回歸分析，篩選出與總生存期(OS)顯著相關的ln cRNA(P<0.01)。然后，建立Lasso回歸模型，并進行10折交叉驗證確定最優(yōu)模型來進一步篩選出關鍵的ln cRNA。最后利用多因素Cox回歸模型進行l(wèi)n cRNA的篩選，并在建模過程中根據(jù)簡單可行的Akaike信息準則(AIC)采用逐步回歸法選擇最佳模型。根據(jù)多因素Cox回歸分析中加權線性組合的回歸系數(shù)β和ln cRNA的表達量構建風險評分預后模型并利用模型來確定每個患者的風險評分。風險評分模型一般具有如下的形式[8,9]

rs=βln cRNA1×exprln cRNA1+βln cRNA2×exprln cRNA2+…+βln cRNAn×exprln cRNAn

根據(jù)風險評分的中位數(shù)將患者劃分為高險組和低風險組。

1.4 預后模型的驗證

采用Kaplan-Meier方法繪制生存曲線，并通過Log-rank檢驗來比較高低風險組的生存是否具有差異；通過單因素和多因素Cox回歸分析進行獨立預后分析；并繪制時間依賴性ROC曲線和計算ROC曲線下的面積(AUC值)以驗證模型的預測效果。

2 實驗結果

2.1 數(shù)據(jù)處理

從UCSC Xena網(wǎng)站下載了612名HNSCC患者的臨床信息，并將臨床信息和基因表達數(shù)據(jù)進行匹配，得到539個匹配樣本；其中包括495個腫瘤樣本和44個正常組織樣本；去除生存信息缺失和生存時間小于30天的樣本，最終得到486個樣本。

根據(jù)基因注釋文件，從HNSCC的全基因轉錄組數(shù)據(jù)中得到15 878個ln cRNA的表達數(shù)據(jù)，通過差異表達分析(FDR<0.05,|log2FC|≥1)得到1 729個差異表達的ln cRNA。根據(jù)ImmPort數(shù)據(jù)庫下載的免疫相關基因信息，從轉錄組數(shù)據(jù)中提取了1 279個免疫相關基因的表達數(shù)據(jù)；采用Pearson相關分析篩選與免疫相關基因相關的lncRNA，結果確定了988個免疫相關ln cRNA(|R|>0.4，P<0.001)。

2.2 關鍵免疫相關ln cRNA的篩選和預后模型的構建

按照3:2的比例將樣本隨機分為訓練集和測試集。在訓練集中，對包含988個ln cRNA的表達數(shù)據(jù)進行了單因素Cox回歸分析，以確定具有預后差異的免疫相關ln cRNA。從分析中發(fā)現(xiàn)，有20個免疫相關ln cRNA與HNSCC患者的OS相關(P<0.01)。再基于這20個ln cRNA建立Lasso回歸模型，通過10折交叉驗證確定最優(yōu)模型，選擇使模型達到最小損失時的懲罰參數(shù)λ(λ=0.019)，得到加入模型的14個ln cRNA。接著對這14個ln cRNA，根據(jù)AIC值進一步建立多因素Cox回歸模型，該模型的C指數(shù)為0.7。最終，將篩選得到的7個免疫相關ln cRNA納入風險評分模型，其相應系數(shù)見表1。風險評分模型如下：

表1 多因素Cox回歸分析篩選出納入預后模型的免疫相關ln cRNA

風險評分=[RP11-91K9.1表達量*(0.100 39)]+[ZFY-AS1表達量*(-0.200 79)]+[PRKG1-AS1表達量*(0.189 03)]+[HOTTIP表達量*(0.190 59)]+[RP11-401O9.3表達量*(0.157 04)]+[RP11-466A19.1表達量*(-0.150 03)]+[RP5-1057J7.7表達量*(-0.208 11)]。

2.3 預后模型的驗證

在訓練集上，根據(jù)上述風險評分模型計算出每個病人的風險得分，并根據(jù)風險評分的中位數(shù)將訓練集中風險評分大于中位數(shù)的患者歸為高風險組，反之則歸為低風險組(圖1(a))；所有病例的生存狀態(tài)分布見圖1(b)。熱圖(圖1(c))顯示了高風險組和低風險組之間7個ln cRNA的表達水平，這7個ln cRNA表達水平的變化趨勢與風險評分計算公式中對應的系數(shù)一致。風險評分隨著RP11-91K9.1、PRKG1-AS1、HOTTIP和RP11-401O9.3的表達水平的增加而增加，隨著ZFY-AS1、RP11-466A19.1和RP5-1057J7.7的表達水平的增加而減少。Kaplan-Meier生存分析結果顯示，低風險組患者3年總生存率比高風險組患者高約30%，5年總生存率高約20%(圖1(d))。此外，繪制了時間依賴性ROC曲線來評估模型的預后價值，計算得到3年和5年總生存率AUC值分別為0.770和0.688(圖1(e))。為了進一步評估生存預后模型的性能，在測試集中進行了相同的分析，并且得到了類似的結果，結果見圖2?？紤]到其他臨床特征變量(如年齡、性別、分期和分級)，進行了單因素和多因素Cox回歸分析來驗證風險評分可作為HNSCC的獨立預后因素。單因素Cox回歸表明，風險評分與較差的生存率相關；多因素Cox回歸分析結果顯示，基于7個關鍵ln cRNA的風險評分是HNSCC患者的獨立預后因素(P<0.05)(圖3)。綜上可知，該模型具有一定的生存預后預測能力。

(a)風險評分的分布

(b)生存狀態(tài)和時間

(c)7個ln cRNA表達水平的熱圖

(d)高風險組和低風險組患者的Kaplan-Meier生存曲線 (e)風險評分的ROC曲線圖1 訓練集預后風險評分模型的構建

(a)風險評分的分布 (b)生存狀態(tài)和時間

(c) 7個ln cRNA表達水平的熱圖

(d) 高風險組和低風險組患者的Kaplan-Meier生存曲線

(e)風險評分的ROC曲線 圖2 測試集預后風險評分模型的驗證

(a)訓練集的單因素Cox回歸分析 (b)訓練集的多因素Cox回歸分析

(c) 測試集的單因素Cox回歸分析 (d) 測試集的多因素Cox回歸分析

(e)整個數(shù)據(jù)集的單因素Cox回歸分析 (f)整個數(shù)據(jù)集的多因素Cox回歸分析圖3 風險評分和臨床特征的Cox回歸分析

3 討 論

近年來，ln cRNA成為腫瘤研究的熱點，人們發(fā)現(xiàn)ln cRNA在腫瘤的發(fā)生和發(fā)展中起著重要作用，ln cRNA的異常表達被認為是影響腫瘤活性的關鍵因素；其中與免疫相關的ln cRNA被證明與惡性腫瘤的不良預后有關。[12,13]因此，我們的研究目的是通過對TCGA數(shù)據(jù)庫中HNSCC患者數(shù)據(jù)的生物信息學分析，確定與HNSCC預后相關的免疫相關ln cRNA。通過對HNSCC患者的相關數(shù)據(jù)進行單因素Cox回歸分析、Lasso回歸分析和多因素Cox回歸分析，我們確定了7個關鍵的免疫相關ln cRNA(RP11-91K9.1、PRKG1-AS1、HOTTIP、RP11-401O9.3、ZFY-AS1、RP11-466A19.1、RP5-1057J7.7)，建立了用于HNSCC生存預后的風險評分模型，并根據(jù)風險評分的中位值將HNSCC患者劃分為高風險組和低風險組。Kaplan-Meier曲線顯示，基于模型劃分的高風險組和低風險組的生存率存在差異且具有統(tǒng)計學意義(P<0.05)。同時在訓練集和測試集中，模型5年生存率的ROC曲線的AUC值分別為0.688和0.635，表明該模型具有較好的預測能力。另外在篩選得到的7個ln cRNA中，ln cRNA HOTTIP被認為是幾乎所有類型癌癥中的致癌ln cRNA，也是人類惡性腫瘤中公認的生物標志物和治療靶點[14]；ln cRNA PRKG1-AS1被證明在口腔鱗癌的細胞生長、遷移和侵襲中起促進作用[15]。而其他ln cRNA目前沒有相關的報道，還需要進行進一步的探索。

本文基于對TCGA數(shù)據(jù)庫中HNSCC患者數(shù)據(jù)的分析，確定了7個與預后相關的免疫相關ln cRNA，并成功構建了一個用于HNSCC患者預后評估的風險評分模型。這個預后模型能在一定程度上反應HNSCC患者的預后情況，并且這7個免疫相關ln cRNA有望成為HNSCC潛在的免疫治療靶標。