結(jié)腸癌免疫相關(guān)基因篩選和預(yù)后模型建立

齊榮暄, 劉格良, 穆俊芳, 李晨龍,李 淵, 賀培鳳, 于 琦

(山西醫(yī)科大學, 1. 管理學院, 3. 基礎(chǔ)醫(yī)學院, 山西 太原, 030001;2. 山西省臨床決策研究大數(shù)據(jù)重點實驗室, 山西 太原, 030001)

結(jié)腸癌(CC)是臨床常見的消化道惡性腫瘤,其發(fā)病率、致死率在癌癥中均位居前列[1]。早發(fā)現(xiàn)、早診斷和早治療可顯著改善CC患者的預(yù)后,但多數(shù)患者確診CC時已達晚期,并不適合接受手術(shù)治療,生存質(zhì)量受到嚴重影響[2]。目前, TNM分期仍然是臨床實踐中預(yù)后預(yù)測的主要依據(jù),但其局限性日益凸顯,TNM分期相同的CC患者在治療方法及預(yù)后方面存在顯著差異[3]。這是因為TNM分期僅局限于腫瘤病理特征這一因素,而忽略了腫瘤微環(huán)境(TME)和免疫成分的影響,預(yù)測結(jié)果可能不夠精確。近年來,免疫微環(huán)境在腫瘤發(fā)生、發(fā)展、侵襲和轉(zhuǎn)移過程中的作用機制不斷被發(fā)現(xiàn)[4]。研究[5-6]表明,在CC的發(fā)生和發(fā)展過程中,免疫相關(guān)基因通過調(diào)節(jié)炎癥反應(yīng)或影響腫瘤免疫監(jiān)視發(fā)揮著重要作用。以免疫相關(guān)基因構(gòu)建的預(yù)后模型已被證明在胃癌、膀胱癌中具有較好的預(yù)后預(yù)測潛力[7-8], 然而目前尚無基于免疫相關(guān)基因的預(yù)后模型用于評估CC的TME并預(yù)測CC患者的預(yù)后。本研究根據(jù)CC的免疫分型,分析腫瘤免疫間的相互作用,利用Lasso-Cox方法篩選出免疫相關(guān)基因作為生物標志物,構(gòu)建免疫相關(guān)預(yù)后模型,并建立可用于預(yù)測CC患者OS的列線圖,以期幫助醫(yī)生對CC患者做出個體化治療決策。

1 材料與方法

1.1 數(shù)據(jù)來源

從癌癥基因組圖譜(TCGA)數(shù)據(jù)庫選擇TCGA-COAD下載標準化基因表達的RNA-seq數(shù)據(jù)、臨床數(shù)據(jù)和體細胞突變數(shù)據(jù)。根據(jù)體細胞突變數(shù)據(jù)計算每個樣本每百萬堿基中被檢測出的體細胞基因編碼錯誤、堿基替換、基因插入或缺失錯誤的總數(shù),定義為腫瘤突變負荷(TMB)。從癌癥免疫組圖譜(TCIA)數(shù)據(jù)庫選擇COAD下載微衛(wèi)星不穩(wěn)定性(MSI)數(shù)據(jù)。從高通量基因表達(GEO)數(shù)據(jù)庫下載歸一化基因表達數(shù)據(jù)和臨床數(shù)據(jù),數(shù)據(jù)集登錄號為GSE39582。

1.2 免疫細胞和基質(zhì)細胞含量、免疫細胞浸潤分析

