王亮 周璇 梁曉杰 何穎芝

作者單位：100730 北京 1首都醫(yī)科大學(xué)附屬北京同仁醫(yī)院血液科；2北京大數(shù)據(jù)精準(zhǔn)醫(yī)療高精尖創(chuàng)新中心，北京航空航天大學(xué)&首都醫(yī)科大學(xué)，北京同仁醫(yī)院；510280 廣州 3南方醫(yī)科大學(xué)第二臨床醫(yī)學(xué)院，南方醫(yī)科大學(xué)珠江醫(yī)院血液科

彌漫大B細(xì)胞淋巴瘤（diffuse large B-cell lymphoma，DLBCL）為最常見的非霍奇金淋巴瘤亞型。隨著利妥昔單抗時(shí)代的到來，DLBCL患者的治療效果獲得極大改善，但仍有部分患者復(fù)發(fā)，預(yù)后較差［1］。研究表明國際預(yù)后指數(shù)（international prognostic index，IPI）、改良國際預(yù)后指數(shù)（revised-IPI，R-IPI）及 NCCN-IPI等在預(yù)后評估中有較好的實(shí)用性［2-3］。近年來反映腫瘤炎癥微環(huán)境的指標(biāo)也被用于DLBCL的預(yù)后評估［4］。然而，對預(yù)后有重要影響的遺傳學(xué)標(biāo)記、細(xì)胞起源等因素并未納入其中，這可能造成部分患者預(yù)后評估的偏差。隨著基因芯片及測序技術(shù)的廣泛應(yīng)用，基因標(biāo)記亦被廣泛應(yīng)用于腫瘤患者的風(fēng)險(xiǎn)分層［5］。本研究基于公開的數(shù)據(jù)集，構(gòu)建可預(yù)測DLBCL患者生存的基因預(yù)測模型，并將其與其他預(yù)后評估指標(biāo)相結(jié)合，以期能夠更好地評估患者的預(yù)后，以進(jìn)一步指導(dǎo)治療方案的選擇。

1 方法

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

分別從 Gene Expression Omnibus（GEO）和The Cancer Genome Atlas（TCGA）下載DLBCL患者的基因表達(dá)數(shù)據(jù)和臨床數(shù)據(jù)，包括GSE10846、GSE32918、NCICCR數(shù)據(jù)集。將隨訪信息不全及隨訪時(shí)間≤0 d的樣本剔除后，納入GSE10846數(shù)據(jù)集 412例患者，GSE32918數(shù)據(jù)集172例患者，NCICCR數(shù)據(jù)集234例患者；剔除GSE32918數(shù)據(jù)集中未采用標(biāo)準(zhǔn)治療的32例患者后，共786例患者符合標(biāo)準(zhǔn)納入分析。GEO數(shù)據(jù)集的原始數(shù)據(jù)處理步驟如下：⑴使用limma包中的normalizeBetweenArrays函數(shù)對GEO樣本進(jìn)行組間矯正。⑵GSE10846、GSE32918數(shù)據(jù)集的探針分別使用GPL570、GPL8423平臺文件進(jìn)行注釋。對于對應(yīng)多個探針的基因，保留最大平均值。NCICCR數(shù)據(jù)集的原始數(shù)據(jù)處理步驟如下：⑴NCICCR數(shù)據(jù)集的IDs使用hg38參考基因組進(jìn)行注釋，對于對應(yīng)多個IDs的基因，保留最大平均值。⑵使用limma包的voom函數(shù)對NCICCR原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理［6］。

1.2 模型的構(gòu)建及驗(yàn)證

采用單因素Cox回歸篩選GSE10846、GSE32918數(shù)據(jù)集中與DLBCL患者總生存期（OS）相關(guān)的基因［7］。風(fēng)險(xiǎn)比（HR）>1為危險(xiǎn)基因，HR<1為保護(hù)基因，將P＜0.01設(shè)定為基因與生存具有相關(guān)性的閾值標(biāo)準(zhǔn)。將兩個數(shù)據(jù)集共有的危險(xiǎn)基因及保護(hù)基因取交集，并與Lasso回歸、Cox回歸相結(jié)合，構(gòu)建最終的基因風(fēng)險(xiǎn)預(yù)測模型［8］。

