基于TCGA 和GEO 數(shù)據(jù)庫探索結(jié)腸癌腫瘤微環(huán)境中的免疫相關(guān)預(yù)后因子*

操利超，巴穎，丁世濤，翁琦，盧曉萍，張核子（深圳市核子基因科技有限公司，廣東深圳 518071）

結(jié)腸癌是世界上死亡率最高的惡性腫瘤之一［1］。盡管最近幾年在診斷和治療方面取得了一定的進展，但結(jié)腸癌患者的總體預(yù)后仍然很差，主要原因之一是缺乏有效的預(yù)后生物標志物［2］。因此，探討結(jié)腸癌有價值的預(yù)后標志物和治療靶點是迫切且必要的。由于腫瘤免疫微環(huán)境（tumor microenvironment，TME）的異質(zhì)性和復(fù)雜性，只有一小部分患者受益于免疫治療。最近，有學者利用生物信息學和機器學習方法發(fā)現(xiàn)多種類型的免疫相關(guān)生物標志物與結(jié)腸癌的預(yù)后相關(guān)［3］。然而，對于TME的分子特征需要進一步研究。在本研究中，筆者利用生物信息學方法篩選出有效的分子標志物，構(gòu)建和驗證免疫相關(guān)的預(yù)后模型，并分析TME 特征，以期為尋找結(jié)腸癌新的治療靶點提供參考依據(jù)。

1 材料和方法

1.1 數(shù)據(jù)下載和獲取 從 UCSC Xena 平臺（https：／ ／xenabrowser.net／datapages／）下 載 結(jié) 腸 癌 的mRNA 表達數(shù)據(jù)和對應(yīng)的臨床信息，選擇隊列為GDC TCGA Colon Cancer（COAD）作為訓(xùn)練數(shù)據(jù)集，具體樣本信息見表1。在GEO數(shù)據(jù)庫中下載數(shù)據(jù)集GSE39582的基因表達譜和臨床信息作為驗證數(shù)據(jù)集，具體信息見表2。從ImmPort 數(shù)據(jù)庫（https：／ ／immport.niaid.nih.gov／）下載免疫相關(guān)基因列表，總計 1 509 個。

表1 TCGA數(shù)據(jù)集的具體樣本信息［n（%）］

表2 數(shù)據(jù)集GSE39582的具體樣本信息［n（%）］

1.2 構(gòu)建預(yù)后風險模型 利用 R 包 limma 對TCGA數(shù)據(jù)集進行差異基因分析，過濾標準為adjusted P 值小于 0.05 和差異倍數(shù)大于 1.5 倍（｜log2FC｜＞0.585），與免疫基因列表取交集，得到免疫相關(guān)的差異表達基因（differentially expressed genes，DEGs）。利用R包Survival對免疫相關(guān)的差異基因與患者總生存時間（overall survival，OS）進行單因子回歸分析，篩選P＜0.05 的基因作為候選預(yù)后因子。通過多因子回歸分析確定每個預(yù)后因子的回歸系數(shù)，建立預(yù)后風險評估模型，預(yù)測患者生存率。

公式：風險分數(shù)＝∑差異基因的回歸系數(shù)χi×歸一化處理后的基因表達量βi

1.3 統(tǒng)計學分析 利用 R 包 SurvivalROC 繪制ROC曲線以評估預(yù)后模型的性能。以真陽性率與假陽性率之間的差值最大處作為最佳風險評分臨界值，高于臨界值的患者歸為高風險評分組，低于臨界值的患者歸為低風險評分組，并使用R包Survminer繪制兩組的生存曲線。采用相同的風險計算公式對驗證數(shù)據(jù)集GSE39582進行統(tǒng)計學分析。

1.4 構(gòu)建和驗證列線圖 使用R包rms，通過整合風險評分模型和臨床信息，包括年齡、性別和腫瘤分期，構(gòu)建列線圖模型，可視化不同患者特征的預(yù)后價值。通過繪制森林圖展示臨床信息與OS之間的關(guān)系，其中，一致性指數(shù)（C-index）表示列線圖的預(yù)測準確性。

1.5 估計及比較腫瘤免疫浸潤細胞類型的構(gòu)成和差異 基于構(gòu)建的預(yù)后風險模型，將腫瘤樣本分為高風險評分組和低風險評分組，采用CIBERSORT算法分析22種腫瘤浸潤免疫細胞的比例［4］。通過非配對t檢驗比較高風險評分組和低風險評分組之間的免疫圖譜。

1.6 分析腫瘤免疫微環(huán)境的特征 使用 R 包maftools分析和比較高風險評分組和低風險評分組的突變譜［5］，并計算TMB值。通過未配對t檢驗統(tǒng)計分析高風險評分組和低風險評分組之間TMB 的差異。使用Kaplan-Meier分析風險評分與OS 之間的關(guān)聯(lián)。利用Wilcoxon 檢驗比較高風險評分組和低風險評分組的免疫檢查點及其配體的mRNA 表達水平的差異。

