基于TCGA和GEO数据库建立了肝内胆管癌的预后风险模型及验证分析

毛 俊，沈秀芬，马 润，何 薇，瞿巧莉，胡 莹（昆明医科大学第二附属医院检验科，昆明 650101）

肝内胆管癌（intrahepatic cholangiocarcinoma，ICCA）是仅次于肝细胞癌的第二常见原发性肝癌，约占20%[1]。ICCA 现在有望治愈的唯一途径仍然是手术治疗，而能进行手术治疗的患者数量不到三分之一[1]，并且目前诊断ICCA 通常采取实验室检查、影像学检查和组织活检等[2]，缺乏早期特异性标志物，往往错过手术的最佳时期。且经手术后的大多数患者，其5年生存率仅20%～35%[3]。因不能手术而去接受放、化疗等治疗的ICCA 患者中位生存期仅12.9 个月[1]。因此，确定更精准的胆管癌预后指标对临床制定个体精准医疗方案及评估患者生存预后具有重要意义。本研究基于TCGA 和GEO 数据库，从TCGA 数据库下载TCGA-CHOL数据作为训练集，GEO 数据库下载GSE107943 数据作为验证集，通过差异表达基因筛选，单因素COX 和LASSO 回归等生信分析建立了ICCA 的预后风险模型，并在验证集中进行验证。以期能够更好地应用于临床，辅助评估ICCA 患者的预后，制定针对性的治疗策略。

1 材料与方法

1.1 研究对象 从TCGA 数据库下载胆管癌（TCGA-CHOL）相关数据，包括表达谱数据、临床病理特征和生存信息。最终获得36 例癌组织和9 例癌旁组织，利用perl 脚本进行整理合并，得到行为样本名称，列为基因名的基因表达矩阵。选择ICCA 作为研究对象，共有31 例ICCA 癌组织及9例癌旁组织的相关数据，其中男性13例，女性18例，平均年龄62.45±13.23 岁；15 例报告死亡，16 例报告存活，平均存活时间2.34±1.45年；26 例白种人，3 例黄种人，2 例黑种人；M 分期：26 例M0期，3 例M1 期，2 例M 分期未知（MX）；N 分期：24 例N0 期，4 例N1 期，3 例N 分期未知(NX)；T 分期：17 例T1 期，9 例T2 期，5 例T3 期。临床病理分期：17 例Stage Ⅰ期，8 例Stage Ⅱ期，1例Stage Ⅲ期，5 例Stage Ⅳ期，将其作为训练集。从GEO 数据库下载GSE107943，获取表达谱数据、临床病理特征和生存信息，共获得30 例ICCA 癌组织及27 例癌旁组织，其中男性24 例，女性6 例，平均年龄65.6±8.59 岁；17 例报告死亡，13 例报告存活，平均存活时间2.32±1.71年；乙型肝炎4例，非乙型肝炎26 例；肿瘤大小：6.13±3.21cm；有血管侵入12 例，血管未侵入18 例，临床病理分期：15 例Stage Ⅰ期，6 例Stage Ⅱ期，1 例Stage Ⅲ期，8 例Stage Ⅳ期，作为验证集。

纳入标准：将肿瘤切除部位为肝内胆管的病例资料作为纳入标准；排除标准：①其他肿瘤切除部位包括肝外胆管、胆囊、肝脏等；②生存时间不足30 天的患者。

1.2 方法

1.2.1 筛选差异表达基因：利用R 软件“DESeq2”包对两组数据的癌组织和癌旁组织进行差异表达分析，校正P值（adjustPvalue，padj）＜0.05，|log2 Fold change|＞2 被认为差异有统计学意义。

1.2.2 单因素COX 回归分析：利用“survival”包[4]将得到的差异表达基因纳入单因素COX回归分析，筛选出预后相关的有统计学差异的表达基因。将两组数据中有预后统计学意义的基因取其交集，得到在GEO,TCGA数据中均有预后统计学意义的基因。

