基于免疫基因的RNA-seq数据构建结直肠癌预后生存预测模型
2021-02-19李伟华赵鹏宇黎鸿坚刘林江李运洁张芳邹海军王玉里
李伟华 赵鹏宇 黎鸿坚 刘林江 李运洁 张芳 邹海军 王玉里
结直肠癌是全球五种最常诊断的癌症之一，占所有癌症的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个预后相关免疫基因的潜在分子机制，以及这些免疫基因在基因组层面的改变对于结直肠癌患者预后影响，这些免疫基因对于结直肠癌可能具有诊断和治疗潜力。
