CN109841280B - 食管癌相关特征通路的识别及早期预测模型的构建方法 - Google Patents

Abstract

本发明属于肿瘤诊断技术领域，具体涉及食管癌相关特征通路的识别及早期诊断模型的构建方法，具体包括表达谱预处理、差异表达基因提取、样本聚类分析、基因聚类分析、特异基因集功能通路分析、通路失常得分比较、功能差异比较分析、食管癌特异性共表达网络构建、对基因进行特征选择和深度学习模型预测等步骤。本发明方法将基因根据其表达相似性和功能一致性划分为不同的分组，以基因集合的形式进行分析，一方面可以避免传统方法中假阳性率高、随机误差大、结果不稳定等缺点，另一方面也可以更加特异性的识别出和食管癌显著相关的功能。

Description

食管癌相关特征通路的识别及早期预测模型的构建方法

技术领域

本发明属于肿瘤预测技术领域，具体涉及食管癌相关特征通路的识别及早期预测模型的构建方法。

背景技术

食管癌是一种非常普遍的消化道癌症，具有较高的发病率和致死率。然而，尽管在临床上食管癌的发病率很高，但食管癌的发病机制仍然不是很明确。临床上常规治疗手段包括手术治疗，放射线治疗，化疗，分子治疗等。然而由于缺乏早期诊断标识，对于进入进展期或晚期的食管癌患者，常规治疗方案往往得不到理想的治疗效果。因此食管癌的发病机制以及对应的有效治疗策略是我们当前迫切需要解决的问题。同时，识别出潜在的诊断标识和新的治疗靶点对于指导食管癌治疗具有重要的启发式意义。

传统分析手段基于单个基因表达水平筛选出和疾病具有统计学显著性的基因，其前提假设是基因特征之间相对独立。然而在生物体内，基因与基因之间并不是相对独立的，而是存在表达相关性和功能一致性，不同基因之间相互作用，共同参与调控下游的生物学作用。这种基因之间的相互作用可分为协同作用，拮抗作用，补偿作用。协同作用是指基因之间相互协调补充，彼此刺激，增强其调控作用；拮抗作用是指基因之间有竞争或抑制关系，通过拮抗作用可以使得基因的调控作用始终处于平衡状态，当其中一方发生异常，这种基因之间的失衡则会引起下游功能水平的异常；补偿作用是指基因之间属于平行关系，参与调控相同或相似的功能。通常情况只有其中一条通路发挥作用，只有当这条通路被阻断，另一条通路才被激活，从而起到功能代偿的作用。

发明内容

本发明主要提供一种食管癌相关特征通路的识别及早期预测模型的构建方法，其将基因根据其表达相似性和功能一致性划分为不同的分组，以基因集合的形式进行分析，一方面可以避免传统方法中假阳性率高、随机误差大、结果不稳定等缺点，另一方面也可以更加特异性的识别出和食管癌显著相关的功能。其技术方案如下：

一种食管癌相关特征通路的识别及早期预测模型的构建方法，包括以下步骤：

（1）表达谱预处理：从GEO数据库https://www.ncbi.nlm.nih.gov/geo/下载人类食管癌表达谱数据GSE100843，GSE100843共包含76个样本，其中40个be segment，36个正常对照组织normal squamous mucosa用control表示，将所有表达值经过z检验标准化校正，校正过程包括计算control组中的表达均值μ和标准差sd，然后利用公式1对所有样本的表达值进行校正，

，

其中，adj.exp为标准化校正后的表达值，exp为原始表达值，为基因表达均值，sd为标准差；

（2）差异表达基因提取：以正常normal squamous mucosa样本作为对照组，besegment作为实验组，通过limma Rpackage计算基因显著性，以显著性P值小于0.01且logFC超出5%和95%分位数的基因作为显著差异表达的基因；

（3）样本聚类分析：利用差异表达基因对所有食管癌样本和正常对照样本进行层次聚类，聚类过程通过R heatmap package实现，相似性矩阵采用pearson相关系数算法，最后通过热图的形式可视化；

（4）基因聚类分析：先利用主成分分析识别出最大的三个主成分，然后利用前三个主成分结合mclust算法（PMID: 27818791）对基因进行聚类，根据BIC准则评价最优的分类个数，以及每个基因聚类对应的基因集；

（5）特异基因集功能通路分析：利用每个子集内的基因进行功能富集分析，分析方法采用fisher exact test，每个基因子集可能富集到多个生物学功能，选择显著性最强的前两个功能来描述这个基因子集的主要功能，取显著性P值小于0.05的功能通路作为显著通路；

（6）通路失常得分比较：采用欧氏距离的量化方法计算通路失常得分，比较通路相对于正常样本在食管癌样本中的动态功能变化；

