一种用于定量分析药物干预前后生物分子网络中模块变化的方法与流程

本发明属于生物信息技术领域。具体而言，本发明涉及复杂生物分子网络如蛋白质相互作用网络、基因表达调控网络、药物代谢和药物靶点等网络中的模块在经药物干预后由一种构象转变为另一种构象的量化方法。

背景技术：

网络结构分析是将一个感兴趣的系统表示成网络，进而采用量化指标来研究该网络的拓扑结构特性的过程。模块作为网络结构层次的中等粒度，代表了网络的局部特征和网络构成。因此，可以将网络拓扑参数度量方法用于对模块拓扑特征的描述。

网络药理学认为当生物系统为稳态和平衡态时，机体处于健康状态。当生物系统(生物网络)的平衡(健康状态)被扰乱或破坏时，则会导致病理或疾病状态。而药物对疾病的治疗作用可以看作生物系统或生物网络对外部扰动信息的应答，即有效药物的作用将使平衡向能够减弱这种改变的方向移动，在于重建生物系统的平衡或减轻平衡被破坏的程度。然而，药物的影响是否能有效地体现在网络的全局改变上还值得商榷，有研究认为聚焦于网络局部及片段(模块)的改变可能会更直接有效地评估药物的影响。

既往研究对模块改变程度的测量多从某一个或几个拓扑参数分别进行量化分析，忽略了对模块整体(多维)构象改变程度的综合定量分析。因此目前本领域仍然需要提供能够综合量化模块整体构象变化程度的新方法。

技术实现要素：

针对上述技术问题，本发明提供一种能够定量分析生物分子网络中模块整体构象变化程度的新方法，该方法整合模块不同维度(即属性)的拓扑参数，构建模块综合评价体系，能够定量分析药物干预前后模块整体构象的变化程度。

定义

本文采用的术语“生物分子网络”是指在生物系统中以不同组织形式存在的生物系统网络，其由代表各种生物分子的节点与代表该生物分子之间的相互作用关系的边组成。常见的生物分子网络包括蛋白质相互作用网络、基因共表达网络、基因转录调控网络、生物代谢网络、表观遗传网络、表型网络、信号传导网络等。

本文采用的术语“节点”是指生物分子网络中不同的个体生物分子。如蛋白质相互作用网络中，网络中的节点表示个体蛋白质。

本文采用的术语“边”是指生物分子网络中不同的个体生物分子(节点)之间具有的某种特定的关系，在生物分子网络中两个个体生物分子之间的关系用边表示。如蛋白质相互作用网络中，蛋白质(节点)之间相互关联、相互影响的关系通过它们之间的边来表示。

本文采用的术语“模块”是指生物分子网络中由至少三个不同的生物分子(节点)组成的，在功能上和形态上相对独立的实体。在网络拓扑结构上，模块内部的节点连接表现为相对稠密，而模块与模块之间的节点连接表现为相对稀疏。模块具有独立的功能，这些功能来源于组成它们的个体生物分子(节点，例如蛋白质、DNA、RNA和小分子)之间的相互作用。

本文采用的术语“模块变化”是指在不同的时空间(如药物干预前后)，模块由于外界扰动(例如药物干预)而进行的改变(如分裂、合并、减缩、增长等)。

本文采用的术语“(模块的)维度”是指模块拓扑结构的不同属性，可以用参数如邻居节点、网络密度、网络中心性、介数中心性等进行描述。在本文中，“维度”与“属性”可互换使用。

本文采用的术语“变化模块对”是指由分别来自药物干预前的生物分子网络的一个模块及药物干预后的生物分子网络的一个模块共同组成的模块对，所述两个模块至少有一个节点(≥1)重叠。

本文采用的术语“有效变化模块对”是指药物干预前后其变化程度超过一定阈值的变化模块对。

本文采用的术语“非有效药物”是指在具体目标病理状态下，应用后未达到预期改善效果的药物。

本发明的具体技术方案如下。

本发明提供了一种用于定量分析药物干预前后生物分子网络中模块变化的方法，所述方法包括以下步骤：

