一种基因网络数量性状定位检测方法和系统与流程

1.本发明涉及数量性状基因座定位技术，具体涉及一种基于网络的数量性状基因座定位方法和系统。


背景技术：

2.生物信息学是生命科学和计算机科学交叉领域的一个研究热点。后基因组时代的一个巨大挑战是如何解释和解决遗传变异对不同生物表型的影响。全基因组关联研究(gwas)是一种从巨大的遗传变异中识别复杂性状潜在遗传位点的方法。gwas可以在全基因组范围内计算出存在的序列变异，即单核苷酸多态性(snp)，从中筛选出与疾病相关的snp位点。但由于gwas研究的设计方法和遗传统计方法本身存在一些无法避免的问题，例如人群混杂、多重比较，这些影响会造成结果的假阳性，因此需要通过重复研究来确保遗传标记与疾病间存在真关联。
3.表达数量性状基因座(eqtl)可以整合基因型数据和基因表达数据，寻找控制表达数量性状的基因在基因组中的位置，这为探究基因调控的分子机制和调控网络开启了新的思路。但传统的方法会存在一些问题，例如会忽略snp位点或基因的系统作用，从而忽略了许多网络相关的表型决定因素。


技术实现要素：

4.以下给出一个或多个方面的简要概述以提供对这些方面的基本理解。此概述不是所有构想到的方面的详尽综览，并且既非旨在指认出所有方面的关键性或决定性要素亦非试图界定任何或所有方面的范围。其唯一的目的是要以简化形式给出一个或多个方面的一些概念以为稍后给出的更加详细的描述之序。
5.本发明的目的在于解决上述问题，提供了一种基因网络数量性状定位检测方法和系统，克服现有的表达数量性状基因座方法存在的忽略基因相互作用的不足，不仅可以识别用于分析复杂生物学过程的边对/网络特征，而且还可以识别表征每种表型及其子类型的边对/网络特征。
6.本发明的技术方案为：本发明揭示了一种基因网络数量性状定位检测方法，方法包括：
7.步骤1：准备网络qtl模型的数据集；
8.步骤2：基于数据集进行网络qtl模型的基因边对矩阵的预构造；
9.步骤3：对构建的基因边对矩阵进行网络qtl模型的关联分析，确定候选网络数量性状定位；
10.步骤4：基于关联分析筛选出网络数量性状定位的关联矩阵，导出网络数量性状定位关联的基因网络模块；
11.步骤5：基于导出的数量性状定位关联的基因网络模块进行生物功能和生物特征分析。
12.根据本发明的基因网络数量性状定位检测方法的一实施例，步骤1中，构建多组模拟数据分别用于表达数量性状基因座、网络数量性状基因座的计算，这些数据集包括基因型数据集、分子表型数据集、生理表型数据集和包括协变量在内的可选数据集，以及包括蛋白质互作网络在内的背景生物网络信息；且所有类型的数据都在相同的样本有对应的数值。
13.根据本发明的基因网络数量性状定位检测方法的一实施例，步骤2进一步包括：
14.利用单样本网络计算基因边对值；
15.再利用蛋白质互作网络对计算出的基因边对进行筛选。
16.根据本发明的基因网络数量性状定位检测方法的一实施例，步骤3中的网络qtl模型的关联分析所需要的数据由两个数据矩阵组成：一个是基因序列变异数据矩阵，其中每一行代表一个snp，每一列代表一个样本，每一个元素代表一个样本中一个snp的基因型；另一个是基因边对矩阵，其中一行表示一个基因边对，一列表示与基因序列变异数据矩阵中相同的匹配样本，每个元素表示一个样本中一个基因边对的共表达水平。
17.根据本发明的基因网络数量性状定位检测方法的一实施例，步骤4进一步包括：
18.基于步骤3的关联分析结果，筛选出热点snp位点和热点基因对：如某一snp位点与基因边对的交互关系比其他snp位点更显著，或者某一基因边对比其他基因边对有更多的snp交互关系，则筛选出这样的snp位点或者基因边对；
19.基于筛选出的snp位点与基因边对的关联，获取边对网络特征，导出网络数量性状定位关联的基因网络模块。
20.根据本发明的基因网络数量性状定位检测方法的一实施例，步骤5进一步包括：
21.通过go，kegg富集分析或ipa，确定每个基因网络模块的生物学功能；通过偏最小二乘判别分析或典型相关分析推断网络特征和生理表型之间的联系。
22.本发明还揭示了一种基因网络数量性状定位检测系统，系统包括：
23.数据集准备模块，准备网络qtl模型的数据集；
24.矩阵预构造模块，基于数据集进行网络qtl模型的基因边对矩阵的预构造；
25.关联分析模块，对构建的基因边对矩阵进行网络qtl模型的关联分析，确定候选网络数量性状定位；
26.关联矩阵筛选模块，基于关联分析筛选出网络数量性状定位的关联矩阵，导出网络数量性状定位关联的基因网络模块；
27.分析模块，基于导出的数量性状定位关联的基因网络模块进行生物功能和生物特征分析。
28.根据本发明的基因网络数量性状定位检测系统的一实施例，数据集准备模块配置为构建多组模拟数据分别用于表达数量性状基因座、网络数量性状基因座的计算，这些数据集包括基因型数据集、分子表型数据集、生理表型数据集和包括协变量在内的可选数据集，以及包括蛋白质互作网络在内的背景生物网络信息；且所有类型的数据都在相同的样本有对应的数值。
29.根据本发明的基因网络数量性状定位检测系统的一实施例，矩阵预构造模块配置为先利用单样本网络计算基因边对值，再利用蛋白质互作网络对计算出的基因边对进行筛选。
30.根据本发明的基因网络数量性状定位检测系统的一实施例，关联分析模块中对网络qtl模型的关联分析所需要的数据由两个数据矩阵组成：一个是基因序列变异数据矩阵，其中每一行代表一个snp，每一列代表一个样本，每一个元素代表一个样本中一个snp的基因型；另一个是基因边对矩阵，其中一行表示一个基因边对，一列表示与基因序列变异数据矩阵中相同的匹配样本，每个元素表示一个样本中一个基因边对的共表达水平。
31.根据本发明的基因网络数量性状定位检测系统的一实施例，关联矩阵筛选模块配置为基于关联分析模块的关联分析结果，筛选出热点snp位点和热点基因对：如某一snp位点与基因边对的交互关系比其他snp位点更显著，或者某一基因边对比其他基因边对有更多的snp交互关系，则筛选出这样的snp位点或者基因边对；基于筛选出的snp位点与基因边对的关联，获取边对网络特征，导出网络数量性状定位关联的基因网络模块。
32.根据本发明的基因网络数量性状定位检测系统的一实施例，分析模块配置为通过go，kegg富集分析或ipa，确定每个基因网络模块的生物学功能；通过偏最小二乘判别分析或典型相关分析推断网络特征和生理表型之间的联系。
33.本发明对比现有技术有如下的有益效果：传统的eqtl方法往往忽略了snps或基因的系统作用，从而忽略了许多基因网络相关的表型决定因素。这个问题催生了本发明的网络qtl模型(network qtl,nqtl)，其目的是检测级联关联为基因型
→
基因网络
→
表型，而非传统eqtl方法的基因型
→
基因表达
→
表型。本发明在单样本网络的基础上开发了一个新的nqtl框架，它不仅可以识别用于分析复杂生物学过程的边对/网络特征，而且还可以识别表征每种表型及其子类型的边对/网络特征。结果表明，与传统的eqtl相比，本发明的nqtl能够更好地计算出各种模拟数据中snps与基因边对之间的关联。此外，在生物医学数据集分析中，可以发现nqtl能够有效地检测具体生物问题下的网络特征，可以发现许多网络特征对应不同的表型，这能够解释网络特征影响疾病子类型化、疾病预后、药物反应等不同因素。所有这些结果都可以说明新的nqtl模型可以整合多层次组学数据，同时检测基因型
→
基因网络
→
表型关联和具有功能重要性的潜在基因网络特征，不仅可以解释生物或生物医学意义,还可以发现nqtl作为网络与表型关联的细节功能。
附图说明
34.在结合以下附图阅读本公开的实施例的详细描述之后，能够更好地理解本发明的上述特征和优点。在附图中，各组件不一定是按比例绘制，并且具有类似的相关特性或特征的组件可能具有相同或相近的附图标记。
35.图1示出了本发明的基因网络数量性状定位检测方法的一实施例的流程图。
36.图2示出了本发明的基因网络数量性状定位检测系统的一实施例的原理图。
具体实施方式
37.以下结合附图和具体实施例对本发明作详细描述。注意，以下结合附图和具体实施例描述的诸方面仅是示例性的，而不应被理解为对本发明的保护范围进行任何限制。
38.图1示出了本发明的基因网络数量性状定位检测方法的一实施例的流程图。下面对本实施例的实现步骤进行详细描述。
39.步骤1：准备网络qtl模型的数据集。
40.构建多组模拟数据分别用于eqtl(表达数量性状基因座)、nqtl(网络数量性状基因座)的计算。这些数据集包括基因型数据集、分子表型数据集、生理表型数据集和一些可选数据集(如协变量)，以及背景生物网络信息(例如蛋白质相互作用网络ppi)。所有类型的数据都必须在相同的样本有对应的数值，例如，样本及其标签在所有数据集中是一一匹配的。
41.举例来说，构建8组模拟数据，对于前4组数据，设置基因和snp位点之间存在相关性，即为表达特征。公式如下：
42.pop＝0.2*rnorm(n)
…
(1)
43.snps＝rnorm(n*nss) pop
…
(2)
44.gene＝rnorm(n*ngs) pop snps[,sample(1:nss,ngs)]*1:ngs)/ngs
9/2
[0045]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ…
(3)
[0046]
上述公式中，n为样本数量，ngs为基因数量，nss为snp数量，pop为所有变量中共同的信号，snps表示模拟生成的基因型数据，gene表示模拟生成的基因表达数据，sample(1:nss,ngs)表示在snp数量范围内随机抽取和基因数量相同的数据，snps[,sample(1:nss,ngs)]则表示随机抽取ngs个基因型数据，rnorm()表示产生一系列的随机数。
[0047]
示例中设定样本数量n为100，基因数量ngs为2000和5000,snp数量nss分别为2000和200000。于是基因个数和snp位点个数有四种数据组合。对于其他四组合成数据，snp和基因数据是随机产生的，将基因对数据与snp，即为边缘性状，按照如下公式联系起来：
[0048]
gene＝rnorm(n*ngs)
…
(4)
[0049]
cordata＝cor pop snps[,sample(1:nss,ncs)]*1:ncs)/ncs
9/2
…
(5)
[0050]
其中cor表示基因边对值，cordata表示与基因型数据相关的基因边对值，ncs表示基因对个数。然后设置样本数量n为100，基因数量ngs和基因对数量cor分别设置为2000、1000*2001和5000、2500*5001，snp位点数量nss设置为2000和20000。同样，有四种数据组合来自基因/基因对和snp。
[0051]
步骤2：基于数据集进行网络qtl模型的基因边对矩阵的预构造。
[0052]
步骤2主要由步骤2-1和步骤2-2这两个处理过程组成，其中：
[0053]
步骤2-1：利用单样本网络计算基因边对值。对于基因边对值的计算，则利用给定原始的基因表达矩阵，群体样本中基因i和j之间的皮尔逊相关系数(pcc)可计算为：
[0054][0055]
上式中，n表示样本数，对于第k个样本，基因i和基因j在该样本上的表达值分别为x
ik
,x
jk
，σ
i
、σ
j
分别表示基因i和基因j的标准差值，e
ijk
表示基因i和基因j在样本k上的单样本基因边对值。对于基因i或者基因j，在所有样本中的表达均值分别为μ
i
，μ
j
。由此值e
ijk
代表样本中一对基因对值，与样本对应的突变数据进一步计算nqtl结果。
[0056]
步骤2-2：再利用蛋白质互作网络(ppi)对计算的基因边对进行筛选，以减少计算时间。ppi筛选的具体过程是现有技术，采用了string数据库中人类蛋白质相互作用信息表中的score值，进行筛选处理，其中留下score>900的基因边对。
[0057]
步骤3：对构建的基因边对矩阵进行网络qtl模型的关联分析，确定候选网络数量性状定位。关联分析是利用现有的matrixeqtl方法对基因序列变异数据矩阵和基因边对矩
阵进行关联分析，并得到关联矩阵。
[0058]
网络qtl模型的关联分析需要的数据由两个数据矩阵组成。一个是基因序列变异数据矩阵，其中每一行代表一个snp，每一列代表一个样本，每一个元素代表一个样本中一个snp的基因型。另一个是基因边对矩阵，其中一行表示一个基因边对，一列表示与基因序列变异数据矩阵中相同的匹配样本，每个元素表示一个样本中一个基因边对的共表达水平。
[0059]
步骤4：基于步骤3的关联分析结果筛选网络数量性状定位的关联矩阵，导出网络数量性状定位关联的基因网络模块。
[0060]
网络数量性状定位关联是指网络数量性状定位中的边对网络特征。
[0061]
基于步骤3的关联分析结果，筛选出热点snp位点和热点基因对，如某个snp位点与基因边对的交互关系比其他snp位点更显著，或者某个基因边对明显比其他基因边对有更多的snp交互关系，则说明筛选过后的基因型和基因边对的网络数量性状定位关联矩阵将会比全部原始数据(基因序列变异数据矩阵或基因表达矩阵)的关联矩阵密度更高。这样的网络数量性状定位关联的snp位点或者基因边对将用于后续分析。
[0062]
基于筛选出的snp位点与基因边对的关联，获取边对网络特征，例如nqtl模块，这是利用对网络qtl关联矩阵上的层次聚类获得的不同(基因边对组成的)基因网络模块。
[0063]
步骤5：基于导出的数量性状定位关联的基因网络模块进行生物功能和生物特征分析。
[0064]
通过go(gene ontology)，kegg(kyoto encyclopedia of genes and genomes)富集分析(或ipa，ingenuity pathway analysis)，确定每个基因网络模块的生物学功能。通过偏最小二乘判别分析(partial least squares discriminant analysis，pls-da)或典型相关分析(canonical correlation analysis，cca)推断网络特征(如基因边对共表达值)和生理表型(如临床指标)之间的联系。
[0065]
图2示出了本发明的基因网络数量性状定位检测系统的一实施例的原理。请参见图2，本实施例的系统包括：数据集准备模块、矩阵预构造模块、关联分析模块、关联矩阵筛选模块、分析模块。
[0066]
数据集准备模块用于准备网络qtl模型的数据集。
[0067]
数据集准备模块配置为构建多组模拟数据分别用于表达数量性状基因座、网络数量性状基因座的计算，这些数据集包括基因型数据集、分子表型数据集、生理表型数据集和包括协变量在内的可选数据集，以及包括蛋白质互作网络在内的背景生物网络信息；且所有类型的数据都在相同的样本有对应的数值，例如，样本及其标签在所有数据集中是一一匹配的。
[0068]
举例来说，构建8组模拟数据，对于前4组数据，设置基因和snp位点之间存在相关性，即为表达特征。公式如下：
[0069]
pop＝0.2*rnorm(n)
…
(1)
[0070]
snps＝rnorm(n*nss) pop
…
(2)
[0071]
gene＝rnorm(n*ngs) pop snps[,sample(1:nss,ngs)]*1:ngs)/ngs
9/2
[0072]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ…
(3)
[0073]
上述公式中，n为样本数量，ngs为基因数量，nss为snp数量，pop为所有变量中共同
的信号，snps表示模拟生成的基因型数据，gene表示模拟生成的基因表达数据，sample(1:nss,ngs)表示在snp数量范围内随机抽取和基因数量相同的数据，snps[,sample(1:nss,ngs)]则表示随机抽取ngs个基因型数据，rnorm()表示产生一系列的随机数。
[0074]
示例中设定样本数量n为100，基因数量ngs为2000和5000,snp数量nss分别为2000和200000。于是基因个数和snp位点个数有四种数据组合。对于其他四组合成数据，snp和基因数据是随机产生的，将基因对数据与snp，即为边缘性状，按照如下公式联系起来：
[0075]
gene＝rnorm(n*ngs)
…
(4)
[0076]
cordata＝cor pop snps[,sample(1:nss,ncs)]*1:ncs)/ncs
9/2
…
(5)
[0077]
其中cor表示基因边对值，cordata表示与基因型数据相关的基因边对值，ncs表示基因对个数。然后设置样本数量n为100，基因数量ngs和基因对数量cor分别设置为2000、1000*2001和5000、2500*5001，snp位点数量nss设置为2000和20000。同样，有四种数据组合来自基因/基因对和snp。
[0078]
矩阵预构造模块用于基于数据集进行网络qtl模型的基因边对矩阵的预构造。
[0079]
矩阵预构造模块配置为先利用单样本网络计算基因边对值，再利用蛋白质互作网络对计算出的基因边对进行筛选。
[0080]
对于基因边对值的计算，则利用给定原始的基因表达矩阵，群体样本中基因i和j之间的皮尔逊相关系数(pcc)可计算为：
[0081][0082]
上式中，n表示样本数，对于第k个样本，基因i和基因j在该样本上的表达值分别为x
ik
,x
jk
，σ
i
、σ
j
分别表示基因i和基因j的标准差值，e
ijk
表示基因i和基因j在样本k上的单样本基因边对值。对于基因i或者基因j，在所有样本中的表达均值分别为μ
i
，μ
j
。由此值e
ijk
代表样本中一对基因对值，与样本对应的突变数据进一步计算nqtl结果。
[0083]
ppi筛选的具体过程是现有技术，采用了string数据库中人类蛋白质相互作用信息表中的score值，进行筛选处理，其中留下score>900的基因边对。
[0084]
关联分析模块用于对构建的基因边对矩阵进行网络qtl模型的关联分析，确定候选网络数量性状定位。
[0085]
关联分析模块中对网络qtl模型的关联分析是利用现有的matrixeqtl方法对基因序列变异数据矩阵和基因边对矩阵进行关联分析，并得到关联矩阵。关联分析模块所需要的数据由两个数据矩阵组成：一个是基因序列变异数据矩阵，其中每一行代表一个snp，每一列代表一个样本，每一个元素代表一个样本中一个snp的基因型；另一个是基因边对矩阵，其中一行表示一个基因边对，一列表示与基因序列变异数据矩阵中相同的匹配样本，每个元素表示一个样本中一个基因边对的共表达水平。
[0086]
关联矩阵筛选模块用于基于关联分析筛选出网络数量性状定位的关联矩阵，导出网络数量性状定位关联的基因网络模块。
[0087]
网络数量性状定位关联是指网络数量性状定位中的边对网络特征。关联矩阵筛选模块配置为基于关联分析模块的关联分析结果，筛选出热点snp位点和热点基因对：如某一snp位点与基因边对的交互关系比其他snp位点更显著，或者某一基因边对明显比其他基因边对有更多的snp交互关系，则说明筛选过后的基因型和基因边对的网络数量性状定位关
联矩阵将会比全部原始数据(基因序列变异数据矩阵或基因表达矩阵)的关联矩阵密度更高，则筛选出这样的snp位点或者基因边对将用于后续分析；基于筛选出的snp位点与基因边对的关联，获取边对网络特征，导出网络数量性状定位关联的基因网络模块。
[0088]
分析模块用于基于导出的数量性状定位关联的基因网络模块进行生物功能和生物特征分析。
[0089]
分析模块配置为通过go(gene ontology)，kegg(kyoto encyclopedia of genes and genomes)富集分析或ipa(ingenuity pathway analysis)，确定每个基因网络模块的生物学功能；通过偏最小二乘判别分析(partial least squares discriminant analysis，pls-da)或典型相关分析(canonical correlation analysis，cca)推断网络特征和生理表型之间的联系。
[0090]
尽管为使解释简单化将上述方法图示并描述为一系列动作，但是应理解并领会，这些方法不受动作的次序所限，因为根据一个或多个实施例，一些动作可按不同次序发生和/或与来自本文中图示和描述或本文中未图示和描述但本领域技术人员可以理解的其他动作并发地发生。
[0091]
本领域技术人员将进一步领会，结合本文中所公开的实施例来描述的各种解说性逻辑板块、模块、电路、和算法步骤可实现为电子硬件、计算机软件、或这两者的组合。为清楚地解说硬件与软件的这一可互换性，各种解说性组件、框、模块、电路、和步骤在上面是以其功能性的形式作一般化描述的。此类功能性是被实现为硬件还是软件取决于具体应用和施加于整体系统的设计约束。技术人员对于每种特定应用可用不同的方式来实现所描述的功能性，但这样的实现决策不应被解读成导致脱离了本发明的范围。
[0092]
结合本文所公开的实施例描述的各种解说性逻辑板块、模块、和电路可用通用处理器、数字信号处理器(dsp)、专用集成电路(asic)、现场可编程门阵列(fpga)或其它可编程逻辑器件、分立的门或晶体管逻辑、分立的硬件组件、或其设计成执行本文所描述功能的任何组合来实现或执行。通用处理器可以是微处理器，但在替换方案中，该处理器可以是任何常规的处理器、控制器、微控制器、或状态机。处理器还可以被实现为计算设备的组合，例如dsp与微处理器的组合、多个微处理器、与dsp核心协作的一个或多个微处理器、或任何其他此类配置。
[0093]
结合本文中公开的实施例描述的方法或算法的步骤可直接在硬件中、在由处理器执行的软件模块中、或在这两者的组合中体现。软件模块可驻留在ram存储器、闪存、rom存储器、eprom存储器、eeprom存储器、寄存器、硬盘、可移动盘、cd-rom、或本领域中所知的任何其他形式的存储介质中。示例性存储介质耦合到处理器以使得该处理器能从/向该存储介质读取和写入信息。在替换方案中，存储介质可以被整合到处理器。处理器和存储介质可驻留在asic中。asic可驻留在用户终端中。在替换方案中，处理器和存储介质可作为分立组件驻留在用户终端中。
[0094]
在一个或多个示例性实施例中，所描述的功能可在硬件、软件、固件或其任何组合中实现。如果在软件中实现为计算机程序产品，则各功能可以作为一条或更多条指令或代码存储在计算机可读介质上或藉其进行传送。计算机可读介质包括计算机存储介质和通信介质两者，其包括促成计算机程序从一地向另一地转移的任何介质。存储介质可以是能被计算机访问的任何可用介质。作为示例而非限定，这样的计算机可读介质可包括ram、rom、
eeprom、cd-rom或其它光盘存储、磁盘存储或其它磁存储设备、或能被用来携带或存储指令或数据结构形式的合意程序代码且能被计算机访问的任何其它介质。任何连接也被正当地称为计算机可读介质。例如，如果软件是使用同轴电缆、光纤电缆、双绞线、数字订户线(dsl)、或诸如红外、无线电、以及微波之类的无线技术从web网站、服务器、或其它远程源传送而来，则该同轴电缆、光纤电缆、双绞线、dsl、或诸如红外、无线电、以及微波之类的无线技术就被包括在介质的定义之中。如本文中所使用的盘(disk)和碟(disc)包括压缩碟(cd)、激光碟、光碟、数字多用碟(dvd)、软盘和蓝光碟，其中盘(disk)往往以磁的方式再现数据，而碟(disc)用激光以光学方式再现数据。上述的组合也应被包括在计算机可读介质的范围内。
[0095]
提供对本公开的先前描述是为使得本领域任何技术人员皆能够制作或使用本公开。对本公开的各种修改对本领域技术人员来说都将是显而易见的，且本文中所定义的普适原理可被应用到其他变体而不会脱离本公开的精神或范围。由此，本公开并非旨在被限定于本文中所描述的示例和设计，而是应被授予与本文中所公开的原理和新颖性特征相一致的最广范围。