一种基于多样本数据的因果关系预测方法与流程

1.本发明涉及计算生物学和生物信息学技术领域，更具体的说，它涉及一种基于多样本数据的因果关系预测方法。


背景技术：

2.从观测数据进行因果推理是生物学、地球科学、经济学、医学、神经科学和机器学习等自然科学和工程各个研究领域的核心问题。格兰杰因果关系(gc)推理作为一种具有代表性的方法，是基于时间序列数据来推断潜在因果关系的一种方法，自1969年提出以来，基于gc的方法在许多领域得到了广泛的应用。然而，大多数生物数据是基于与时间无关的数据，如表型或阶段，而不是时间序列，因此这种方法并不适合。
3.分子之间的因果关系通常可以用因果图来表示，节点代表不同的分子，有向边表征分子之间的直接因果关系。现有的预测因果关系的主流算法包括著名的granger因果关系(gc)，收敛交叉映射算法(ccm)，传递熵(te)和贝叶斯理论。具体地说，gc基于时滞信息，并用过去信息解释的当前状态，例如基于自回归模型。作为gc方法的一种非线性版本，te考虑了信息在时间上的不对称性来确定因果关系，广泛应用于神经科学和生理学等生物系统。gc和te都是基于原始状态空间推测因果。另一方面，作为gc理论的补充，ccm是在延迟嵌入空间中测量的，反映了当前状态下两个变量之间的非线性关系。所有上述算法都是在需要时间相关数据或时间序列数据的情况下度量因果关系。相比之下，贝叶斯网络和结构因果模型(scm)能够基于统计独立性和干预处理与时间无关的数据，并能够识别分子之间的直接因果关系。然而，这些方法都依赖于有向无环图的结构来推断因果关系，这限制了对通常具有反馈环或环状相互作用的生物分子网络的应用。
4.因此，基于时间无关的数据而不是依赖于时间数据的数据，来高精度地推断真实生物/分子系统中的因果关系，从而更好地提供数据间的关系，为实际决策提供有力的数据，仍然是一个待提高和攻克的技术问题。


技术实现要素：

