結直腸癌樣本免疫細胞浸潤及腫瘤突變負荷分析

王 娜，王 平，何儀佳，曾 晶，韓聰玲，王思睿，軒小燕，

(1．鄭州大學基礎醫學院，河南 鄭州 450001；2．黃河科技學院，河南 鄭州 450006)

結直腸癌(colorectal cancer，CRC)是全球發病率第3、死亡率第4 的腫瘤，每年有超過120 萬新增病例，60多萬人死于CRC[1-3]。手術是目前治愈早期CRC的主要方法，但由于起病隱匿，多數患者發現時已有轉移。化療、放療等延長了晚期CRC患者生存期，但治療效果并不持久，耐藥性不可避免地通過各種機制產生[4]。基因突變尤其是驅動基因如APC、TP53、KRAS的突變在不同腫瘤中突變頻率不同，是影響腫瘤發生發展的重要因素[5-7]，同時也影響腫瘤靶向治療效果。以免疫檢查點抑制劑尤其是PD-1/PD-L1 為代表的腫瘤免疫治療是近年發展起來的非常有前景的治療方法，在一些晚期CRC患者中具有顯著療效，但研究發現僅對微衛星不穩定性轉移性CRC 患者有效[8]，獲益人群有限，仍急需開發新的治療方法。研究CRC的腫瘤免疫細胞浸潤微環境、基因突變圖譜可為深入探討CRC發生發展的分子機制以及為CRC的靶向治療或免疫治療提供理論基礎。

轉錄組測序技術以及生物信息技術的發展，可以幫助我們多維度分析免疫細胞浸潤、腫瘤細胞突變等影響腫瘤發生、發展、轉移及治療等重要因素。本文分析了來自癌癥基因組圖譜集(The Cancer Genome Atlas，TCGA)和基因表達綜合數據庫(Gene Expression Omnibus，GEO)的CRC 樣本的免疫細胞浸潤、基因集富集通路、腫瘤突變負荷等信息，為進一步研究CRC發生和轉移的分子機制提供生物信息學數據支撐。

1 材料與方法

1.1 數據來源