（7）功能差异比较分析：利用所有正常样本的失常得分计算均值μ和标准差sd，然后对于每个样本都做Z检验校正，若某个样本失常得分显著高，则信号在Z检验下被放大，反之则信号被削弱，对正常样本和食管癌两组失常得分做wilcox检验，根据显著性阈值p<0.05识别出食管癌相关的功能，这些筛选出的功能自身在食管癌发生过程中发生明显功能水平的偏差，或发生功能亢进，或被抑制，同时参与调控这些功能的基因也发生明显差异表达；

（8）食管癌特异性共表达网络构建：将差异表达基因随机两两组合并计算在所有样本中的Pearson相关系数，根据所有基因对的相关系数分布设定阈值，超过阈值的基因对认为存在共表达相关性，根据任意两基因之间的共表达相关性构建系统网络，通过对网络拓扑性质的分析，识别网络中hub基因，这些基因在网络中具有较高的度或介数；

（9）对基因进行特征选择：利用遗传算法模拟生物进化过程，通过亲本染色体重组过程，淘汰适应度低的子代，扩增适应度强的子代，从而优化出最适合的遗传信息组合，将所有特征基因随机组合形成特征链，初始化特征链长度为特征总数的50%，然后抽取一对特征链进行重组，使得子代特征链同时包含了来自亲本的特征信息，并计算子代适应度，用高适应度的子代特征链淘汰低适应度的子代特征链，此过程循环直至达到最大进化代数或模型收敛找到最优特征链，最后获得最优特征基因组合；

（10）深度学习模型预测：利用遗传算法对特征基因进行进化筛选，获得食管癌显著相关的基因组合，利用这些特征基因在两组样本中的表达值作为特征值，结合神经网络深度学习模型进行训练并预测，对分析数据随机排序，取50%作为训练集，余下50%作为测试集，训练过程利用gridsearch算法进行参数优化，优化参数包括激活函数，隐含层层数，学习速率，最后通过ROC曲线评价模型分类预测效能。

优选的，步骤（2）中数据共包含32321个基因探针，每个探针的表达值都经过标准化预处理，最后根据所有探针logFC的分布，取5%和95%的分位数作为阈值，显著性p值的阈值为0.01，共筛选出1616个上调基因探针和1616个下调基因探针。

优选的，步骤（6）中采用公式2进行通路失常得分计算，

，

其中，A（P）为功能的失常得分，m为通路中上调差异表达基因的个数，n为下调差异表达基因的个数，Xi为基因i的表达值，Xj是基因j的表达值，为基因在正常样本中的表达均值，若A(P)=0，则说明功能中上调基因和下调基因平衡；若A(P)>0，则说明上调基因占优势，功能发生上调偏差；若A(P)<0，则说明通路中下调基因占优势，功能发生下调偏差。

优选的，步骤（6）中利用公式1对所有富集到的通路进行打分，若多个基因子集富集到同一个功能，则合并取均值处理，最后共获得68条功能term，并利用功能得分构建68x76的矩阵，对于每一条功能均利用Wilcox检验验证其显著性，最后识别出9条生物学功能，均在两组样本间存在显著性。

优选的，步骤（8）中食管癌特异性网络构建的结果为，利用差异表达基因之间的表达值计算相似性矩阵，根据所有基因对之间的Pearson相关系数绝对值的95%分位数为阈值，95%的分位数对应相关系数为0.88，即筛选所有相关系数大于0.88的基因对作为存在显著关联的基因对，利用这些显著相关的基因对构建相似性共表达网络，网络中包含节点个数481，边个数1256，根据网络中每个节点所连结临近节点的个数，统计了每个节点的度分布，并进行2为底的对数转换，log2(degree)=7.3为阈值，所有大于阈值的节点作为hub节点。

优选的，步骤（9）中最优特征基因为HNF1B、 CREB3L1、IL35三个解释基因。

优选的，步骤（10）中深度学习模型预测的具体方法为，利用三个解释基因作为特征训练模型，首先利用gridsearch对模型参数进行优化，神经网络初始化参数为激活函数：sigmoid，学习速率0.01，输入层5个单元，隐藏层10个单元，输出层1个单元，经过gridsearch参数优化，优化后模型参数变更为激活函数relu，学习速率0.1，两个隐藏层，分别由10个，5个单元构成，样本随机重排后，在其中50%用于训练模型后，在剩下50%上进行测试，结果用ROC曲线表示，经过特征选择和参数优化后，训练集和测试集精度达到93%和87%，说明模型没有发生过拟合。

采用上述方案，本发明具有以下优点：

（1）本发明方法创新性在于克服传统分析方法依赖于单个基因的表达水平，从而引入大量假阳性结果，造成结论的可重复性差。我们通过基因表达相似性和功能一致性两个约束条件，限制了基因的随机波动，并有效识别出9个具有特异性的基因子集，每个基因子集均富集到显著特异的功能通路，提示了在食管癌样本中仍存在更精细的分子机制，结合共表达网络挖掘hub节点以及遗传算法进行特征选择，最后我们筛选出三个新的食管癌相关基因；