(1)以节点数目≥3为标准，分别对药物干预前和药物干预后的生物分子网络进行模块识别，以获得药物干预前的生物分子网络的一组模块和药物干预后的生物分子网络的一组模块；

优选地，进行模块识别的方法选自下述中的一种或多种：MCODE、MCL、CFinder、CPM、SPC、G-N algorithm、ModuLand、WGCNA、DME、MINE、SVD等。根据本发明的具体实施方式，在本发明的方法中进行模块识别的方法为MCODE或WGCNA。

优选地，所述生物分子网络为蛋白质相互作用网络、基因共表达网络、基因转录调控网络、生物代谢网络、表观遗传网络、表型网络、信号传导网络等。更优选地，所述生物分子网络为蛋白质相互作用网络或基因共表达网络。

(2)将经步骤(1)识别的药物干预后的生物分子网络的一组模块与药物干预前的生物分子网络的一组模块相匹配，获得一个或多个变化模块对。

(3)对于经步骤(2)获得的一个或多个变化模块对，采用整体综合指标k值对模块变化程度进行量化，其中，整体综合指标k值由描述不同维度的拓扑参数及其权值计算得到，计算公式如下：

所述公式中n为选择的拓扑参数的个数，i为1至n的整数，a_i为模块组中药物干预前的模块的第i个拓扑参数的空间向量特征，b_i为模块组中药物干预后的模块的第i个拓扑参数的空间向量特征，di为第i个拓扑参数的空间向量特征在药物干预前后的距离，w_i为每个拓扑参数的权值。

其中，本发明方法的步骤(3)中构建量化模块变化的整体综合指标的不同维度拓扑参数可以选自表1：

表1：模块不同维度拓扑参数

优选地，n≥5，优选n≥8，更优选n≥10。

其中，可以用层次分析法、特尔斐法法、因子分析权数法、信息量权数法、主成分分析法、熵权法、优序图法、标准离差法等确定权值。根据本发明的具体实施方式，用分层次分析法确定权值。

上述方法的技术路线图见图8。

另一方面，本发明提供一种用于鉴定药物干预前后生物分子网络中有效变化模块对的方法，所述方法包括以下步骤：

(A)执行上述方法的步骤(1)至步骤(3)，以获得表示模块变化程度的整体综合指标k值；

(B)采用非有效药物的物质替换步骤(1)中的药物，重复上述方法的步骤(1)至步骤(3)，以获得一个或多个变化模块对中表示模块变化程度的整体综合指标k值，选择最大k值作为阈值；

(C)将步骤(A)中获得的k值与步骤(B)中获得的阈值相比较，如果k值≥阈值，则将具有所述k值的变化模块对鉴定为有效变化模块对。

本发明的模块变化量化方法可以用于生物网络分析、新药研发设计、药理机制研究等领域。例如根据本发明的具体实施方式，采用精制清开灵的五个有效组分干预小鼠脑缺血模型，可以获得疾病状态下与药物干预后的分别的生物分子网络，通过本发明的方法进行模块识别、匹配、模块变化的量化，以及与非有效药物珍珠母进行比较，可以获得药物干预前后的变化模块对与有效变化模块对，为药物开发与药理机制研究提供相关信息。

附图说明

以下，结合附图来详细说明本发明的实施方案，其中：

图1为实施例1中各组网络构建结果。网络中节点表示基因或蛋白质，节点之间的连线表示蛋白质之间的相互作用。a.Vehicle组；b.BJ组；c.JU组。

图2为实施例1中采用MCODE方法对各个网络进行模块识别的结果。a.Vehicle组；b.BJ组；c.JU组。

图3为实施例1中药物干预前(A状态)到药物干预后(B状态)模块变化的整体综合指标k值的计算公式。

图4为实施例1中Vehicle-BJ组有效变化模块对。

图5为实施例1中Vehicle-JU组有效变化模块对。

图6为实施例4中Vehicle-black与JA-cyan变化模块对的k值。

图7为实施例4中Vehicle-violet与JA-lightyellow、Vehicle-violet与JA-red、Vehicle-violet与JA-yellow变化模块对的k值。

