免疫基因和代謝基因在胃癌患者預(yù)后分析中的比較研究

周 霞,江紹萍,葉明亮,喻文霞,鐘 琦

(云南民族大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院, 云南 昆明 650500)

在惡性腫瘤中, 胃癌發(fā)病率與死亡率分別排行第4、第3[1], 其診治和預(yù)后評(píng)估一直是專家學(xué)者和臨床醫(yī)生關(guān)注的重點(diǎn)之一. 胃癌起源于胃黏膜上皮, 其癌變的發(fā)展由多種因素, 多個(gè)階段造成, 是一種目前還不能完全治愈的疾病, 且患者多數(shù)為50歲以上. 我國(guó)胃癌發(fā)病率極高, 胃癌患者約占世界胃癌患者一半, 縱然近些年來(lái)胃癌的診治取得了一些良好的進(jìn)展, 但胃癌患者的預(yù)后效果普遍較差[2].

現(xiàn)在越來(lái)越多的證據(jù)表明, 免疫基因在胃癌的發(fā)生和發(fā)展中起著關(guān)鍵作用[3], 如章婧文[2]運(yùn)用CIBERSORT算法解析胃癌中22種免疫細(xì)胞的構(gòu)成, 并篩選免疫細(xì)胞去建立免疫風(fēng)險(xiǎn)評(píng)分模型和個(gè)體化列線圖模型以預(yù)測(cè)胃癌患者的死亡風(fēng)險(xiǎn). 王堯等[4]采用ESTIMATE及CIBERSORT分別評(píng)估胃癌樣本中免疫微環(huán)境及免疫細(xì)胞的組成, 總結(jié)出胃癌免疫分型有助于判斷患者的預(yù)后, 且高免疫組患者具有較高的免疫原性及抗腫瘤免疫活性. 張靜等[5]通過(guò)生物信息學(xué)方法分析發(fā)現(xiàn)胃癌Her-2基因高表達(dá)可能會(huì)對(duì)免疫系統(tǒng)產(chǎn)生一定的抑制作用, 從而影響免疫系統(tǒng)發(fā)揮抗腫瘤的作用. 這些都表明了免疫基因能夠?qū)ξ赴┗颊弋a(chǎn)生直接或間接的影響.

胃癌的發(fā)生發(fā)展也與代謝基因的變化密切相關(guān). 陶柱萍等[6]基于代謝基因發(fā)現(xiàn)了多種胃癌的內(nèi)源性代謝生物標(biāo)志物, 總結(jié)歸納出代謝組學(xué)技術(shù)在胃癌早期診斷及預(yù)后方面的研究進(jìn)展，為胃癌的研究提供了1個(gè)參考. 向麗娟等[7]根據(jù)轉(zhuǎn)錄組測(cè)序, 篩選出胃癌與相鄰非癌樣本之間的所有基因表達(dá)量, 進(jìn)行富集分析, 認(rèn)為脂代謝中的某些基因與胃癌的相關(guān)通路有關(guān)聯(lián), 得到胃癌脂代謝中關(guān)鍵基因的異常表達(dá)可為尋找胃癌的診斷和預(yù)后指標(biāo)提供線索. 趙麗沙等[8]從數(shù)據(jù)處理方式、取得的差異代謝物及其調(diào)節(jié)的代謝通路等方面研究了代謝組學(xué)技術(shù)在胃癌早期診斷與預(yù)后判斷的臨床應(yīng)用, 為代謝基因應(yīng)用于胃癌早期診斷及預(yù)后判斷提供了參考的依據(jù). 這些都說(shuō)明了代謝基因能夠?qū)ξ赴┗颊叩脑\斷和預(yù)后分析提供1個(gè)參考價(jià)值.

綜上可以看出免疫基因和代謝基因都可能影響胃癌患者的預(yù)后分析. 因而本研究將基于LASSO回歸構(gòu)建預(yù)后風(fēng)險(xiǎn)模型[9], 并根據(jù)特異性與靈敏度之和最大的點(diǎn)來(lái)劃分患者的風(fēng)險(xiǎn), 并從生存曲線, 生存風(fēng)險(xiǎn)分布等圖表分別對(duì)免疫基因和代謝基因進(jìn)行比較, 進(jìn)一步了解哪個(gè)基因能更好的對(duì)胃癌患者進(jìn)行預(yù)后及生存狀況分析.

1 數(shù)據(jù)來(lái)源與預(yù)處理