2 結(jié)果

2.1 基于訓(xùn)練數(shù)據(jù)集的免疫相關(guān)預(yù)后模型構(gòu)建對訓(xùn)練數(shù)據(jù)集中的腫瘤樣本和正常樣本的表達譜進行差異性分析，得到571個DEGs。其中，275 個基因上調(diào)，296 個下調(diào)（圖1A）。與免疫基因列表取交集后，得到102 個免疫相關(guān)的DEGs，其中19個基因上調(diào)，83個基因下調(diào)。

通過單因子回歸和多因子回歸分析表明，有6個免疫相關(guān)的 DEGs與OS 相關(guān)（P 值＜0.05），見表3。其中，回歸系數(shù)見圖1B、1D。根據(jù)逐步回歸模型，Akaike信息標準（AIC）為983.38，C 指數(shù)為0.63，見圖1E。

表3 與結(jié)腸癌預(yù)后相關(guān)的基因信息

圖1 鑒定預(yù)后相關(guān)的差異表達基因

由圖2A～2F 可知，這6 個免疫相關(guān)的差異基因可作為獨立的預(yù)后因子，均與OS 顯著相關(guān)（P＜0.05）。其中，BMP5的高表達與良好的預(yù)后顯著相關(guān)，而其他5個基因則相反。

圖2 基因表達量高低與患者生存狀態(tài)的關(guān)系

2.2 預(yù)后模型的性能評估 基于已構(gòu)建的預(yù)后風險模型，將結(jié)腸癌患者分為高風險評分組和低風險評分組，其中，cut-off值設(shè)為－0.19。由圖 3A、3B 可知，隨著風險評分的增加，生存時間呈縮短的趨勢，且高風險評分組的死亡比例（圖3B 右側(cè)的紅點和藍點比例）比低風險評分組高，這與Kaplan-Meier生存分析結(jié)果一致。由圖3D 可知，高風險評分患者的OS比低風險評分患者預(yù)后較差（P＜0.000 1）。由圖3C可知，基因BMP5在低風險組表達量高，在高風險組表達量低，而其他5個基因趨勢相反。

圖3 訓(xùn)練數(shù)據(jù)集預(yù)后風險模型的構(gòu)建

2.3 預(yù)后模型的統(tǒng)計分析 通過繪制ROC曲線和腫瘤分層分析來進一步評估預(yù)后風險模型的性能。由圖4A可知，ROC曲線下面積（AUCROC）在3年時為 0.668，4 年時為 0.699，5 年時為 0.696。Wilcoxon檢驗表明，較高的風險評分與較高的病理分期（P ＝0.000 29）、T 分期（P ＝ 0.000 05）、M 分期（P ＝0.023）、N 分期（P ＝0.001 3）顯著相關(guān)。

圖4 訓(xùn)練數(shù)據(jù)集預(yù)后風險模型的統(tǒng)計分析

2.4 預(yù)后模型的驗證 基于構(gòu)建的預(yù)后風險模型，在驗證數(shù)據(jù)集GSE39582 中進一步驗證。低風險評分與良好的預(yù)后顯著相關(guān)（P ＝0.008 1，圖5A），且較高的風險評分與較高的病理分期（P ＝0.003 7）、T 分期（P ＝0.047）、M 分期（P ＝0.047）、N分期（P ＝0.032）顯著相關(guān)（圖 5B～5E）。

圖5 驗證數(shù)據(jù)集預(yù)后風險模型的統(tǒng)計分析

2.5 構(gòu)建和評估列線圖模型 在列線圖中，每個變量的得分映射到分數(shù)軸上，通過計算總分來估算3年、4年和5年的生存概率（圖6A）。由森林圖可以看出，患者特征，包括年齡（＞60 歲）、腫瘤分期（Ⅲ和Ⅳ）和風險評分與 OS 顯著相關(guān)（P 值＜0.05，圖6B）。此外，通過繪制校準曲線以驗證列線圖的性能，可看出預(yù)測曲線接近理想曲線（圖6C ～6E），這表明構(gòu)建的列線圖模型性能良好。此外，該列線圖（C-index：0.75）的預(yù)測準確性高于風險評分模型（C-index：0.63）。

圖6 列線圖模型的構(gòu)建與驗證