图8为本发明的技术路线图。

具体实施方式

以下参照具体的实施例来说明本发明。本领域技术人员能够理解，这些实施例仅用于说明本发明，其不以任何方式限制本发明的范围。

本发明的目的是量化不同状态下(如药物干预前(有可能为疾病状态)和药物干预后的生物分子网络)模块的变化程度，从而识别出对药物干预起到贡献性响应的有效变化模块，为指导疾病治疗和药物研发提供依据。以下的实施例证明了本发明方法的有效性和可行性。这些实施例是非限制性的，本发明的方法还可以应用其他类型的网络。

下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的药材原料、试剂材料等，如无特殊说明，均为市售购买产品。

实施例1采用清开灵有效组分干预小鼠脑缺血模型，量化蛋白质相互作用网络中模块的变化，以及识别有效变化模块对(10个拓扑参数融合)

数据来源

本实施例数据来源于：利用精制清开灵五个有效组分即黄芩苷(BA)、栀子苷(JA)、胆酸(UA)、黄芩苷+栀子苷(BJ；BA+JA)、栀子苷+胆酸(JU；JA+UA)及无效组分珍珠母(CM)对小鼠脑缺血再灌注损伤模型进行干预后，采用Ingenuity Pathway Analysis(IPA)对其显著差异表达基因进行分析，提取IPA生物功能富集结果中与生物功能注释相关的有统计学意义的基因集。其中模型组(Vehicle)包含149个有统计学意义的基因，BA组74个，JA组121个，UA组104个，BJ组70个，JU组107个，CM组40个。将这7个组的基因作为目标基因分别映射于全局背景(全局的小鼠基因和蛋白质相互作用数据)。

本实施例用模型组作为药物干预前的数据，具有加合效应的BJ组与具有协同效应的JU组作为药物干预后的数据。

参见图1，脑缺血模型组(Vehicle组)的蛋白质相互作用网络(图1-a)，由3750个节点和9162条边组成；

BJ组的蛋白质相互作用网络(图1-b)，由2968个节点和6273条边组成；

JU组的蛋白质相互作用网络(图1-c)，由3429个节点和8111条边组成。

药物干预前后模块变化的定量分析过程如下：

步骤1，采用MCODE方法(模块节点数目≥3)对各个网络进行模块识别，结果如图2所示：由脑缺血模型组(Vehicle)的蛋白质相互作用网络识别的模块见图2-a。由BJ组的蛋白质相互作用网络识别的模块见图2-b，由JU组的蛋白质相互作用网络识别的模块见图2-c。

步骤2，追踪药物干预前后网络模块的变化情况，在本实施例中，将药物组模块(BJ、JU)分别与Vehicle组模块相匹配，至少有一个节点的重叠定义为一个变化模块对(modular reconstructional pairs)。最终，BJ、JU与Vehicle组变化模块对匹配结果如表2、3所示：

表2-Vehicle与BJ变化模块对匹配

表3-Vehicle与JU变化模块对匹配

步骤3，应用量化模块变化的整体综合指标(多个指标的融合，同时组合表示具体拓扑参数的重要程度的权重，形成综合的模块拓扑参数指标)对模块的变化程度进行量化分析。在本实施例中，采用非重叠节点、非重叠边、重叠节点、重叠边、平均邻居节点、网络密度、网络中心性、平均介数中心性、网络平均权重、最短路径这10个代表模块不同维度的拓扑参数进行融合(在排除不同参数间共线性的基础上还可以选择其他参数)。将在药物干预前(Vehicle组)，模块A的结构设置为a，对于它的拓扑状态，应用这10个参数(不同属性变量)将它的空间向量特征描述为a₁、a₂、a₃……a₁₀；药物干预后(BJ或JU组)，模块A的拓扑结构发生变化，转变为模块B，相对应的拓扑结构b描述为b₁、b₂、b₃……b₁₀。同时组合表示不同参数的重要程度的相应的权值(层次分析法)，如表4所示。将模块从A到B的综合拓扑参数变化用k值表示，如图3所示。k取值在0-1之间，越小证明两个模块的整体差异越小；反之，k值越大证明两个模块间的差异越大。

