李夢涵，肖瓊，高鵬，付昱，孫晨蕊，宋永喜

0 引言

我國是消化系統(tǒng)腫瘤高發(fā)國家。2020年全球癌癥統(tǒng)計(jì)報(bào)告顯示，結(jié)腸癌、胃癌、直腸癌、食管癌均位于世界新發(fā)與死亡癌癥的前十位，占新發(fā)病例總數(shù)的18.5%、新增死亡例數(shù)的22.4%[1]。得益于免疫治療、靶向治療等新技術(shù)的不斷發(fā)展[2-3]，消化道腫瘤的發(fā)病率、死亡率持續(xù)下降，但患者的總生存期（overall survi val,OS）依然較低，預(yù)后不良[4]。因此，我們急需行之有效的預(yù)后預(yù)測模型來指導(dǎo)臨床治療。

腫瘤分子標(biāo)志物在基因組學(xué)的研究不斷深入，使得腫瘤生物學(xué)研究逐漸步入泛癌研究階段。泛癌同時(shí)對多種不同類型的腫瘤基因組進(jìn)行分析，能夠?qū)ふ也煌[瘤的共同特征，尋找多種腫瘤的診斷、治療及預(yù)后廣譜靶點(diǎn)[5]。由于消化道解剖結(jié)構(gòu)與生理功能相似，食管癌、胃癌、結(jié)腸癌與直腸癌數(shù)據(jù)合并分析并建立消化道泛癌的預(yù)后模型有助于尋找四種腫瘤在發(fā)生發(fā)展過程中的相同點(diǎn)。

腫瘤的精準(zhǔn)治療要求根據(jù)患者特異性指標(biāo)尋找精準(zhǔn)的生物標(biāo)志物。長非編碼RNA是一類特殊的的基因轉(zhuǎn)錄本，雖然絕大多數(shù)不參與蛋白質(zhì)的編碼[6-7]，但與轉(zhuǎn)錄及轉(zhuǎn)錄后調(diào)控密切相關(guān)[8]，廣泛參與了腫瘤的發(fā)生發(fā)展[9]，能夠成為腫瘤診斷和預(yù)后的潛在標(biāo)志物[10]。

癌癥基因組圖譜TCGA數(shù)據(jù)庫共收集了33種癌癥的11000多例患者的相關(guān)數(shù)據(jù)，極大地提高了研究者對腫瘤生物學(xué)的理解[11]。本研究收集整理了來自TCGA數(shù)據(jù)庫的消化道腫瘤相關(guān)數(shù)據(jù)，運(yùn)用Cox回歸和Lasso回歸分析進(jìn)行預(yù)后相關(guān)基因篩選和模型構(gòu)建，以期能夠?qū)εR床治療進(jìn)行輔助和指導(dǎo)。

1 資料與方法

1.1 資料來源

從UCSC-Xena（http://xena.ucsc.edu）數(shù)據(jù)庫下載經(jīng)過歸一化處理后的TCGA數(shù)據(jù)庫中的消化道腫瘤患者（共1117例，其中食管癌172例、胃癌381例、結(jié)腸癌487例、直腸癌167例）的RNA測序（RNA Seq）數(shù)據(jù)和臨床相關(guān)數(shù)據(jù)，整合數(shù)據(jù)，用于后續(xù)分析。

1.2 方法

1.2.1 消化道腫瘤lncRNA表達(dá)數(shù)據(jù)及臨床資料收集與處理 首先對病例樣本進(jìn)行篩選：（1）排除無完整臨床病理資料的病例；（2）排除預(yù)后相關(guān)信息不全的病例。納入有完整臨床病理及預(yù)后資料（包括性別、年齡、病理分期和生存狀況）的病例，Ensemble數(shù)據(jù)庫進(jìn)行基因注釋得到納入研究的lncRNA表達(dá)譜。

1.2.2 消化道腫瘤lncRNA差異分析及單因素Cox分析 使用配對t檢驗(yàn)對94例患者的腫瘤樣本及配對正常樣本進(jìn)行統(tǒng)計(jì)學(xué)分析，計(jì)算logFC值并初步篩選差異表達(dá)的lncRNA。對結(jié)果進(jìn)行單因素Cox回歸分析篩選預(yù)后相關(guān)的差異lncRNA。