1.2.3 构建预后模型：在TCGA数据中利用“survival”包glmnet 函数将在两组数据均有预后意义的差异表达基因进行10 折交叉验证的LASSO 回归分析，构建ICCA 的预后风险模型。并计算每个样本的风险分数，公式如下：

风险评分=∑βx×Expx

βx表示LASSO 回归分析筛选出的各个基因的系数，Expx表示这些基因的表达水平。以风险分数中位数作为截断值将训练集分成高、低风险组。

1.2.4 绘制预后模型的Kaplan-Meier 曲线和受试者工作特征（receiver operating characteristic,ROC）曲线：对高低风险组利用“survival”包绘制Kaplan-Meier 曲线进行生存分析，并利用“timeROC”包绘制时间依赖的ROC 曲线评估模型1，3，5年生存的预测效能。

1.2.5 利用GEO 数据集验证：按上述公式以相同的系数计算验证集各个样本的风险分数并同样以中位数为截断值将验证集划分为高、低风险组，并绘制Kaplan-Meier 曲线和ROC 曲线进行验证。

1.2.6 独立预后因子的筛选：将各个样本病理资料与风险分数进行单、多因素COX 回归分析，以明确影响ICCA 总体生存率的独立危险因素有哪些，并绘制列线图展示。

1.2.7 基因本体论（gene ontology，GO）和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)分析：高、低风险组间差异表达基因将通过“DESeq2”包来获取，padj＜0.05 且|log2 Fold change|＞2 被认为差异有统计学意义，并进行GO 和KEGG 富集分析（“clusterProfiler”包），以P＜0.05 为差异有统计学意义。

1.2.8 基因集富集分析（Gene Set Enrichment Analysis，GSEA）：利用“clusterProfiler”包对高低风险表达组进行GSEA 富集分析，同时下载MSigDB 数据库中的“c2.cp.kegg.v7.5.1.entrez.gmt”基因集。以此基因集中的通路和基因作为参考，显著富集的通路以P＜0.05,FDR＜0.25 为判断标志。

1.2.9 单样本基因集富集分析（Single Sample Gene Set Enrichment Analysis,ssGSEA）：读入文献定义的28 种免疫细胞的参考基因集[5]，利用“GSVA”包ssGSEA[6]对高低风险组进行免疫浸润分析，分析高低风险组之间的免疫浸润差异。以P＜0.05为差异有统计学意义。

1.3 统计学分析 采用RStudio 4.1.2 对数据进行统计分析，利用单因素COX 回归筛选具有预后意义的差异表达基因及预后相关病理信息，多因素COX 回归用以确定影响总体生存率的独立危险因素。LASSO 回归分析选择变量和调整模型复杂度，从而避免过度拟合，并降维构建预后风险模型。Kaplan-Meier 曲线用以评价高低风险组的总体生存率，ROC 曲线评价模型预测效能。以P＜0.05 为差异有统计学意义。

2 结果

2.1 TCGA 与GEO 的差异表达基因 TCGA 数据集癌组织和癌旁组织进行差异表达分析后，共获得2 922 个基因，其中1 753 个基因表达量升高，1 169个基因表达量降低。而GEO 数据集差异表达分析共筛选出3 075 个差异表达基因，其中1 920 个表达上调基因，1 155 个表达下调基因（均P＜0.05）。

2.2 预后相关基因 TCGA 的差异表达基因经单因素COX 回归分析共筛选出68 个预后相关基因（HR=0.13～7.2，均P＜0.05），GEO 共筛选出413个预后相关基因（HR=0.17～215.1，均P＜0.05）。取两者重叠部分，共获得9 个基因，见图1。其中GOLGA6B,PIWIL4,PRICKLE1,FUT4 和COL4A3 风险比（HR）＜1，可能是影响ICCA 患者预后的保护因素。MTFR2，TPM2，EPHX4 和DIO2 的HR＞1，可能是影响ICCA 患者预后的危险因素。

