李偉華 趙鵬宇 黎鴻堅(jiān) 劉林江 李運(yùn)潔 張芳 鄒海軍 王玉里

結(jié)直腸癌是全球五種最常診斷的癌癥之一，占所有癌癥的6%，是癌癥死亡的第三大原因［1］。在中國，CRC的年發(fā)病率排名第三，是癌癥相關(guān)死亡的第五大原因［2］。雖然近年來結(jié)腸癌的臨床治療模式不斷進(jìn)步［3］，患者治愈率逐漸升高，但是結(jié)直腸癌患者的預(yù)后生存情況一直不容樂觀，5年生存率僅超過50%［4］，迫切需要尋找準(zhǔn)確的預(yù)后標(biāo)志物。

隨著新一代測序技術(shù)的發(fā)展及分子生物學(xué)的進(jìn)步，從分子層面解析疾病發(fā)生與發(fā)展機(jī)理、藥物反應(yīng)差異以及預(yù)后差異成為可能。近年來，免疫微環(huán)境的影響逐漸成為結(jié)直腸癌研究的熱點(diǎn)［5］。有證據(jù)表明，免疫相關(guān)基因通過介導(dǎo)炎癥或免疫監(jiān)視逃避在結(jié)直腸癌發(fā)生與發(fā)展中起著重要作用［6-8］，并對(duì)結(jié)直腸癌患者的預(yù)后生存具有顯著的影響［4，9］。

本文基于公開已發(fā)表的數(shù)據(jù)集，通過多種生物信息學(xué)分析方法對(duì)免疫相關(guān)基因在結(jié)直腸癌患者中的表達(dá)以及功能狀態(tài)進(jìn)行全面刻畫，并篩選得到對(duì)結(jié)腸癌患者預(yù)后生存影響最顯著的5個(gè)免疫相關(guān)基因構(gòu)建了Risk Score模型。構(gòu)建的Risk Score在獨(dú)立驗(yàn)證集中顯示出較高的準(zhǔn)確率，證明了該模型的有效性與普遍性。我們的研究結(jié)果為預(yù)測結(jié)直腸癌的預(yù)后生存提供了有效的生物標(biāo)志物模型。

材料與方法

一、數(shù)據(jù)來源

從 The Cancer Genome Atlas［10］（TCGA，https://www.cancer.gov/tcga）數(shù)據(jù)庫中獲取結(jié)腸癌（Colon adenocarcinoma，COAD）組織與正常組織樣本的mRNA表達(dá)譜數(shù)據(jù)及相應(yīng)的臨床信息，包括478個(gè)癌癥樣本和41個(gè)正常樣本。從Immport［11］（https://www.immport.org/shared/home）數(shù)據(jù)庫中下載得到2013個(gè)免疫相關(guān)基因。從Gene Expression Omnibus（GEO，https://www.ncbi.nlm.nih.gov/GEO）數(shù)據(jù)庫中下載GSE39582表達(dá)譜及其臨床信息作為驗(yàn)證集，包括586個(gè)癌癥樣本和19個(gè)正常樣本。

二、數(shù)據(jù)預(yù)處理與差異表達(dá)分析

將TCGA下載的結(jié)腸癌與正常樣本的mRNA表達(dá)譜進(jìn)行預(yù)處理。然后，使用“Deseq2”包［12］對(duì)mRNA表達(dá)譜進(jìn)行差異分析，其中|log2FC|＞1且FDR＜0.05的基因被認(rèn)為是具有統(tǒng)計(jì)學(xué)意義的差異基因。

三、構(gòu)建免疫預(yù)后特征模型

首先，我們將差異分析得到的差異基因與免疫相關(guān)基因取交集得到免疫差異基因集。然后，結(jié)合相應(yīng)的臨床信息進(jìn)行單因素COX回歸分析，定義P＜0.05的mRNA為候選的免疫相關(guān)mRNA。接下來，對(duì)其進(jìn)行“向前”多因素COX回歸分析，評(píng)估它們作為獨(dú)立預(yù)后因素對(duì)患者生存的貢獻(xiàn)。因此，我們確定了五個(gè)免疫相關(guān)的mRNA作為預(yù)后特征，并基于其表達(dá)水平和多因素COX回歸模型的回歸系數(shù)進(jìn)行線性組合構(gòu)建預(yù)后模型。具體公式如下：