1.2.3 Lasso及多因素Cox回歸分析 最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator,Lasso）算法，可調(diào)整自變量的回歸系數(shù)，首先將對模型貢獻(xiàn)小的變量的系數(shù)壓縮為零，再進(jìn)行剔除，降低模型過度擬合的程度，使模型中均為相對更加重要的變量。

將1117例樣本隨機(jī)分成訓(xùn)練集（894例）和驗(yàn)證集（223例）兩部分。R語言中的survival包和glmnet包進(jìn)行Lasso回歸，選擇最優(yōu)的λ納入模型中，得到與預(yù)后顯著相關(guān)的lncRNA，再利用多因素Cox回歸分析計(jì)算每個(gè)lncRNA的回歸系數(shù)，根據(jù)每個(gè)lncRNA的表達(dá)值及回歸系數(shù)按照如下公式構(gòu)建消化道腫瘤預(yù)后風(fēng)險(xiǎn)模型：

Risk Score(RS)=∑explncRNAn*αlncRNAn

該公式中，Risk Score代表風(fēng)險(xiǎn)值，exp代表每個(gè)lncRNA在樣本中的表達(dá)值，α代表多因素Cox回歸分析中計(jì)算得到的每個(gè)lncRNA的回歸系數(shù)。

1.2.4 模型效能評估 基于預(yù)后風(fēng)險(xiǎn)模型，得出每例患者的評分，將患者進(jìn)一步分為高、低風(fēng)險(xiǎn)組。同時(shí)將模型中的各個(gè)基因在不同腫瘤高低風(fēng)險(xiǎn)組中的表達(dá)量進(jìn)行分析。

采用時(shí)間依賴的受試者工作特征曲線（receiver operating characteristic curve,ROC）評估模型預(yù)測性能，Kaplan-Meier生存曲線和Log rank檢驗(yàn)分析患者生存差異，Cox回歸綜合評價(jià)包括模型風(fēng)險(xiǎn)評分在內(nèi)的各個(gè)臨床病理特征對預(yù)后的影響。

1.2.5 功能富集分析 根據(jù)每例消化道腫瘤樣本的高低風(fēng)險(xiǎn)狀態(tài)劃分，R包“clusterProfiler”對高風(fēng)險(xiǎn)組與低風(fēng)險(xiǎn)組mRNA表達(dá)數(shù)據(jù)進(jìn)行差異分析后分別對上調(diào)基因與下調(diào)基因進(jìn)行基因富集分析。

1.3 統(tǒng)計(jì)學(xué)方法

數(shù)據(jù)預(yù)處理、模型的構(gòu)建、驗(yàn)證及評價(jià)均基于R語言（Version 4.0.3）及相關(guān)的R程序包。所有統(tǒng)計(jì)學(xué)檢驗(yàn)均為雙側(cè)檢驗(yàn)，P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 差異表達(dá)lncRNA篩選

本研究以校正后的P<0.05且差異表達(dá)倍數(shù)>1（FDR<0.05且|logFC|>1）為篩選標(biāo)準(zhǔn)，最終得到了1426個(gè)差異表達(dá)的lncRNA，其中1280個(gè)lncRNA表達(dá)上調(diào)，146個(gè)lncRNA表達(dá)下調(diào)。

2.2 構(gòu)建模型并評估效能

采用單因素Cox回歸在全體差異表達(dá)的lncRNA中篩選與預(yù)后顯著相關(guān)的lncRNA，以P<0.001為標(biāo)準(zhǔn)選擇出94個(gè)lncRNA。為避免模型過度擬合并篩選出與消化道腫瘤患者預(yù)后更加相關(guān)的lncRNA，對單因素分析得到的lncRNA進(jìn)行Lasso回歸分析，λ取最優(yōu)值時(shí)，得到36個(gè)lncRNA（圖請掃描本文OSID碼），對結(jié)果繼續(xù)進(jìn)行多因素Cox回歸分析，得到最相關(guān)的13 個(gè)lncRNA，見圖1A。計(jì)算每個(gè)lncRNA的危險(xiǎn)比及95%CI，通過將每個(gè)lncRNA對應(yīng)的系數(shù)代入公式得到消化道腫瘤預(yù)后風(fēng)險(xiǎn)評分模型，具體如下：Riskscore=0.5117*expAC006357.1+0.9657* expAC010998.1+0.5926*expAC011586.2+(-0.4004)*expAC099792.1+(-0.3537)*expAC110772.1 +0.4013*expAC12493.1+0.3443*expAL513123.1+(-0.6706)* expAL590483.1+0.8417*expDSCR10+1.0330*expLINC02382+0.7515*expLINC02866+(-0.2830)* expMIR3142H G+0.3331*expRHOXF1P1。