在癌癥基因組圖譜(the cancer genome atlas, TCGA)數(shù)據(jù)庫(kù)中獲得胃癌完整的信使RNA表達(dá)數(shù)據(jù)和臨床信息. 本次研究包含375個(gè)胃癌樣本和32個(gè)相鄰的非癌性樣本. 研究中先從這些數(shù)據(jù)提取相應(yīng)的臨床信息, 包括患者的基本資料, 腫瘤分期, 腫瘤分級(jí), TNM分類(見(jiàn)表1). 然后從Immport數(shù)據(jù)庫(kù)中下載免疫基因相關(guān)信息. 使用R語(yǔ)言的limma包提取差異表達(dá)的免疫相關(guān)基因, 于是在胃癌的 56 753 個(gè)基因中找到 1 651 個(gè)免疫相關(guān)基因, 并提取所有免疫基因的表達(dá)量. 最后, 利用基因表達(dá)匯編(gene expression omnibus, GEO)數(shù)據(jù)庫(kù)獲得具有詳細(xì)生存信息的433個(gè)胃癌樣本及提取免疫基因表達(dá)量. 把從TCGA數(shù)據(jù)庫(kù)和GEO數(shù)據(jù)庫(kù)中獲得的免疫基因矩陣取交集, 得到 23 355 組免疫基因?qū)?

表1 部分臨床樣本數(shù)據(jù)表

同樣的, 我們從下載的信使RNA表達(dá)數(shù)據(jù)和臨床信息中獲取代謝相關(guān)的基因, 最終在胃癌的 56 753 個(gè)基因中找到944個(gè)代謝相關(guān)基因, 并且把代謝基因的表達(dá)量提取出來(lái). 最后把胃癌患者的生存信息與免疫基因和代謝基因合并.

2 實(shí)驗(yàn)設(shè)計(jì)

2.1 方法

用單因素Cox分析篩選出與胃癌患者總體生存率顯著相關(guān)的免疫基因和代謝基因, 并計(jì)算出每個(gè)基因的風(fēng)險(xiǎn)比(hazard ratio, HR). 基于套索算法(least absolute shrinkage and selection operator, LASSO)構(gòu)建預(yù)后模型. 依據(jù)風(fēng)險(xiǎn)值=∑(基因系數(shù)×基因表達(dá)量)計(jì)算出每個(gè)患者的風(fēng)險(xiǎn)評(píng)分, 與獲得的最佳截?cái)嘀?cutoff)進(jìn)行比較, 繪出接收者操作特征曲線(receiver operating characteristic curve, ROC), 將患者區(qū)分為高低風(fēng)險(xiǎn)兩組. 最后, 通過(guò)乘積極限法(kaplan-meier survival estimate, K-M)對(duì)胃癌患者進(jìn)一步分析, 并采用對(duì)數(shù)秩檢驗(yàn). 繪制相關(guān)的風(fēng)險(xiǎn)曲線去分析預(yù)測(cè)胃癌患者的總體生存有效性, 計(jì)算ROC曲線下面積(area under curve, AUC), 用于分析模型的準(zhǔn)確率.

2.2 實(shí)驗(yàn)環(huán)境

所有的統(tǒng)計(jì)分析均使用R4.0.1(www.r-project.org)和Strawberry Perl (www.perl.org/)進(jìn)行,P值小于0.05為差異顯著且具有統(tǒng)計(jì)意義.

2.3 預(yù)后基因篩選

通過(guò)單因素Cox分析方法篩選出與375例胃癌患者總體生存率顯著相關(guān)的免疫基因和代謝基因, 得到與預(yù)后相關(guān)的37對(duì)免疫基因和16對(duì)代謝基因以及它們每個(gè)基因的HR值(見(jiàn)表2, 表3).

表2 部分免疫基因的HR值

表3 部分代謝基因的HR值

表3中, 若HR值大于1表示該基因參與預(yù)后的風(fēng)險(xiǎn)較高, HR值小于1表示它參與預(yù)后的風(fēng)險(xiǎn)較低; 若coxPvalue值大于0.05說(shuō)明這個(gè)基因與生存不相關(guān), coxPvalue值小于0.05說(shuō)明該基因與生存相關(guān), 可以作為預(yù)后分析基因.

2.4 基于LASSO回歸的預(yù)后模型建立與驗(yàn)證

2.4.1 LASSO回歸模型

構(gòu)建風(fēng)險(xiǎn)評(píng)分模型的最佳方法是LASSO回歸, 它經(jīng)過(guò)構(gòu)造懲罰參數(shù), 將變量的系數(shù)實(shí)行壓縮, 使某些系數(shù)變?yōu)?, 從而取得一個(gè)更為精煉的模型. 同時(shí)為了避免出現(xiàn)過(guò)度擬合現(xiàn)象, 運(yùn)用同一體系的參數(shù)來(lái)控制該模型的復(fù)雜性. LASSO回歸可通過(guò)普通線性模型增加L1懲罰項(xiàng)得到. 普通線性模型為:

Y=Xβ+ε,

(1)

(2)

當(dāng)矩陣X不滿足列滿秩時(shí), 回歸系數(shù)將不能使用OLS去求解, 在這種情況下, 引入懲罰方法, 即將最小懲罰似然函數(shù)的值作為回歸系數(shù)的估計(jì)值, 從而可以獲得LASSO回歸的模型：

(3)

2.4.2 基于免疫基因和代謝基因的胃癌預(yù)后分析的LASSO回歸模型

對(duì)預(yù)后相關(guān)的37對(duì)免疫基因使用LASSO回歸, 最終選擇了19對(duì)免疫基因用于構(gòu)建模型(見(jiàn)表4): HLA-DQB1|CXCR4, HLA-F|HSPA1A, MICB|ESM1, MICB|NR2F1, MICB|NRP1, CXCL13|CD79A, S100A9|AGT, APOD|AHNAK, APOD|DMBT1, APOD|BMP4, ITGAV|MAPK3, TLR1|CTLA4, ARRB1|STC1, TNFSF4|CXCR6, VAV3|SCG2, PTPN6|PLAUR, EDN3|DKK1, SBDS|PLXNB2, FAM3C|ESRRA.

同樣的，對(duì)預(yù)后相關(guān)的16個(gè)代謝基因使用LASSO回歸, 最終選擇了13個(gè)代謝基因用于構(gòu)建模型(見(jiàn)表5): GSTA2, POLD3, GLA, GGT5, DCK, CKMT2, ASAH1, OPLAH, ME1, ACYP1, NNMT, POLR1A, RDH12.

表4 免疫基因?qū)跋鄳?yīng)系數(shù)

表5 代謝基因及相應(yīng)系數(shù)

使用隨時(shí)間變化的ROC曲線進(jìn)行分析, 分別得到了1, 3, 5年的最佳截?cái)嘀?cutoff),即特異性與靈敏度之和最大的點(diǎn)(見(jiàn)圖1). 它們以假陽(yáng)性率(1-特異性)為橫軸, 真陽(yáng)性率即靈敏度為縱軸. 特異性(Specificity)和靈敏度(Sensitivity)的公式[10]為:

(4)

(5)

其中TP為真陽(yáng)性,Fp為假陽(yáng)性,TN為真陰性,FN為假陰性.

按照風(fēng)險(xiǎn)值=∑(基因系數(shù)×基因表達(dá)量), 分別計(jì)算每個(gè)患者的風(fēng)險(xiǎn)值, 再讓其與cutoff值進(jìn)行比較, 將患者劃分為高危組和低危組患者. 由圖可以直觀的看出, 免疫基因下1年的ROC曲線圖最平緩, 最接近左上角的坐標(biāo)軸, 這表明它有很高的診斷準(zhǔn)確性, 同時(shí)對(duì)胃癌患者的預(yù)后研究具備較高的價(jià)值.

圖1 免疫基因和代謝基因下患者在不同年份的ROC曲線

2.4.3 LASSO回歸模型的準(zhǔn)確性驗(yàn)證

計(jì)算出在免疫基因和代謝基因下患者的ROC曲線在不同時(shí)間下的AUC值(如表6). 用AUC值來(lái)驗(yàn)證模型的準(zhǔn)確性, 其值越接近1表示準(zhǔn)確性越高.發(fā)現(xiàn)在相同的年份里, 免疫基因的AUC值都比代謝基因的AUC值高, 且免疫基因的AUC值隨著年份的增加而不斷增加, 但是代謝基因的AUC值卻隨著年份的增加而減少, 說(shuō)明利用免疫基因來(lái)構(gòu)建LASSO回歸模型具有較高的準(zhǔn)確性.

表6 免疫基因和代謝基因在不同時(shí)間下的AUC值

3 結(jié)果及分析

本文采用K-M法對(duì)已經(jīng)分好的胃癌高低風(fēng)險(xiǎn)患者進(jìn)行生存分析[11]. K-M法是1種非參數(shù)方法, 可以通過(guò)觀察胃癌患者的生存時(shí)間來(lái)對(duì)患者的生存概率進(jìn)行估計(jì). 生存概率是指某段時(shí)間開(kāi)始時(shí)存活的樣本至該時(shí)間結(jié)束時(shí)仍然存活的可能性大小[12]. 生存概率可以表示為:

(6)

其中,S表示生存概率,U表示活過(guò)某段時(shí)間的人數(shù),V表示該時(shí)段初期觀察人數(shù).

通過(guò)生存概率能夠得到高低風(fēng)險(xiǎn)組隨著時(shí)間的變化而變化的生存曲線. 在生存曲線上半部分,橫軸為時(shí)間(年), 縱軸為生存率. 下半部分為在每個(gè)時(shí)間點(diǎn)上高低風(fēng)險(xiǎn)組患者存活的數(shù)目(見(jiàn)圖2和圖3).

由圖可知,隨著隨訪時(shí)間的增加, 兩種基因的高低風(fēng)險(xiǎn)患者生存率都有所下降. 在免疫基因中, 高低風(fēng)險(xiǎn)組曲線之間存在的差異較大, 高?；颊叩纳媛手饾u趨于0, 低危患者的生存率大概到達(dá)55%趨于穩(wěn)定, 且隨著時(shí)間的增加生存曲線比較平緩; 在代謝基因中, 高低風(fēng)險(xiǎn)組曲線之間存在交叉, 即差異較小, 生存率分別到達(dá)13%和35%就趨于穩(wěn)定. 5年后, 免疫基因中高低風(fēng)險(xiǎn)患者存活的數(shù)目分別為1和6; 代謝基因中高低風(fēng)險(xiǎn)患者存活的數(shù)目分別為1和2.

免疫基因的生存曲線P值為0, 代謝基因的生存曲線P值為4.927e-06, 它們都小于0.05, 經(jīng)過(guò)Log-rank檢驗(yàn), 它們的整體生存曲線間有明顯的差異, 具有統(tǒng)計(jì)學(xué)意義.

同時(shí)我們還考察了患者的風(fēng)險(xiǎn)值分布和生存狀況圖(見(jiàn)圖4). 把患者按照風(fēng)險(xiǎn)值從左到右由低到高的順序進(jìn)行排序, 并把這些患者作為橫軸, G和I以風(fēng)險(xiǎn)值為縱軸. 隨著胃癌患者往右增加, 它們的風(fēng)險(xiǎn)值也增加. 由于免疫基因的高低風(fēng)險(xiǎn)組之間差異較大, 風(fēng)險(xiǎn)值分布圖呈現(xiàn)逐步上升趨勢(shì). 而代謝基因的風(fēng)險(xiǎn)值分布圖上升較緩慢. H和J以生存時(shí)間為縱軸. 隨著胃癌患者風(fēng)險(xiǎn)值的增加, 死亡的患者越來(lái)越多. 從免疫基因的生存狀況分布圖看出，表示死亡的紅點(diǎn)由少變多較為明顯, 能更好呈現(xiàn)出胃癌患者的生存時(shí)間.

圖4 免疫基因和代謝基因的風(fēng)險(xiǎn)分布圖和生存狀況圖

總結(jié)出利用免疫基因能更好的分析胃癌患者的生存狀態(tài). 且隨著時(shí)間的增加, 高低風(fēng)險(xiǎn)組之間的差異越來(lái)越明顯. 患者風(fēng)險(xiǎn)值增加, 其風(fēng)險(xiǎn)評(píng)分增加, 死亡的病人也逐漸增多, 這與預(yù)期的結(jié)果一致.

4 結(jié)語(yǔ)

本研究通過(guò)單因素Cox分析篩選出與預(yù)后相關(guān)的免疫基因和代謝基因, 并利用LASSO回歸構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型, 找到cutoff值進(jìn)行比較, 將胃癌患者劃分為高低風(fēng)險(xiǎn)組, 進(jìn)行生存分析的一系列圖表對(duì)比研究, 結(jié)果顯示利用免疫基因能更好的分析胃癌患者的生存狀況, 即依據(jù)免疫基因在時(shí)間依賴的ROC曲線中找到了最佳截?cái)帱c(diǎn), 它能更好的將胃癌患者分為高危組和低危組, 且準(zhǔn)確性較高. 由免疫基因計(jì)算得到的生存曲線顯示, 高低風(fēng)險(xiǎn)組患者生存時(shí)間有明顯的差異, 具有統(tǒng)計(jì)學(xué)意義. 綜上所述, 免疫基因或許可成為預(yù)測(cè)胃癌患者生存時(shí)間的一個(gè)新指標(biāo), 為進(jìn)一步研究胃癌的預(yù)后和臨床輔助治療提供一個(gè)參考價(jià)值.