基于生物信息学胰腺腺癌关键基因的筛选及支持向量机诊断模型的构建

张波，徐涛，徐浩，夏雨，周文策，

（1.兰州大学第一临床医学院，甘肃 兰州 730000；2.兰州大学第一医院 普通外科，甘肃 兰州 730000）

胰腺腺癌（pancre atic a denocarcinoma，PAAD）是最常见的胰腺肿瘤，占所有胰腺癌的绝大多数；由于其早期诊断困难，超过52%的患者发现时已有远处转移，大多数患者确诊在晚期，丧失了手术机会，对放化疗敏感度又差，5年生存率不足10%[1-3]。因此，寻找PAAD早期诊断和治疗的新靶点是近年来的研究热点。

基因改变对PAAD 的发病至关重要，研究PAAD发生的分子机制是延长患者总生存期的关键。然而，目前尚未发现对胰腺癌敏感度强、特异度高的肿瘤标记物，也没有找到有效的治疗靶点[4]。单基因对肿瘤的诊断效果多不理想，局限性较大，多基因联合检测有望在诊断领域开辟新方向。支持向量机（support vector machine，SVM）是机器学习领域中一种重要的分类算法，它通过估计每个样本的内部联系和规则来判别样本类型，具有较高的分类精度[5]。随着生物信息技术的发展，基因芯片可快捷、高效地分析多个基因联合检测对肿瘤的诊断效能。为了探究PAAD的发病机理，本研究从基因表达数据库（Gene Expression Omnibus，GEO）中下载了3个芯片，筛选PAAD组织和正常组织之间的差异表达基因（differentially expression genes，DEGs），对这些DEGs 进行基因本体论（Gene Ontology，GO）富集和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）信号通路富集分析，并通过蛋白互作网络（proteinprotein interaction network，PPI）分析和关键子网络分析筛选出候选基因，探讨可能与PAAD诊断和预后相关的关键基因，并构建分类PAAD和正常样本的SVM模型，为以后研究PAAD的诊治和分子机制提供理论依据。

1 资料与方法

1.1 原始数据下载及预处理

从GEO数据库中下载3个PAAD基因表达谱芯片（GSE28735、GSE62165、GSE62452）；数据集GSE28735和GSE62452基于GPL6244 Affymetrix Human Gene 1.0 ST Array，而GSE62165芯片数据在GPL13667 Affymetrix Human Genome U219 Array平台上注释。其中GSE28735包括45例PAAD样本和45例正常组织样本；GSE62452含有69例PAAD样本和61例正常组织样本；GSE62165包括118例PAAD样本和13例正常组织样本。整理数据集中的样本分组及临床相关信息，并对表达谱按芯片注释信息进行重注释。从UCSC Xena下载得到癌症基因组图谱（The Cancer Genome Atlas，TCGA）数据库中胰腺癌的生存数据和基因表达谱，并对表达谱进行Log2（count+1）的标准化处理，以消除量纲差异。

1.2 DEGs 的筛选

芯片数据预处理后，采用R语言Limma包[6]进行PAAD组织和正常组织的差异基因分析，筛选标准为|log FC|＞1和P＜0.05；通过Venn图得到3个芯片交集的DEGs，其中在3个芯片中均上调的基因作为上调DEGs，均下调的基因作为下调DEGs。

1.3 DEGs 的GO 分析和KEGG 分析

将得到的DEGs上传至STRING数据库（http://www.string-db.org）中进行GO分析（包括生物学过程BP、细胞组分CC、分子功能MF）和KEGG通路富集分析，分析参数采用数据库默认，并将TOP10的显著富集结果以气泡图进行可视化。

1.4 DEGs 的PPI 分析及关键子网络分析

利用STRING数据库对DEGs进行PPI分析[7]，蛋白互作的combined score＞0.4作为阈值条件，再通过Cytoscape软件（http://cytoscape.org）对该网络进行可视化[8]。利用Cytoscape软件中的MCODE插件识别该PPI网络中的关键子网络，参数设置为：Degree Cutoff=2，Node Score Cutoff=0.2，Max.Depth=100，K-Core=2，MCODE score≥5。

1.5 关键节点的筛选及其功能富集分析

