一种胆道闭锁潜在分子亚型及其核心基因的识别方法

技术领域

本申请涉及胆道闭锁基因识别技术领域，特别是涉及一种胆道闭锁潜在分子亚型及其核心基因的识别方法。

背景技术

胆道闭锁(biliary atresia，BA)是新生儿常见的一种进行性的阻塞性胆道疾病，临床表现为梗阻性黄疸，如果BA患儿得不到及时治疗，易发展为胆汁淤积性肝硬化、门静脉高压，最终导致肝衰竭而死亡。目前，研究者提出了许多理论来解释BA的病因，包括基因缺陷、病毒感染、免疫形态异常、基因的多态性等。BA的复杂致病因素，暗示了BA中可能存在不同的分子亚型。

目前，基于高通量数据对疾病进行分子亚型的研究主要是根据样本基因表达水平的相似性。通过无监督聚类的方法，例如，非负矩阵分解(NMF)、层次聚类和潜在因子分析(LF)等挖掘潜在的分子亚型，但这类方法往往基于较为复杂的数量统计模型，同时由于基因芯片的检测受实验条件、试剂批号、实验操作人员及操作流程的影响较大，导致由不同实验室检测的样本及由同一个实验室不同批次检测的样本存在系统性偏差，当从一组样本中识别了预后标志并估计了预测预后的阈值后，该阈值并不能直接应用于由其它实验室采用同样的基因芯片平台检测的样本中，需要被重新训练，因此难以利用不同实验室检测的数据直接验证一组标志，更难以应用于单个样本的个体化判断。

发明内容

基于此，有必要针对上述技术问题，提供一种新的胆道闭锁潜在分子亚型及其核心基因的识别方法。

一种胆道闭锁潜在分子亚型及其核心基因的识别方法，所述方法包括：

步骤A：从基因表达综合数据库获取胆道闭锁组织原始表达谱数据，对所述表达谱数据进行预处理得到包含正常组织样本、疾病组织样本的背景基因；

步骤B：将背景基因两两组合，比较所述正常组织样本、所述疾病组织样本中基因对的相对表达秩次关系差异，确定具有潜在疾病分子亚型鉴别能力的基因对，根据所述具有潜在疾病分子亚型鉴别能力的基因对识别疾病分子亚型；

步骤C：基于预先设计的加权基因共表达网络提取所述疾病分子亚型相关的基因模块作为候选模块；

步骤D：对所述候选模块内的基因进行蛋白与蛋白之间的互作分析，筛选得到核心基因。

在其中一个实施例中，所述步骤A：从基因表达综合数据库获取胆道闭锁组织原始表达谱数据，对所述表达谱数据进行预处理得到包含正常组织样本、疾病组织样本的背景基因，具体包括：

从GEO下载原始胆道闭锁组织表达谱数据，利用R软件Affy软件包提供的RMA算法对其进行背景校正，并利用平台注释文件将探针ID对应到基因ID上。

在其中一个实施例中，所述步骤B：将背景基因两两组合，比较所述正常组织样本、所述疾病组织样本中基因对的相对表达秩次关系差异，确定具有潜在疾病分子亚型鉴别能力的基因对，根据所述具有潜在疾病分子亚型鉴别能力的基因对识别疾病分子亚型，具体包括：

步骤B1：基于基因对中两基因的基因表达水平大小关系，筛选所述正常组织样本中基因间表达稳定的基因对作为稳定基因对；

步骤B2：基于所述稳定基因对在疾病组织样本中相对表达秩次关系的改变情况，筛选所述疾病组织样本中的逆转基因对；

步骤B3：基于所述逆转基因对在疾病样本中的相对表达秩次关系，筛选得到具有潜在疾病分子亚型鉴别能力的基因对；所述具有潜在疾病分子亚型鉴别能力的基因对满足逆转对的秩次关系的改变只在部分疾病患者之中发生，同时逆转对将疾病样本分成两组，且两组间识别的差异基因存在生物学功能的扰动；

步骤B4：根据所述潜在疾病分子亚型鉴别能力的基因对识别疾病分子亚型。

在其中一个实施例中，所述步骤B1：基于基因对中两基因的基因表达水平大小关系，筛选所述正常组织样本中基因间表达稳定的基因对作为稳定基因对，具体包括：