图1 预后相关基因的筛选

2.3 预后风险模型及生存分析 将9 个交集基因纳入LASSO 回归分析，进一步降维筛选，最后在训练集TCGA 数据集中建立了一个6 个基因（MTFR2，TPM2，PIWIL4，PRICKLE1，DIO2，COL4A3）的预后风险模型。风险评分（Risk Score）=0.464×表达量MTFR2 ＋0.550×表达量TPM2-0.511×表达量PIWIL4-0.097×表达量PRICKLE1 ＋0.215×表达量DIO2-0.313×表达量COL4A3，以此计算公式得出训练集中每个样本的风险评分，以得分的中位数1.43 将训练集分成高、低风险组。分析发现，随着得分的增加，ICCA 患者的总体生存期逐渐减少，死亡人数逐渐增加。且MTFR2，TPM2 和DIO2 随着风险分数的增加表达量也逐渐增加，而PIWIL4，PRICKLE1 和COL4A3 表达量则显示出逐渐降低的趋势。表明PIWIL4，PRICKLE1 和COL4A3 可能是影响ICCA 患者预后的保护因素，MTFR2，TPM2和DIO2 可能是影响ICCA 患者预后的危险因素。Kaplan-Meier 曲线表明高风险组总体存活时间少于低风险组（P＜0.001），见图2A；ROC 曲线中预测ICCA 患者1，3，5年总生存期的AUC 值分别为0.971（cutoff=0.22）,0.921（cutoff=2.33） 和0.701（cutoff=1.52），见图2B，表明构建的6 基因模型预测能力良好。

图2 训练集组ICCA 患者生存分析

2.4 预后风险模型的验证 利用与TCGA 训练集一样的评分计算公式得出GEO 验证组各个样本风险分数，并采用中值2.48 为界点使验证集分成高、低风险组，与训练集结果一致，随着风险分数的增加，ICCA 病患的总体生存时间下调，相应的死亡人数上升；且同样能够表明PIWIL4，PRICKLE1和COL4A3 可能是影响ICCA 患者预后的保护因素，MTFR2，TPM2 和DIO2 可能是影响ICCA 患者预后的危险因素。Kaplan-Meier 曲线同样证明两组总体生存时间差异有统计学意义(P=0.004)，高风险组少于低风险组，见图3A。ROC 曲线中预测ICCA 患者1，3，5年总生存期的AUC 值分别为0.908（cutoff=3.23）,0.851（cutoff=1.02） 和0.752（cutoff=2.70），见图3B。表明构建的6 基因模型在验证集中预测能力同样良好。

图3 验证集组ICCA 患者生存分析

2.5 ICCA 患者的独立预后因子 为了进一步评价所建立的预后风险模型的预测价值，将TCGA 训练集风险分数与临床病理特征一同纳入单因素COX 回归分析和多因素COX 回归分析，单因素COX 回归分析结果显示风险分数HR=5.18(95%CI:2.15～12.49,P＜0.001)，多因素COX 回归分析结果显示风险分数HR=72.5(95%CI:4.52～1162.9,P=0.002)，表明预后模型风险分数是影响ICCA 预后的独立风险因素；将可能影响患者预后的临床病理特征与风险分数绘制列线图展示，见图4。

图4 训练集组风险分数及临床病理特征列线图

在GEO 验证集中，风险分数与临床病理特征的单因素COX 回归分析结果显示风险分数HR=2.76(95%CI:1.65～4.60),P＜0.001，多因素COX 回归分析结果显示风险分数HR=4.68(95%CI:2.13～10.3),P＜0.001，证明风险分数在验证集中同样是独立危险因子。绘制列线图见图5。

图5 验证集组风险分数及临床病理特征列线图