其中，Coefi代表估計(jì)的回歸系數(shù)，xi代表免疫相關(guān)的mRNA表達(dá)值。

四、功能富集分析

將識(shí)別得到的免疫預(yù)后相關(guān)的風(fēng)險(xiǎn)mRNA使用R包“ClusterProfiler”［13］分別對(duì)其進(jìn)行基因本體（GO）［14］富集分析，以確定與免疫預(yù)后關(guān)鍵基因相關(guān)的生物過程、分子功能、細(xì)胞組分和信號(hào)通路，研究其潛在的生物學(xué)功能。

五、免疫預(yù)后模型的評(píng)估與驗(yàn)證

為評(píng)估模型預(yù)后能力，根據(jù)COAD患者免疫預(yù)后風(fēng)險(xiǎn)評(píng)分，使用“survival”包［15］比較總體生存時(shí)間，繪制Kaplan-Meier（KM）生存曲線，并進(jìn)行Log-rank檢驗(yàn)；使用“timeROC”包［16］進(jìn)行繪制時(shí)間相關(guān)的受試者工作特征（ROC）曲線，以評(píng)估免疫預(yù)后模型的有效性。然后，從GEO數(shù)據(jù)庫中下載獨(dú)立數(shù)據(jù)集GSE39582繪制KM生存曲線和時(shí)間相關(guān)的ROC曲線來驗(yàn)證模型的魯棒性，檢驗(yàn)?zāi)Ｐ偷念A(yù)后價(jià)值。

六、統(tǒng)計(jì)學(xué)分析

所有數(shù)據(jù)的統(tǒng)計(jì)分析均使用R軟件（版本4.0.2，https://www.r-project.org/） 進(jìn) 行 。Log-rank檢驗(yàn)和秩和檢驗(yàn)用于計(jì)算兩組數(shù)據(jù)間的差異。通過Wald檢驗(yàn)或Fisher精確檢驗(yàn)分析分類數(shù)據(jù)。P值多重檢驗(yàn)校正采用Benjamini&Hochberg（BH）方法。P＜0.05認(rèn)為差異具有統(tǒng)計(jì)學(xué)意義。

結(jié) 果

一、結(jié)腸癌差異表達(dá)mRNA

將從TCGA下載的結(jié)腸癌的癌癥與正常樣本進(jìn)行差異分析，共識(shí)別出4 456個(gè)差異基因，這些基因由2 264個(gè)上調(diào)基因和2 192個(gè)下調(diào)基因組成。根據(jù)差異結(jié)果做火山圖（圖1A），紅色代表上調(diào)基因，藍(lán)色代表下調(diào)基因，灰色代表差異不顯著的基因。從Immpot中下載免疫相關(guān)基因，與差異表達(dá)基因取交集，最終獲得362個(gè)結(jié)腸癌免疫相關(guān)基因（圖1B）。

圖1 TCGA中結(jié)腸癌免疫差異基因的鑒定。1A：結(jié)腸癌與正常樣本之間差異表達(dá)mRNA的火山圖；1B：識(shí)別COAD免疫差異基因的韋恩圖

二、關(guān)鍵mRNA的篩選

