基于TCGA 和GEO 数据库构建结肠癌预后模型

操利超，巴 颖，卢晓萍，张核子

（深圳市核子基因科技有限公司，广东 深圳 518071）

结肠癌（colorectal cancer，CRC）是一种常见的恶性肿瘤，是世界上第二大致死原因[1]。尽管结肠癌的诊断和治疗已经取得了很大的进展，但结肠癌患者通常会出现复发和转移，导致5 年生存率显著下降[2]。因此，迫切需要改善结肠癌患者的诊断、治疗和预后。近些年来，分子诊断技术已广泛应用于肿瘤的治疗、预后领域[3-5]。生物信息学和机器学习技术已广泛应用于肿瘤诊断或预后分子标志物的识别，这种分子标志物类型多种多样，如microRNAs[6]、长链非编码RNA[7]、差异表达基因[8]、DNA 甲基化[9]等。其中，差异表达基因作为潜在的肿瘤诊断或预后标志物应用最为广泛。为得到广泛验证的结肠癌相关的差异表达基因，本文利用生物信息学方法，从多个数据集、不同的数据库中寻找共同的结肠癌相关的差异表达基因，并进一步利用机器学习的方法，从这些差异基因中挑选出结肠癌预后相关的预测因子，并建立预后风险评估模型。

1 材料与方法

1.1 数据下载和获取 通过GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）下载基因芯片表达数据集GSE44076、GSE28000 和GSE39582，每个参考数据集的正常和肿瘤样本情况见表1。TCGA 中mRNA表达数据集和对应的临床信息从UCSC Xena 平台（https://xenabrowser.net/datapages/）下载，选择队列为GDC TCGA Colon Cancer（COAD），样本信息见表2。

表1 3 个GEO 数据集的样本量情况

表2 TCGA 数据集的样本信息[n（%）]

1.2 差异基因分析和统计分析 利用R 包分别对3个GEO 数据集和TCGA 数据集进行差异基因分析，过滤标准为adjustedP-value＜0.05 和差异倍数1.5倍（|log2FC|＞0.585），然后取交集，得到共同的上调差异基因和下调差异基因。

1.3 构建和评估预后风险评分模型 为了确定与生存相关的差异表达基因，使用R 包Survival 进行单变量Cox 比例风险回归模型（P＜0.05）。接着，使用LASSO 回归分析进一步缩减预后因子数量，通过多因子回归分析确定每个预后因子的回归系数，建立预后风险评估模型，预测患者生存率。公式为：

1.4 绘制生存曲线和ROC 曲线 根据风险评分预后模型，计算每个肿瘤样本的预后因子风险评分。利用R 包survivalROC 绘制ROC 曲线，用以展示构建的风险评估模型的敏感性和特异性。在ROC 曲线的转折点选择最佳风险评分临界值，转折点处真阳性和假阳性之间的差异最大。高于临界值的患者属于高危组，低于临界值的患者属于低危组。使用未配对t检验估计两组正态分布变量的统计显著性，并使用R 包Survminer 绘制两组的生存曲线。

1.5 构建和验证列线图 为了提高预后模型的性能，通过整合风险评分模型和临床信息，包括年龄、性别和肿瘤分期，可视化不同患者特征的预后价值，构建列线图。该分析使用R 软件包rms 绘制校准曲线，以评估预测概率，并与理想预测线进行比较。此外，基于单因子回归分析的森林图说明了临床信息与OS 之间的关系。其中，一致性指数（C-index）表明了列线图的预测准确性。

2 结果

2.1 差异基因分析 通过比较测试肿瘤样本组和正常样本组，TCGA 数据集和三个GEO 数据集的差异基因数量分布见表3 和图1，可以看到共同的上调差异基因为48 个，共同的下调差异基因数为77 个。

表3 差异表达基因的统计信息

图1 三个GEO 数据集和TCGA 数据集的差异基因数量情况

2.2 结肠癌预后模型的建立 通过单因子回归分析表明，有14 个DEGs 与总生存期（OS）有关（P＜0.05），见表4。进一步LASSO 回归将基因数量缩减为10 个，见图2。根据逐步回归模型，Akaike 信息标准（AIC）为995.94，C 指数为0.63。

图2 LASSO 回归分析结果

表4 与结肠癌预后相关的基因信息

2.3 结肠癌预后模型的性能评估 使用预后模型公式计算每个结肠癌患者的风险评分，然后根据R软件包survminer 中预后因子相关风险评分的最佳临界值（cut-off 为0.16），将结肠癌患者分为高评分组和低评分组。结果显示，随着风险得分的增加，生存时间呈现缩短的趋势，并且高危组的死亡比例比低危组高，高风险评分患者的OS 比低评分患者预后更差，其中基因CILP 和C7 在低风险组表达量低，在高风险组表达量低，而其余8 个基因趋势相反，见图3。