（2）本研究基于无监督算法识别出9个特异基因子集，通过量化分析识别出9个食管癌中发生显著失衡的功能，这些功能一方面解释了特异的食管癌发生机制，另一方面也提示了新的个性化治疗靶点。利用共表达网络和遗传算法筛选出3个食管癌相关基因。最后基于这三个基因为特征构建深度学习模型，可有效预测食管癌患者，对于临床上实现早期诊断具有重要意义。

附图说明

图1为差异表达基因分布图；

图2为差异表达基因的样本聚类热图；

图3为Mclust结合主成分分析聚类热图；

图4为9个基因子集之间的相关性矩阵图；

图5为9个基因子集的功能富集图；

图6为6个基因子集重要功能箱形图；

图7为3个基因子集重要功能箱形图；

图8为差异表达基因之间的相关分布图；

图9为用显著相关的基因对构建相似性共表达网络图；

图10为节点度分布图；

图11为基因对特征适应度图；

图12为样本在三个解释基因水平上的分布图；

图13为模型预测的ROC曲线图。

具体实施方式

以下实施例中的实验方法如无特殊规定，均为常规方法，所涉及的实验试剂及材料如无特殊规定均为常规生化试剂和材料。

1.表达谱预处理

从GEO数据库下载人类食管癌表达谱数据GSE100843（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE100843）。GSE100843共包含76个样本，其中40个be segment，36个正常对照组织normal squamous mucosa用control表示。由于基因之间的固有表达水平不同，有些基因对刺激敏感，表达波动大，而有些保守基因则相对波动较小，因此为了让所有基因之间具有可比性，同时便于作为深度学习模型的输入数据，我们将所有表达值经过z检验标准化校正，校正过程包括计算control组中的表达均值μ和标准差sd，然后利用公式1对所有样本的表达值进行校正。

，

其中：adj.exp为标准化校正后的表达值，exp为原始表达值，为基因表达均值，sd为标准差。

2.差异表达基因提取

以正常normal squamous mucosa样本作为对照组，be segment作为实验组，通过limma Rpackage计算基因显著性。为了严格控制基因显著性，避免因基因随机波动而引入假阳性结果，我们最后以显著性P值小于0.01且logFC超出5%和95%分位数的基因作为显著差异表达的基因。差异表达基因在背景基因中的分布利用火山图可视化。

数据共包含32321个基因探针，每个探针的表达值都经过标准化预处理，最后根据所有探针logFC的分布，取5%和95%的分位数作为阈值，显著性p值的阈值为0.01，最后我们共筛选出1616个上调基因探针和1616个下调基因探针，差异表达基因分布情况见图1所示。图1中右边深色圆圈为上调探针，左边深色三角为下调探针，中间灰色节点为背景非差异表达基因。横线对应p值为0.01的cutoff，竖线分别对应logfc为-1.4和1.24的cutoff。我们发现基因的显著性p值和logfc之间存在一定的关联。拥有较大logfc的基因，也倾向于伴随更显著的p值。同时整体基因中，上调基因和下调基因的数量比较平衡，没有发生明显的偏差。对于差异表达的探针，我们根据平台数据对基因名进行转换，探针id均被转换为genesymbol。对于多个探针对应同一个gene symbol的情况，我们合并后取均值处理，探针合并后共获得2598个差异表达基因，这些基因在食管癌样本和正常对照样本中存在明显差异表达，logFC和显著性p值均超过阈值标准。

3.对样本聚类分析

为了进一步观察食管癌相关基因在疾病状态和正常状态下的差异，我们利用差异表达基因对所有食管癌样本和正常对照样本进行层次聚类。聚类过程通过R heatmappackage实现。相似性矩阵采用pearson相关系数算法，最后通过热图的形式可视化。通过层次聚类一方面可以验证差异表达基因对于区分两组样本的效果，另一方面可以观察样本内部是否存在异常样本或特异的局部特征。聚类结果用热图可视化，如图2所示。

图2中每一行为一个基因，每一列为一个样本，横坐标左半部分标记normal对照组，横坐标右半部分标记be_segment食管癌组。图中浅色色块代表低表达，深色色块代表高表达。从热图中可以明显观察出：1. 我们识别出的差异表达基因可以有效区分食管癌样本和正常样本；2. 差异表达基因之间存在差异的表达模式，例如有些基因在食管癌中高表达，而在正常组中低表达，反之有些基因在食管癌中低表达，而在正常组中高表达；3. 正常样本所在聚类明显均一性更高，即上调基因和下调基因存在非常稳定的分布，而在食管癌组中则体现出明显的混乱度，原本在正常组中的上调基因和下调基因，在食管癌样本中发生了转换，进一步提示了在食管癌发生过程中发生了某些基因的表达模式转换以及功能失调。

4.对基因聚类分析

由于基因之间也存在功能相关性，因此我们同时也对基因进行聚类，聚类在一起的基因被认为具有一定的功能相似性。为了识别合理的基因聚类，我们先利用主成分分析识别出最大的三个主成分，然后利用前三个主成分结合mclust算法对基因进行聚类，根据BIC准则评价最优的分类个数，以及每个基因聚类对应的基因集。

