吳嘉雯，唐加山

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

0 前 言

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

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

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

1 材料與方法

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

HNSCC患者的轉(zhuǎn)錄組數(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ù)進行整理和合并，使得整理后的樣本均有對應(yīng)的表達數(shù)據(jù)和臨床數(shù)據(jù)，并刪除生存信息缺失和生存時間小于30天的樣本。

1.2 免疫相關(guān)lncRNA的提取

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

1.3 預(yù)后模型的構(gòu)建

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

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

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

1.4 預(yù)后模型的驗證

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

2 實驗結(jié)果

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

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

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

2.2 關(guān)鍵免疫相關(guān)ln cRNA的篩選和預(yù)后模型的構(gòu)建

按照3:2的比例將樣本隨機分為訓(xùn)練集和測試集。在訓(xùn)練集中，對包含988個ln cRNA的表達數(shù)據(jù)進行了單因素Cox回歸分析，以確定具有預(yù)后差異的免疫相關(guān)ln cRNA。從分析中發(fā)現(xiàn)，有20個免疫相關(guān)ln cRNA與HNSCC患者的OS相關(guān)(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個免疫相關(guān)ln cRNA納入風(fēng)險評分模型，其相應(yīng)系數(shù)見表1。風(fēng)險評分模型如下：

表1 多因素Cox回歸分析篩選出納入預(yù)后模型的免疫相關(guān)ln cRNA

風(fēng)險評分=[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 預(yù)后模型的驗證

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

(a)風(fēng)險評分的分布

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

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

(d)高風(fēng)險組和低風(fēng)險組患者的Kaplan-Meier生存曲線 (e)風(fēng)險評分的ROC曲線圖1 訓(xùn)練集預(yù)后風(fēng)險評分模型的構(gòu)建

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

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

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

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

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

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

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

3 討 論

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

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