表4各变量代表参数及权值(10个参数)

应用图3描述的k值方法，对模块从药物干预前的疾病状态到药物干预后状态的变化程度表述如下，见表5(Vehicle--BJ)、表6(Vehicle--JU)

表5Vehicle--BJ变化模块对k值(10个参数)

可以看出，Vehicle—BJ各变化模块对中最大k值为0.587(Vehicle4-BJ7)，最小k值为0.180(Vehicle2--BJ2)，平均0.355。k值大部分在0.3-0.5之间。

表6Vehicle--JU变化模块对k值(10个参数)

Vehicle--JU各变化模块对中最大k值为0.556(Vehicle12--JU15)，最小k值为0.137(Vehicle16--JU4)，平均0.326。

步骤4，本实施例中，在经步骤3对Vehicle--BJ、Vehicle--JU模块变化程度量化分析的基础上，利用阈值k₀来识别出能够对药物干预疾病网络产生贡献性响应的有效变化模块。其中阈值k₀通过如下方式获得：

利用无效药物珍珠母(CM)干预脑缺血网络后，与模型组相比，珍珠母组模块变化的最大k值为0.456，将其作为阈值，Vehicle--BJ、Vehicle--JU中任意变化模块对ki≥0.456，则认为是有效变化模块对。

表7Vehicle--CM变化模块对k值(10个参数)

图4中为Vehicle--BJ组确定的3对有效变化模块对，分别是:a.Vehicle4--BJ2、Vehicle4--BJ7；b.Vehicle11--BJ16。图5中为Vehicle--JU组确定的3对有效变化模块对，分别是：a.Vehicle12--JU15；b.Vehicle16--JU19；c.Vehicle20--JU13。

实施例2采用清开灵有效组分干预小鼠脑缺血模型，量化蛋白质相互作用网络中模块的变化，以及识别有效变化模块对(8个拓扑参数融合)

采用与实施例1相同的网络数据。在本实施例中，步骤1和步骤2与实施例1相同。

步骤3，采用8个参数，即非重叠节点、非重叠边、重叠节点、重叠边、网络密度、网络中心性、网络平均权重、最短路径这8个代表模块不同维度的拓扑参数进行融合，同时各参数的权值见表8。应用实施例1中描述的k值方法，对模块从疾病状态到药物干预后JU组状态的变化程度表述如下：

表8各变量代表参数及权值(8个参数)

表9Vehicle--JU变化模块对k值(8个参数)

发现Vehicle-JU各变化模块对中拥有最大k值的仍为Vehicle12--JU15；拥有最小k值仍为Vehicle16--JU4，与实施例1相同，平均0.355。

本实施例中，采用相同方法利用清开灵组分无效药物组-珍珠母(CM)得到阈值k₀，识别出能够对药物干预疾病网络产生贡献性响应的有效变化模块对。在8个拓扑参数标融合的量化指标中，CM组干预脑缺血网络后模块变化的最大k值为0.535(表10)，Vehicle--JU中任意变化模块对ki≥0.535，则认为是有效的变化模块对。

表10Vehicle--CM变化模块对k值(8个参数)

Vehicle--JU组确定的3对有效变化模块对分别是Vehicle12--JU15(k＝0.640)、Vehicle16--JU19(k＝0.543)、Vehicle20--JU13(k＝0.573)，与实施例1确定的有效变化模块对相同。

实施例3采用清开灵有效组分干预小鼠脑缺血模型，量化蛋白质相互作用网络中模块的变化，以及识别有效变化模块对(5个拓扑参数融合)

采用与实施例1相同的网络数据。在本实施例中，步骤1和步骤2与实施例1相同。

步骤3，采用5个参数，即非重叠节点、非重叠边、重叠节点、重叠边、网络中心性这5个代表模块不同维度的拓扑参数进行融合，同时各参数的权值见表11。应用实施例1中描述的k值方法，对模块从疾病状态到药物干预后BJ组状态的变化程度表述如下：

表11各变量代表参数及权值(5个参数)

