結腸癌免疫相關基因篩選和預后模型建立

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

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

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

1 材料與方法

1.1 數據來源

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

1.2 免疫細胞和基質細胞含量、免疫細胞浸潤分析

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

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

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

1.4 預后相關免疫基因(PIGs)篩選

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

1.5 富集分析

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

1.6 預后模型的構建與驗證

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

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

2 結 果

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

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

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

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

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

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

2.3 構建和驗證免疫基因相關預后模型

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

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

表1 GEO GSE39582數據集556例CC患者的臨床特征

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

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

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

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

2.5 列線圖的構建與驗證

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

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

3 討 論

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

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

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

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

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

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