数据集GSE28735和GSE62452的临床信息中包含生存信息，将所有关键子网络中包含的基因以及PPI分析结果中度值（degree）≥15的基因，分别在两个数据集中使用R 语言的survival 包进行生存分析，将在两个数据集中都具有显著生存意义（P＜0.05）的基因作为关键节点。将得到的关键节点上传至Metascape数据库[9]进行功能富集 分析。

1.6 最优特征基因的筛选和及其表达验证

递归特征消除（recursive feature elimination，RFE）是一种贪婪的算法，通过重复构建模型筛选出分类的最佳基因组合[10]。为了进一步缩小候选基因的范围，准确识别最优特征基因，将数据集GSE28735作为训练集，其他两个数据集作为验证集。在训练集中，利用R语言caret包中的RFE算法从关键节点里筛选出最优特征基因组合。在10倍交叉验证中，以均方根误差（RMSE）最小、准确率最高的基因组合为最佳基因组合。由于TCGA 数据库中正常胰腺样本较少，故在GEPIA[11]数据库对最优特征基因的差异表达进行验证，筛选条件为|log2FC|＞1和P＜0.01。

1.7 SVM 模型的构建及验证

为探索8 个最优特征基因联合检测在分类PAAD和正常样本中的作用，使用 R 语言e1071 包[12]构建基于这些基因的SVM模型，最终选择参数为SVM-Type：eps-regression；SVM-Kernel：radial；cost：1；gamma：0.125；epsilon：0.1，并在训练集GSE28735 和验证集（GSE6216 5 和GSE62452）中采用R 语言 pROC包[13]绘制受试者工作特征（receiver operating characteristic，ROC）曲线对该模型进行验证。

1.8 最优特征基因的预后分析

将TCGA数据库中的PAAD的样本挑选出来，剔除正常样本，对缺乏生存信息的样本删除，使用R语言survminer包[14]计算每个最优特征基因的最佳截断值（optimal cutoff），对于mRNA表达值＞optimal cutoff视为高表达，≤optimal cutoff作为低表达，结合生存时间和生存状态进行生存分析，以P＜0.05为差异有统计学意义。

2 结 果

2.1 DEGs 的筛选

通过对GSE28735、GSE62452、GSE62165数据集的分析，分别鉴定出388个（243个上调，145个下调）、283个（185个上调，98个下调）和3063个（1993个上调，1070个下调）差异基因。韦恩图分析发现257个DEGs均在3个数据集中差异表达，其中168个为上调DEGs，89个是下调DEGs（图1）。

图1 DEGs 的Venn 图 A：上调DEGs；B：下调DEGsFigure 1 Venn diagram of DEGs A: Up-regulated DEGs; B: Down-regulated DEGs

2.2 DEGs 的功能和通路富集分析

GO富集分析结果显示：257个DEGs主要参与BP的细胞外基质的组成、细胞外结构组成、细胞黏附、组织发育等过程；CC主要聚集于细胞外区、细胞外区域部分、细胞外间隙等；MF主要与丝氨酸肽酶活性、丝氨酸内肽酶活性、肽酶活性等有关。KEGG通路分析发现，DEGs主要参与蛋白质的消化和吸收、胰腺的分泌、细胞外基质受体相互作用、黏着斑、PI3K-Akt信号通路等（均P＜0.05）（图2）。

2.3 PPI 及关键子网络分析

通过STRING数据库构建257个DEGs的PPI网络，使用Cytoscape进行可视化，该PPI网络共有210个节点和822条边（图3）。在PPI 网络中筛选出29个度值≥15的重要基因（表1）。利用MCODE插件筛选出4个关键子网络（图4）。结合关键子网络，发现一些新的在PAAD的进展中起调控作用的基因，如COL5A2、OAS2、DDX60、CELA2A，这为以后研究PAAD的分子机制提供更多的依据。

图2 257 个DEGs 的功能和通路富集气泡图Figure 2 Function and pathway enrichment bubble graph of 257 DEGs

图3 PPI 图（红色为上调基因，绿色为下调基因）Figure 3 PPI (red color showing the up-regulated genes, and green color showing the down-regulated genes)

表1 筛选出度值≥15 的29 个基因Table 1 29 screened genes with a degree ≥ 15

