1. 引言
食管癌(Esophageal cancer, ESCA)是一种常见的消化系统恶性肿瘤，在中国因癌症死亡的病例中标化死亡率排第五位，在我国超过全球死亡总例数30%的癌种中位列第一[1]。由于食管癌的早期临床症状并不明显，大于50%的患者确诊时已经处于疾病发展的中晚期，导致其治疗效果和预后水平较差。此外，食管癌患者的死亡率逐年升高，患者的5年生存率不足10%，明显低于其他恶性肿瘤。目前根治性切除和纵隔淋巴结清扫术依旧是该病临床治疗的首选方法，但部分患者术后淋巴结病理仍呈阳性，需进行术后辅助化疗，预后效果并不理想。对于食管癌患者，早期确诊和及时手术利于预后，而淋巴结转移不利于预后，有研究表明，当食管癌发展至粘膜下浸润 > 1000 μm时，淋巴结转移的风险会显著增加。因此，早期诊断对于提高食管癌患者的生存率十分重要[2]-[4]。本研究基于基因表达数据库(Gene Expression Omnibus, GEO)和癌症基因组图谱(The Cancer Genome Atlas, TCGA)数据库，查找出与食管癌发病相关的基因后，进行多种机器学习算法(Boruta、Lasso、XGBoost、随机森林和支持向量机)后进行多组学分析，评估了ADAM12、CTHR1、IBSP、OLR1四个基因对食管癌发生发展的机制，旨在寻找食管癌的免疫学标志物，为其早期治疗提供帮助，为后续的研究提供思路。
2. 资料与方法
2.1. 数据集的处理和差异表达基因分析
2.1.1. GEO、TCGA数据的获取
从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)和TCGA数据库(https://portal.gdc.cancer.gov)中以“Esophageal cancer”为关键词检索，获得食管癌组织基因表达数据和正常食管组织基因表达数据。
2.1.2. 差异基因筛选
选用基因芯片GSE149609，使用R语言再次标准化数据，并采用limma包分别对数据集进行差异分析，按照p < 0.05的标准筛选上调和下调的差异基因，最后进行可视化作图展示差异基因。
2.2. 机器学习进行基因进阶筛选
2.2.1. Boruta机器学习
(1) 创建差异基因的副本后，打乱重复副本的值，用以消除其与目标变量的相关性，即为影子特征。
(2) 将原副本与洗过的副本合并后，使用随机森林方法，比较差异基因每个变量的重要性，进一步消除不相关的特征。其中，重要性显著高于影子属性的属性会被标记为“重要”的，显著低于影子属性的属性会被标记为“不重要”。
2.2.2. Lasso机器学习
采用glmnet包来执行Lasso回归，指定family参数为binomial后，设置弹性网络中的混合参数alpha参数是1，即在损失函数中加入一个L1正则化项。选用cv.glmnet函数，deviance为度量指标执行十折交叉验证，从而实现特征选择。
2.2.3. XGBoost机器学习
借助xgb.DMatrix函数将数据转换为XGBoost的内部数据结构。在进行模型训练时，选定xgboost函数，将目标函数确定为“binary:logistic”，用于二分类问题并采用逻辑回归损失，同时将迭代次数设置为25。随后，运用xgb.importance函数进行特征选择。最后，选用xgb.ggplot.importance函数绘制特征重要性图。
2.2.4. 随机森林
(1) 对原始数据进行预处理后，通过自主采样的方式抽出多个数据集后进行特征随机选择，进一步增加多样性。
(2) 使用数据集构建决策树，Out-of-Bag (OOB)误差率进行交叉验证
(3) 选择用佳树的数量训练随机森林模型后，使用importance函数计算特征重要性，得出特征在模型中的重要性得分，得分越高该特征对模型的影响越大。选取前20%的特征基因。
2.2.5. 向量机递归特征消除(SVM-RFE)算法
基于支持向量机进行机器学习，设定交叉验证的折数为10，对数据集进行10折划分后，对每一折应用SVM-RFE算法，获取特征的重要性排名，对所有折的排名结果进行整合后，计算平均排名，根据平均排名对特征进行排序。
2.3. 通路分析
采用limma包执行差异分析，clusterProfiler包执行基因集富集分析后，选取分别在高/低表达组显著富集的top5通路进行可视化。
2.4. 受试者工作特征曲线(Receiver Operating Characteristic Curves Analysis, ROC)分析
使用pROC包对差异基因进行ROC分析，计算95%置信区间，获得曲线下总面积，使用pROC的ci函数评估AUC，绘制平滑ROC曲线。
2.5. 生存分析
2.5.1. 生存风险比Meta分析
采用逆方差法，选择对数HR的值为主要指标，对单因素cox生存分析的结果进行Meta分析，标准误差使用置信区间95% CI计算后，使用R (4.3.2版本0)的“Meta”包进行统计分析和可视化。
2.5.2. Kaplan-Meier生存分析
采用survival包进行Kaplan-Meier生存分析，P < 0.05为有统计学意义。使用survfit函数进行log rank test评估高、低表达组的显著性。
2.6. 单细胞测序
对食管癌的单细胞测序数据进行整合，分析目标基因在食管癌中单个细胞水平层面的表达。
3. 结果
3.1. 基因的差异分析
对GSE149609中的基因进行分析，其中红色为上调，蓝色为下调(图1)。以20个基因为一组进行筛选，最终选定一基因组中的四个基因：ADAM12、CTHRC1、IBSP和OLR1。
Figure 1. Differential expression analysis of the GSE149609 dataset
图1. GSE149609数据集的差异表达分析
使用GTEx的正常样本TPM表达量与ADAM12、CTHRC1、IBSP和OLR1四个基因的TCGA肿瘤TPM表达量进行配对后，通过(x − μ)/σ将数据转化为无单位的Z-Score分值，使得数据标准统一化。当z-score大于3.0或小于−3.0时，可归类为离群值予以去除。采用Wilcoxon Rank Sum Tests比较四个DGEs与GTEx正常组织之间的表达量统计学差异，并绘制表达箱线图(图2)，其中四个基因的差异均具有统计学意义(p < 0.05)。
Figure 2. Box plots of differential analysis between DGEs and normal tissues. A: Box plot of ADAM12; B: Box plot of CTHRC1; C: Box plot of IBSP; D: Box plot of OLR1
图2. DGEs与正常组织的差异分析箱线图。A：ADAM12的箱线图；B：CTHRC1的箱线图C：IBSP的箱线图；D：OLR1的箱线图
3.2. 机器学习
3.2.1. Boruta
Figure 3. Boruta machine learning
图3. Boruta机器学习
图中y轴表示每个特征的重要性分数，分数越高表明该特征在模型中的重要性就越高。其中ADAM12、CTHRC1、IBSP和OLR1均为青色，表示已被Boruta算法确认为与预测变量显著相关的“已确认”特征(图3)。
3.2.2. Lasso
将选定基因组中的20个DGEs纳入Lasso回归模型后，使用plot函数绘制十折交叉验证结果，其中第一条线对应于交叉验证误差最小的λ值，为最优的λ值(图4(A))。在最优的λ值下，ADAM12、CTHRC1、IBSP和OLR1在模型中有较大影响(图4(B))。
Figure 4. A: Cross-validation curve of Lasso regression; B: Coefficient path of Lasso regression
图4. A：Lasso回归交叉验证曲线；B：Lasso回归系数路径
3.2.3. XGBoost
Figure 5. Evaluation of XGBoost variable importance
图5. XGBoost变量重要性评价
将基因组中的基因纳入xgboost模型，使用xgb.importance函数获取特征的重要性评估，xgb.ggplot.importance函数绘制XGBoost模型的特征重要性图(图5)，获得了12个重要特征，包括：IBSP、CTHRC1、ADAM12、OLR1、APOC1、NETO1、FNDC1、PLPP4、CHI3L1、POSTN、SHISA2、INHBA，其中IBSP得分最高。
3.2.4. 随机森林
绘制随机森林的Out-of-Bag误差率随树数量变化的图形，直观地展示模型误差率的变化趋势(图6(A))。选取最优随机树量后，对影响因素的重要性进行排序，鉴定了重要性前20%的特征，包括ADAM12、CTHRC1、IBSP和OLR1 (图6(B))。
Figure 6. A: Variation plot of random forest error rate; B: Importance ranking plot
图6. A：随机森林误差率变化图；B：重要性排序图
3.2.5. SVM-RFE
Figure 7. Average ranking plot
图7. 平均排序图
采用支持向量机递归特征消除(SVM-RFE)算法在基因表达数据中识别并排序最重要的特征(图7)，重要性top10的特征基因如下：LRRC15、IBSP、ADAM12、APOC1、OLR1、FBN2、CTHRC1、POSTN、SHISA2、SFRP4。
3.3. 通路分析
将表达最高的30%的样本定义为高表达组，表达最低的30%的样本定义为低表达组，对选中的四个基因进行通路分析(图8~9)。不同的颜色代表不同的基因集，柱状图方向为左，意味着在低表达组显著富集，向右则意味在高表达组显著富集(图10)。结果提示ADAM12与CTHRC1的Angiogenesis、EMT、Invasion通路，IBSP的EMT、Invasion通路，OLR1的Angiogenesis、Apoptosis、Differentiation、EMT、Hypoxia、Inflammation、Invasion、Metastasis、Proliferation、Quiescence、Stemness通路呈显著高表达富集趋势。
Figure 8. Pathway analysis plots of ADAM12 and CTHRC1
图8. ADAM12通路分析图与CTHRC1通路分析图
Figure 9. Pathway analysis plots of IBSP and OLR1
图9. IBSP通路分析图与OLR1通路分析图
Figure 10. Multiple gene set enrichment analysis of high/low expression groups performed by clusterProfiler package. A. Enrichment analysis of ADAM12; B. Enrichment analysis of CTHRC1; C. Enrichment analysis of IBSP; D. Enrichment analysis of OLR1
图10. clusterProfiler包执行高/低表达组的多个基因集富集分析。A. ADAM12富集分析；B. CTHRC1富集分析；C. IBSP富集分析；D. OLR1富集分析
3.4. ROC分析
ROC曲线分析提示，采用ADAM12、CTHR1、IBSP、OLR1的表达预测食管癌疾病组与正常组均具有很高的准确性，AUC值分别为0.971、0.974、0.938、0.967，95%置信区间分别为0.958~0.983、0.963~0.983、0.916~0.957、0.947~0.982 (图11)。
Figure 11. ROC analysis evaluating the diagnostic efficacy of gene expression for distinguishing tumor groups from normal groups
图11. ROC评估基因表达对肿瘤组与正常组的诊断效能
3.5. 生存分析
Figure 12. Kaplan-Meier survival analysis plot of ADAM12
图12. ADAM12Kaplan-Meier生存分析图
对ADAM12、CTHR1、IBSP和OLR1进行3个生存期(DSS，PFI和DFI)的Kaplan-Meier生存分析(图12~15)，红色为基因高表达组，蓝色为低表达组，p < 0.05即为显著。四组中p均小于0.05。结果表明，食管癌患者中，ADAM12、CTHR1、IBSP和OLR1低表达组相比于高表达组预后更加良好。
Figure 13. Kaplan-Meier survival analysis plot of CTHRC1
图13. CTHRC1 Kaplan-Meier生存分析图
Figure 14. Kaplan-Meier survival analysis plot of IBSP
图14. IBSP Kaplan-Meier生存分析图
Figure 15. Kaplan-Meier survival analysis plot of OLR1
图15. OLR1 Kaplan-Meier生存分析图
采用逆方差法对上述生存分析的结果进行Meta分析，其中对数HR值为主要测量指标，标准误差使用95% CI (置信区间)计算(图16)。结果表明ADAM12、CTHR1、IBSP、OLR1的OS、DSS和DFI相关(p < 0.05)，且HR均>1，均属危险因素，荟萃分析的结果与之一致，不同的生存期的结果存在较低的异质性。
Figure 16. Meta-analysis plot of survival hazard ratios. A. ADAM12 analysis plot; B. CTHRC1 analysis plot; C. IBSP analysis plot; D. OLR1 analysis plot
图16. 生存风险比Meta分析图。A. ADAM12分析图；B. CTHR1分析图；C. IBSP分析图；D. OLR1分析图
3.6. 单细胞分析基因的表达
Figure 17. Single-Cell analysis plot of ADAM12 (Top Left), Single-Cell analysis plot of CTHRC1 (Top Right), Single-Cell analysis plot of IBSP (Bottom Left), Single-Cell analysis plot of OLR1 (Bottom Right)
图17. ADAM12的单细胞分析图(左上)，CTHR1的单细胞分析图(右上)，IBSP的单细胞分析图(左下)，OLR1的单细胞分析图(右下)
对数据集GSE149609进行单细胞分析。在分析之前首先对细胞进行质量控制，如消除批次效应，以获得质量较好的细胞，便于后续的分析。三个细胞系免疫细胞、恶性细胞和基质细胞中包括13种细胞群：B细胞、CD4常规T细胞、CD8+T细胞、树突状细胞、内皮细胞、胚胎成纤维细胞、恶性肿瘤细胞、肥大细胞、中性粒细胞、周细胞、浆细胞、T细胞和调节性T细胞。
对ADAM12、CTHR1、IBSP、OLR1进行单基因表达UMAP定位后，对其在不同细胞中的表达和细胞比例进行差异分析(图17)。结果显示，ADAM12、CTHR1、IBSP和OLR1在基质细胞、恶性细胞和胚胎成纤维细胞中均表达显著，同时ADAM12在免疫细胞与调节性T细胞表达较为显著，CTHR1在内皮细胞中表达较为显著，IBSP与OLR1在周细胞中的表达最为显著。
结果表明四个基因均在恶性细胞与胚胎成纤维细胞中显著表达，当四个基因中任一基因的表达呈现升高趋势时，提示可能有食管癌的发生。
4. 讨论
食管癌(Esophageal cancer, ESCA)是最常见的恶性肿瘤之一，在2020年的相关调查中，ESCA的发病率于所有癌症中位列第七，总死亡率居第六，这表明每18例因癌症死亡的病例中，就有1例是死于食管癌的，此外，发展中国家人群的发病率显著高于发达国家，特别是中国，根据2018年的统计，我国食管癌的发病率位居榜首，全世界半数的食管癌患者均在我国[5]。食管癌分为食管鳞状细胞癌(Esophageal squamous cell carcinoma, ESCC)和食管腺癌(Esophageal adenocarcinoma, EAC)两种亚型，食管鳞状细胞癌是我国食管癌的主要病理类型，占所有食管癌患者中的90%，这可能与我国部分地区的饮食习惯如大量咀嚼槟榔、食用腌制食物和较烫的食物，遗传变异、基因和环境的相互作用等有关[6]。
食管癌的临床症状包括进食后梗噎感、胸骨后疼痛、进行性吞咽困难、背痛等，早期症状隐匿，通常并不明显，当出现明显的不适时病情往往已进展至中晚期，肿瘤狭窄引起的吞咽困难是晚期食管癌患者的主要症状，吞咽困难程度与食管癌根治术后患者的生存率息息相关，因此早期诊断与治疗尤为重要[7]。目前临床上根治性切除术和淋巴结清扫术是治疗食管癌的主要方式，在尚无用于临床早期诊断ESCA的生物标志物之时，准确评估患者的淋巴结状态对手术预后有重要意义[8] [9]。AJCC/UICC发布的第8版TNM分期系统仍然是全球公认的食管癌分期标准，但是由于缺乏明确的食管癌手术指南，这种分期方法容易出现偏倚，对患者的生存期造成影响。其他的淋巴结分期系统还包括N分期、淋巴结率(LNR)和淋巴结阳性对数几率(LODDS)等，其中LODDS系统不仅关注受累淋巴结的数目，还可以分析阴性淋巴结的数目，在多种肿瘤中均显示出优秀的预测价值，但其作用主要是对于食管癌患者在淋巴结转移数低于12枚的情况下进行分期，在淋巴结并未受侵袭的早期患者身上并不适用，并且对其的研究还较为局限，没有确定的分类标准，在临床上的适用性仍存在争议[10]-[12]。在过去的临床食管癌诊断中常常用到CT和X线，它们不仅可以显示食管腔内的情况，还可以直接地观察到病变组织的黏膜转变情况和是否存在充盈和缺损的问题，但其诊断的准确性相对较低，日常实践中食管癌诊断的假阳性率和假阴性率高，并不可靠。早期食管癌浸润黏膜层或黏膜下层时，可通过内镜下黏膜切除术(endoscopic mucosal resection, EMR)和内镜黏膜下剥离术(endoscopic submucosal dissection, ESD)进行根治[13]，尤其是ESD，可以内镜下完全剥离病灶，对周边组织和血管的损伤少，不仅能降低对机体造成的损伤，还能降低因为病灶的残留导致复发的可能性，有研究表明接受ESD治疗的患者中引发并发症和复发的患者仅为5.36%，但当疾病进展到中晚期甚至发生远处转移时，单纯手术效果不佳，患者甚至需要进行食管切除术与全身放化疗，许多患者由于化疗药物的毒性作用，为了延长生存期只能减少化疗药物剂量或提前终止化疗，但同时可能导致治疗效果不佳，反而增加了风险[14] [15]。因此，探究可用于食管癌早期诊断的标志物成为当下研究的重点方向[16]。
本研究基于生物信息学，通过GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)和TCGA数据库
(https://portal.gdc.cancer.gov)获得了食管癌组织基因表达数据和正常食管组织基因表达数据，筛选出癌组织与正常组织表达明显差异的基因组。进一步采用Boruta机器学习、Lasso机器学习、极端梯度提升(eXtreme Gradient Boost)、随机森林(RF)和向量机递归特征消除(SVM-RFE)算法五种机器学习算法(ML)进行筛选，其中XGBoost模型是基于一组决策树的ML算法，使用梯度提升框架，通过梯度下降算法最大限度地减少错误，是硬件和软件优化技术的理想组合，可使用最少的计算资源获得卓越的结果；RF 模型可以减轻训练变化并增强模型泛化和集成，是一种ML分类器，可通过多个树进行训练和样本预测，得到最重要的影响因素；SVM是分类和回归问题的线性模型，可以解决线性和非线性问题，并将数据分成几类，最终得到特征的重要性排名[17]。进一步筛选出了ADAM12、CTHR1、IBSP、OLR1四个基因，进行通路分析、受试者工作特征曲线(receiver operating characteristic curves analysis, ROC)分析、生存分析与单细胞分析来对基因进行验证，最终得到了运用ADAM12、CTHR1、IBSP和OLR1的表达预测食管癌疾病组与正常组均具有很高的准确性，当四个基因表达显著升高时，可能提示有食管癌的发生。
综上所述，本研究指出了ADAM12、CTHR1、IBSP、OLR1在食管癌的早期诊断和预后方面有重要意义，可能成为食管癌诊断的新型标志物。但本研究的结论还需要实验与检测临床样本来进行进一步的验证。
基金项目
河北北方学院大学生创新创业训练计划(项目编号：XJ10092027)。
NOTES
*共同第一作者。
#通讯作者。