利用該模型對訓(xùn)練集894例樣本進(jìn)行時(shí)間依賴ROC分析，結(jié)果顯示，三年曲線下面積（AUC）=0.746，提示該模型對患者預(yù)后具有良好的預(yù)測能力，見圖1B。

根據(jù)模型，計(jì)算訓(xùn)練集患者的風(fēng)險(xiǎn)評分，采用評分中位數(shù)（1.279）作為cut-off值將所有樣本分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。生存時(shí)間分布圖顯示，隨著風(fēng)險(xiǎn)評分的升高，消化道腫瘤患者死亡例數(shù)逐漸增多，見圖1C。K-M生存曲線的Log rank檢驗(yàn)顯示，比較五年總生存率，低分險(xiǎn)組顯著高于高風(fēng)險(xiǎn)組（P<0.01），見圖1D。

分析模型中各個(gè)基因在食管癌、胃腺癌、結(jié)腸癌和直腸癌高低風(fēng)險(xiǎn)組中的表達(dá)，結(jié)果顯示，在絕大多數(shù)分組中，高低風(fēng)險(xiǎn)組之間差異顯著，見圖2。

2.3 lncRNA預(yù)后模型的驗(yàn)證與獨(dú)立性檢驗(yàn)

利用223例驗(yàn)證集樣本數(shù)據(jù)進(jìn)行內(nèi)部驗(yàn)證分析，結(jié)果顯示，該模型的三年AUC為0.704，提示該模型在驗(yàn)證集中預(yù)測能力良好，具有較好的外推性，見圖3A。生存分析顯示，低風(fēng)險(xiǎn)組總生存率依然顯著高于高風(fēng)險(xiǎn)組（P<0.01），見圖3B。與其他臨床性狀評價(jià)模型效果相比，風(fēng)險(xiǎn)評分作為變量時(shí)AUC最大，見圖3C。

由于本研究需要評價(jià)該13-lncRNA模型是否能獨(dú)立于其他臨床病理特征進(jìn)行預(yù)后預(yù)測，將894例樣本數(shù)據(jù)進(jìn)行Cox回歸分析，提示該預(yù)后風(fēng)險(xiǎn)評分模型可以獨(dú)立于其他臨床病理特征進(jìn)行預(yù)后預(yù)測，見圖3D。同時(shí)，熱圖提示風(fēng)險(xiǎn)評分與年齡、性別、TNM分期、T分級和N分級顯著相關(guān)，見圖3E。

為證明模型在不同類型腫瘤中的預(yù)測價(jià)值，本研究將該預(yù)后模型在172例食管癌，381例胃癌，487例結(jié)腸癌與167例直腸癌數(shù)據(jù)中獨(dú)立進(jìn)行驗(yàn)證。結(jié)果顯示，四種腫瘤的高低風(fēng)險(xiǎn)組均有顯著的生存差異，見圖4。

2.4 功能富集分析

為進(jìn)一步確認(rèn)消化道腫瘤發(fā)生發(fā)展過程中涉及到的基因功能、信號通路，將高低風(fēng)險(xiǎn)組中具有顯著差異的基因進(jìn)行GO、KEGG功能富集分析。通路主要富集在細(xì)胞頂端成分的組成和檢測知覺等化學(xué)刺激的生物過程中，見圖5。

3 討論