2.6 GO 和KEGG 富集分析 为了探究高、低风险组预后差异的原因，分析了高、低风险组差异基因的功能和可能参与的通路，将筛选到的98 个在高、低风险组表达有差异的基因进行GO,KEGG富集分析。GO 富集分析提示高低风险差异基因功能主要富集在机体免疫反应及调节相关方面，如“体液免疫反应”（P＜0.001,校正P=0.005，qvalue=0.004）、“体液免疫调节”（P＜0.001,校正P=0.02，qvalue=0.02）、“免疫反应与激活细胞表面受体信号通路及信号传导”（P＜0.001,校正P=0.02,qvalue=0.02）等。KEGG 结果与GO 富集结果遥相呼应，都与机体免疫相关，包括“B细胞受体信号通路”（P=0.006,校正P=0.24，qvalue=0.23）、“原发性免疫缺陷”（P=0.01,校正P=0.33,qvalue=0.32）等。结合GO，KEGG富集结果，表明高低风险组之间的差异可能是由于机体的免疫反应及免疫调节导致的。

2.7 GSEA 富集分析 为了进一步验证高、低分险组是否通过影响免疫反应及免疫调节来导致两组间的预后差异，进行了GSEA 富集分析，结果显示高低风险组基因富集于“B 细胞受体信号通路”（富集分数=-0.57，NES=-2.10，P＜0.001,校正P值＜0.001，FDR＜0.001）、“原发性免疫缺陷”（富集分数=-0.77，NES=-2.44，P＜0.001,校正P值＜0.001，FDR＜0.001）等与KEGG 结果一致，此外还富集在“趋化因子信号通路”（富集分数=-0.52，NES=-2.29，P＜0.001，校正P值＜0.001，FDR＜0.001）、“T 细胞受体信号通路”（富集分数=-0.46，NES=-1.81，P＜0.001,校正P值=0.001，FDR=0.001）等免疫相关信号通路，并且这些免疫信号通路在高风险组均处于抑制状态（富集分数＜0）。这提示高分险组可能通过抑制机体的免疫反应来促进肿瘤的恶性进展，进而导致高风险组的预后不良。

2.8 ssGSEA 免疫浸润分析 在确定高、低风险组的预后差异主要由于影响机体免疫反应及免疫调节的前提下，对高、低风险组进行了ssGSEA 免疫浸润分析，以下载的28 种免疫细胞基因集为参照对象，评估这些细胞在高低风险组的分布及差异情况。依据CHAROENTONG 等[5]人定义的基因集将28种免疫细胞分为三类，分别为抗肿瘤免疫细胞、促肿瘤免疫细胞和功能不明确的免疫细胞。ssGSEA结果表明促肿瘤免疫细胞中CD56-自然杀伤细胞在高风险组表达量相对于低风险组更高（富集分数：高风险组0.31～0.46，低风险组0.31～0.42，P＜0.05），而抗肿瘤细胞中CD4 中心记忆型T 细胞则表现出相反的结果（富集分数：高风险组0.51～0.59，低风险组0.44～0.62，P＜0.05），此外，激活的B 细胞，激活的CD8+T 细胞也表现出在高风险组表达抑制的状态，见图6。这进一步验证了富集分析的高风险免疫抑制状态。

图6 ssGSEA 免疫浸润分析

3 讨论

本研究基于TCGA 和GEO 数据库构建了一个6 基因预后风险模型，包括MTFR2，TPM2，DIO2，PIWIL4，PRICKLE1 和COL4A3，模型的质量也在GEO 数据集中得到了验证。其中MTFR2和TPM2 的高表达与多种肿瘤的不良预后有关，包括肝细胞癌[7]、肺腺癌[8]、乳腺肿瘤[9-11]和口腔癌[12]等，DIO2 在浆液性卵巢癌肠道转移中过表达，表明DIO2 可能参与肿瘤的恶性转移[19]。此外，已有研究表明，PIWIL4 能作为ICCA 的预后生物标志物[13]，且其低表达与肾透明细胞癌患者的不良预后有关[14]。COL4A3 的表达与鼻咽癌[15]和胃癌[16]的预后呈负相关关系。PRICKLE1 的高表达减少了神经母细胞瘤的生长，并与患者的预后呈负相关关系[17]。以上研究表明这6 个基因在ICCA 中发挥的作用和其他肿瘤大致相同，表明研究建立的预后风险模型有望成为ICCA 患者预测性能良好的预后生物标志物。

