基于免疫基因的RNA-seq數據構建結直腸癌預后生存預測模型

李偉華 趙鵬宇 黎鴻堅 劉林江 李運潔 張芳 鄒海軍 王玉里

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

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

本文基于公開已發表的數據集，通過多種生物信息學分析方法對免疫相關基因在結直腸癌患者中的表達以及功能狀態進行全面刻畫，并篩選得到對結腸癌患者預后生存影響最顯著的5個免疫相關基因構建了Risk Score模型。構建的Risk Score在獨立驗證集中顯示出較高的準確率，證明了該模型的有效性與普遍性。我們的研究結果為預測結直腸癌的預后生存提供了有效的生物標志物模型。

材料與方法

一、數據來源

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

二、數據預處理與差異表達分析

將TCGA下載的結腸癌與正常樣本的mRNA表達譜進行預處理。然后，使用“Deseq2”包［12］對mRNA表達譜進行差異分析，其中|log2FC|＞1且FDR＜0.05的基因被認為是具有統計學意義的差異基因。

三、構建免疫預后特征模型

首先，我們將差異分析得到的差異基因與免疫相關基因取交集得到免疫差異基因集。然后，結合相應的臨床信息進行單因素COX回歸分析，定義P＜0.05的mRNA為候選的免疫相關mRNA。接下來，對其進行“向前”多因素COX回歸分析，評估它們作為獨立預后因素對患者生存的貢獻。因此，我們確定了五個免疫相關的mRNA作為預后特征，并基于其表達水平和多因素COX回歸模型的回歸系數進行線性組合構建預后模型。具體公式如下：

其中，Coefi代表估計的回歸系數，xi代表免疫相關的mRNA表達值。

四、功能富集分析

將識別得到的免疫預后相關的風險mRNA使用R包“ClusterProfiler”［13］分別對其進行基因本體（GO）［14］富集分析，以確定與免疫預后關鍵基因相關的生物過程、分子功能、細胞組分和信號通路，研究其潛在的生物學功能。

五、免疫預后模型的評估與驗證

為評估模型預后能力，根據COAD患者免疫預后風險評分，使用“survival”包［15］比較總體生存時間，繪制Kaplan-Meier（KM）生存曲線，并進行Log-rank檢驗；使用“timeROC”包［16］進行繪制時間相關的受試者工作特征（ROC）曲線，以評估免疫預后模型的有效性。然后，從GEO數據庫中下載獨立數據集GSE39582繪制KM生存曲線和時間相關的ROC曲線來驗證模型的魯棒性，檢驗模型的預后價值。

六、統計學分析

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

結 果

一、結腸癌差異表達mRNA

將從TCGA下載的結腸癌的癌癥與正常樣本進行差異分析，共識別出4 456個差異基因，這些基因由2 264個上調基因和2 192個下調基因組成。根據差異結果做火山圖（圖1A），紅色代表上調基因，藍色代表下調基因，灰色代表差異不顯著的基因。從Immpot中下載免疫相關基因，與差異表達基因取交集，最終獲得362個結腸癌免疫相關基因（圖1B）。

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

二、關鍵mRNA的篩選

為了挖掘出與疾病更為相關的mRNA，結合臨床信息，使用單因素COX對免疫相關基因進行篩選。獲得與患者生存更為相關的關鍵mRNA。在預后mRNA中截取P＜0.05的為關鍵mRNA（表1）。關鍵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個。這32個mRNA都顯示出與結腸癌有著很強的預后相關性。對上述32個mRNA進行功能富集分析（圖2），預后相關基因調控T細胞介導的免疫、T細胞介導的細胞毒性的正調等相關的生物學功能，進而影響患者生存。

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

圖2 預后風險mRNA富集分析

為了找出這32個mRNA中與結腸癌最為相關的免疫預后mRNA，對其進行多因素COX回歸分析。篩選出其中預后更為相關的mRNA以構建模型。選取P＜0.05的mRNA為預后風險mRNA。篩選出的預后風險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個基因，構建風險得分模型，即：風險評分=（0.191*SCTR的表達值）+（0.151*XCL1） +（0.231*NGF） +（ -0.241*CD1B）+（-0.126*EREG）。然后，通過survival包surv_cutpoint函數確定最佳風險評分閾值將患者分為高低風險組，繪制K-M生存曲線（圖3A），其中低分組的總生存率更高。同時，用AUC評估風險模型的預測能力，AUC越大，模型預測能力越好（圖3B）。其中一年生存預測效能最好（AUC=0.743），三年的預測效能次之（AUC=0.73），五年生存預后效能最低（AUC=0.633）。

圖3 TCGA-COAD中5個基因特征模型的預后分析。3A：TCGA患者高低風險組OS的K-M曲線；3B：OS的時間相關ROC曲線

預后風險免疫mRNA在癌癥和正常樣本中呈現不同的表達模式，其中CD1B、EREG在癌癥樣本中傾向于高表達，而NGF、SCTR、XCL1在正常樣本中傾向于高表達（圖4）。

圖4 TCGA-COAD中預后風險免疫mRNA的表達對比。

三、獨立數據集驗證

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

圖5 GSE39582-COAD的驗證。5A：GSE3952患者高低風險組OS的K-M曲線；5B：OS的時間相關ROC曲線

討 論

本研究整理了免疫基因數據和結直腸癌表達譜數據，通過免疫基因在結腸癌患者于正常樣本間的差異表達分析篩選識別出了疾病相關的差異表達免疫基因，進一步的分析得到與結腸癌患者預后生存相關的風險免疫基因，構建得到結腸癌預后風險模型，并在獨立數據集中得到較為準確的預測率。此外，我們系統了對免疫基因表達情況對于患者生存影響的潛在機制進行了刻畫，并評估了其預后價值。

考慮到多個免疫基因對于患者生存均存在顯著性，我們基于多因素COX回歸分析，將多個基因整合為一個打分公式，對結腸癌患者計算預后風險評分。我們對所有的統計計算結果P值均進行了Bonferroni嚴格校正，以確保結果的顯著性。

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

由于數據的影響，本文僅對結直腸癌免疫基因的表達情況進行了闡述，并構建了結直腸癌預后風險預測模型，對影響結直腸癌患者預后生存的biomarker的篩選提供了一個基礎的工作框架，以期為結直腸癌患者預防和預后提供參考價值。需要進一步研究以揭示這5個預后相關免疫基因的潛在分子機制，以及這些免疫基因在基因組層面的改變對于結直腸癌患者預后影響，這些免疫基因對于結直腸癌可能具有診斷和治療潛力。