图4 关键子网络分析（包括4 个模块）Figure 4 Key sub-module analysis of the PPI (including 4 modules)

2.4 关键节点及其功能富集分析

4 个关键子网络包含的所有基因及PPI网络中degree≥15的基因共有52个。将这52 个基因分别在数据集 GSE62452和GSE28735中进行生存分析，分别得到26 个和22 个与PAAD 预后相关的基因（P＜0.05），其中14个基因（ANLN、BGN、CENPF、DDX60、FN1、IFI44L、ITGA3、LAMA3、MET、MKI67、MMP14、OAS2、PLAU、TOP2A）在两个数据集中均与预后相关（图5A）。利用Metascape对这14个关键节点进行功能富集分析，发现这些基因在肿瘤侵犯和肿瘤发生中发挥一定作用（图5B）。

图5 关键节点的筛选及其功能富集分析 A：Venn 示GSE62452 和 GSE28735 数据集中与预后相关基因的交叉验证结果； B：Metascape 对14 个关键节点进行的功能和通路富集分析结果Figure 5 Screening of key nodes and function enrichment analysis A: The Venn diagram showing the result of cross-validation of prognostic related genes in the GSE62452 and GSE28735 datasets; B: Function and pathway enrichment analysis of 14 key nodeSBy using Metascape

2.5 最优特征基因及其表达验证

在最优参数（最小RMSE=0.3429，最大准确度=0.8694）下，用RFE算法筛选出8个最优特征基因：LAMA3、FN1、ITGA3、MET、PLAU、CENPF、MMP14、OAS2（图6）。GEPIA数据库验证发现这8 个最优特征基因在PAAD 组织中的表达均高于正常组织，差异有统计学意义（均P＜0.01）。

2.6 SVM 模型及其验证

通过R语言e1071包构建8个最优特征基因诊断PAAD的SVM建模如图7所示，该模型在训练集GSE28735中的曲线下面积（AUC）为0.898，灵敏度和特异度均为0.911；在验证集GSE62452中的AUC为0.905，敏感度为0.855，特异度为0.918；在验证集GSE62165中的AUC为1，敏感度和特异度均为1；表明该SVM模型可以较好地区分PAAD和正常样本。

图6 最优特征基因组合的准确度曲线（横轴代表基因变量的数量，纵轴代表交叉验证的准确性，点标记是最优特征基因组合对应的基因数）Figure 6 Accuracy curve for screening the optimal feature gene combination (the horizontal axis representing the number of gene variables, the vertical axis representing the cross-validation accuracy, and the marked content standing for the number of genes corresponding to the optimal gene combination)

图7 训练集（GSE28735）及验证集（GSE62165 和GSE62452）对SVM 模型的ROC 曲线Figure 7 ROC curves of the training set (GSE28735) and the validation sets (GSE62165 and GSE62452) on the SVM classifier

2.7 最优特征基因的生存分析

从TCGA数据库中筛选出178例PAAD样本，Kaplan-Meier生存曲线发现，与高表达LAMA3、ITGA3、MET、PLAU、CENPF及OAS2组的PAAD患者相比，这些基因的低表达组总生存率更长（均P＜0.05）（图8），而FN1、MMP14上调对PAAD的生存率无明显影响（均P＞0.05）。

图8 基于最佳截断值进行的TCGA 中最优特征基因与PAAD 关系的生存分析Figure 8 Survival analysis of the relationship between optimal feature genes and PAAD in TCGA database based on the optimal cutoff

3 讨 论