功能异常以及识别重要的疾病相关基因具有重要意义。我们首先对基因进行pca主成分分析，利用前三个主成分结合mclust算法对基因进行聚类分析，聚类结果如图3所示。图3中pc1-pc3分别对应前三个主成分，根据BIC准则最后基因被分成9个基因子集，用9种深度不同的颜色表示，可以看出在前三个主成分的水平上，9个基因子集在分布上存在明显区别。根据给定9个基因子集，计算基因之间的相关性矩阵，如图4所示。图4中用9种深度不同的颜色标记9个基因子集，色块代表子集内部基因相关性强弱，颜色越深说明内部相关性越强，因此可以观察出任意一个基因子集内部都体现出极强的相关性，证明了这些基因子集内部存在表达相关性，进一步提示了每个基因子集的功能一致性。

5. 特异基因集功能通路分析

基于基因在不同疾病状态下的表达水平，我们对基因进行聚类，每个聚类中的基因集具有相似的表达模式，提示了在功能上这些基因存在紧密的相关性。同时某些基因集在食管癌样本以及正常对照组中存在显著差异，因此我们利用这些特异基因集进行功能富集分析，富集方法采用DAVID软件实现，分析方法采用fisher exact test（费舍尔精确检验），取显著性P值小于0.05的功能通路作为显著通路。每个基因子集可能富集到多个生物学功能，我们选择显著性最强的前两个功能来描述这个基因子集的主要功能，如图5所示。图5中行为负对数转换的p值，纵轴为生物学功能term，term1为Phenylalaninemetabolism，term2为Sphingolipid metabolism，term3为TGF-beta signaling pathway，term4为Tightjunction，term5为regulation of cilium assembly，term6为protein phosphorylation，term7为regulation of protein kinase C signaling，term8为mitral valvemorphogenesis，term9为regulation small GTPase mediated signal transduction，term10为Pathogenic Escherichia coil infection，term11为microtubule-basedprocess，term12为Hippo signaling pathway，term13为Glutathione metabolism，term14为Fructose and mannose metabolism，term15为Starch and sucrose metabolism，term16为Gastric acid secretion，term17为Pathways in cancer，我们用9个深度不同的颜色标记9个基因子集。可以观察出每个功能子集都有其特异性的生物学作用，例如geneset1集中在氨基酸和鞘脂代谢；geneset2集中在信号通路和细胞连接；geneset3集中在蛋白磷酸化；geneset4集中在蛋白激酶c信号通路；geneset5集中在GTP酶街道的信号转换；geneset6集中在大肠杆菌感染和微管发生过程；geneset7集中在谷胱甘肽代谢；geneset8集中在糖类代谢；geneset9集中在胃酸分泌和癌症通路。

6.通路失常得分比较

利用不同聚类中的基因进行功能富集，这些功能的显著性一方面体现在差异基因高度富集在相应功能term中，另一方面体现在功能量化水平上的显著改变。我们采用欧氏距离的量化方法按照公式2计算通路相对于正常样本，在食管癌样本中的动态功能变化。

，

对于功能term P，A（P）为功能的失常得分，m为通路中上调差异表达基因的个数，n为下调差异表达基因的个数，Xi为基因i的表达值，Xj是基因j的表达值，为基因在正常样本中的表达均值，最后取2为底的对数转换。因此若A(P)=0，则说明功能中上调基因和下调基因平衡；若A(P)>0，则说明上调基因占优势，功能发生上调偏差；若A(P)<0，则说明通路中下调基因占优势，功能发生下调偏差。通过公式2计算通路P偏离正常状态的程度。

通过富集分析，我们识别出每个基因子集所特异的生物学功能，然而为了量化这些功能从而评价哪些功能是和食管癌发生显著相关的，我们利用公式1对所有富集到的通路进行打分。若多个基因子集富集到同一个功能，则合并取均值处理，最后共获得68条功能term，并利用功能得分构建68x76的矩阵。对于每一条功能我们均利用Wilcox检验验证其显著性，最后识别出9条生物学功能，均在两组样本间存在显著性，如图6、图7所示。图6、7中横轴为两组样本，左边的为食管癌，右边的为normal对照，纵轴为每个功能的打分。可以看出9个功能在两组样本中存在显著的差异分布。同时也发现除了detection of chemicalstimulus involved in sensory perception和Pathogenic Escherichia coliinfection两条功能是食管癌相对于正常样本体现出功能下调，其他功能都在食管癌样本中体现出明显的上调趋势。

7.功能差异分析比较