表12Vehicle--BJ变化模块对k值(5个参数)

发现Vehicle-BJ各变化模块对最大k值为0.858(Vehicle11--BJ16)，最小k值为0.255(Vehicle2--BJ7)，平均0.496。

本实施例中，采用相同方法利用清开灵组分无效药物组-珍珠母(CM)得到阈值k₀，识别出能够对药物干预疾病网络产生贡献性响应的有效变化模块对。在5个拓扑参数标融合的量化指标中，CM组干预脑缺血网络后模块变化的最大k值为0.722(表13)，Vehicle-JU中任意变化模块对ki≥0.722，则认为是有效的变化模块对。

表13Vehicle--CM变化模块对k值(5个参数)

Vehicle-BJ组确定的3对有效变化模块对分别是Vehicle4--BJ2(k＝0.733)、Vehicle4--BJ7(k＝0.850)、Vehicle11--BJ16(k＝0.858)。本实施例确定的有效变化模块对与实施例1确定的有效变化模块对相同，只是三个有效模块对的k值大小排序稍有差别：在实施例1中，三个有效变化模块k值从大到小依次为Vehicle4--BJ2、Vehicle11--BJ16、Vehicle4--BJ7；而本实施例中Vehicle11--BJ16模块对的k值大于Vehicle4--BJ7。

实施例4采用清开灵有效组分干预小鼠脑缺血模型，量化基因共表达网络中模块的变化

本实施例以精制清开灵有效组分栀子苷(JA)干预小鼠脑缺血模型的基因表达谱数据为例，具体实施过程如下：

1.数据来源：