尽管PAAD的治疗已经取得了一些进展，但由于缺乏准确诊断PAAD的特异度生物标志物及个体化治疗的有效靶点，PAAD仍然是一种难治性癌 症[15]。探索PAAD的分子机制，找出潜在的诊断及预后的生物标志物，具有重要意义。在本研究中，从GSE28735、GSE62165及GSE62452数据集中筛选出257个在胰腺癌和正常组织之间的DEGs，其中168个上调DEGs和89个下调DEGs。通过GO和KEGG通路分析，发现这些DEGs主要参与的细胞外基质[16]、细胞黏附[17]、细胞迁移[18]、胰腺的分泌[19]、PI3K-Akt信号通路[20]已被证明与PAAD的发生和发展密切相关，表明本研究筛选的DEGs可靠性高。PPI预测和关键子网络分析发现一些新的在PAAD进展中起调控作用的基因，如COL5A2、OAS2、DDX60、CELA2A。既往研究发现这些基因与疾病的发生或治疗密切相关[21-24]，如DDX60的表达可调节乳腺癌患者对放疗的敏感度[23]，CELA2A突变与代谢综合征的发生相关[24]等；但他们在PAAD中尚无相关研究，这为以后研究PAAD的分子机制提供了参考。数据集中的生存分析筛选到14个与生存相关的关键节点，这些基因参与肿瘤的发生和侵犯，值得进一步研究。由递归特征消除算法筛选到8个最优特征基因，经GEPIA数据库证实这8个最优特征基因在PAAD组织中的表达高于正常组织，与芯片结果一致。鉴于单个靶点对肿瘤诊断价值有限，本研究通过R语言的e1071包对8个最优特征基因构建SVM模型，ROC曲线验证发现3个数据集的AUC均在0.85以上，敏感度和特异度高。据我们所知，用于PAAD诊断的SVM模型以前很少有报道，本研究构建的SVM模型可有效区分PAAD样本和正常样本，有望用于临床实践。

在TCGA 数据库中验证8 个最优基因的预后价值，其中CENPF、ITGA3、LAMA3、MET、OAS2、PLAU与PAAD患者的预后有关，即基因高表达者预后差、基因低表达者总生存率更长；该6个基因被视为关键基因。既往研究表明，ITGA3、LAMA3、CENPF在胰腺癌组织中高表达，且这些基因高表达的胰腺癌患者预后更差[25-27]，这与本研究结果一致。Jiao等[25]证实，ITGA3对胰腺癌的早期诊断有一定价值，其表达水平与组织学类型、生存状态以及复发相关。肿瘤间质相互作用是肿瘤治疗的重要靶点。研究[28]表明ITGA3是PAAD肿瘤间质相互作用的靶点之一，靶向ITGA3可能是PAAD治疗的潜在方法。Yang等[26]指出LAMA3不仅与PAAD 患者的预后有关，还有诊断PAAD的潜力。研究[29]发现LAMA3的沉默可抑制胰腺癌细胞的增殖、迁移和侵袭，并促进细胞凋亡，是PAAD潜在的治疗靶点。Chen等[27]发现CENPF的表达可促进胰腺癌细胞的增殖、迁移、上皮间质转化（EMT）及有丝分裂G2/M的转化，这种调节与TNF信号通路有关。PLAU可促进细胞外基质降解，参与细胞的侵袭和迁移。雷公藤甲素通过下调PLAU抑制胰腺癌细胞增殖和迁移，在此过程中，EMT信号通路被激活[30]。MET是受体酪氨酸激酶家族成员，包含MET的9个基因组合可有效预测PAAD患者的预后[31]。然而，这些关键基因在PAAD中的具体作用机制尚不清楚，仍需进一步研究。OAS2是2'，5'-寡腺苷酸合成酶中一员，在口腔鳞状细胞癌中高表达，与患者预后呈负相关[32]；而在乳腺癌中，OAS2高表达的患者预后较好[22]，这与本研究OAS2对PAAD的预后存异。可见，OAS2在不同肿瘤中的作用并不一致，目前缺乏OAS2在PAAD中的研究，需要进一步实验来研究OAS2在PAAD中的具体作用。

本研究通过生物信息学分析寻找可能参与PAAD发病的DEGs，通过对这些DEGs进行分析，以更好地了解PAAD致癌的机制。研究结果显示，COL5A2、OAS2、DDX60、CELA2A为探究PAAD的分子机制提供新路经；关键基因LAMA3、ITGA3、MET、PLAU、CENPF、OAS2可能成为PAAD诊断或治疗的新靶点；基于8个最优特征基因构建的SVM模型可有效诊断PAAD。本研究为今后研究PAAD的分子机制、生物标志物及治疗靶点提供了新的思路。然而，该研究的主要局限性是缺乏基础实验来证实这些结果的真实性，分析结果的可靠性严重依托于本报告中涉及数据集的准确性。未来的研究应通过基础实验和临床实践来验证本研究的结果。