我们利用公式2计算出功能通路在每个样本中的失常得分，为了识别出食管癌相关的功能我们做了如下分析。首先利用所有正常样本的失常得分计算均值μ和标准差sd，然后对于每个样本都做Z检验校正，若某个样本失常得分显著高，则信号在Z检验下被放大，反之则信号被削弱。我们对正常样本和食管癌两组失常得分做wilcox检验，根据显著性阈值p<0.05识别出食管癌相关的功能。这些筛选出的功能自身在食管癌发生过程中发生明显功能水平的偏差，或发生功能亢进，或被抑制，同时参与调控这些功能的基因也发生明显差异表达，因此这些功能对于解释食管癌发生机制具有重要意义。

8.食管癌特异性共表达网络构建

为了从系统的角度分析食管癌相关的差异表达基因，我们将差异表达基因随机两两组合并计算在所有样本中的Pearson相关系数，根据所有基因对的相关系数分布设定阈值，超过阈值的基因对认为存在共表达相关性。通过对网络拓扑性质的分析，识别网络中hub基因，这些基因在网络中具有较高的度或介数，即一个hub节点与周围多个临近节点有边。当hub基因发生表达异常时，往往可以影响周围多个基因同时发生差异表达，进而引起下游功能异常。因此这些网络中具有较高度分布的hub基因很可能是重要的疾病靶点或诊断标记物。为了筛选出具有功能调控作用的特征基因，我们利用功能比较分析中wilcox显著的功能的基因取并集作为特征基因。

我们利用差异表达基因之间的表达值计算相似性矩阵，根据所有基因对之间的Pearson相关系数绝对值的95%分位数为阈值，如图8所示。横轴为相关系数，纵轴为基因对的频数，95%的分位数对应相关系数为0.88，即我们筛选所有相关系数大于0.88的基因对作为存在显著关联的基因对。我们利用这些显著相关的基因对构建相似性共表达网络，网络中包含节点个数481，边个数1256，如图9所示。图9中六边形节点为上调基因，菱形节点为下调基因，根据网络中每个节点所连结临近节点的个数，我们统计了每个节点的度分布，并进行2为底的对数转换，节点度分布如图10所示。横轴为2为底对数转换的节点度，纵轴为density分布，可以看出网络中的节点度存在明显的两个峰，peak1对应的是节点度较少的点（log2(degree)<1.2），而在log2(degree)从5开始则为平台期并逐渐进入度增加的阶段，另一个峰是对应的节点度最高的点（log(degree)>7.3）。因此我们卡log2(degree)=7.3为阈值，所有大于阈值的节点作为hub节点。

9.利用遗传算法对基因进行特征选择

遗传算法模拟生物进化过程，通过亲本染色体重组过程，淘汰适应度低的子代，扩增适应度强的子代，从而优化出最适合的遗传信息组合。我们将所有特征基因随机组合形成“特征链”，初始化特征链长度为特征总数的50%。然后抽取一对特征链进行重组，使得子代特征链同时包含了来自亲本的特征信息，并计算子代适应度（本研究中适应度为特征链的预测精度），用高适应度的子代特征链淘汰低适应度的子代特征链，此过程循环直至达到最大进化代数或模型收敛找到最优特征链。我们设置最大进化代数为100代，容忍度为0.001。之后逐步缩小初始化特征链长度，重复以上过程，最后获得最优特征基因组合。当相邻两代的适应度相差不足0.001时模型提前终止，否则说明模型未收敛，则一直进化直至达到100代。

我们利用富集到显著差异的功能中的基因为特征，在100代进化过程中模型精度逐渐收敛，最后获得近似最优解，即最优的特征基因组合，如图11所示。图11中展示了在100次进化后适应度超过0.5的基因特征，纵轴为对应的适应度，筛选的8个基因中有四个基因具有明显较强的适应度，分别是APPC2、CREB3L1、HNF1B、IL35。同时我们还比较了这些基因的节点度信息，如表1所示。

表1 主要基因信息

注：第一列为基因名，第二列为适应度，第三列为节点在网络中经对数转化的的度。

可以看出fitness最高的四个基因有三个基因同时也是网络中的hub节点，即CREB3L1、HNF1B、IL35。提示了这三个基因很可能是食管癌相关的重要的解释基因。所有样本在这三个基因水平上的分布如图12所示。图12中x,y,z三个轴分别对应三个解释基因，深色圆形样本为normal组，浅色三角形样本为食管癌组，可以看出在这三个解释基因的作用下，两组样本可以完全被区分开。同时也可以观察出，正常样本在这三个基因的水平上比较收敛，彼此距离相对紧凑。然而在食管癌样本中则体现的更加离散。并且可以直观的发现食管癌样本并不总是聚集在一起，有些食管癌样本彼此接近，可以聚类成簇。这也再一次验证了我们的推测，食管癌样本可以进一步分成不同的分子亚型。从整体趋势上看，相对于正常样本，在食管癌样本中，HNF1B, CREB3L1, IL35三个基因具有不同水平的上调。

10.深度学习模型预测