將預(yù)后基因的表達(dá)值與其回歸系數(shù)加權(quán)后計(jì)算風(fēng)險(xiǎn)評分［9］。以各隊(duì)列的中位風(fēng)險(xiǎn)評分作為分界點(diǎn)，將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。采用Kaplan-Meier法評估高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組DLBCL患者的生存率，組間比較采用Log-rank檢驗(yàn)；采用多因素Cox回歸分析及分層分析評估風(fēng)險(xiǎn)評分與OS的關(guān)系；采用受試者工作特征曲線（ROC曲線）的下面積評價(jià)該模型的預(yù)測價(jià)值并與其他臨床特征進(jìn)行比較［10］；在多因素Cox回歸分析的基礎(chǔ)上構(gòu)建列線圖模型，根據(jù)模型中各影響因素對結(jié)局變量的貢獻(xiàn)程度，給影響因素的每個取值水平進(jìn)行賦值，再將各個評分相加得到總評分，計(jì)算個體結(jié)局事件的預(yù)測值，并評估列線圖的預(yù)測能力［11］；采用Pearson相關(guān)性檢驗(yàn)識別與風(fēng)險(xiǎn)評分有關(guān)的基因（風(fēng)險(xiǎn)評分相關(guān)基因定義為r>0.4且P<0.001），并對這些基因進(jìn)行GO與KEGG分析，探討相關(guān)基因的生物學(xué)功能［12］。詳細(xì)流程見圖1。

采用R（版本3.5.2）軟件進(jìn)行數(shù)據(jù)分析。采用“glmnet”包行 Lasso 回歸分析，“survival”包行單因素、多因素Cox回歸分析及繪制生存曲線，“survival ROC”包繪制ROC曲線，“pheatmap”包繪制風(fēng)險(xiǎn)相關(guān)性熱圖，“clusterProfiler”包行GO及KEGG分析；根據(jù)多因素分析結(jié)果，用“rms”包繪制列線圖及校準(zhǔn)圖。觀察的主要結(jié)局是OS，定義為從確診到任何原因死亡的時(shí)間或末次隨訪時(shí)間。

圖1 流程圖Fig.1 Study flow diagram

2 結(jié)果

2.1 基因預(yù)測模型的構(gòu)建

對GSE10846、GSE32918數(shù)據(jù)集進(jìn)行單因素Cox回歸分析后，在P<0.01和HR<1篩選條件下，GSE10846和GSE32918數(shù)據(jù)集中分別有2 011個基因和126個基因被確定為候選保護(hù)基因。在P<0.01和HR>1篩選條件下，GSE10846數(shù)據(jù)集有 2 352個基因、GSE32918數(shù)據(jù)集有102個基因被確定為候選危險(xiǎn)基因。對候選保護(hù)基因及候選危險(xiǎn)基因取交集后共保留27個基因。以GSE10846數(shù)據(jù)集為訓(xùn)練集，通過Lasso回歸對變量進(jìn)行篩選，保留13個基因，見圖2。通過多因素Cox回歸識別出 6個基因（LNPEP、SNX20、GTPBP10、CALR、BDH1、C5orf30）作為 DLBCL 患者的預(yù)后標(biāo)志物，見表1。

圖2 Lasso回歸確定模型所需基因數(shù)量Fig.2 Determination of gene numbers in the model through Lasso regression analysis

表1 基于GSE10846數(shù)據(jù)集構(gòu)建的6基因風(fēng)險(xiǎn)模型Tab.1 Construction of 6-gene based risk model from GSE10846 dataset

2.2 基因預(yù)測模型的驗(yàn)證

將預(yù)后基因的表達(dá)值與其回歸系數(shù)加權(quán)得出風(fēng)險(xiǎn)評分，根據(jù)風(fēng)險(xiǎn)評分的中位數(shù)將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。生存曲線結(jié)果顯示，高風(fēng)險(xiǎn)組患者的OS短于低風(fēng)險(xiǎn)組。在訓(xùn)練集GSE10846中，6基因模型預(yù)測3年總生存率的AUC為0.722，在驗(yàn)證集GSE32918、NCICCR中分別為0.758、0.693。見圖3。

圖3 訓(xùn)練集及驗(yàn)證集中6基因模型的Kaplan-Meier曲線和ROC曲線Fig.3 Kaplan-Meier and ROC curves for the 6-gene based model in the training and testing cohorts

在該基因模型與臨床相關(guān)性熱圖中（去除臨床信息不完整的患者后，共保留344例患者），發(fā)現(xiàn)保護(hù)基因LNPEP、SNX20、GTPBP10在低風(fēng)險(xiǎn)組中高表達(dá)，在高風(fēng)險(xiǎn)組中低表達(dá)。危險(xiǎn)基因CALR、BDH1、C5orf30在低風(fēng)險(xiǎn)組中低表達(dá)，在高風(fēng)險(xiǎn)組中高表達(dá)；驗(yàn)證隊(duì)列中亦觀察到相似的結(jié)果。風(fēng)險(xiǎn)評分與分期、治療方案、細(xì)胞來源、年齡及生存狀態(tài)有關(guān)，見圖4。ROC曲線結(jié)果顯示，與其他臨床指標(biāo)相比，風(fēng)險(xiǎn)評分的AUC最高，且在驗(yàn)證集中AUC高于IPI評分（0.706 vs 0.674），見圖5。