图3 预后风险评分分组和评估

2.4 预后模型的统计分析 为了进一步评估预后风险评分模型的性能，绘制ROC 曲线和肿瘤分层分析。通过将风险评分的预后准确性作为一个连续变量进行研究，OS 预后模型的ROC 曲线下面积（AUC）在3 年时为0.628，4 年时为0.678，5 年时为0.730，见图4。Wilcoxon 检验表明，较高的风险评分与较高的病理分期（P=0.0019）、T 分期（P=0.049）、M分期（P=0.003）、N 分期（P=0.0015）相关。

图4 预后风险评估模型性能评价

2.5 列线图模型的构建与验证 在列线图中，每个变量的得分可以在分数表上找到，然后通过计算总分来估计3 年、4 年和5 年的生存概率，见图5A。森林图显示患者特征，包括年龄（＞60）、肿瘤分期（Ⅲ和Ⅳ）和风险评分与OS 相关（P＜0.05），见图5B。为了验证列线图的性能，绘制校准曲线，可观察到预测曲线接近理想曲线，性能良好，见图5C～图5E。此外，该列线图（C-index：0.74）的预测准确性高于风险评分模型（C-index：0.63）。

图5 列线图模型的构建与验证

3 讨论

本研究中利用3 个GEO 数据集和TCGA 数据集来挖掘共同的差异表达基因。其中，GEO 数据集是基于基因芯片平台，TCGA 数据集是基于二代测序平台，不同数据集交叉验证，使得得到的差异表达基因具有相对广泛适用性。回归分析研究表明，有10 个差异表达基因与结肠癌预后显著相关。其中，SLC4A4 全称Solute Carrier Family 4 Member 4，已有研究报道该基因在结肠癌患者中低表达，与结肠癌癌较差的预后相关，也与淋巴结浸润和远处转移有关[10，11]。CD177 被认为是一种干细胞因子受体，有实验证明CD177 可作为结直肠癌患者对含贝伐单抗的抗癌治疗反应的潜在预测生物标记物[12]。另有研究[13]提出CD177 调节胃癌中的肿瘤细胞粘附和迁移，是生存预后的因素。有报道表明[14]，CD177 的异常表达与结肠癌的发生发展相关。C7 作为一种潜在的肿瘤抑制因子，被报道与前列腺癌的免疫相关预后生物标志物[15]。此外，细胞膜上表达的C7 是过度促炎反应的调节因子[16]，而非小细胞肺癌（NSCLC）中C7 的低表达可能是肿瘤抑制剂，与肿瘤进展和预后相关[16]。UGT2A3 是葡萄糖醛酸转移酶（UGT）家族成员之一，UDP-葡萄糖醛酸转移酶负责外源和内源性化合物的葡萄糖醛酸化，包括药物、环境麻醉剂、类固醇、神经递质、胆汁酸和其他激素。据报道，UGT2A3 主要在与药物清除有关的组织中表达水平最高，其中肝脏是表达最多的器官，其次是胃肠道和肾脏。研究表明[17]，原发性结肠癌肝转移患者的UGT2A3 水平明显高于无肝转移患者。DNASE1L3 与自身免疫性疾病相关，它可以消化凋亡细胞释放的微粒中的染色质，这是一种潜在的自身抗原，它的过度积累将导致身体产生自身免疫反应[18，19]。有研究表明，DNASE1L3 可能是结肠癌免疫浸润的重要生物标志物，并将为结肠癌免疫治疗靶点的选择提供理论依据[20]。HEPACAM2 是粘附基因免疫球蛋白家族的成员，已有研究表明该基因与结肠癌的发生发展有关[21]。ITLN1 可作为多种癌症的肿瘤抑制因子，如胃癌[22]、卵巢癌[23]、神经母细胞瘤[24]和结肠癌[25]。有研究通过免疫组织化学发现[25]，148例大肠癌中87 例（59%）ITLN1 蛋白低表达，ITLN1表达低的结肠癌患者的M 分级高于ITLN1 表达高的结肠癌患者（P=0.0017），且ITLN1 表达高的患者比ITLN1 表达低的患者预后更为良好。MMP3 和MMP10 均属于基质金属蛋白酶（matrix metalloproteinase，MMP）家族，该家族成员参与正常生理过程中细胞外基质的分解，如胚胎发育、生殖和组织重构，以及疾病的发生发展。研究表明[26，27]，MMP3 和MMP10 均可作为结肠癌的预后相关。CILP 基因编码软骨中间层蛋白，在早期骨关节病软骨中增加，目前还未见该基因与肿瘤的相关性报道。

综上所述，本研究基于以上10 个预后相关的基因，构建了结肠癌风险评分模型和列线图模型，结果表明，构建的模型在结肠癌预后和肿瘤分层中表现良好，具有一定的应用价值。