我们利用遗传算法对特征基因进行进化筛选，获得食管癌显著相关的基因组合，利用这些特征基因在两组样本中的表达值作为特征值，结合神经网络深度学习模型进行训练并预测。我们对分析数据随机排序，取50%作为训练集，余下50%作为测试集。训练过程利用gridsearch算法进行参数优化，优化参数包括激活函数，隐含层层数，学习速率。最后通过ROC曲线评价模型分类预测效能。

我们利用三个解释基因作为特征训练模型，首先利用gridsearch对模型参数进行优化，神经网络初始化参数为激活函数：sigmoid，学习速率0.01，输入层5个单元，隐藏层10个单元，输出层1个单元。经过gridsearch参数优化，优化后模型参数变更为激活函数relu，学习速率0.1，两个隐藏层，分别由10个，5个单元构成。样本随机重排后，在其中50%用于训练模型后，在剩下50%上进行测试，结果用ROC曲线表示，如图13所示。图13中横轴为特异性，纵轴为敏感性，三角线条为训练集精度，矩形线条为测试集精度，经过特征选择和参数优化后，训练集和测试集精度达到93%和87%，一方面模型精度达到较高水平，另一方面训练集和测试集结果比较接近，说明模型没有发生过拟合。因此我们识别出的3个解释基因能够准确的识别出食管癌患者和正常对照，利用神经网络诊断模型，基于这三个基因的表达水平可以有效实现早期预测食管癌患者，同时也对研究食管癌致病机制提供了新的思路。

综上分析，本方法严格控制筛选差异表达基因的阈值，根据limma显著性p值小于0.01，且logfc超出5%和95%分位数筛选出显著差异表达的基因。在食管癌样本和正常对照组中，这些基因均体现出显著差异的表达模式，有些在正常对照组中表达上调的基因，在食管癌样本中变为下调，而另一部分在正常组织中下调的基因，则在食管癌中体现为上调。通过无监督的聚类算法我们发现基于这种差异的表达模式，两组样本可以被完全区分开。我们发现在正常组中，差异表达基因根据他们的表达模式可以分成两部分，即上调基因和下调基因。同时上调基因或下调基因内部相对均一，不存在较大的波动。而在食管癌样本中，不仅差异表达基因的表达方向发生了转换，同时也体现出明显的不均一性和波动性。可以观察出某些特定基因在有些食管癌样本中转换的幅度大，而在另一些食管癌样本中转换的幅度小。这种不均一的波动一方面来自食管癌本身的异质性：在食管癌进展的不同阶段，肿瘤细胞所处的分期，亚克隆，细胞核型，以及周围内环境均有所差异，因此导致同一个基因在不同的肿瘤组织中体现出差异的波动。另一方面来自于癌症样本的不均一性：食管癌患者根据癌症进展程度，分期，以及个体化遗传差异可以被进一步分层。同一个基因在处于不同水平的个体中，表达状态也有所差异。因此，在食管癌样本中的不均一波动提示我们，基于分子水平差异，食管癌样本也许可以被进一步分成不同的分子亚型。由于基因集在不同亚型中的表达模式不同，因此每个分子亚型都可能存在特异的发病机制，同时也对个性化治疗提供了新的思路。

通过层次聚类分析我们发现，在所有差异表达基因中仍存在具有特异性的基因子集，这些基因子集在所有食管癌样本中体现出差异的表达模式，但却在特定的某些样本中体现出一致性。为了识别出这些具有特异性关联的基因子集，我们利用主成分分析提取了前三个主成分用于对基因进行进一步分类。根据BIC准则，我们最后识别出9个基因子集。每个基因子集都进行独立的功能富集分析。由于基因子集在样本中体现出明显的差异性和特异性，因此这些功能也体现了个性化特征。geneset1集中在氨基酸和鞘脂代谢；geneset2集中在信号通路和细胞连接；geneset3集中在蛋白磷酸化；geneset4集中在蛋白激酶c信号通路；geneset5集中在GTP酶街道的信号转换；geneset6集中在大肠杆菌感染和微管发生过程；geneset7集中在谷胱甘肽代谢；geneset8集中在糖类代谢；geneset9集中在胃酸分泌和癌症通路。我们怀疑在食管癌样本中仍存在更精细的分子机制。尽管这些功能的异常最终都会引起食管癌发生，但根据其发生机理的不同，可以有效制定个性化治疗方案。

我们利用欧式距离计算每个功能中上调基因和下调基因的偏差，从而评价功能的失衡程度。若一个功能中大部分基因都发生上调，则认为在食管癌发生过程中，这个功能的水平收到了激活或亢进，反之则认为功能受到了抑制或阻断。通过对功能水平进行量化分析，结合Wilcox差异分析，我们可以有效识别出在食管癌发生或进展过程中，功能水平明显发生失衡的显著功能。这些发生明显失衡的功能包括Retinol metabolism, Tightjunction, protein phosphorylation, regulation of Rho protein signaltransduction, detection of chemical stimulus involved in sensory perception,Pathogenic Escherichia coli infection, Maturity onset diabetes of the young,Amino sugar and nucleotide sugar metabolism, Thyroid hormone synthesis。其中detection of chemical stimulus involved in sensory perception和PathogenicEscherichia coli infection功能水平下调，其他7个功能水平上调。这是一个非常有趣的结论，我们可以联想当化学刺激信号的感知被阻断，外源性刺激和内源性胃酸腐蚀作用，将会持续的损伤食道。慢性的炎症刺激是食管癌发生的重要诱因。抑制大肠杆菌感染通路可能会引起线粒体功能异常和细胞骨架失控，从而诱导细胞凋亡。