2.6 探索結(jié)腸癌TME 基于CIBERSORT算法，估算每例結(jié)腸癌患者中22 種免疫細胞的比例，并比較高風險評分組和低風險評分組間的免疫細胞比例，發(fā)現(xiàn)有6種免疫細胞類型間的差異有統(tǒng)計學意義（P＜0.05）。其中，CD4 記憶活化 T 細胞（T cells CD4 memory activated）和Tregs調(diào)控T細胞（T cells regulatory，Tregs）最為顯著（P＜0.000 1，圖 7A）。由繪制的結(jié)腸癌患者的突變譜可以看出，前20 位顯著突變的基因分別為 APC、TP53、TTN、KRAS、SYNE1、MUC16、PIK3CA、FAT4、RYR2、OBSCN、ZFHX4、DNAH5、PCLO、CSMD3、ABCA13、DNAH11、LRP1B、FAT3、USH2A 和 CSMD1（圖 7B）。通過計算并比較每個樣本的TMB，發(fā)現(xiàn)高風險評分組和低風險評分組間無顯著性差異（P ＝0.25，圖 7C），高TMB 與較差患者的 OS 顯著相關(guān)（P ＝ 0.028，圖7D）。利用Wilcoxon 檢驗比較高風險評分組和低風險評分組間免疫檢查點及其配體的表達水平，結(jié)果表明，高風險評分組 PD-L1（P ＝ 0.001）、PD-1（P ＝0.000 22）和 CTLA-4（P ＝0.000 21）的表達水平較高（圖7E～7G）。

圖7 結(jié)腸癌腫瘤免疫微環(huán)境特征

3 討論

TME 中的免疫細胞在腫瘤的進展中起著重要作用［3］。研究表明，免疫檢查點抑制劑（immune checkpoint inhibitors，ICIs）在結(jié)腸癌免疫治療中具有很大的潛力［6］。然而，由于對TME 特性的理解不足，在一定程度上阻礙了免疫治療的廣泛應(yīng)用。近年來，大量研究鑒定了與結(jié)腸癌診斷和預(yù)后相關(guān)的生物標志物［7-8］。然而，為了最大限度地發(fā)揮免疫治療的作用，還需要探索更多可靠的生物標志物。

本研究構(gòu)建1 個基于6 個免疫相關(guān)差異基因的結(jié)腸癌預(yù)后風險評估模型，該6 個基因分別為BMP5、OXTR、STC2、UCN、GRP 和 AMH。研究表明，GMP5、OXTR、STC2 和 UCN 基因與結(jié)腸癌的預(yù)后顯著相關(guān)［9-12］。BMP5 基因參與 TGF-β ／Smad 信號通路，其表達與上皮標志物呈正相關(guān)，但與間充質(zhì)標志物呈負相關(guān)［13］。同時，BMP5 基因被證實與PI3K-AKT 和 MAPKs 信 號 通 路 相 互 作 用［14］。OXTR基因及其配體催產(chǎn)素（OXT）存在于胃腸系統(tǒng)中，并參與腫瘤的發(fā)生、發(fā)展［15］。AMH基因被證實與結(jié)腸癌的發(fā)生、發(fā)展有關(guān)，并可作為槲皮素的靶點，降低結(jié)腸癌引起的死亡率［16］。對于STC2、UCN和GRP基因，還未見其與腫瘤進展的相關(guān)性研究。

為了進一步探索建立的預(yù)后模型的潛在應(yīng)用價值，筆者分析了訓(xùn)練數(shù)據(jù)集和驗證數(shù)據(jù)集中不同病理分期患者的風險評分，結(jié)果表明隨著風險評分的增加，其腫瘤分期等級越高，這預(yù)示著建立的預(yù)后模型在腫瘤分期方面有一定的預(yù)測能力。此外，通過系統(tǒng)分析結(jié)腸癌TME 的特征，發(fā)現(xiàn)CD4 記憶活化T 細胞、Tregs 調(diào)控T 細胞的免疫浸潤程度和TMB可能是結(jié)腸癌的獨立預(yù)后標志。此外，高風險評分患者的免疫檢查點及其配體（PD-L1、PD-1 和CTLA-4）的表達水平較高，表明構(gòu)建的預(yù)后模型能夠為免疫治療提供一定的參考價值。

綜上所述，本研究構(gòu)建了1個性能良好的免疫相關(guān)的結(jié)腸癌風險評估預(yù)后模型，并探索其TME特征，這可能有助于結(jié)腸癌患者的預(yù)后風險預(yù)測、腫瘤分期預(yù)測和免疫治療。為進一步提升預(yù)后模型的預(yù)測性能，本研究基于風險評分預(yù)后模型和臨床指標（包括年齡、性別、病理分期）構(gòu)建了列線圖模型，使預(yù)后模型的準確性得到較大提升。然而，本研究構(gòu)建的預(yù)后模型是利用生物信息學和生物統(tǒng)計學分析的技術(shù)手段，基于公開的數(shù)據(jù)集進行分析，篩選出的預(yù)測因子的生物學功能還需要進一步的實驗驗證。期望本研究的研究成果能為結(jié)腸癌患者的預(yù)后評估、腫瘤分期預(yù)測和免疫治療提供新的研究方向。