将正常组织样本中的背景基因两两组合，根据公式(1)，筛选满足P(G

其中n表示所有的正常样本数，t代表第t个样本，G

在其中一个实施例中，所述步骤B2：基于所述稳定基因对在疾病组织样本中相对表达秩次关系的改变情况，筛选所述疾病组织样本中的逆转基因对，具体包括：

基于正常组织样本中得到的稳定基因对G

在其中一个实施例中，所述步骤B3：基于所述逆转基因对在疾病样本中的相对表达秩次关系，筛选得到具有潜在疾病分子亚型鉴别能力的基因对，具体包括：

获取在疾病样本中筛选得到的任意逆转基因对G

根据所述逆转基因对G

利用t检验，控制假阳性发现率FDR为5％，识别两组间的差异基因，对两组间识别的差异基因进行功能富集分析；

如果一个逆转基因对按其在疾病样本中的相对表达秩次关系分成的两组样本间存在差异表达基因，同时存在生物学功能的扰动，则该基因对为具有潜在疾病分子亚型鉴别能力的基因对。

在其中一个实施例中，所述步骤B4：根据所述潜在疾病分子亚型鉴别能力的基因对识别疾病分子亚型，具体包括：

根据所述具有潜在疾病分子亚型鉴别能力的基因对扰动通路的数目；

将所述基因对进行降序排列，构建一个具有潜在疾病分子亚型鉴别能力的基因对×疾病样本的矩阵，所述基因对×疾病样本的矩阵是由基因对数作为行、疾病样本数作为列的矩阵；

利用欧式距离对所述基因对×疾病样本矩阵进行聚类分析，识别疾病分子亚型。

在其中一个实施例中，所述步骤C：基于预先设计的加权基因共表达网络提取所述疾病分子亚型相关的基因模块作为候选模块，具体包括：

基于预先设计的加权基因共表达网络WGCNA分析提取与疾病亚型显著相关的基因模块作为候选模块；

对所述候选模块中背景基因的基因表达矩阵相关系数进行加权，使基因间的相互作用关系符合无标度分布；

对所述候选模块的参数进行设置，将具有相似表达模式的基因划分为一个模块，最小模块大小设置为100，其他参数设置为默认值。

在其中一个实施例中，在执行所述步骤D之前，还可以对所述候选模块内的基因功能进行富集分析，具体为：

基于京都基因和基因组百科全书KEGG数据库对模块基因的功能注释和通路富集析，用R包clusterProfile实现。

在其中一个实施例中，所述步骤D：对所述候选模块内的基因进行蛋白与蛋白之间的互作分析，筛选得到核心基因，具体包括：

使用STRING在线网站对所述候选模块内的基因进行蛋白与蛋白之间的互作分析，由Cytoscape显示；其中，Cytoscape中的cytoHubba插件用于实现网络最大集团度评分算法，如公式(2)；

通过所述cytoHubba插件识别PPI网络中的关键基因，

MCC(ν)＝∑

其中，集合S表示集合中元素的数量，S(ν)是包含ν的最大集团的集合，(|C|-1)！是所有小于|C|的正整数的乘积。

上述一种胆道闭锁潜在分子亚型及其核心基因的识别方法，基于样本内基因间相对表达秩次关系识别潜在疾病分子的亚型，摆脱了基因表达绝对值的限制，相当程度地克服了不同实验室间表达值不可比的弱点。大部分片间标准化方法对表达值的处理是单调变化的，因而不影响基因表达秩次关系的变化，而片间标准化方法对基因值的处理参照群体表达水平，因此，基于基因表达秩次关系的方法还能够最大程度地消除群体差异对结果造成的影响。除此之外，这类方法相对于基于绝对检测值的方法对个体间的生物学变异更加不敏感，更易得到有生物学意义的结果，通过这种方法识别的胆道闭锁亚型关键基因可以应用于单个个体样本。

附图说明

图1为一个实施例中一种胆道闭锁潜在分子亚型及其核心基因的识别方法的流程示意图；

图2为一个实施例中GSE46960和GSE15235聚类分析结果示意图；

图3为一个实施例中无标度拓扑拟合指数与平均连通性的示意图；

图4为一个实施例中模块与表型相关分析示意图；

图5为一个实施例中KEGG富集分析图；

图6为一个实施例中PPI网络示意图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本申请提供的一种胆道闭锁潜在分子亚型及其核心基因的识别方法，基于样本内基因间相对表达秩次关系，分析胆道闭锁基因表达谱数据集(GSE46960，GSE15235)，挖掘潜在的具有亚型识别能力的基因对，再通过聚类分析识别分子亚型并验证，然后利用加权基因共表达网络分析识别胆道闭锁亚型相关的基因模块筛选亚型相关的核心基因。

在一个实施例中，如图1所示，提供了一种胆道闭锁潜在分子亚型及其核心基因的识别方法，具体包括：

步骤102，从基因表达综合数据库获取胆道闭锁组织原始表达谱数据，对表达谱数据进行预处理得到包含正常组织样本、疾病组织样本的背景基因；

从基因表达综合数据库GEO获取基因表达谱数据集GSE46960和GSE15235，利用软件对基因表达谱数据集进行背景校正，并将探针ID与基因ID对应，得到包含正常组织样本、疾病组织样本的背景基因。

步骤104，将背景基因两两组合，比较正常组织样本、疾病组织样本中基因对的相对表达秩次关系差异，确定具有潜在疾病分子亚型鉴别能力的基因对，根据具有潜在疾病分子亚型鉴别能力的基因对识别疾病分子亚型；

将背景基因i和j组合得到基因对i，j，根据i和j的表达水平大小关系G

步骤106，基于预先设计的加权基因共表达网络提取疾病分子亚型相关的基因模块作为候选模块；

步骤108，对候选模块内的基因进行蛋白与蛋白之间的互作分析，筛选得到核心基因。

上述一种胆道闭锁潜在分子亚型及其核心基因的识别方法中，样本内基因间相对表达秩次关系识别潜在疾病分子的亚型，摆脱了基因表达绝对值的限制，可以整合来自不同实验室的样本进行生物信息学分析，对基因表达水平的个体间生物变异具有鲁棒性。

在一个实施例中，步骤：基于基因间相对表达秩次关系识别BA分子亚型，具体包括：

(1)筛选正常组织样本中基因间相对表达稳定的的基因对

将正常组织样本中的基因两两组合，根据公式(1)，满足P(G

其中n表示所有的正常样本数，t代表第t个样本，G

(2)基于稳定基因对筛选疾病组织样本中的逆转基因对

基于正常样本中得到的稳定基因对，评估其在疾病中相对表达秩次关系改变情况。对于稳定对G

(3)筛选具有潜在疾病分子亚型鉴别能力的基因对

对于在疾病中筛选得到的逆转基因对，按如下步骤进一步筛选具有潜在疾病分子亚型鉴别能力的基因对：(1)逆转对的秩次关系的改变只发生在20-80％的疾病患者之中发生，确保基因对具有鉴别疾病分子亚型的能力。(2)对于任意逆转对G

(4)识别分子亚型

首先，按照分子亚型候选基因对对扰动通路的数目，将基因对进行降序排列。然后，构建一个分子亚型候选基因对×疾病样本的矩阵，其中元素r

随机实验表明，对64例BA样本分别进行聚类分析，前3000对，4000对，5000对时的聚类结果与随机情况一致。因此，将前3000对基因对的聚类结果作为BA分子亚型识别依据(图2A)。聚类结果显示64例BA样本明显聚集为两类：一类样本个数为54个，另一类样本个数为10个。

为了进一步验证BA中存在不同的分子亚型，我们把得到的3000个亚型基因对应用于GSE15235数据集中，对其聚类结果(图2B)所示，GSE15235数据集中43例BA样本中依然明显的聚集为两大类。通过数据集GSE46960，GSE15235间两类样本间的差异基因比较发现：在GSE46960两类间识别了436个差异基因，GSE15235两类间识别个584差异基因，共交叠了75差异基因(超几何分布p＝4.8*10

在一个实施例中，在执行步骤：对候选模块内的基因进行蛋白与蛋白之间的互作分析以及核心基因筛选之前，还可以对候选模块内的基因功能进行富集分析，具体为：基于加权基因共表达网络WGCNA分析提取与疾病亚型显着相关的模块。首先，对背景基因基因表达矩阵相关系数进行加权，使基因间的相互作用关系符合无标度分布。然后对基因进行分类，并将具有相似表达模式的基因分成一个模块，最小模块大小(min-ModuleSize)设置为100，其他参数设置为默认值。同一模块的基因往往表现出相似的表达模式和功能。基于京都基因和基因组百科全书KEGG数据库对模块基因的功能注释和通路富集析，用R包clusterProfile实现。

通过WGCNA包分析数据集GSE46960构建加权共表达网络。使用无标度拓扑标准选择β＝8(图3)，采用动态剪切法划分模块，合并相似度大于75％的模块，最终构建了17个共表达模块(图4)。进一步分析了BA不同亚型与各种基因模块的相关性，发现洋红色(Magenta)模块与BA亚型最相关(|r|＝0.58，p＝6*10

基于KEGG数据库，对Magenta模块中涉及的基因和富集分析结果显示(图5)：这些基因主要涉及到PI3K-Akt信号通路，粘附斑激酶通路(FAK)和ECM受体相互作用通路等。有研究表明PI3K-Akt信号通路可能抑制NF-κB和NLRP3炎症通路的蛋白表达，从而减少炎症因子的分泌，FAK和ECM受体相互作用通路抑制TGF-β1信号传导，降低成纤维细胞分化。

在一个实施例中，步骤：对候选模块内的基因进行蛋白与蛋白之间的互作分析以及核心基因筛选，具体包括：使用STRING在线网站对候选模块内的基因进行蛋白与蛋白之间的互作分析，由Cytoscape显示，Cytoscape中的cytoHubba插件用于网络最大集团度评分算法，如公式(2)，识别PPI网络中的关键基因，

MCC(ν)＝∑

其中，集合S表示集合中元素的数量，S(ν)是包含ν的最大集团的集合，(|C|-1)！是所有小于|C|的正整数的乘积。

将所有Magenta模块内的基因导入STRING数据库以构建PPI网络，通过MCC算法发现了10个核心基因(LUM，COL6A3，FBN1，SPARC，DCN，LAMA4，FAP，ANTXR1，LAMA2，COL1A2)可能是BA亚型的关键核心基因，红色显示的基因表示与其他蛋白质的相互作用更频繁(图6)。

在一个实施例中，步骤：从基因表达综合数据库获取胆道闭锁组织表达谱数据，对表达谱数据进行预处理，具体包括：从GEO下载原始表达谱数据，利用R软件Affy软件包提供的RMA算法对其进行背景校正，并利用相应的平台注释文件将探针ID对应到基因ID上。

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。