一个功能通常受到多个基因的调控作用，这些基因的贡献并不相同。有些基因可通过自身表达变化同时影响多个下游基因，具有重要的“枢纽”作用。因此往往只通过少数几个基因就可以代表整个功能的变化水平，同时这些具有重要贡献的基因也往往更可能成为多效靶点。共表达网络是一个非常有效的手段来研究基因之间的相互关联。若基因之间存在相互调控作用，则他们的表达水平也呈现一定的相关性。若两个基因之间是激活刺激关系，则体现为正相关，若两个基因之间是抑制反馈关系，则体现为负相关。通过对共表达网络的节点度分析，我们识别出显著的hub gene。

最后我们利用遗传算法对所有参与调控失衡通路的基因进行特征选择。通过最多100代进化最后筛选出显著特征组合作为近似最优解。筛选出的基因在食管癌样本中均发生显著差异表达，同时参与调控发生明显失衡的功能通路。我们还从特征基因中筛选出三个hub gene，CREB3L1,HNF1B,IL35。这三个基因均在共表达网络中具有较高的贡献度，其自身差异表达可同时影响多个下游基因发生异常。由基因creb3l1编码的蛋白主要被发现于内质网膜。然而由于内质网外压力作用下，编码蛋白被清除，同事释放的细胞质转录因子转移到细胞核。在核内通过结合box-b元件激活了靶基因的转录。Creb3l1已经被证实是许多癌症的重要诱因，包括直肠癌，前列腺癌，乳腺癌等。Hnf1b编码了转录因子超家族成员蛋白。这个蛋白结合dna，和相关蛋白1-alpha形成二聚体。这个基因的表达水平在某些癌症中发生明显改变。已经有很多报道证明这个基因的转录本存在突变。Il-35编码了一种分泌性糖蛋白，并且和23 kda蛋白形成异质二聚体，构成IL27。Il27通过参与激活jak/stat通路调控t细胞和炎症应答。最后以这三个基因为特征我们构建了神经网络分类起作为诊断预测模型，训练集和测试集的精确度分别为93％，和87％。 本研究的创新性在于克服传统分析方法依赖于单个基因的表达水平，从而引入大量假阳性结果，造成结论的可重复性差。我们通过基因表达相似性和功能一致性两个约束条件，限制了基因的随机波动，并有效识别出9个具有特异性的基因子集。每个基因子集均富集到显著特异的功能通路，提示了在食管癌样本中仍存在更精细的分子机制。结合共表达网络挖掘hub节点以及遗传算法进行特征选择，最后我们筛选出三个新的食管癌相关基因。

最后本研究基于无监督算法识别出9个特异基因子集。通过量化分析识别出9个食管癌中发生显著失衡的功能。这些功能一方面解释了特异的食管癌发生机制，另一方面也提示了新的个性化治疗靶点。利用共表达网络和遗传算法筛选出3个食管癌相关基因。最后基于这三个基因为特征构建深度学习模型，可有效预测食管癌患者。对于临床上实现早期诊断具有重要意义。

对本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及形变，而所有的这些改变以及形变都应该属于本发明权利要求的保护范围之内。

Claims (3)

1.一种食管癌相关特征通路的识别及早期预测模型的构建方法，其特征在于：包括以下步骤：

(1)表达谱预处理：从GEO数据库https://www.ncbi.nlm.nih.gov/geo/下载人类食管癌表达谱数据GSE100843，GSE100843共包含76个样本，其中40个be segment，36个正常normal squamous mucosa样本用control表示，将所有表达值经过z检验标准化校正，校正过程包括计算control中的均值μ和标准差sd，然后利用公式1对所有样本的表达值进行校正，

，

其中，adj.exp为标准化校正后的表达值，exp为原始表达值，μ为基因表达均值，sd为标准差；

(2)差异表达基因提取：以36个正常normal squamous mucosa样本作为对照组，40个besegment作为实验组，通过limma Rpackage计算基因显著性，以显著性P值小于0.01且logFC超出5％ 和95％分位数的基因作为显著差异表达的基因；

(3)样本聚类分析：利用差异表达基因对所有食管癌样本和正常对照样本进行层次聚类，聚类过程通过R heatmap package实现，相似性矩阵采用pearson相关系数算法，聚类结果用热图可视化；

(4)基因聚类分析：先利用主成分分析识别出最大的三个主成分，然后利用前三个主成分结合mclust算法对基因进行聚类，根据BIC准则评价最优的分类个数，以及每个基因聚类对应的基因子集；