消化道腫瘤早期癥狀隱匿，并且篩查在大部分地區(qū)未得到普及，患者往往出現(xiàn)癥狀后到醫(yī)療機(jī)構(gòu)就診，診斷結(jié)果很多已處于晚期，手術(shù)治療后容易出現(xiàn)轉(zhuǎn)移復(fù)發(fā)，嚴(yán)重影響患者的預(yù)后[12-14]。本研究通過制定風(fēng)險(xiǎn)評分模型輔助臨床評估患者的預(yù)后。

既往研究大多采用血清學(xué)標(biāo)志物等對消化道腫瘤的預(yù)后進(jìn)行預(yù)測，Zheng等對血清中外源性miR-590-5p作為胃癌診斷和預(yù)后生物標(biāo)志物的有效性進(jìn)行了研究驗(yàn)證[15]，但這些標(biāo)志物單獨(dú)應(yīng)用的預(yù)測效能不高，且不同研究沒有一致的判斷標(biāo)準(zhǔn)。本研究使用lncRNA測序數(shù)據(jù)構(gòu)建模型，隨著基因測序技術(shù)的不斷突破為腫瘤研究提供大量數(shù)據(jù)，曾被認(rèn)為不發(fā)揮生物學(xué)功能的lncRNA被證實(shí)通過染色體修飾、轉(zhuǎn)錄激活和干擾等方式，在轉(zhuǎn)錄和轉(zhuǎn)錄后及表觀遺傳水平調(diào)控相關(guān)靶基因，從而在多種惡性腫瘤的發(fā)生發(fā)展過程中發(fā)揮作用[16]。因此使用lncRNA構(gòu)建的預(yù)后風(fēng)險(xiǎn)評分模型是綜合考慮每一患者基因?qū)W的疾病預(yù)后評價(jià)體系。不同部位的腫瘤細(xì)胞的形態(tài)、增殖和轉(zhuǎn)移等顯著不同造成的腫瘤異質(zhì)性是腫瘤治療面臨的重要難題[17]，但是這些腫瘤細(xì)胞在分子水平上的相似性超過同部位腫瘤內(nèi)部的相似性，因此根據(jù)患者的基因信息等可以進(jìn)行準(zhǔn)確治療與預(yù)后評估，從而實(shí)現(xiàn)精準(zhǔn)治療[18]。

消化道腫瘤泛癌分析通過聯(lián)合食管癌、胃癌、結(jié)腸癌與直腸癌數(shù)據(jù)分析腫瘤之間的相似與差異，從分子水平上認(rèn)識腫瘤發(fā)展導(dǎo)致不同患者出現(xiàn)不同預(yù)后狀況的機(jī)制，最終指導(dǎo)預(yù)后判斷。目前已有的預(yù)后模型大多為單個(gè)腫瘤的研究，Liu等對食管癌預(yù)后相關(guān)lncRNA進(jìn)行研究并構(gòu)建6-lncRNA風(fēng)險(xiǎn)評分模型[19]，該模型AUC為0.785，提示lncRNA模型對食管癌患者預(yù)后風(fēng)險(xiǎn)具有較好的預(yù)測效能。但是不能從整體上對消化道腫瘤的發(fā)展過程做出分析，從而無法基于消化道的相似性做出整體預(yù)后評價(jià)。本文是針對消化道泛癌的分析，彌補(bǔ)了當(dāng)下科學(xué)研究中消化道泛癌預(yù)后狀態(tài)評價(jià)方面的空白。

本研究尚具有一定的局限性：（1）高通量測序數(shù)據(jù)具有一定的背景噪音，即使在分析前進(jìn)行了批次矯正，仍需要臨床試驗(yàn)來證明；（2）通過生物信息學(xué)及統(tǒng)計(jì)學(xué)方法獲得13基因模型，分子機(jī)制仍需進(jìn)一步通過基礎(chǔ)實(shí)驗(yàn)驗(yàn)證。

綜上所述，本研究基于TCGA數(shù)據(jù)庫中的消化道腫瘤患者的lncRNA表達(dá)和預(yù)后數(shù)據(jù)，使用Lasso和Cox回歸，構(gòu)建了消化道腫瘤lncRNA預(yù)后風(fēng)險(xiǎn)評分模型。結(jié)果顯示：模型對消化道腫瘤患者預(yù)后具有較好的預(yù)測能力，且模型所得評分為獨(dú)立預(yù)后因素。