将模型组基因表达谱数据作为药物干预前的数据，JA组基因表达谱数据作为药物干预后的数据。每组基因表达谱数据都由12个样本的374个基因组成(Tbp、Zeb1、Pou2f1、Foxb1、Creb1、Camk2g、Csf1、F5、Hspd1、Matn2、Mt1、Adamts1、Klf6、Dffa、Rgs18、Rhoa、Kcnmb1、Pdcd11、Pdpk1、Casp8ap2、Mogat1、Rps26、Ak1、Csnk2a2、Dkk2、Ppm1e、Tnfrsf22、Trp53i11、Smpd3、Grin1、Cdk5、Jund、E2f1、Apoe、Il1b、Prkar1b、Il7r、Ngfb、Rela、Ifnar1、Adcy6、Bak1、Fzd6、Prkch、Rgs4、Actg1、Gck、Rgs9、Sox9、Rgs1、Dgke、Rgs20、Map2k2、Pin1、Prkcn、Dgkz、Csnk1g1、Dusp4、Il11、Grb2、Shc1、Syk、Sim2、Ywhah、Fgf13、Bid、Gstm2、Rarg、Pou3f1、Camk2b、Mapkapk2、Tcf4、Sos1、Stat5a、Vegfb、Bad、Etv3、Id1、Lcat、Nf1、Gsn、Bbc3、Clu、Capn9、Ercc5、Comt、Ctsl、Amph、Vegfc、Bax、Cyp51、Sox10、Nfyc、Gata2、Id3、Lef1、Pou6f1、6330503C03Rik、Ech1、Ccl4、Itm2a、Hspa1a、Cbx3、Klf10、Idh3g、Gpx2、Map2k5、Daxx、E2f3、Fgf12、Ikbkg、Btrc、Ikbkap、Ifnar2、Cdk5、Psmb1、Sufu、Gab1、Sox30、Pxn、Pygo2、Ctnnb1、Grin2a、Il5ra、Cdk4、Bcl2l1、Actb、Myb、Prkca、Csf2rb2、Gnaq、B-raf、Wnt6、Adcy7、Cacna1b、Fzd7、Prkcm、Rock1、Adcy8、Prkcc、Sub1、Tuba1b、Rgs6、Plcb1、Mknk1、Diablo、Mef2c、Lrp1b、Dgkg、Rgs12、Serpina5、Hspb1、Ppm1b、Dlk1、Cdc42、Fadd、Mdfi、Fgf11、Map3k4、Klk1b3、Il6ra、Tgfb2、Wnt11、Ccna1、Map2k6、Htr1f、Zmat3、Bnip3、Tsg101、Vim、Srf、D14Abb1e、Cdh11、Vdac2、Tfdp1、Gak、Ccna2、Vegfa、Vegfa、Hdac1、Srebf1、Stch、E2f1、Nfatc1、Gna12、Gna13、Cacnb3、Zic1、Pou4f3、Tcf12、Ldb1、Capns1、Fxyd2、Gcgr、LOC100304588、Syt11、Gadd45a、Pbx2、Ier3、Mapk9、Ctnnbip1、Fgf15、Smad3、Nlk、Mecp2、Sigirr、Rgs18、Ptk2b、Sap30bp、Pcmt1、Tcf3、Braf、Ankrd6、Rgs5、Rap1gap、Adcy1、Grin2b、Gap43、Map2k1、Mapk10、Tgfb1、Lta、Rps6ka1、Wnt3、Rara、Prkcd、Atf4、Adcyap1r1、Cycs、Hint1、Rdx、Src、Adcy9、Prkce、Shcbp1、Elk3、Rgs14、Rgs17、Dusp10、Tubb3、Cyc1、Dusp16、Plcg2、Fzd10、Dgkd、Stat3、Mapk14、Map2k4、Htr1a、Map3k2、Frat1、Casp7、Eef2k、Thbd、Rarb、Camk4、Htr2c、E2f5、Met、Htr7、Camk2b、Stat6、Sod1、Efna4、Vdac3、Adora1、Bmp1、Vdac1、Grb2、Igfbp2、Top2b、Rpl35、Bdnf、Ppp3cb、Raf1、Cpe、Cacnb3、0610007C21Rik、Gna14、Gna11、Tuba1a、Zic3、Mlx、Id4、Ldb2、Sepp1、Prodh、S100a9、Pgam2、Rcan1、Abcc5、Ccr5、Ap1m1、Map3k5、Csnk1e、Axin1、Freq、Sh2b1、Rps6ka4、Wif1、Nkd1、Pam、Crem、Tgm2、Barhl1、Tradd、Plcd4、Ppp2r4、Otud7b、Rgs7、Casp2、Junb、Il2rg、Bad、Il1a、Egr1、Pdgfa、Gapdh、Eif4e、Apc、Prkcz、Parp1、Egfr、Prkcb1、Rgs2、Traf2、Ccr3、Rgs16、Smpd1、Tbp、Dgka、Mos、B230120H23Rik、Eif4e2、Rgs19、Adcy3、Creb5、Taf7、Pik3ca、Stat1、Il15、Atf3、Dvl3、Map3k3、Casp4、Kcnq1、Ptp4a3、fosB、Wnt3a、Calm1、Htr3a、Crkl、Casp3、Lhx1、Camk4、Selenbp2、Tcfe2a、Scg5、Pold3、Mmp2、Farp2、Pold2、Pold1、Gpx4、App、Mlh3、Rbl2、Tpp2、Cdh3、Fmo2、Pold4、Arf1、Sox1、Arhgef1)

2.网络构建和模块划分

用加权共表达网络分析(WGCNA)工具分别构建模型组(Vehicle)及JA组数据的基因共表达网络(节点数目均为374个)并划分模块(将最小模块设定为三个节点)。Vehicle组得到48个模块，JA组得到42个模块。选择模型组和JA组的四个变化模块对进行k值分析，分别是：Vehicle-black与JA-cyan；Vehicle-violet与JA-lightyellow；Vehicle-violet与JA-red；Vehicle-violet与JA-yellow。应用实施例1中10个拓扑参数的融合进行模块变化的量化追踪，四个变化模块对相应的k值如表14、图7和图8所示(虚线圆形强调的节点为变化模块对重叠的节点)。

表14Vehicle-JA部分变化模块对k值(基因共表达网络)

以上对本发明具体实施方式的描述并不限制本发明，本领域技术人员可以根据本发明作出各种改变或变形，只要不脱离本发明的精神，均应属于本发明所附权利要求的范围。