分子遗传学与基因治疗德赢娱乐国际

全文

研究文章
基于gaw18数据集卡尔曼滤波模型的动态关联映射

Burak Karacaoren

土耳其安塔利亚Akdeniz大学农学院动物科学系

*通讯作者:Burak Karacaören, Akdeniz大学农学院动物科学系,土耳其安塔利亚,E-mail: burakkaracaoren@akdeniz.edu.tr

摘要
背景

在横断面基因组学研究中,遗传关系的检测和校正常采用带有系谱信息的线性混合模型。本研究的主要目的是利用随机游走-卡尔曼滤波方法对GAW18数据集进行动态关联映射分析。我们还扩展了模型,将动态基因和环境影响的自回归结构纳入平稳过程。

方法

我们使用随机游走模型,如下图所示

{a_{t + 1}} = {a_t} + {eta _t},{eta _t} proto {mkern 1mu} N左({O,A\sigma _n^2}右)$$ In(1)第一个方程称为观测方程,第二个方程称为状态方程。我们假设观察结果,yt,取决于不可观测的数量,αt,我们的目的是做统计推断αt(州)

结果

用基因组关系矩阵预测舒张压的误差、遗传和永久性环境方差分量分别为17.3(0.0006)、10.9(0.0007)和8.2(0.0006),用家系关系矩阵预测舒张压的误差分量为18.0(0.0007)、9.2(0.0008)和8.4(0.0006)。时间点1,2,3的11711953卢比被发现与MAP4基因相关。

结论

可能由于时间点少,模型没有检测到所有真实的基因组信号。基因组关系矩阵给出了较好的膨胀因子。随机游动是一个非平稳过程,本文通过调整∆参数对平稳情况下的随机游动模型进行了扩展。在基因组研究中,如果不考虑纵向基因和环境的影响,随着时间的推移,可能会导致无法检测到真实信号和/或也可能由于随机误差导致假阳性。

关键字

协会的映射;随机漫步;卡尔曼滤波器;吉布斯抽样

背景

在横断面基因组学研究中,结合家系信息的混合模型方法常用来检测和校正遗传关系。由于基因表达可能随时间变化,重复测量对于检测相关基因组信号[2]更有用。动态关联研究通常采用功能映射方法。然而,从生物学上解释(非)参数函数的回归系数的结果可能是困难的。动态关联映射[3]采用随机回归系数模型。然而,随机回归模型在获得轨迹[4]开始和结束时的准确估计方面有局限性。

此外,上述两种模型(以及动态关联映射文献中的其他模型)都需要一系列的观察才能进行预测:这也可能会造成等待数月(如果不是数年)才能获得预测的问题。在这项研究中,我们假设基因组信号随时间的变化可以通过状态空间形式的随机漫步-卡尔曼滤波模型来获得纵向残差。由于卡尔曼滤波,我们不必等待收集整个数据集来进行模型评估,因此,一旦测量进行,就可以得到估计。由于在关联映射中使用了纵向残差,生物推理也可以很容易地推导出来,只要信号是真实的。

最近,我们在贝叶斯语境[5]中扩展了[1]的语法模型。本文采用[6]模型进行动态关联映射,采用随机游走-卡尔曼滤波方法对GAW18数据集进行动态关联映射分析。我们还扩展了模型,将动态基因和环境影响的自回归结构纳入平稳过程。

方法

GAW 18提供了3个时间点849个个体200个重复的模拟表型。我们使用舒张压(DBP)表型进行关联映射。在第一次复制的3个时间点,利用相关的849个个体分析了第3染色体的65519个SNPs。

质量控制