為了挖掘出與疾病更為相關(guān)的mRNA，結(jié)合臨床信息，使用單因素COX對(duì)免疫相關(guān)基因進(jìn)行篩選。獲得與患者生存更為相關(guān)的關(guān)鍵mRNA。在預(yù)后mRNA中截取P＜0.05的為關(guān)鍵mRNA（表1）。關(guān)鍵mRNA包括TPM2、BMP5、MAPT、SCTR、PTH1R、NGFR、NRG1、XCL1、NGF、CD1A、CD1B、 PLXNA3、 IL13RA2、 OXT、 PGF、TNFRSF19、 MC1R、 LTB4R、 HAMP、 JAG2、LHB、 NMB、 VGF、 GRP、 INHBB、 UCN、SLC11A1、 OXTR、 CXCL1、 TDGF1、 EREG、STC2總共32個(gè)。這32個(gè)mRNA都顯示出與結(jié)腸癌有著很強(qiáng)的預(yù)后相關(guān)性。對(duì)上述32個(gè)mRNA進(jìn)行功能富集分析（圖2），預(yù)后相關(guān)基因調(diào)控T細(xì)胞介導(dǎo)的免疫、T細(xì)胞介導(dǎo)的細(xì)胞毒性的正調(diào)等相關(guān)的生物學(xué)功能，進(jìn)而影響患者生存。

表1 TCGA總體生存率的單因素cox回歸分析

圖2 預(yù)后風(fēng)險(xiǎn)mRNA富集分析

為了找出這32個(gè)mRNA中與結(jié)腸癌最為相關(guān)的免疫預(yù)后mRNA，對(duì)其進(jìn)行多因素COX回歸分析。篩選出其中預(yù)后更為相關(guān)的mRNA以構(gòu)建模型。選取P＜0.05的mRNA為預(yù)后風(fēng)險(xiǎn)mRNA。篩選出的預(yù)后風(fēng)險(xiǎn)mRNA包括SCTR（HR=1.21，P=0.02）、XCL1（HR=1.16，P=0.04）、NGF（HR=1.26，P=0.047）、CD1B （HR=0.79，P=0.013）、EREG（HR=0.88，P=0.019）。

利用獲得的5個(gè)基因，構(gòu)建風(fēng)險(xiǎn)得分模型，即：風(fēng)險(xiǎn)評(píng)分=（0.191*SCTR的表達(dá)值）+（0.151*XCL1） +（0.231*NGF） +（ -0.241*CD1B）+（-0.126*EREG）。然后，通過survival包surv_cutpoint函數(shù)確定最佳風(fēng)險(xiǎn)評(píng)分閾值將患者分為高低風(fēng)險(xiǎn)組，繪制K-M生存曲線（圖3A），其中低分組的總生存率更高。同時(shí)，用AUC評(píng)估風(fēng)險(xiǎn)模型的預(yù)測能力，AUC越大，模型預(yù)測能力越好（圖3B）。其中一年生存預(yù)測效能最好（AUC=0.743），三年的預(yù)測效能次之（AUC=0.73），五年生存預(yù)后效能最低（AUC=0.633）。

圖3 TCGA-COAD中5個(gè)基因特征模型的預(yù)后分析。3A：TCGA患者高低風(fēng)險(xiǎn)組OS的K-M曲線；3B：OS的時(shí)間相關(guān)ROC曲線

預(yù)后風(fēng)險(xiǎn)免疫mRNA在癌癥和正常樣本中呈現(xiàn)不同的表達(dá)模式，其中CD1B、EREG在癌癥樣本中傾向于高表達(dá)，而NGF、SCTR、XCL1在正常樣本中傾向于高表達(dá)（圖4）。

圖4 TCGA-COAD中預(yù)后風(fēng)險(xiǎn)免疫mRNA的表達(dá)對(duì)比。

三、獨(dú)立數(shù)據(jù)集驗(yàn)證

在GEO數(shù)據(jù)庫下載了結(jié)腸癌的RNA表達(dá)譜GSE39582。同樣，通過surv_cutpoint函數(shù)確定其最佳風(fēng)險(xiǎn)評(píng)分閾值，將GSE39582中的患者分為高低風(fēng)險(xiǎn)組，然后用Log-rank檢驗(yàn)比較兩組生存差異（圖5A），進(jìn)而驗(yàn)證風(fēng)險(xiǎn)模型的魯棒性。風(fēng)險(xiǎn)模型的預(yù)測能力和TCGA呈現(xiàn)相同的趨勢（圖5B），一年生存預(yù)后最好，三年生存和五年生存次之。由此驗(yàn)證我們的風(fēng)險(xiǎn)模型可以對(duì)結(jié)腸癌患者進(jìn)行較好的風(fēng)險(xiǎn)評(píng)估。