(5)特异基因集功能通路分析：利用每个基因子集内的基因进行功能富集分析，分析方法采用fisher exact test，每个基因子集富集到多个生物学功能，选择显著性最强的前两个功能来描述这个基因子集的主要功能，取显著性P值小于0.05的功能通路作为显著通路；

(6)通路失常得分比较：采用欧氏距离的量化方法计算通路失常得分，比较通路相对于正常normal squamous mucosa样本在食管癌样本中的动态功能变化；

(7)功能差异比较分析：利用所有正常normal squamous mucosa样本的失常得分计算均值μ和标准差sd，然后对于每个正常normal squamous mucosa样本都做Z检验校正，若某个样本失常得分显著高，则信号在Z检验下被放大，反之则信号被削弱，对正常normalsquamous mucosa样本和食管癌样本两组失常得分做wilcox检验，根据显著性阈值p<0.05识别出食管癌相关的功能；

(8)食管癌特异性共表达网络构建：将差异表达基因随机两两组合并计算在所有样本中的Pearson相关系数，根据所有组合相关系数分布设定阈值，超过阈值的组合认为存在共表达相关性，根据任意两基因之间的共表达相关性构建系统网络，通过对网络拓扑性质的分析，识别网络中hub节点；

(9)对基因进行特征选择：利用遗传算法模拟生物进化过程，通过亲本染色体重组过程，优化出最适合的遗传信息组合，将所有特征基因随机组合形成特征链，初始化特征链长度为特征总数的50％，然后抽取一对特征链进行重组，使得子代特征链同时包含了来自亲本的特征信息，并计算子代适应度，适应度为特征链

的预测精度，此过程循环直至达到最大进化代数或模型收敛找到最优特征链，最后获得最优特征基因组合；

(10)深度学习模型预测：利用遗传算法对特征基因进行进化筛选，获得食管癌显著相关的基因组合，利用这些特征基因在两组样本中的表达值作为特征值，结合神经网络深度学习模型进行训练并预测，对分析数据随机排序，取50％作为训练集，余下50％作为测试集，训练过程利用gridsearch算法进行参数优化，优化参数包括激活函数，隐含层层数，学习速率，最后通过ROC曲线评价模型分类预测效能；

步骤(6)中利用公式1对所有通路进行打分，若多个基因子集富集到同一个功能，则合并取均值处理，最后共获得68条功能term，并利用功能得分构建68x76的矩阵，对于每一条功能均利用Wilcox检验验证其显著性，最后识别出9条生物学功能，均在两组样本间存在显著性；

步骤(8)中食管癌特异性网络构建的结果为，利用差异表达基因之间的表达值计

算相似性矩阵，根据所有组合之间的Pearson相关系数绝对值的95％分位数为阈值，95％的分位数对应相关系数为0.88，即筛选所有相关系数大于0.88的组合作为存在显著关联的组合，利用这些显著相关的组合构建相似性共表达网络，网络中包含节点个数481，边个数1256，根据网络中每个节点所连结临近节点的个数，统计了每个节点的度分布，并进行2为底的对数转换，log2(degree)＝7.3为阈值，将所有大于阈值的节点作为hub节点；

步骤(9)中最优特征基因为HNF1B、CREB3L1、IL35三个解释基因；

步骤(10)中深度学习模型预测的具体方法为，利用所述三个解释基因作为特征训练模型，首先利用gridsearch算法对模型参数进行优化，神经网络初始化参数为激活函数：sigmoid，学习速率0.01，输入层5个单元，隐藏层10个单元，输出层1个单元，经过gridsearch算法对模型参数进行优化，优化后模型参数变更为激活函数relu，学习速率0.1，两个隐藏层，分别由10个，5个单元构成，样本随机重排后，在其中50％用于训练模型后，在剩下50％上进行测试，结果用ROC曲线评价模型分类预测效能。

2.根据权利要求1所述的食管癌相关特征通路的识别及早期预测模型的构建方法，其特征在于：步骤(2)中数据共包含32321个基因探针，每个基因探针的表达值都经过标准化预处理，根据所有基因探针logFC的分布，取5％和95％的分位数作为阈值，显著性p值的阈值为0.01，共筛选出1616个上调基因探针和1616个下调基因探针。

3.根据权利要求1所述的食管癌相关特征通路的识别及早期预测模型的构建方法，其特征在于：步骤(6)中采用公式2进行通路失常得分计算，

，

其中，A(P)为功能的失常得分，m为通路中上调差异表达基因的个数，n为下调差异表达基因的个数，X_i为基因i的表达值，X_j是基因j的表达值，μ为基因在正常样本中的表达均值，若A(P)＝0，则说明功能中上调基因和下调基因平衡；若A(P)>0，则说明上调基因占优势，功能发生上调偏差；若A(P)<0，则说明通路中下调基因占优势，功能发生下调偏差。