5.本发明克服了现有技术的不足，提供一种基于多样本数据的因果关系预测方法。
6.本发明的技术方案如下：
7.一种基于多样本数据的因果关系预测方法，不考虑时间的依赖性推断因果，具体包括如下步骤：
8.101)构建初始网络步骤：利用特定表象的表达数据作为初始输入数据，计算初始输入数据之间的相关系数并利用初始输入数据之间的相关系数构建初始网络；
9.102)判断因果关系步骤：在初始相关性网络的基础上，对每一对边进行h0假设和h1假设；h0这个边代表因果关系不存在，h1这个边代表因果关系存在；在不同的假设基础上利用部分表达数据进行不同的回归得到每对边的线性回归方程；
10.103)构建因果网络步骤：对于两边假设下不同的回归方程，带入相同的余下的表
达数据得到两个回归的误差；通过比较两者之间的拟合误差的大小，即可判断该因果关系是否存在，从而构建整个因果网络。
11.进一步的，步骤101)中利用偏相关系数评估两个变量之间的相关性，以消除基因的间接影响，从而获取基因之间的直接相关性，构建相关性网络。
12.进一步的，pearson相关系数的计算公式如下：
[0013][0014]
式中x、y为两个变量；
[0015]
一阶偏相关系数的计算公式如下：
[0016][0017]
式中r
12
表示变量1与变量2之间的相关系数，r
13
表示变量1与变量3之间的相关系数，r
23
表示变量2与变量3之间的相关系数。
[0018]
进一步的，步骤102)中h0假设中利用不包含因变量的协变量对果变量进行回归拟合；h1假设中利用因变量和协变量对果变量进行回归拟合；具体的回归公式分别如下：
[0019]
h0:
[0020]
h1:y＝f(x,z1,z2,
…
,z
n-2
) ε
[0021]
其中，y代表果变量；x代表因变量；z为其它不包含x，y变量的其他变量；ε和为噪音。
[0022]
进一步的，把特定表象的表达数据划分为训练集和测试集，训练集用作回归，测试集用作判断误差；划分方式包括k折交叉验证，留一法。
[0023]
进一步的，在不同的回归方程中，带入训练集数据，进行方程的拟合，得到方总误差，具体方式如下：
[0024]
先通过拟合方程的输出与真实值的距离作为误差，多次交叉检验得到的多次误差进行加和得到总误差，公式如下：
[0025][0026]
其中ei表示第i次交叉检验中得到的误差，m代表交叉验证的总次数。
[0027]
进一步的，步骤103)中的误差可由距离度量手段计算；比较的具体方式由因果强度的比较实现；其中，因果强度的计算公式为：
[0028][0029]
其中e代表h1假设中，包含因变量和其他随机变量对果变量的回归，得到的交叉验证的测试集的平均残差平方和，代表h0假设中，包含除因变量外的所有随机变量对果变量的回归，得到的交叉验证的测试集的平均残差平方和。
[0030]
进一步的，特定表象的表达数据包括基因表达数据、生物链数据、疾病传播模型数据。
[0031]
进一步的，构建装置包括初始网络构建模块、数据切分模块、回归拟合模块、拟合误差模块、因果判断模块和网络判断模块；
[0032]
初始网络构建模块，利用特定表象的表达数据作为原始输入数据，计算输入数据之间的相关系数并利用数据之间的相关系数构建初始网络；
[0033]
数据切分模块，将原始输入数据切分为训练集和测试集；训练集用于回归拟合模块的训练，测试集用于拟合误差模块误差的计算；
[0034]
回归拟合模块，在初始相关性网络的基础上，对每一对边进行h0假设和h1假设；h0这个边代表因果关系不存在，h1这个边代表因果关系存在；在不同的假设基础上利用部分表达数据进行不同的回归得到每对边的线性回归方程；
[0035]
拟合误差模块，用测试集输入带入回归拟合模块得到的回归方程，计算得到的拟合误差；
[0036]
因果判断模块，根据数据切分模块得到多次数据集的切分后，在拟合误差模块得到的多次误差，利用多次误差对因果强度进行判断，得到真正的因果关系；
[0037]
网络判断模块，利用以上模块方式，对初始网络中边对应的不同结点对之间的遍历，最终得到真正的因果关系网络。
[0038]
本发明相比现有技术优点在于：本发明完全提出了一个新的因果概念和建立方法，简称为cvp方法，这是一个数据驱动的无模型算法，用于处理与时间无关的数据。cvp方法通过交叉验证的可预测性和观测变量在时间独立数据上的统计独立性来量化，类似于但不同于时间相关数据上的gc。cvp方法不需要多余的先验信息，适用于多种背景的真实数据，包括基因调控网络和其他带有反馈环的因果网络。cvp方法能够有效地识别复杂生物系统中的因果关系，可以更好地揭示调控机制和解释生物功能。特别是，推导出特定的疾病因果/调控网络可以揭示分子效应的基本机制，为进一步的精准治疗提供精准数据和决策指导。
附图说明
[0039]
图1为本发明的整体框架原理图；
具体实施方式
[0040]
下面详细描述本发明的实施方式，其中自始至终相同或类似的标号表示相同或类似的元件或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明而不能作为对本发明的限制。
[0041]
本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科技术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样的定义，不会用理想化或过于正式的含义来解释。
[0042]
各实施方式中提到的有关于步骤的标号，仅仅是为了描述的方便，而没有实质上先后顺序的联系。各具体实施方式中的不同步骤，可以进行不同先后顺序的组合，实现本发明的发明目的。
[0043]
下面结合附图和具体实施方式对本发明进一步说明。
[0044]
实施例：
[0045]
如图1所示，一种基于多样本数据的因果关系预测方法，具体包括如下步骤：
[0046]
101)构建初始网络步骤：利用特定表象的表达数据作为初始输入数据，计算初始输入数据之间的相关系数并利用初始输入数据之间的相关系数构建初始网络。特定表象的表达数据包括基因表达数据、生物链数据、疾病传播模型数据等。比如正常情况下的表达数据作为疾病研究的参考数据，此类数据主要是基因表达数据，主要以rna-seq数据为主，对于其他符合标准的数据，如生物链数据，疾病传播模型数据亦可。
[0047]
具体利用皮尔逊相关系数或偏相关系数评估两个变量之间的相关性，偏相关系数去除基因的间接影响，从而获取基因之间的直接相关性，构建相关性网络。
[0048]
具体皮尔逊pearson相关系数的计算公式如下：
[0049][0050]
式中x、y为两个变量。
[0051]
即用皮尔逊相关系数来衡量两个变量之间的相关性强度，当使用pcc/ptcc构建相关性网络，网络节点表示基因，边表示基因之间的关联强度。显然，相关性的大小受许多因素的影响，如样本量、测序深度等。为了减少技术因素对不同数据集的影响，会使用自适应阈值而不是固定阈值来进行限定。如特征维度低于1000，则相关网络选择保留所有边，否则，将保留平均度数为100的高度相关的边。
[0052]
一阶偏相关系数的计算公式如下：
[0053][0054]
式中r
12
表示变量1与变量2之间的相关系数，r
13
表示变量1与变量3之间的相关系数，r
23
表示变量2与变量3之间的相关系数。
[0055]
即用统计学中偏相关系数来评价两个变量之间的直接相关性时，偏相关系数是剔除一组随机变量后得到相关网络，通过计算两基因间的ptcc，可以得到对应的假设检验p值，筛选出ptcc的p值小于0.1的边缘，从而得到偏相关网络。
[0056]
102)判断因果关系步骤：在初始相关性网络的基础上，对每一对边进行h0假设和h1假设。h0这个边代表因果关系不存在，h1这个边代表因果关系存在。在不同的假设基础上利用部分表达数据进行不同的回归得到每对边的线性回归方程。即从相关性网络中分离出的每个基因具有一阶邻域的子网络，通过本方案量化相邻基因，对该基因的因果强度(cs)。
[0057]
其中，h0假设中利用不包含因变量的其它变量对果变量进行回归拟合；h1假设中利用因变量和其它变量对果变量进行回归拟合；具体的回归公式分别如下：
[0058]
h0:
[0059]
h1:y＝f(x,z1,z2,
…
,z
n-2
) ε
[0060]
式中，y代表果变量；x代表因变量；z为其它不包含x，y变量的协变量；ε和为噪音。
[0061]
103)构建因果网络步骤：对于两边假设下不同的回归方程，带入相同的余下的表达数据得到两个回归的误差。通过比较两者之间的拟合误差的大小，即可判断该因果关系是否存在，从而构建整个因果网络。整体上可简述为基于因果强度(cs)显著性检验，剔除冗
余的因果关系，得到每个基因的最优因果子网络。整合所有基因的子网，得到最终的因果网络。
[0062]
具体需要把特定表象的表达数据划分为训练集和测试集，训练集用作回归，测试集用作判断误差；划分方式包括k折交叉验证，留一法。
[0063]
在不同的回归方程中，带入训练集数据，进行方程的拟合，得到方总误差，具体方式如下：
[0064]
先通过拟合方程的输出与真实值的距离作为误差，多次交叉检验得到的多次误差进行加和得到总误差，公式如下：
[0065][0066]
其中ei表示第i次交叉检验中得到的误差，m代表交叉检验的总次数。
[0067]
而误差可由距离度量手段计算；比较的具体方式由因果强度的比较实现；其中，因果强度的计算公式为：
[0068][0069]
其中e代表h1假设中，包含因变量和其他随机变量对果变量的回归，得到的交叉验证的测试集的平均残差平方和，代表h0假设中，包含除因变量外的所有随机变量对果变量的回归，得到的交叉验证的测试集的平均残差平方和。
[0070]
因果强度cs
x
→y的显著性检验可以定义为不包含变量x和包含变量x对y的解释能力之间的差距，即也可以用来检验e和之间的差异，即显著性检验。但是，统计检验对于推断一个大规模的网络来说太费时间了。因此利用因果强度ω
x
→y有效地计算了测试样本误差项的差异，且更适用于大规模网络。
[0071]
具体过程如下：
[0072]
1)根据初始网络，选取目标变量对，做出两种不同的假设，h0和h1假设。
[0073]
在本步骤中h0假设：其中一个变量(因变量)对另一个变量(果变量)存在因果效应。
[0074]
在本步骤中h1假设：其中一共变量(因变量)对另一个变量(果变量)不存在因果效应。
[0075]
2)利用交叉验证的方式划分数据集，把数据集划分为训练集和测试集，并根据不同的假设选择性掩盖其中因变量的数据。
[0076]
在本步骤中的选择性掩盖表示为在h0假设中，掩盖因变量的信息。在h1假设中，加入因变量的信息。
[0077]
在本步骤中的交叉验证划分数据方式包括k折交叉验证，留一法。
[0078]
3)在不同的假设中，利用训练集的数据，针对因变量建立不同的回归方程。
[0079]
在本步骤中建立的回归方程如下：
[0080]
在h0假设中的回归方程如下：
[0081][0082]
在h1假设中的回归方程如下：
[0083]
y＝f(x,z1,z2,
…
,z
n-2
) ε
[0084]
在不同的回归方程中，带入训练集数据，进行方程的拟合，得到方总误差。总误差先通过拟合方程的输出与真实值的距离作为误差，多次交叉检验得到的多次误差进行加和得到总误差。具体公式如下：
[0085][0086]
比较两种假设下得到的误差，即可确定该对因果变量间的因果关系是否存在。误差的比较方式是通过因果强度定义的。因果强度为正，认为该因果关系存在；反之，因果强度为负，认为该因果关系不存在。因果强度的具体公式如下：
[0087][0088]
对初始网络中存在连边的变量对，重复上述步骤，即可得到初始网络中所有的因果关系，即得到了最终的因果关系网络(目标网络)。
[0089]
尽管为使解释简单化将上述方法图示并描述为一系列动作，但是应理解并领会，这些方法不受动作的次序所限，因为根据一个或多个实施例，一些动作可按不同次序发生和/或与来自本文中图示和描述或本文中未图示和描述但本领域技术人员可以理解的其他动作并发地发生。
[0090]
本领域技术人员将进一步领会，结合本文中所公开的实施例来描述的各种解说性逻辑板块、模块、电路、和算法步骤可实现为电子硬件、计算机软件、或这两者的组合。为清楚地解说硬件与软件的这一可互换性，各种解说性组件、框、模块、电路、和步骤在上面是以其功能性的形式作一般化描述的。此类功能性是被实现为硬件还是软件取决于具体应用和施加于整体系统的设计约束。技术人员对于每种特定应用可用不同的方式来实现所描述的功能性，但这样的实现决策不应被解读成导致脱离了本发明的范围。
[0091]
结合本文所公开的实施例描述的各种解说性逻辑板块、模块、和电路可用通用处理器、数字信号处理器(dsp)、专用集成电路(asic)、现场可编程门阵列(fpga)或其它可编程逻辑器件、分立的门或晶体管逻辑、分立的硬件组件、或其设计成执行本文所描述功能的任何组合来实现或执行。通用处理器可以是微处理器，但在替换方案中，该处理器可以是任何常规的处理器、控制器、微控制器、或状态机。处理器还可以被实现为计算设备的组合，例如dsp与微处理器的组合、多个微处理器、与dsp核心协作的一个或多个微处理器、或任何其他此类配置。
[0092]
结合本文中公开的实施例描述的方法或算法的步骤可直接在硬件中、在由处理器执行的软件模块中、或在这两者的组合中体现。软件模块可驻留在ram存储器、闪存、rom存储器、eprom存储器、eeprom存储器、寄存器、硬盘、可移动盘、cd-rom、或本领域中所知的任何其他形式的存储介质中。示例性存储介质耦合到处理器以使得该处理器能从/向该存储介质读取和写入信息。在替换方案中，存储介质可以被整合到处理器。处理器和存储介质可驻留在asic中。asic可驻留在用户终端中。在替换方案中，处理器和存储介质可作为分立组件驻留在用户终端中。
[0093]
在一个或多个示例性实施例中，所描述的功能可在硬件、软件、固件或其任何组合中实现。如果在软件中实现为计算机程序产品，则各功能可以作为一条或更多条指令或代
码存储在计算机可读介质上或藉其进行传送。计算机可读介质包括计算机存储介质和通信介质两者，其包括促成计算机程序从一地向另一地转移的任何介质。存储介质可以是能被计算机访问的任何可用介质。作为示例而非限定，这样的计算机可读介质可包括ram、rom、eeprom、cd-rom或其它光盘存储、磁盘存储或其它磁存储设备、或能被用来携带或存储指令或数据结构形式的合意程序代码且能被计算机访问的任何其它介质。任何连接也被正当地称为计算机可读介质。例如，如果软件是使用同轴电缆、光纤电缆、双绞线、数字订户线(dsl)、或诸如红外、无线电、以及微波之类的无线技术从web网站、服务器、或其它远程源传送而来，则该同轴电缆、光纤电缆、双绞线、dsl、或诸如红外、无线电、以及微波之类的无线技术就被包括在介质的定义之中。如本文中所使用的盘(disk)和碟(disc)包括压缩碟(cd)、激光碟、光碟、数字多用碟(dvd)、软盘和蓝光碟，其中盘(disk)往往以磁的方式再现数据，而碟(disc)用激光以光学方式再现数据。上述的组合也应被包括在计算机可读介质的范围内。
[0094]
提供对本公开的先前描述是为使得本领域任何技术人员皆能够制作或使用本公开。对本公开的各种修改对本领域技术人员来说都将是显而易见的，且本文中所定义的普适原理可被应用到其他变体而不会脱离本公开的精神或范围。由此，本公开并非旨在被限定于本文中所描述的示例和设计，而是应被授予与本文中所公开的原理和新颖性特征相一致的最广范围。