我们使用来自3号染色体的849个家系个体,共65519个snp进行关联作图。我们排除了729,208个SNPs(由于小等位基因频率<1%),208个SNPs(由于Hardy Weinberg检验(p<0.001), 2个SNPs(由于缺失检验(p>0.1)),在分析[7]中留下58080个SNPs。我们排除了44个基因分型过低的个体,在数据集中留下805个个体。Kolmogrow-Smirnow检验用于评估反应变量的正态性。时间、性别、吸烟状况、年龄和家系数作为固定效应被纳入后续分析,这些分析基于预测和观察之间的相关性进行初步分析。

随机游走模型

我们使用随机游走模型,如下图所示

$ $ {y_t} ={\α_t} + {\ varepsilon _t}, {\ varepsilon _t} \ propto N(0 \σ_e ^ 2) $ $

$ ${\α_ {t + 1}} ={\α_t} +{\埃塔_t},{\埃塔_t} \ propto N(0 \σ_n”^ 2 )............( 1) $ $

(1)第一个方程称为观测方程,第二个方程称为状态方程。我们$ $ \ eqalign{&假定\ \ {\ rm {}}, {\ rm{}} \观察,{y_t}, \, rm{}}{\ \,取决于\ {\ rm {}}, {\ rm{}},难以察觉的\ \ {\ rm{}}数量,{\α_t}, {\ rm {}} \,, \ \ cr & \ {\ rm {}}, {\ rm {}} {\ rm{}} \目标,是\ \ {\ rm {}}, {\ rm{}} \,统计\ {\ rm {}}, {\ rm{}},推理\ \ {\ rm{}}, \,α_t}{\ \({州}\右)。\,我们\ {\ rm {}}, {\ rm {}} \, {\ rm{}}常数\ \ cr & {\ rm{}}方差\ \ {\ rm {}}, {\ varepsilon _t} \,{\埃塔_t} \ \,, \ \σ_e ^ 2和\ \σ_n”^ 2 \ \,分别与\ {\ rm {}}, {\ rm{}}独立\ cr} $ $相同和正态分布随机变量与零的意思。我们假设既有基因影响,也有永久的环境影响。

性状遗传分析采用混合模型;

$ $ y = X \β+ {Z_a} + {Z_p} p + e ................( 2) $ $

其中y为观测向量,β为固定效应向量,a为随机效应向量,p为随机永久环境效应向量,X, Z一个, Zp为设计矩阵,e为随机$$\eqalign{& residual\,{\rm{}}效应的向量。\ \,σ_a ^ 2 \ \,, \σ_p ^ 2 \ {\ rm {}} \, {\ rm{}}, \ \σ_e ^ 2 \; \, \ {\ rm {}}, {\ rm{}}遗传、\,{\ rm{}}永久\环境\ {\ rm {}}, {\ rm {}} \, {\ rm{}}错误\ \ cr & {\ rm{}}差异。\,A,是个体的\,{rm{}},{rm{}}加性\cr} $$遗传关系矩阵;是单位矩阵。A是利用基因型和家系的共祖矩阵系数得到的。

下面,我们展示了基于贝叶斯原理的KF-RW方法的一般假设。基于以下递归关系[8],使用(2)的(3)中不含常数项的联合后验比例分布;

$ $ p \左({{\θ_t} \左|{{\θ_{\离开({- t} \右)}},{最大}}\。} \) \ propto p \左({{\θ_t} \左|{{\θ_{\离开({t - 1} \右)}},{Y_ {t - 1}}} \。左({}\右)p \{\θ_ {t + 1}} \左|{{\θ_t}, {Y_ {t - 1}}} \。左({}\右)p \ {{\ rm {Y}} _t} \左|{{\θ_t}, {Y_ {t - 1}}} \。} \右),$ $

$ $ f \离开({y / b、p \σ_a ^ 2 \σ_p ^ 2 \σ_e ^ 2} \) \ propto{\离开({\σ_e ^ 2} \右)^ {- {1 \ / 2}N}} {\ mkern 1μ}\ exp \离开({-{1 \ / 2}{{\离开({y - Xb - {Z_a} - {Z_p} p} \右)}^ \ '}\离开({y - Xb - {Z_a} - {Z_p} p} \右)/ \σ_e ^ 2} \右)$ $

$ $ {\ rm {X}} \三角洲{\离开({\σ_a ^ 2} \右)^ -}^ {{1 \ / 2}{N_a}} \ exp \离开({{1 \ / 2}{{{rm \{‘}}}_1}{{\ rm{一}}^ {- 1}}{a_1} / \σ_a ^ 2}左\)\[{\刺激\ limits_ {t = 2} ^ t{{{\离开({\σ_a ^ 2} \右)}^ {- {1 \ / 2}{N_a}}}} \ exp \离开({{1 \ / 2}{{\离开({{{rm \{一}}_ {\ rm {t}}}, {{\ rm{一}}_ {{\ rm {t}} - 1}}} \右)}^ \ '} {{\ rm{一}}^{- 1}}\离开({{{rm \{一}}_ {\ rm {t}}}, {{\ rm{一}}_ {{\ rm {t}} - 1}}} \右)/ \σ_a ^ 2} \右)}\]$ $

$ $ {\ rm {X}} \三角洲{\离开({\σ_p ^ 2} \右)^ {- {1 \ / 2}{N_p}}} \ exp \离开({{1 \ / 2}{{{rm \ {p '}}} _1} {{rm \ p{}} _1} / \σ_e ^ 2}左\)\[{\刺激\ limits_ {t = 2} ^ t{{{\离开({\σ_p ^ 2} \右)}^ {- {1 \ / 2}{N_p}}}} \ exp \离开({{1 \ / 2}{{\离开({{{rm \ p {}} _ {\ rm {t}}}, {{rm \ p {}} _ {{\ rm {t}} - 1}}} \右)}^ \ '}{{\ rm{我}}^ {-1}} \离开({{{rm \ p {}} _ {\ rm {t}}}, {{rm \ p {}} _ {{\ rm {t}} - 1}}} \右)/ \σ_p ^ 2} \右)}\]$ $

$ $ {\ rm {X}}{\离开({\σ_e ^ 2} \右)^{- \离开({{{{\ν_e}} \ / 2} + 1} \右)}}\ exp \离开了 [ { - {{{\ ν_e} {S_e}} \ /{2 \σ_e ^ 2}}} \右]{\离开({\σ_a ^ 2} \右)^ -}^{\离开({{{{\ν_a}} \ / 2} + 1} \右)}\ exp \离开了 [ { - {{{\ ν_a} {S_a}} \ /{2 \σ_a ^ 2}}} \]{\左右({\σ_p ^ 2} \右)^{\离开({{{\νp} \ / 2}+ 1} \right)}}\exp \left[ { - {{{\nu _p}{S_p}} \over {2\sigma _p^2}}} \right]...................(3)$$

(3)的最后一行是假设先验方差参数的比例倒卡方分布密度的乘积。假设随机行走模型的∆为1。经过代数处理,条件分布可以写成:

$ $ {b_t} \左|{\σ_a ^ 2 \σ_p ^ 2 \σ_e ^ 2, {a_t}, {p_t}, {y_t} \ sim N \离开({{{\离开({X 'X} \右)}^ {- 1}}X ' \离开({y - {Z_a} - {Z_p} {p_t}} \右),{{\离开({X 'X} \右)}^{- 1}}\σ_e ^ 2} \右)}\ $ $

$ $ {a_t} \左|{\σ_a ^ 2 \σ_p ^ 2 \σ_e ^ 2, {b_t}, {p_t}, {y_t}} \。\ sim $ $

$ $ N \左({{{\左({{1 \ /{\σ_e ^ 2}} {{{rm \ {Z '}}} _ {{a_t}}} {{\ rm {Z}} _ {{a_t}}} +{2 \ /{\σ_a ^ 2}} {{\ rm{一}}^{- 1}}}\右)}^{- 1}}\离开({{1 \ /{\σ_e ^ 2}} {\ rm {Z '}} \离开({{{rm \ {y}} _t}, {{\ rm {X}} _ {\ rm{一}}}{{\ rm {b}} _t}, {{\ rm {Z}} _ {{{rm \ p {}} _ {\ rm {t}}}}} {{rm \ p {}} _ {\ rm {t}}}} \右)+{1 \ /{\σ_a ^ 2}} {{\ rm{一}}^{- 1}}{现代{t+ 1}}} \右),{{\离开({{1 \ /{\σ_e ^ 2}} {{{rm \ {Z '}}} _ {{{rm \{一}}_ {\ rm {t}}}}} {{\ rm {Z}} _ {{{rm \{一}}_ {\ rm {t}}}}} +{2 \ /{\σ_a ^ 2}} {{\ rm{一}}^{- 1}}}\右)}^{- 1}}}\右)$ $

$ $ {p_t} \左|{\σ_a ^ 2 \σ_p ^ 2 \σ_e ^ 2, {a_t}, {b_t}, {y_t}} \。\ sim $ $

$ $ N \左({{{\左({{1 \ /{\σ_e ^ 2}} {{{rm \ {Z '}}} _ {{p_t}}} {{\ rm {Z}} _ {{p_t}}} +{2 \ /{\σ_p ^ 2}} {{\ rm{我}}^{- 1}}}\右)}^{- 1}}\离开({{1 \ /{\σ_e ^ 2}} {\ rm {Z '}} \离开({{{rm \ {y}} _t}, {{\ rm {X}} _ {\ rm{一}}}{{\ rm {b}} _t}, {{\ rm {Z}} _ {{{rm \ p {}} _ {\ rm {t}}}}} {{\ rm{一}}_ {\ rm {t}}}} \右)+{1 \ /{\σ_p ^ 2}} {{\ rm{我}}^{- 1}}{现代{t+ 1}}} \右),{{\离开({{1 \ /{\σ_e ^ 2}} {{{rm \ {Z '}}} _ {{{rm \ p {}} _ {\ rm {t}}}}} {{\ rm {Z}} _ {{{rm \ p {}} _ {\ rm {t}}}}} +{2 \ /{\σ_p ^ 2}} {{\ rm{我}}^{- 1}}}\右)}^{- 1}}}\右)$ $

$ $ \σ_a ^ 2 \左|{\σ_p ^ 2 \σ_e ^ 2, {p_t}, {a_t},} \。{b_t}, {y_t} \ sim{{\离开({{Q_a} +{\ν_a} {S_a}} \右)}\ /{\气_ {DF} ^ 2}} $ $

$ $ \σ_p ^ 2 \左|{\σ_a ^ 2 \σ_e ^ 2, {p_t}, {a_t},} \。{b_t}, {y_t} \ sim{{\离开({{Q_p} +{\ν_p} {S_p}} \右)}\ /{\气_ {DF} ^ 2}} $ $

$ $ \σ_e ^ 2 \左|{\σ_a ^ 2 \σ_p ^ 2, {p_t}, {a_t},} \。{b_t}, {y_t} \ sim{{\离开({{Q_e} +{\ν_e} {S_e}} \右)}\ /{\气_ {DF} ^ 2}} $ $

其中最后两行Qa, Qp和Qe分别表示误差项和DF自由度的二次形式。我们使用DBP的5000次迭代老化周期运行该模型,迭代次数为10,000次。为了减少自相关性,我们每10次迭代采样一次。我们尝试了不同参数的逆Wishart先验分布来获得残差。

我们使用混合模型进行全基因组关联分析[9,7],使用R软件[10]:

$ $ {\ rm {y =}} \, {\ rm {Xb}} \, {\ rm { + }}\,{\ rm {e }}..........{\ rm {(4)}} $ $

其中y包含(3)的残差或随机效应,b指定固定效应(SNP), X和是关联$$矩阵,而\,e\,是\,一个\,向量,包含\,残差\,\,和,假设\,正态\,分布\,带有\,I\,\sigma _e^2。I是一个单位矩阵,∑_e^2$$是剩余方差。在关联作图中,我们使用了5%的错误发现阈值来检测基因组信号。为了便于比较,我们还对每个时间点使用了横截面GRAMMAR[1]方法。利用基因组共祖矩阵[9]和系谱信息估计DBP的遗传率分别为0.299和0.259。

结果

在不了解底层仿真模型的情况下进行分析。然而,我们在讨论结果时使用了GAW18的答案。

我们使用Kolmogrow-Smirnov检验确认了正态性。然而,由于我们采用了贝叶斯残差:所有响应变量转换为正态分布(P > 0.01)。时间、性别、吸烟状况、年龄和家系数作为固定效应被纳入后续分析,这些分析基于预测和观察之间的相关性进行初步分析。我们发现预测和观测之间的相关性在0.15时最高。利用基因组关系矩阵预测DBP的误差、遗传和永久环境方差分别为17.3(0.0006)、10.9(0.0007)和8.2(0.0006),用系谱关系矩阵预测DBP的误差、遗传和永久环境方差分别为18.0(0.0007)、9.2(0.0008)和8.4(0.0006)。DBP模拟遗传率为0.317,而基因组亲缘关系估计更接近其真实值(表1-3)。

讨论

关于基因进化和长期环境影响的假设可能很重要。一定程度的自回归结构可能比随机游走模型更现实。我们简单地根据(3)中参数空间的限制来调整模型。使用∆=0.1和∆=0.9,我们考虑了偏离随机行走的两种极端情况。随机漫步假设基因效应可以随着时间的推移在上下两个方向上缓慢变化,∆=1.0,通过调整∆,可以引入基因和环境效应的自回归结构,得到平稳分布。在这里,时间序列将沿着平均轨迹分布。

使用∆=0.1和∆=0.9预测DBP的误差、遗传和永久环境方差分量分别为10.26(0.0006)、1.98(0.0001)和1.98(0.0003)。使用∆=1.0、∆=0.1和∆=0.9预测遗传度分别为0.299、0.139和0.246。与自回归结构相比,随机游走的结果更好(DBP模拟结果为0.317)。我们假设:随着时间的推移,随着相关性和子结构信息的积累,增加时间点应该会降低基因组膨胀因子[7]。我们在混合模型中使用了基于基因组和系谱的关系矩阵(3)。基因组关系矩阵显示,与基于系谱的关系矩阵1.57、1.83和1.59相比,在三个时间点上,DBP的基因组膨胀因子为1.40、1.33和1.63。由于时间点少(t=3),我们仍然获得了高水平的基因组膨胀因子(λ > 1)。从表1和表2可以看出,基因组关系和系谱关系矩阵在不同时间点检测到的snp多为不同集。

但是,由于样本量小,时间点少,可能会导致假阳性和假阴性。特别是第一次可以如此观点:基因关系矩阵检测到154个snp错误发现率为5%(罗斯福)(134和56个snp检测时间点2和3分别在5%罗斯福)和血统关系矩阵检测到216个snp在5%罗斯福(300和96个snp检测时间点2和3分别在5%罗斯福)。由于较小的基因组膨胀因子,我们研究了因果单核苷酸多态性的基因组关系矩阵结果。在时间点1,2,3上发现rs11711953与MAP4基因相关。

我们使用GRAMMAR方法对每个时间点(以及它们的平均值)进行横断面分析(表3)。然而,在多次假设修正后,我们没有检测到任何基因组信号。虽然ARHGEF3附近存在rs1948722从时间点2发出的信号(p<0.00012),但FDR多次假设校正后SNP变得不显著。与随机漫步模型(表1,2)的p值相比,GRAMMAR p值的大小(表3)更大。这清楚地表明,需要通过适当的方法来考虑基因和环境随时间的纵向影响。否则,由于基因组信号将被随机错误污染,这可能导致信号未被检测或也可能导致假阳性。随机游动是一个非平稳过程,本文通过调整∆参数对平稳情况下的随机游动模型进行了扩展。然而,如果非平稳假设对基因的动力学和永久性环境影响有用或没用,则需要理论和实证的动态关联研究。

结论

与系谱信息相比,基因组关系矩阵提供了更好的膨胀因子和遗传力估计。基于卡尔曼滤波器递归结构的随机游动模型在实际应用中具有一定的应用价值。当纵向观测可用时(例如每日或每月),模型可以根据卡尔曼滤波顺序预测在线基因组信号。在基因组研究中,如果不考虑纵向基因和环境的影响,随着时间的推移,可能会导致无法检测到真实信号和/或也可能由于随机误差导致假阳性。

确认

基因分析研讨会得到了美国国立卫生研究院(NIH) R01 GM031575国家普通医学科学研究所的资助。本研究得到了Akdeniz大学106项目研究基金的资助。作者希望感谢与Luc Janss博士关于自回归结构的有益讨论。

表1:利用基因组关系矩阵随机游走模型获得DBP首次复制的前10个单核苷酸多态性及其对应的原始p值

表2:利用系谱关系矩阵随机游走模型获得DBP首次复制的前10个单核苷酸多态性及其对应的原始p值

表3:利用GRAMMAR对DBP首次复制得到的前10个SNPs及其对应的原始p值

相互竞争的利益

没有相互竞争的利益。

参考文献
  1. auulchenko SY, de Koning DJ, Haley C(2007)基于混合模型和回归的全基因组快速关联:基于全基因组系谱的数量性状位点关联分析的一种快速和简单的方法。遗传学177:577 - 585。[Ref。
  2. [吴锐,林敏(2006)功能作图——如何作图和研究动态复杂性状的遗传结构。]Nat Rev Genet 7: 229- 237。[Ref。
  3. Macgregor S, Knott SA, White I, Visscher PM(2005)复杂系谱中纵向数量性状的数量位点分析。遗传学171:1365 - 1376。[Ref。
  4. Faro EL, Cardoso VL, Albuquerque LG(2008)应用随机回归模型对Caracu小母牛(牛科偶蹄目)测试日产奶量的方差分量估计。Genet Mol Biol 31。[Ref。
  5. Karacaören B(2014)使用祖先信息标记对F2小鼠数据集生长相关性状的混合定位。J Bioinform Comput Biol 12。[Ref。
  6. Karacaören B(2015)一种映射动态数量性状的贝叶斯随机游走方法。J应用非线性动力学。
  7. Purcell S, Neale B, tod - brown K, Thomas L, Ferreira MA, et al. (2007) PLINK:全基因组关联和基于群体连锁分析的工具集。Am J Hum Genet 81: 559-575。[Ref。
  8. 贝叶斯预测与动态模型。施普林格。[Ref。
  9. Endelman JB(2011)利用R包rrBLUP进行基因组选择的Ridge回归和其他核。植物Gen 4: 250-255。[Ref。
  10. R Development Core Team(2014)统计计算的语言和环境。统计计算基础,维也纳,奥地利。[Ref。

在此下载临时PDF

PDF

条信息

Aritcle类型:研究文章

引用:Karacaören B(2015)基于GAW 18数据集Kalman滤波模型的动态关联映射。Int J Mol Genet Gene Ther 1 (1): http://dx.doi.org/10.16966/2471-4968.101

版权:©2015 Karacaören b。这是一篇开放获取的文章,在知识共享署名许可协议的条款下发布,该协议允许在任何媒体上无限制地使用、发布和复制,前提是注明原作者和来源。

出版的历史:

  • 收到日期:2015年5月27日

  • 接受日期:2015年7月6日

  • 发表日期:2015年7月10