根據(jù)基因表達數(shù)據(jù)并利用R包"estimate"估計樣本免疫細胞、基質(zhì)細胞的含量水平分值,以及樣本純度數(shù)值。利用CIBERSORTx(https://cibersortx.stanford.edu/)工具,基于每個樣本基因表達數(shù)據(jù)估算免疫細胞浸潤的豐度。

1.3 單樣本基因集富集分析和層次聚類分析

本研究采用R包(GSVA、GSEABase和limma)的ssGSEA算法,基于人類分子特征數(shù)據(jù)庫(MSigDB)中29個免疫基因集,綜合評估本研究中每個樣本的免疫學特征。將原始值x通過min-max標準化映射成在區(qū)間[0, 1]中的值x′。然后,利用R軟件中基于dist方法的歐氏距離和基于hclust方法的Ward′s linkage進行層次聚類分析,確定CC的亞型。通過R包"Rtsne"的tSNE算法來驗證CC亞型的準確性。

1.4 預(yù)后相關(guān)免疫基因(PIGs)篩選

基于免疫分型將樣本分為高免疫組和低免疫組,采用R包"limma"篩選差異表達基因(DEGs), 篩選條件為|log2FC|>0.585(FC為差異倍數(shù))、錯誤發(fā)現(xiàn)率(FDR)<0.05。從ImmPort數(shù)據(jù)庫下載免疫相關(guān)基因列表,將列表基因與DEGs取交集,得到差異表達免疫基因(DEIGs)。通過對DEIGs進行單因素Cox回歸分析,得到PIGs(P<0.05)。

1.5 富集分析

利用基因集富集分析(GSEA)探討涉及DEGs的關(guān)鍵信號通路,該分析可分別展示京都基因與基因組百科全書(KEGG)和Hallmark基因組中顯著富集的通路(FDR<0.01)。

1.6 預(yù)后模型的構(gòu)建與驗證

利用Lasso懲罰的Cox回歸分析篩選PIGs, 構(gòu)建CC的免疫相關(guān)預(yù)后模型,參與預(yù)后模型構(gòu)建的PIGs被稱為關(guān)鍵預(yù)后相關(guān)免疫基因(KIGs), 風險評分的計算公式如下:

式中,Expri表示基因i的表達水平,Coefi表示模型中基因i的回歸系數(shù)。基于得到的風險評分計算中位數(shù),根據(jù)中位數(shù)將所有患者樣本分為高風險組、低風險組,并采用R包"survival" "survminer"進行Kaplan-Meier生存分析。采用R包"timeROC"生成具有時間依賴性的受試者工作特征(ROC)曲線,并計算1、3、5年總生存率(OS)的曲線下面積(AUC)來驗證預(yù)后模型的預(yù)測能力。風險評分與臨床特征、免疫浸潤細胞、TMB、MSI的相關(guān)性采用Spearman相關(guān)分析法進行分析,2組比較采用Wilcoxon秩和檢驗,多組比較采用Kruskal-Wallis秩和檢驗。使用R包"maftools"工具獲取、分析和可視化CC樣本基因突變瀑布圖。對風險評分和其他臨床因素進行單因素和多因素Cox回歸分析,確定獨立預(yù)后因素,基于上述協(xié)變量通過R包"rms" "nomogramEx" "regplot"建立列線圖,并采用ROC曲線和校準曲線驗證建立的列線圖是否適用于臨床。P<0.05為差異有統(tǒng)計學意義。

2 結(jié) 果

2.1 基于免疫基因組分析鑒定CC的2種亞型

基于29個免疫基因集,采用ssGSEA算法和層次聚類算法分析來自TCGA-COAD隊列的479個腫瘤樣本,將樣本分為高免疫組和低免疫組(圖1A、圖1B)。R包"estimate"評估結(jié)果顯示,與低免疫組相比,高免疫組基質(zhì)評分、免疫評分和綜合評分水平較高,差異有統(tǒng)計學意義(P<0.001), 見圖1C。免疫細胞浸潤分析結(jié)果顯示,高免疫組與低免疫組的初始B細胞、CD8+T細胞、靜止的CD4+記憶T細胞、記憶激活CD4 T細胞、濾泡輔助性T細胞、單核細胞、M1型巨噬細胞、M2型巨噬細胞、中性粒細胞存在顯著差異(圖1D)。進一步采用tSNE算法對CC患者進行免疫水平聚類驗證,結(jié)果證明了聚類算法的可靠性(圖1E)。本研究比較高免疫和低免疫這2種亞型中的HLA基因表達情況,發(fā)現(xiàn)24種HLA基因在高免疫亞型中的表達水平均高于低免疫亞型,差異有統(tǒng)計學意義(P<0.05), 見圖1F。

A: 樣本分為高免疫和低免疫這2種亞型; B: TCGA-COAD中的腫瘤微環(huán)境熱圖和免疫特征; C: 高免疫組和低免疫組免疫評分、基質(zhì)評分、綜合評分比較; D: 高免疫組和低免疫組的免疫細胞浸潤分析; E: tSNE算法驗證分型結(jié)果; F: 高免疫與低免疫亞型的HLA基因表達差異。圖1 層次聚類分析結(jié)果

2.2 DEGs、DEIGs和PIGs的鑒定及功能富集分析

對高免疫組和低免疫組進行基因差異表達分析,得到1 439個DEGs, 其中1 074個基因表達上調(diào),365個基因表達下調(diào)。將ImmPort數(shù)據(jù)庫中免疫基因與DEGs取交集,得到379個DEIGs, 其中365個基因表達上調(diào), 14個基因表達下調(diào)(圖2A～圖2D)。采用單因素Cox回歸分析篩選出15個PIGs, 見圖2E。為了深入了解免疫基因表達影響CC生物學過程的整體模式,本研究通過GSEA確定DEGs涉及的通路。在Hallmark基因集中, DEGs主要富集在免疫反應(yīng)、信號傳導、細胞生長與死亡等相關(guān)通路,見圖2F。KEGG通路中, DEGs主要富集在免疫系統(tǒng)、信號分子與相互作用信號傳導等相關(guān)通路,見圖2G。

A: DEGs與免疫基因的交集; B: DEGs熱圖; C: DEIGs熱圖; D: DEGs火山圖; E: 基于單因素Cox風險回歸分析的PIGs及其風險比; F: DEGs的Hallmark基因集富集分析; G: DEGs的KEGG基因集富集分析。圖2 DEGs、DEIGs和PIGs的鑒定及富集分析

2.3 構(gòu)建和驗證免疫基因相關(guān)預(yù)后模型

將PIGs作為預(yù)選基因進行Lasso-Cox回歸分析,篩選出12個KIGs用于構(gòu)建預(yù)后模型,見圖3A～B。下載GEO數(shù)據(jù)庫中556例CC患者的樣本數(shù)據(jù)并篩選相應(yīng)臨床數(shù)據(jù)(表1), 作為驗證集驗證預(yù)后模型的準確性。分別計算TCGA-COAD、GEO GSE39582數(shù)據(jù)集中每例CC患者的風險評分,并根據(jù)風險評分中位數(shù)將患者分為高風險組和低風險組(圖3C～D)。與低風險組相比,高風險組患者的生存時間較短,病死率較高,見圖3E～F。Kaplan-Meier曲線分析顯示,低風險組的OS高于高風險組,差異有統(tǒng)計學意義(P<0.05), 見圖3G～H。ROC曲線顯示,預(yù)后模型預(yù)測TCGA-COAD中患者1、3、5年OS的AUC分別為0.735、0.725、0.699, 預(yù)后模型預(yù)測驗證集中患者1、3、5年OS的AUC分別為0.615、0.607、0.584, 表明預(yù)后模型有著良好的預(yù)測能力和準確性,見圖3I～J。

A: 在Lasso模型中調(diào)整參數(shù)進行10倍交叉驗證(2條虛線之間的區(qū)域表示合理值); B: 12個KIGs在Lasso模型中的coef系數(shù)剖面圖; C、D: TCGA-COAD、GEO GSE39582數(shù)據(jù)集中CC患者風險評分的分布和中位值; E、F: TCGA-COAD、GEO GSE39582數(shù)據(jù)集中CC患者生存狀態(tài)、生存時間和風險評分的分布情況; G、H: TCGA-COAD、GEO GSE39582數(shù)據(jù)集中不同生存風險患者的Kaplan-Meier生存曲線; I、J: 預(yù)后模型在TCGA-COAD、GEO GSE39582數(shù)據(jù)集中預(yù)測1、3、5年OS的ROC曲線。圖3 預(yù)后模型的構(gòu)建和驗證

表1 GEO GSE39582數(shù)據(jù)集556例CC患者的臨床特征

2.4 風險評分與臨床特征、免疫浸潤細胞、TMB、MSI的相關(guān)性

基于TCGA-COAD數(shù)據(jù),對風險評分與臨床特征、免疫浸潤細胞、TMB、MSI的相關(guān)性進行檢驗,患者的臨床特征見表2。不同疾病分期、不同TNM分期患者的風險評分比較,差異有統(tǒng)計學意義(P<0.05), 見圖4A～D。KIGs與眾多免疫細胞均存在相關(guān)性(P<0.05), 見圖4E。風險評分與TMB呈正相關(guān)(r=0.16,P=0.002 1), 見圖4F; 高風險組和低風險組患者的TMB比較,差異有統(tǒng)計學意義(P<0.05), 見圖4G; 低TMB組OS高于高TMB組,差異有統(tǒng)計學意義(P<0.05), 見圖4H; 瀑布圖顯示,腫瘤驅(qū)動突變基因腺瘤性息肉樣腺癌基因(APC)突變頻率最高,見圖4I。高風險組較低風險組患者MSI更強,見圖4J; 進一步比較微衛(wèi)星穩(wěn)定狀態(tài)(MSS)、微衛(wèi)星低不穩(wěn)定性(MSI-L)和微衛(wèi)星高不穩(wěn)定性(MSI-H)的風險評分差異發(fā)現(xiàn), MSS與MSI-H(P=2e-05)、MSI-L與MSI-H(P=0.049)的風險評分差異均有統(tǒng)計學意義,見圖4K。

表2 TCGA-COAD中446例CC患者的臨床特征

A～D: 風險評分與疾病分期、TNM分期的相關(guān)性(不同分期比較所得數(shù)據(jù)為P值); E: KIGs與免疫浸潤細胞的相關(guān)性; F: 風險評分與TMB的相關(guān)性; G: 高風險組與低風險組TMB比較; H: 高風險組與低風險組Kaplan-Meier生存曲線分析; I: 前30個突變基因瀑布圖; J: 高風險組與低風險組MSI比較; K: 不同微衛(wèi)星穩(wěn)定性狀態(tài)的風險評分比較(不同狀態(tài)比較所得數(shù)據(jù)為P值)。圖4 風險評分與臨床特征、KIGs、TMB、MSI的相關(guān)性分析

2.5 列線圖的構(gòu)建與驗證

將風險評分、年齡、性別和疾病分期作為協(xié)變量進行單因素和多因素Cox回歸分析,結(jié)果顯示,風險評分是CC患者的獨立預(yù)后因素(P<0.001), 見圖5A～B。結(jié)合上述因素,構(gòu)建列線圖,進一步發(fā)揮預(yù)后模型的臨床可用性,見圖5C。ROC曲線顯示,列線圖對生存狀態(tài)具有良好的預(yù)測能力和準確性,預(yù)測1、3、5年OS的AUC分別為0.821、0.817、0.778, 見圖5D。校準圖顯示,列線圖的性能與理想模型相似,見圖5E。

A: 單因素Cox分析; B: 多因素Cox分析; C: 基于風險評分和臨床特征構(gòu)建列線圖; D: 列線圖預(yù)測1、3、5年OS的ROC曲線; E: 列線圖的校準圖。圖5 列線圖的建立和驗證

3 討 論

根據(jù)國際癌癥研究機構(gòu)2020年發(fā)布的數(shù)據(jù), CC是世界范圍內(nèi)第3位常見的惡性腫瘤,且是惡性腫瘤相關(guān)死亡的第2位原因,已成為嚴重危害人類身體健康的重要公共衛(wèi)生問題[1]。研究[9]表明, TME中的免疫細胞在腫瘤進展中起著重要的作用。即使腫瘤侵襲性較小,如果獲得性免疫反應(yīng)較弱,預(yù)后也會變差; 相反,無論腫瘤的局部范圍和侵襲程度如何,高密度的獲得性免疫細胞都預(yù)示著良好的預(yù)后[10]。因此,醫(yī)生對患者進行診斷時必須對TME細胞的數(shù)量和類型進行明確,并在治療后盡可能量化,以有效地指導治療。轉(zhuǎn)錄組學分析可利用少量的RNA樣本進行研究,近年來研究者們開始關(guān)注基于轉(zhuǎn)錄組學的基因表達在腫瘤免疫治療中的影響。本研究篩選出CC的KIGs用于構(gòu)建預(yù)后模型,并建立可預(yù)測CC患者OS的列線圖,這可能有助于CC患者的預(yù)后風險預(yù)測、腫瘤分期預(yù)測和免疫治療。

本研究利用ssGSEA算法將CC分為高免疫和低免疫亞型,高免疫亞型的免疫細胞浸潤程度更高,HLA基因表達水平更高,表明其相較于低免疫亞型具有更強的免疫原性,與既往研究[11]結(jié)果一致。DEGs不僅在免疫相關(guān)信號通路中富集,還在許多癌相關(guān)通路中富集,包括凋亡、胰腺癌、JAK/STAT信號通路、MAPK信號通路和p53信號通路,與既往研究[12]結(jié)論(不同的免疫信號與JAK/STAT信號通路呈正相關(guān))相符。值得注意的是,本研究結(jié)果揭示了CC中通路活性與免疫活性之間可能存在聯(lián)系。

本研究基于篩選出的12個KIGs構(gòu)建預(yù)后模型,結(jié)果顯示,高風險評分患者和低風險評分患者的OS曲線明顯分離,且隨著風險評分的升高,疾病分期、TNM分期越高,表明該預(yù)后模型對腫瘤分期具有一定預(yù)測能力。本研究篩選出的12個KIGs中,大多數(shù)與腫瘤的發(fā)生發(fā)展存在關(guān)系,且在CC中也發(fā)揮著非常重要的作用。ULBP2編碼一種主要的組織相容性復合體Ⅰ類相關(guān)分子,該分子與自然殺傷細胞(NK細胞)上的NKG2D受體結(jié)合,觸發(fā)多種細胞因子和趨化因子的釋放,進而促進NK細胞的募集和活化,并破壞靶細胞[13]。ULBP2是維持B7-H3賦予CC細胞對V-δ-2T細胞毒性抗性的必要基因[14]。ANGPTL4可通過增強葡萄糖攝取促進具核梭桿菌定植,進而促進結(jié)直腸癌細胞的糖酵解活性,對結(jié)直腸癌的增殖和侵襲具有促進作用[15]。ANGPTL4可作為前列腺素E2的下游效應(yīng)分子,可促進腫瘤細胞對低氧的適應(yīng),進而促進結(jié)直腸癌發(fā)展[16]。CX3CL1配體與其受體CX3CR1通過將NK細胞和T細胞等抗腫瘤免疫細胞招募到TME而控制腫瘤的生長,從而發(fā)揮抗腫瘤作用。但另一方面,CX3CL1-CX3CR1軸也激活了促腫瘤反應(yīng)[17]。多項證據(jù)[18-19]表明,CX3CL1與結(jié)直腸癌顯著相關(guān),可導致預(yù)后不良。FABP4主要在脂肪組織和巨噬細胞中表達,在肥胖相關(guān)疾病中功能顯著,其通過與多種經(jīng)典途徑(如PI3K/AKT和STAT3/ALDH1信號傳導)相互作用,在腫瘤發(fā)生和侵襲中起著關(guān)鍵作用[20]。FABP4在CC中過表達會增加脂肪酸轉(zhuǎn)運,增強能量和脂質(zhì)代謝,激活A(yù)KT通路和上皮間質(zhì)轉(zhuǎn)化,進而促進CC細胞的遷移和侵襲[21]。

HE X F等[22]發(fā)現(xiàn),BST2通過募集腫瘤相關(guān)巨噬細胞(TAMs)并將其誘導為M2表型而增加TAMs的浸潤率,這有助于形成免疫抑制的TME, 進而促進結(jié)直腸癌發(fā)展。RBP7編碼的蛋白質(zhì)是細胞視黃醇結(jié)合蛋白家族的成員[23],RBP7可作為腫瘤微環(huán)境調(diào)節(jié)因子誘導5-氟尿嘧啶耐藥,從而影響結(jié)直腸癌患者的預(yù)后[24]。另有研究[25]表明,RBP7高表達是早期和晚期CC患者特異性生存率低的獨立生物標志物,在功能上有助于CC細胞的惡性表型。TGFb1屬于轉(zhuǎn)化生長因子-β(TGF-β)家族中的一員, TGF-β1信號在腫瘤細胞中具有雙重性作用,既在癌癥前期細胞中發(fā)揮有效的細胞抑制和促凋亡活性,也有利于惡性轉(zhuǎn)化后期的上皮間質(zhì)轉(zhuǎn)化和轉(zhuǎn)移。此外,TGFb1還可通過各種通路影響癌癥相關(guān)成纖維細胞、內(nèi)皮細胞以及腫瘤免疫浸潤細胞的活性,進而影響腫瘤的發(fā)展[26], 這些通路在結(jié)直腸癌中同樣被驗證[27-29]。APOBEC3F是胞苷脫氨酶基因家族的成員[30], 是在基因簇中發(fā)現(xiàn)的7個相關(guān)基因或假基因之一,這些基因編碼的蛋白可能是RNA編輯酶,在生長或細胞周期控制中起作用。研究[31]表明,APOBEC3F存在意味著成纖維細胞具有強烈的促炎活性,可促進結(jié)直腸癌的發(fā)展。目前,IL1RL2、CHGB在癌癥領(lǐng)域中的相關(guān)研究尚較少,未來可進一步深入分析。總之,這些KIGs從免疫應(yīng)答、腫瘤相關(guān)免疫細胞活性、免疫細胞浸潤、能量和脂質(zhì)代謝、腫瘤細胞的耐氧特性等方面影響腫瘤細胞的增殖、遷移和浸潤。

微衛(wèi)星通常由1～6個短串聯(lián)重復DNA序列組成,MSI被認為是結(jié)直腸癌的主要致癌途徑之一,這種不穩(wěn)定性是由缺陷DNA錯配修復機制引起的,可導致體細胞突變,增加腫瘤突變負擔[32]。研究[33]表明,相較于MSS、MSI-L的CC患者, MSI-H的CC患者接受免疫檢查點抑制劑治療可獲得更好的治療效果。本研究發(fā)現(xiàn),高風險組患者的MSI狀態(tài)更高,表明免疫檢查點抑制劑治療在高風險組患者中可能獲益更大。本研究中,低TMB組CC患者的預(yù)后與高TMB組更好,與既往研究[34]結(jié)論一致。由此提示, TMB或可作為一個獨立的生物標志物,指導更有效的免疫治療策略,并改善預(yù)后。本研究還發(fā)現(xiàn),高風險組與低風險組的TMB存在顯著性差異, TMB與風險評分呈正相關(guān),由于高體細胞突變意味著MSI高,具有更強的免疫治療反應(yīng),更高的TMB與更好的免疫療法療效之間的相關(guān)性已得到證實[35]。

綜上所述,本研究構(gòu)建了基于12個KIGs的免疫相關(guān)預(yù)后模型,并建立了可用于預(yù)測CC患者OS的列線圖,這可能有助于CC患者總體生存情況的預(yù)測和個性化治療方案的優(yōu)化。本研究亦存在一些局限性: 本研究僅基于TCGA公開數(shù)據(jù)集進行分析,未來可結(jié)合其他數(shù)據(jù)庫的數(shù)據(jù)進一步深入研究; 本研究篩選出的12個KIGs與CC患者預(yù)后顯著相關(guān),但僅通過數(shù)據(jù)挖掘方法進行分析,這些基因的功能和機制有待進一步開展實驗研究加以闡明。