2.3 6基因風(fēng)險(xiǎn)評分是DLBCL患者的獨(dú)立預(yù)后因素

多因素Cox回歸分析結(jié)果顯示，在3個獨(dú)立的數(shù)據(jù)集中，6基因風(fēng)險(xiǎn)評分是預(yù)測患者總生存的獨(dú)立預(yù)后因素，見表2。

圖4 GSE10846數(shù)據(jù)集的風(fēng)險(xiǎn)評分與臨床特征的相關(guān)性熱圖Fig.4 Heat map of risk score and clinical relevance in the GSE10846 dataset

2.4 分層分析

對年齡、性別、分期和細(xì)胞起源、ECOG等臨床特征進(jìn)行分層分析，結(jié)果顯示，高風(fēng)險(xiǎn)組患者OS均較低風(fēng)險(xiǎn)組短（均P<0.05），見圖6。

2.5 列線圖模型的構(gòu)建及評估

在多因素Cox回歸的基礎(chǔ)上，將風(fēng)險(xiǎn)評分與臨床因素相結(jié)合構(gòu)建列線圖模型。結(jié)果顯示列線圖模型在GSE10846數(shù)據(jù)集中的AUC為0.796，均高于任何單一因素，見圖7。校準(zhǔn)曲線顯示，在列線圖預(yù)測生存概率和實(shí)際觀測的生存概率之間具有較好的一致性，見圖8。

2.6 功能分析

Pearson相關(guān)性檢驗(yàn)識別出與風(fēng)險(xiǎn)評分有關(guān)的1 123個基因（r>0.4且P<0.001），其中629個基因與風(fēng)險(xiǎn)評分呈正相關(guān)，494個基因與風(fēng)險(xiǎn)評分呈負(fù)相關(guān)。GO與KEGG分析結(jié)果顯示，以上基因主要富集于DNA復(fù)制和修復(fù)、蛋白加工、細(xì)胞周期、病毒致癌機(jī)制相關(guān)等生物功能及通路上，見圖9。

3 討論

DLBCL作為高度異質(zhì)性腫瘤，在治療早期發(fā)現(xiàn)高危患者并及時(shí)調(diào)整治療策略對延長生存期至關(guān)重要［13-14］，然而現(xiàn)有的預(yù)后評分系統(tǒng)盡管具有很好的臨床價(jià)值，但對部分患者仍不能發(fā)揮很好的識別作用。近年來，隨著生物信息學(xué)的發(fā)展，DLBCL的基因表達(dá)譜也已被應(yīng)用于預(yù)測腫瘤特征和預(yù)后。ZHOU等［15］研究確認(rèn)了一組由17個lncRNA組成的可區(qū)別GCB和ABC亞型的生物標(biāo)志物。然而，目前直接用于預(yù)測DLBCL患者總體生存率的基因預(yù)測模型的研究較少。

圖5 多指標(biāo)ROC曲線Fig.5 Multi-index ROC curves

表2 影響DLBCL患者預(yù)后因素的單因素和多因素Cox回歸分析Tab.2 Univariable and multivariable Cox regression analyses of prognostic factors affecting DLBCL patients

圖6 不同細(xì)胞起源、年齡、性別、ECOG評分、分期患者的Kaplan-Meier生存曲線Fig.6 Kaplan-Meier survival curves for patients in different subgroups,stratified by cell origin,age,gender,ECOG scores and stage

圖7 預(yù)測DLBCL患者生存的列線圖Fig.7 Nomogram of survival prediction in patients with DLBCL

圖8 訓(xùn)練隊(duì)列患者3年和5年生存的校準(zhǔn)曲線Fig.8 The calibration curves for predicting survival at 3 years and 5 years in the training cohort

圖9 與風(fēng)險(xiǎn)評分相關(guān)基因的GO及KEGG富集分析Fig.9 Analysis of GO and KEGG enrichment of genes related to risk score