通过GO，KEGG，GSEA 和ssGSEA 分析提示ICCA 患者的预后可能与机体免疫反应有关，这与黄建斌[18]的结论是一致的。而MTFR2，TPM2，DIO2，PIWIL4 和COL4A3 均与机体免疫密切相关，例如：MTFR2 的表达在胃癌中与B 细胞、CD8+T 细胞、CD4+T 细胞、巨噬细胞、中性粒细胞和树突状细胞的渗透水平显著相关[26]，TPM2 与膀胱癌巨噬细胞和NK 细胞浸润呈正相关[27]，DIO2在炎症诱导的巨噬细胞中表达升高，并可能成为参与COPD 机制过程的标志物[28]，PIWIL4 在ICCA中与静息自然杀伤细胞和激活的记忆CD4+T 细胞的富集呈正相关[20]，COL4A 家族基因的表达在肾透明细胞癌中与免疫细胞、肿瘤浸润淋巴细胞、MHC 分子、趋化因子和受体的浸润显著相关[29]。另一方面，机体免疫与肿瘤的关系主要包括免疫监视和肿瘤免疫逃逸。ssGSEA 免疫浸润分析部分指出，肿瘤可以驱动免疫抑制或调节性免疫细胞亚型的产生，也可以招募大量促肿瘤的髓系细胞来建立肿瘤微环境，从而促进肿瘤的进展。肿瘤恶性进展的关键步骤是逃避免疫破坏和启动肿瘤细胞转移，这些步骤可以通过抑制宿主的免疫系统来实现，特别是通过诱导、扩增和重新招募免疫抑制细胞[30]。由此推测ICCA 可能通过免疫逃逸机制来实现肿瘤恶性进展，从而导致ICCA 患者预后不良。早期研究已经表明，ICCA 能够通过B7-H1/PD-1 通路促进CD8+肿瘤浸润性淋巴细胞（CD8+TILs）的凋亡而参与肿瘤免疫逃逸[19]。SURIYO 等[20]的研究也表明ICCA 能够通过PD-L1/PD-1 轴抑制T 细胞介导的免疫反应，使CD8+T 细胞凋亡增加,实现肿瘤细胞的免疫逃逸。而CD8+T 细胞数量的减少与ICCA患者的总体生存时间短密不可分[21]。此外，有学者发现sPDL1 相对于其他类型胆管癌，在ICCA 中表达最高，sPDL1 低水平患者表现出较长的生存时间，表明sPDL1 高水平会增加患者死亡风险[22]。目前临床上针对免疫检查点，已经开发出免疫检查点抑制剂，如PD-1 抑制剂和PD-L1 抑制剂。且在ICCA的治疗中显示出较高的有效率[23]。因此，从肿瘤免疫出发，针对性进行治疗能使ICCA 患者总体生存期升高，这具有重要的临床意义。

综上所述，本研究建立的风险预后模型对评估ICCA 患者预后具有较高的临床应用价值，且对于高风险组患者的针对性疗法具有指导意义，为利用免疫检查点抑制剂来提高高风险组患者预后开创了新平台，但这仍需要大量的临床研究去论证。在未来，有望能够利用免疫抑制剂联合靶向药物进行综合治疗，这或许在胆管癌中能获得更好的治疗效果。

本研究同样存在一定的局限性。首先，本研究是通过在线公共数据库中的数据构建的预后模型，未使用本院ICCA 患者病理标本及病理资料进行模型验证。其次，对组成预后模型的6 个基因的作用及机制还需进行进一步的研究。