從美國國家癌癥研究所TCGA數據庫(https://portal.gdc.cancer.gov/)下載473 例CRC 樣本轉錄組表達數據、452 例CRC 患者臨床信息數據以及399 例CRC 患者突變數據。從美國國家生物信息技術中心GEO 數據庫(https://www.ncbi.nlm.nih.gov/geo/)下載GSE40967 數據集，含585 例CRC 患者樣本芯片表達數據及臨床信息，所有數據下載日期為2022年3月8 日。依據文獻方法[9]將來自TCGA 數據庫中的CRC 樣本基因表達數據由FPKM 轉化為TPM 后，與GSE40967 數據集的樣本表達數據合并作為本文基因表達分析文件，采用ComBat 算法對數據進行批次校正。整理合并TCGA 及GEO 臨床數據作為后續生存分析文件，在R4.2.0版本下運行數據。

1.2 CRC免疫細胞浸潤分析

采用CIBERSORT 和ESTIMATE 算法分析CRC 樣本基因表達數據，獲得樣本免疫評分、基質評分以及22種免疫細胞在樣本中的比例[10]。采用Spearman檢驗分析樣本中免疫細胞浸潤程度、免疫評分和基質評分之間的相關性。根據CRC樣本各免疫細胞比例的相似性，用ConsensusClusterPlus 將樣本聚類為免疫細胞聚類A、B 兩組，采用pheatmap 包繪制免疫細胞浸潤熱圖。采用Kaplan-Meier 分析免疫細胞聚類A、B 兩組間樣本總生存期的差異。采用limma 包分析PD-L1基因表達在免疫細胞聚類A、B兩組間的差異。

1.3 免疫細胞評分模型的構建及樣本的生存分析和基因集富集分析

采用limma包分析免疫細胞聚類A、B兩組間樣本的所有差異表達基因(differentially expressed genes，DEG)，以|Log2(FC)|＞1，P＜0.05 為過濾條件，FC 指基因表達的差異倍數(fold change)。采用Boruta 包查找DEG 的特征基因，用主成分分析法(principal component analysis，PCA)對特征基因降維分析獲得每個CRC 樣本的特征基因的特征值，采用類似基因表達指數的方法[11]利用公式“免疫細胞評分=∑PC1A-∑PC1B”(即每個樣本的所有A 組特征基因的特征值減去所有B 組特征基因的特征值)建立CRC 樣本的免疫細胞評分模型，對每個樣本進行評分。采用surv_cutpoint 函數查找評分的最佳臨界值，根據最佳臨界值將CRC樣本分為免疫細胞評分高、低兩組。采用Kaplan-Meier 分析組間樣本總生存期差異。對樣本進行基因集富集分析(gene set enrichment analysis，GSEA)，闡明免疫細胞評分高、低兩組樣本的基因富集通路。

1.4 腫瘤突變負荷計算及生存分析

腫瘤突變負荷(tumor mutation burden，TMB)指每百萬堿基中發生置換、插入/缺失突變的總數，用來評估腫瘤的基因突變頻率。從TCGA 下載CRC 突變數據，統計非同義突變數目，計算樣本TMB。用surv_cutpoint 函數查找TMB 最佳臨界值，根據最佳臨界值將CRC樣本分成腫瘤突變負荷高、低兩組，采用Kaplan-Meier 分析組間樣本總生存期差異，采用Maftool繪制樣本基因突變瀑布圖。

1.5 統計分析

采用Wilcoxon 檢驗進行兩組之間比較，Kaplan-Meier(log-rank 檢驗)進行生存分析，Spearman 進行相關性分析。P＜0.05為差異具有統計學意義。

2 結 果

2.1 CRC樣本免疫細胞浸潤分析

免疫細胞浸潤分析發現，CRC 腫瘤組織內CD8+T細胞和初始記憶性CD4+T 細胞(r=-0.34，P＜0.05)、活化的記憶性CD4+T 細胞(r=0.37，P＜0.05)、M0 巨噬細胞(r=-0.47，P＜0.05)具有相關性。中性粒細胞和活化的漿細胞具相關性(r=0.41，P＜0.05)。在22種免疫細胞中，免疫評分僅與M1 巨噬細胞呈弱相關(r=0.32)，而與其他21 種細胞均無明顯相關(|r|＜0.3，P＞0.05，圖1A)。根據CRC 樣本中浸潤的免疫細胞比例的相似性，采用ConsensusClusterPlus 對樣本進行聚類為免疫細胞聚類A和B兩群。CD8+T細胞、M0巨噬細胞、活化的記憶性CD4+T細胞等呈現明顯聚類(圖1B)。生存分析提示免疫細胞聚類A、B 兩組患者總生存期無顯著差異(P=0.928，圖1C)。PD-L1基因表達水平在兩組之間差異無統計學意義(P＞0.05，圖1D)。

圖1 CRC樣本免疫細胞浸潤分析

2.2 CRC 樣本免疫細胞評分模型建立及基因集富集分析

采用limma包分析免疫細胞聚類A、B組間樣本的DEG，采用Boruta 包查找DEG 的特征基因，依據文獻方法將特征基因分為A、B 兩組，若特征基因的相對表達量在免疫細胞聚類A組樣本呈高表達而在免疫細胞聚類B 組樣本低表達，該基因定義為A 組特征基因，反之定義為B 組特征基因[12]。經計算，共獲得28個特征基因，其中A 組特征基因18 個，分別為CLCA1、IGHM、IGLJ3、ZG16、ITLN1、ADH1C、CLCA4、CA2、CA4、PIGR、UGT2B17、MS4A12、HEPACAM2、FCGBP、PLA2G2A、CEACAM7、SPINK4、REG4；B 組特征基因10 個，分別為MMP9、SPP1、COL10A1、THBS2、COL11A1、FNDC1、CXCL5、COMP、SFRP4、CXCL8。依據CRC 樣本免疫細胞評分模型，計算每個CRC 樣本的免疫細胞評分，采用surv_cutpoint 函數查找最佳臨界值，依據臨界值將CRC樣本分為免疫細胞評分高、低兩組。生存分析提示免疫細胞評分高組的患者總生存期長于免疫細胞評分低組的患者(P=0.005，圖2A)。

圖2 CRC樣本生存分析及基因集富集分析

GSEA 分析顯示，免疫細胞評分高組的基因主要富集在KEGG_FATTY_ACID_METABOLISM，NITROGEN_METABOLISM、KEGG_CITRATE_CYCLE_TCA_CYCLE、KEGG_PYRUVATE_METABOLISM、KEGG_RETINOL_METABOLISM 等信號通路，這些通路在碳水化合物、脂肪酸代謝中發揮重要作用，與腸道的正常生理功能密切相關。值得注意的是，未見免疫相關基因在免疫細胞評分高組樣本中富集。免疫細胞評分低組的基因主要富集在KEGG_COLORECTAL_CANCER，ECM_RECEPTOR_INTERACTION， KEGG_REGULATION_OF_ACTIN_CYTOSKELETON 等信號通路，這些通路與癌癥的形成、發展和轉移密切相關，這也跟免疫細胞評分低組患者的低生存期相一致(圖2B)。

2.3 腫瘤突變負荷分析

TMB 是腫瘤對免疫治療反應的潛在生物標志物[13]，較高的TMB增加了腫瘤新抗原產生的概率。本文分析了免疫細胞評分高、低兩組樣本的基因突變數據，用maftool繪制突變瀑布圖，列出變異頻率最高的前20個基因。值得注意的是，兩組之間突變頻率最高的20個基因一致，這一結果和其他課題組分析的結果一致[14-16]。CRC 腫瘤驅動突變基因APC突變頻率在免疫細胞評分高、低組分別為77%和73%，TP53分別為47%和57%，KRAS分別為45%和41%。APC基因以移碼插入突變和移碼缺失突變為主，TP53、KRAS基因以錯義突變為主(圖3A、3B)。TMB 在免疫細胞評分高、低組間的差異無統計學意義(P=0.17，圖3C)。生存分析提示腫瘤突變負荷低組患者總生存期長于腫瘤突變負荷高組患者(P=0.035，圖3D)。

圖3 CRC腫瘤突變負荷分析

3 討論

CRC是常見的惡性腫瘤，手術、放療、化療以及免疫治療在一定程度上延長了腫瘤患者的生存時間，但目前的治療方法均不理想。深入分析腫瘤的免疫細胞浸潤、腫瘤突變負荷等信息，可為研究影響CRC發生發展的關鍵因素以及為CRC的治療提供理論基礎。

為了降低不同測序平臺對結果的影響，本研究將來自TCGA 數據庫的CRC 樣本轉錄組數據如前所述標準化后和來自GEO 數據庫的GSE40967 數據集進行合并。通過對免疫細胞浸潤分析，未發現在CRC樣本中具有強相關性的免疫細胞，免疫細胞聚類A、B 兩群患者之間總生存期差異無統計學意義，這些結果提示在CRC 樣本中抗腫瘤免疫未充分活化發揮抗腫瘤作用。PD-L1的表達是免疫檢查點抑制劑是否有效的主要指標[17]，本文分析發現PD-L1表達水平在免疫細胞聚類A、B 兩組之間差異無統計學意義，這也解釋了臨床觀察到的大多數CRC患者對免疫檢查點抑制劑靶向治療無反應。

為進一步探索影響CRC 患者生存差異的關鍵因素，本文構建了CRC樣本免疫細胞評分模型，依據免疫細胞評分將樣本分成免疫細胞評分高、低兩組。免疫細胞評分高組患者總生存期長于免疫細胞評分低組患者，提示該模型具有潛在預測CRC 患者預后的價值。GSEA 分析發現脂肪酸代謝、氮代謝、碳代謝等主要的物質代謝信號通路在免疫細胞評分高組富集，提示免疫細胞評分高組樣本的活躍基因以參與生理代謝為主，推測免疫細胞評分高組樣本腫瘤細胞分化程度高，執行生理功能相關基因具有表達優勢。有趣的是，一些跟腫瘤發生、發展密切相關的通路在免疫細胞評分低組樣本中富集，提示免疫細胞評分低組樣本分化程度低，惡性程度高，這也跟免疫細胞評分低組患者的低生存率相一致。據此，我們推測免疫細胞評分高、低組患者之間的生存差異主要由與腫瘤發生發展相關的一些信號通路決定。GSEA 分析發現在免疫細胞高、低組均未見有免疫活化相關通路富集，這也和上述免疫細胞浸潤觀察到的結果一致，提示在CRC樣本中，腫瘤免疫應答未充分活化而發揮抗腫瘤作用。

腫瘤的發生由基因突變積累引起，這些基因所編碼的蛋白質往往是細胞周期信號通路中的重要蛋白質，控制細胞分裂和生長的驅動基因發生突變，在腫瘤的發生、發展、轉移、免疫逃逸以及治療耐藥中發揮關鍵作用[18-19]。CRC 的基因組已經進行了廣泛深入的研究[20-21]，APC、KRAS、TP53等基因為CRC 的驅動基因[22]，在CRC 發生發展和轉移中發揮重要作用。APC突變導致APC/β-catenin 通路的過度活化[23-24]，KRAS突變使其下游信號通路異常[25]，抑癌基因TP53突變后產生的蛋白導致其生理功能喪失，這些驅動基因以不同的方式參與細胞增殖、侵襲和轉移[26-27]。TMB 分析發現APC、KRAS、TP53的突變在免疫細胞評分高、低組樣本中均占有較大的比例，結合上述免疫細胞評分高、低組GSEA 分析結果，我們認為，在CRC 樣本中，與機體產生的抗腫瘤免疫作用相比較，驅動基因的突變是決定腫瘤惡性程度以及腫瘤發展和轉移的主要因素。本文數據為進一步深入研究CRC發生和轉移機制提供了必要的生信數據支撐。