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

(1．鄭州大學(xué)基礎(chǔ)醫(yī)學(xué)院，河南 鄭州 450001；2．黃河科技學(xué)院，河南 鄭州 450006)

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

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

1 材料與方法

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

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

1.2 CRC免疫細(xì)胞浸潤(rùn)分析

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

1.3 免疫細(xì)胞評(píng)分模型的構(gòu)建及樣本的生存分析和基因集富集分析

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

1.4 腫瘤突變負(fù)荷計(jì)算及生存分析

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

1.5 統(tǒng)計(jì)分析

采用Wilcoxon 檢驗(yàn)進(jìn)行兩組之間比較，Kaplan-Meier(log-rank 檢驗(yàn))進(jìn)行生存分析，Spearman 進(jìn)行相關(guān)性分析。P＜0.05為差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1 CRC樣本免疫細(xì)胞浸潤(rùn)分析

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

圖1 CRC樣本免疫細(xì)胞浸潤(rùn)分析

2.2 CRC 樣本免疫細(xì)胞評(píng)分模型建立及基因集富集分析

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

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

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

2.3 腫瘤突變負(fù)荷分析

TMB 是腫瘤對(duì)免疫治療反應(yīng)的潛在生物標(biāo)志物[13]，較高的TMB增加了腫瘤新抗原產(chǎn)生的概率。本文分析了免疫細(xì)胞評(píng)分高、低兩組樣本的基因突變數(shù)據(jù)，用maftool繪制突變瀑布圖，列出變異頻率最高的前20個(gè)基因。值得注意的是，兩組之間突變頻率最高的20個(gè)基因一致，這一結(jié)果和其他課題組分析的結(jié)果一致[14-16]。CRC 腫瘤驅(qū)動(dòng)突變基因APC突變頻率在免疫細(xì)胞評(píng)分高、低組分別為77%和73%，TP53分別為47%和57%，KRAS分別為45%和41%。APC基因以移碼插入突變和移碼缺失突變?yōu)橹?，TP53、KRAS基因以錯(cuò)義突變?yōu)橹?圖3A、3B)。TMB 在免疫細(xì)胞評(píng)分高、低組間的差異無統(tǒng)計(jì)學(xué)意義(P=0.17，圖3C)。生存分析提示腫瘤突變負(fù)荷低組患者總生存期長(zhǎng)于腫瘤突變負(fù)荷高組患者(P=0.035，圖3D)。

圖3 CRC腫瘤突變負(fù)荷分析

3 討論

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

為了降低不同測(cè)序平臺(tái)對(duì)結(jié)果的影響，本研究將來自TCGA 數(shù)據(jù)庫(kù)的CRC 樣本轉(zhuǎn)錄組數(shù)據(jù)如前所述標(biāo)準(zhǔn)化后和來自GEO 數(shù)據(jù)庫(kù)的GSE40967 數(shù)據(jù)集進(jìn)行合并。通過對(duì)免疫細(xì)胞浸潤(rùn)分析，未發(fā)現(xiàn)在CRC樣本中具有強(qiáng)相關(guān)性的免疫細(xì)胞，免疫細(xì)胞聚類A、B 兩群患者之間總生存期差異無統(tǒng)計(jì)學(xué)意義，這些結(jié)果提示在CRC 樣本中抗腫瘤免疫未充分活化發(fā)揮抗腫瘤作用。PD-L1的表達(dá)是免疫檢查點(diǎn)抑制劑是否有效的主要指標(biāo)[17]，本文分析發(fā)現(xiàn)PD-L1表達(dá)水平在免疫細(xì)胞聚類A、B 兩組之間差異無統(tǒng)計(jì)學(xué)意義，這也解釋了臨床觀察到的大多數(shù)CRC患者對(duì)免疫檢查點(diǎn)抑制劑靶向治療無反應(yīng)。

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

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