本研究利用DLBCL患者的基因表達(dá)譜數(shù)據(jù)構(gòu)建并驗(yàn)證了一種新的6基因預(yù)后模型?；驍?shù)據(jù)是來源于2個獨(dú)立的數(shù)據(jù)集，較以單個數(shù)據(jù)集為基礎(chǔ)所構(gòu)建的預(yù)后模型更穩(wěn)健。此外，本次構(gòu)建的模型僅以6個基因?yàn)榛A(chǔ)，可減少檢測工作量及費(fèi)用，并從分子生物學(xué)層面彌補(bǔ)了臨床指標(biāo)的不足。將其與年齡、亞型、治療方案、ECOG、分期、結(jié)外部位數(shù)量等臨床因素結(jié)合后發(fā)現(xiàn)能更準(zhǔn)確地評估患者的預(yù)后，進(jìn)一步識別出被遺漏的高?；颊咭约澳承┛赡鼙桓吖赖牡臀；颊撸瑫r(shí)列線圖的AUC及校正曲線說明其具有較好的預(yù)測能力。但由于數(shù)據(jù)集缺乏完整的臨床數(shù)據(jù)，尚不能完成外部驗(yàn)證，仍需在今后的臨床工作中進(jìn)一步驗(yàn)證其價(jià)值。

本研究通過多因素Cox回歸識別出可作為DLBCL患者的預(yù)后標(biāo)志物 6個基因?yàn)?LNPEP、SNX20、GTPBP10、CALR、BDH1、C5orf30，其影響患者預(yù)后的可能原因如下：研究表明生酮療法通過減少葡萄糖供給，選擇性切斷腫瘤細(xì)胞的能量供應(yīng)，并提高脂肪供能比例以抑制腫瘤細(xì)胞增殖，而BDH1正是酮體代謝的關(guān)鍵酶，其與BDH2基因相鄰并相互調(diào)節(jié)，而BDH2被發(fā)現(xiàn)可作為抗凋亡因子參與急性白血病的發(fā)生、發(fā)展［16-17］。CALR基因編碼的鈣網(wǎng)蛋白是具有多種生物學(xué)功能的Ca2+結(jié)合蛋白，可參與腫瘤免疫逃逸、未折疊蛋白反應(yīng)（unfolded protein response，UPR）和Ca2+信號傳遞，在人類腫瘤發(fā)生發(fā)展中發(fā)揮重要作用［18-19］。C5orf30被證實(shí)與類風(fēng)濕性關(guān)節(jié)炎的組織破壞程度以及TNF、IL-1、IL-10等因子的表達(dá)相關(guān)，目前認(rèn)為炎癥因子與腫瘤發(fā)生發(fā)展有關(guān)，但暫無C5orf30與腫瘤的相關(guān)報(bào)道［20］。Genecards數(shù)據(jù)庫檢索結(jié)果表明LNPEP與胎盤滋養(yǎng)細(xì)胞腫瘤相關(guān)，其相關(guān)途徑包括囊泡介導(dǎo)的轉(zhuǎn)運(yùn)和Ⅰ類MHC介導(dǎo)的抗原加工和呈遞。GTPBP10在基本細(xì)胞過程，如蛋白質(zhì)合成，核轉(zhuǎn)運(yùn)，膜運(yùn)輸和信號轉(zhuǎn)導(dǎo)的調(diào)節(jié)中起關(guān)鍵作用［21］。本研究通過Pearson相關(guān)性檢驗(yàn)識別出與風(fēng)險(xiǎn)評分相關(guān)的1 123個基因，發(fā)現(xiàn)這些基因主要富集在DNA復(fù)制和修復(fù)、蛋白加工、細(xì)胞周期、病毒致癌等生物功能及通路上，這從側(cè)面反映出模型基因的功能作用。

本研究存在以下局限性：⑴用于構(gòu)建模型的GSE10846及GSE32918數(shù)據(jù)集獲得基因表達(dá)的測序平臺不同，分別是Affymetrix Human Genome U133 Plus 2.0 Array及 Illumina HumanRef-8 WG-DASL v3.0，其探針數(shù)目及所對應(yīng)的基因表達(dá)水平有一定區(qū)別，可能存在一定偏倚。⑵雖然本研究發(fā)現(xiàn)在驗(yàn)證隊(duì)列中，風(fēng)險(xiǎn)評分模型的預(yù)測效能高于IPI評分，但由于各個數(shù)據(jù)集中包含的臨床信息不完整，無法在所有隊(duì)列中進(jìn)行驗(yàn)證。⑶本研究缺少在病理組織標(biāo)本中進(jìn)行相關(guān)基因表達(dá)情況的檢測。

綜上所述，本研究成功構(gòu)建了一個具有較好預(yù)測能力的6基因風(fēng)險(xiǎn)評分預(yù)測模型，與年齡、亞型、治療方案、ECOG、分期、結(jié)外部位數(shù)量等臨床因素構(gòu)建的列線圖模型，有助于更全面地評估患者的預(yù)后，各基因可能是預(yù)后預(yù)測的有效指標(biāo)。