吳嘉雯，唐加山

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

0 前 言

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

長鏈非編碼RNA(ln cRNA)被定義為大于等于200個核苷酸的RNA，目前還沒有證據(jù)表明它們可以被翻譯成肽。隨著研究的深入，ln cRNA已經(jīng)被證明在腫瘤的發(fā)生和發(fā)展中起著不可或缺的作用。[6,7]研究表明，其中與免疫相關(guān)的ln cRNA在癌癥的預(yù)后中具有獨(dú)特的價值，免疫相關(guān)ln cRNA作為新興的癌癥生物標(biāo)志物已被用于多種癌癥的診斷和生存預(yù)測。例如，Shen等人[8]在乳腺癌中識別了11種免疫相關(guān)ln cRNA并構(gòu)建預(yù)后標(biāo)志物，并且證明了這個標(biāo)志物和免疫細(xì)胞亞型的浸潤有關(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ù)后模型，但在這項(xiàng)研究中僅基于訓(xùn)練集數(shù)據(jù)進(jìn)行了模型的構(gòu)建和驗(yàn)證。因此，本文旨在提出一種新的基于免疫相關(guān)ln cRNA的生存預(yù)后模型來預(yù)測HNSCC的預(yù)后，并在訓(xùn)練集和測試集上對構(gòu)建的模型進(jìn)行驗(yàn)證。

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

1 材料與方法

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)結(jié)果

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

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

根據(jù)基因注釋文件，從HNSCC的全基因轉(zhuǎn)錄組數(shù)據(jù)中得到15 878個ln cRNA的表達(dá)數(shù)據(jù)，通過差異表達(dá)分析(FDR<0.05,|log2FC|≥1)得到1 729個差異表達(dá)的ln cRNA。根據(jù)ImmPort數(shù)據(jù)庫下載的免疫相關(guān)基因信息，從轉(zhuǎn)錄組數(shù)據(jù)中提取了1 279個免疫相關(guān)基因的表達(dá)數(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的比例將樣本隨機(jī)分為訓(xùn)練集和測試集。在訓(xùn)練集中，對包含988個ln cRNA的表達(dá)數(shù)據(jù)進(jìn)行了單因素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àn)證確定最優(yōu)模型，選擇使模型達(dá)到最小損失時的懲罰參數(shù)λ(λ=0.019)，得到加入模型的14個ln cRNA。接著對這14個ln cRNA，根據(jù)AIC值進(jìn)一步建立多因素Cox回歸模型，該模型的C指數(shù)為0.7。最終，將篩選得到的7個免疫相關(guān)ln cRNA納入風(fēng)險(xiǎn)評分模型，其相應(yīng)系數(shù)見表1。風(fēng)險(xiǎn)評分模型如下：

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

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

2.3 預(yù)后模型的驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

3 討 論

近年來，ln cRNA成為腫瘤研究的熱點(diǎn)，人們發(fā)現(xiàn)ln cRNA在腫瘤的發(fā)生和發(fā)展中起著重要作用，ln cRNA的異常表達(dá)被認(rèn)為是影響腫瘤活性的關(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ù)進(jìn)行單因素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)險(xiǎn)評分模型，并根據(jù)風(fēng)險(xiǎn)評分的中位值將HNSCC患者劃分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。Kaplan-Meier曲線顯示，基于模型劃分的高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組的生存率存在差異且具有統(tǒng)計(jì)學(xué)意義(P<0.05)。同時在訓(xùn)練集和測試集中，模型5年生存率的ROC曲線的AUC值分別為0.688和0.635，表明該模型具有較好的預(yù)測能力。另外在篩選得到的7個ln cRNA中，ln cRNA HOTTIP被認(rèn)為是幾乎所有類型癌癥中的致癌ln cRNA，也是人類惡性腫瘤中公認(rèn)的生物標(biāo)志物和治療靶點(diǎn)[14]；ln cRNA PRKG1-AS1被證明在口腔鱗癌的細(xì)胞生長、遷移和侵襲中起促進(jìn)作用[15]。而其他ln cRNA目前沒有相關(guān)的報(bào)道，還需要進(jìn)行進(jìn)一步的探索。

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