圖5 GSE39582-COAD的驗(yàn)證。5A：GSE3952患者高低風(fēng)險(xiǎn)組OS的K-M曲線；5B：OS的時(shí)間相關(guān)ROC曲線

討 論

本研究整理了免疫基因數(shù)據(jù)和結(jié)直腸癌表達(dá)譜數(shù)據(jù)，通過免疫基因在結(jié)腸癌患者于正常樣本間的差異表達(dá)分析篩選識(shí)別出了疾病相關(guān)的差異表達(dá)免疫基因，進(jìn)一步的分析得到與結(jié)腸癌患者預(yù)后生存相關(guān)的風(fēng)險(xiǎn)免疫基因，構(gòu)建得到結(jié)腸癌預(yù)后風(fēng)險(xiǎn)模型，并在獨(dú)立數(shù)據(jù)集中得到較為準(zhǔn)確的預(yù)測率。此外，我們系統(tǒng)了對(duì)免疫基因表達(dá)情況對(duì)于患者生存影響的潛在機(jī)制進(jìn)行了刻畫，并評(píng)估了其預(yù)后價(jià)值。

考慮到多個(gè)免疫基因?qū)τ诨颊呱婢嬖陲@著性，我們基于多因素COX回歸分析，將多個(gè)基因整合為一個(gè)打分公式，對(duì)結(jié)腸癌患者計(jì)算預(yù)后風(fēng)險(xiǎn)評(píng)分。我們對(duì)所有的統(tǒng)計(jì)計(jì)算結(jié)果P值均進(jìn)行了Bonferroni嚴(yán)格校正，以確保結(jié)果的顯著性。

另外，已發(fā)表的文獻(xiàn)表明CD1B影響前列腺癌進(jìn)展進(jìn)而影響預(yù)后，CD1B在前列腺癌中低表達(dá)與較差無復(fù)發(fā)生存相關(guān)，而在結(jié)腸癌中更傾向于高表達(dá)［17］，體現(xiàn)了癌癥異質(zhì)性［18］。在頭頸部鱗狀細(xì)胞癌（HNSCC）中，EREG表達(dá)上調(diào)預(yù)示預(yù)后不良，并通過激活表皮生長因子受體（EGFR）信號(hào)通路觸發(fā)HNSCC致癌轉(zhuǎn)化［19］。XCL1是一種C類趨化因子，也稱為淋巴趨化素，在感染和炎癥反應(yīng)期間由T、NK和NKT細(xì)胞產(chǎn)生，而XCL1受體XCR1則由樹突狀細(xì)胞亞群表達(dá)。XCL1-XCR1軸在樹突狀細(xì)胞介導(dǎo)的細(xì)胞毒性免疫反應(yīng)中起重要作用。另有研究證實(shí)，XCL1和XCR1在胸腺中組成性表達(dá)，并調(diào)節(jié)胸腺自我耐受的建立和調(diào)節(jié)性T細(xì)胞的生成［20］。

由于數(shù)據(jù)的影響，本文僅對(duì)結(jié)直腸癌免疫基因的表達(dá)情況進(jìn)行了闡述，并構(gòu)建了結(jié)直腸癌預(yù)后風(fēng)險(xiǎn)預(yù)測模型，對(duì)影響結(jié)直腸癌患者預(yù)后生存的biomarker的篩選提供了一個(gè)基礎(chǔ)的工作框架，以期為結(jié)直腸癌患者預(yù)防和預(yù)后提供參考價(jià)值。需要進(jìn)一步研究以揭示這5個(gè)預(yù)后相關(guān)免疫基因的潛在分子機(jī)制，以及這些免疫基因在基因組層面的改變對(duì)于結(jié)直腸癌患者預(yù)后影響，這些免疫基因?qū)τ诮Y(jié)直腸癌可能具有診斷和治療潛力。