CN103699815B - 一种同源四倍体自然群体的连锁不平衡分析模型的构建方法 - Google Patents

Abstract

本发明提供一种同源四倍体自然群体的连锁不平衡分析模型的构建方法，包括：对来自自然群体中一组同源四倍体植物，设置一个位点A和一个位点B；两个位点产生16单倍型，且分别都具有35种基因型,但在遗传上每个位点能够识别的四倍体基因型有只有15种，因此两个位点共产生225种基因型；从所述自然群体中随机抽取n个同源四倍体个体，得到不同基因型个体数量；根据单体型频率计算相应的配子频率，并根据配子频率计算计算相应的基因型频率；计算基因型频率；计算单倍型基因频率；通过单倍型基因频率获取两个位点的等位基因的估算频率；通过单倍型基因频率计算每两个位点间的连锁不平衡系数；然后进行验证。

Description

一种同源四倍体自然群体的连锁不平衡分析模型的构建方法

技术领域

本发明涉及生物建模技术领域，特别是指一种同源四倍体自然群体的连锁不平衡分析模型的构建方法。

背景技术

半数以上的被子植物是多倍体，或在其进化史上出现过多倍体现象，故此多倍体对于植物的遗传与进化研究具有重要的价值。同源四倍体作为自然界中常见的多倍体，是研究多倍体时常用的材料，分析多位点群体遗传研究的重要性已经越来越被认可。连锁不平衡分析是分析自然群体遗传结构、分化、多样性的重要工具。连锁不平衡分析是基于多个位点，可以通过基因组估算程度和非随机组合的分布提供不同物种结构的额外信息。对于一个随机组合群体，两标记间连锁不平衡世代的比例的衰减取决于两标记间的重组率。连锁不平衡的衰减速率为物种复杂性状的高通量作图建立了一个基本原则，因此通过比较LD遗传距离的衰减率，可以推断一个物种的进化历史，但是这一方面目前仅适用于二倍体自然群体，至今仍没有一种同源四倍体两位点的连锁不平衡分析模型。

发明内容

本发明要解决的技术问题是提供一种能够构建同源四倍体进行两位点连锁不平衡分析模型的方法。

为了解决上述问题，本发明实施例提出了一种同源四倍体自然群体的连锁不平衡分析模型的构建方法，包括：

步骤1、对于来自自然群体中一组同源四倍体植物，设置一个位点A和一个位点B；

其中该位点A具有4个等位基因记为A₁，A₂，A₃，A₄；该位点B具有4个等位基因记为B₁，B₂，B₃，B₄；其中p₁，p₂，p₃，p₄为位点A的基因频率，且p₁+p₂+p₃+p₄=1；其中q₁，q₂，q₃，q₄为位点B的基因频率，且q₁+q₂+q₃+q₄=1；

两个位点产生16单倍型表示为A₁B₁，A₁B₂，A₁B₃，A₁B₄，A₂B₁，A₂B₂，A₂B₃，A₂B₄，A₃B₁，A₃B₂，A₃B₃，A₃B₄，A₄B₁，A₄B₂，A₄B₃，A₄B₄；基因连锁不平衡系数表示为D_ij，其代表A_i和B_j的连锁不平衡系数（i=1，2，3，4；j=1，2，3，4，下同）；单倍型A_iB_j频率表示为p_ij；

步骤2、从所述自然群体中随机抽取n个同源四倍体个体，得到不同基因型个体数量n_u/v；

步骤3、根据单体型频率计算相应的配子频率，并根据配子频率计算相应的基因型频率为p_u/v；

如基因型为A₃/B₁₂的同源四倍体，其包含的基因型有A₃A₃A₃A₃B₁B₁B₁B₂，A₃A₃A₃A₃B₁B₁B₂B₂，A₃A₃A₃A₃B₁B₂B₂B₂，其基因型频率p_3/12＝2p³ ₃₁+3p₃₁ ²p₃₂ ²+2p₃₁p₃₂ ³；

步骤4、

估算参数，未知单倍型频率Ω=（p_ij），根据如下似然函数值计算：logL(Ω)＝Σn_ijlog(p_ij) （1）

由公式（1）获得单倍型频率的最大似然估计

让n_u/v表示含有单倍型A_iB_j的基因型的个体数量，p_u/v表示含有单倍型A_iB_j的基因型频率，p_ij表示p_u/v中含有需要估算的单倍型A_iB_j的基因型频率；基因型中含有需要估算的单倍型并且只含有此单一的单倍型的同源四倍体的数量由m_u/v表示；

需要估算的单倍型为A_iB_j的个体数量为

其中

其中k是与u，v相关的自然数；

根据公式（2）和公式（4）组成EM算法的一个循环；将单倍型频率的初始值来计算在E步骤的公式（4）中φ_k的估算值；然后使用所计算出的φ_k的估算值带入M步骤的公式（2）来估算单倍型频率；重复进行E和M步骤，直至单体型频率估计收敛，得到单倍型基因频率的估算值；

根据以下公式，通过单倍型基因频率获取两个位点的等位基因的估算频率；

根据以下公式，通过单倍型基因频率计算每两个位点间的连锁不平衡系数；

;

步骤5、根据以下公式验证估计分布与实际分布是否存在显著性差异；

H₀：D₁₁=D₁₂=D₁₃=D₂₁=D₂₂=D₂₃=D₃₁=D₃₂=D₃₃=0；

H₁：Not H₀ （5）

其中H₀为同源四倍体各位点不存在连锁不平衡系数；

通过最大似然估计，分别估算和波浪线对应的零假设的最大似然估计，抑扬符号对应的备择假设的最大似然估计；

通过以下公式计算对数似然比的统计学检验

LR是遵循χ²分布，自由度是9；

将计算得到的LR进行卡方检验，判断得到的LR值经卡方检验是否呈显著差异，如果呈显著差异则最终结果为H₀，如果不呈显著差异则最终结果为H₁。

其中，所述方法还包括：

步骤6、将9个连锁不平衡系数分别进行检验。

其中，所述16单倍型的频率包含以下等位基因频率以及九个独立的连锁不平衡系数：

p₁₁ p₁q₁+D₁₁

p₁₂ p₁q₂+D₁₂

p₁₃ p₁q₃+D₁₃

p₁₄ p₁q₄-D₁₁-D₁₂-D₁₃

p₂₁ p₂q₁+D₂₁

p₂₂ p₂q₂+D₂₂

p₂₃ p₂q₃+D₂₃

p₂₄ p₂q₄-D₂₁-D₂₂-D₂₃

p₃₁ p₃q₁+D₃₁

p₃₂ p₃q₂+D₃₂

p₃₃ p₃q₃+D₃₃

p₃₄ p₃q₄-D₃₁-D₃₂-D₃₃

p₄₁ p₄q₁-D₁₁-D₂₁-D₃₁

p₄₂ p₄q₂-D₁₂-D₂₂-D₃₂

p₄₃ p₄q₃-D₁₃-D₂₃-D₃₃

p₄₄ p₄q₄+D₁₁+D₁₂+D₁₃+D₂₁+D₂₂+D₂₃+D₃₁+D₃₂+D₃₃

其中，所述位点A和位点B分别都具有35种基因型，在遗传上能够识别基因型共225种；

位点A有，4种A_iA_iA_iA_i，其中i=1，2，3，4；12种其中i₁≠i₂=1，2，3，4；6种其中i₁≠i₂=1，2，3，4；12种其中i₁≠i₂≠i₃=1，2，3，4；1种A₁A₂A₃A₄；

位点B有，4种B_jB_jB_jB_j，其中j=1，2，3，4；12种B_j1B_j1B_j1B_j2，其中j₁≠j₂=1，2，3，4；6种B_j1B_j1B_j2B_j2，其中j₁≠j₂=1，2，3，4；12种B_j1B_j1B_j2B_j3，其中j₁≠j₂≠j₃=1，2，3，4；1种B₁B₂B₃B₄；

两个位点，在遗传上能够识别的四倍体基因型共225种，表示为Au/Bv（u=1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234；v=1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234）。

本发明的上述技术方案的有益效果如下：

本发明实施例提供了一种同源四倍体的连锁不平衡分析方法，提出了一种根据同源四倍体减数分裂的特点推导出计算分子标记之间连锁不平衡的统计模型，并提供了检测各种连锁不平衡系数的方法；通过计算机模拟试验验证了本发明实施例的方法的有效性，指出了模型的适用条件。利用本研究发展的模型，我们可以对同源四倍体进行连锁不平衡检验，观测其连锁不平衡系数，了解其群体的结构。

附图说明

图1为本发明实施例的方法的流程示意图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明实施例提出了一种同源四倍体自然群体的连锁不平衡分析模型的构建方法，其流程如图1所示的，包括：

步骤1、对于本来自自然群体中一组同源四倍体植物，设置有一个位点A和一个位点B；

其中，该位点A具有4个等位基因，记为A₁，A₂，A₃，A₄；该位点B具有4个等位基因，记为B₁，B₂，B₃，B₄；其中p₁，p₂，p₃，p₄（p₁+p₂+p₃+p₄=1）为位点A的基因频率和q₁，q₂，q₃，q₄（q₁+q₂+q₃+q₄=1）为位点B的基因频率；

两个位点产生16个单倍型（不同位点的等位基因在同一染色体上的线性组合）表示为A₁B₁，A₁B₂，A₁B₃，A₁B₄，A₂B₁，A₂B₂，A₂B₃，A₂B₄，A₃B₁，A₃B₂，A₃B₃，A₃B₄，A₄B₁，A₄B₂，A₄B₃，A₄B₄；基因连锁不平衡系数表示为D_ij（i=1，2，3，4；j=1，2，3，4），其代表A_i和B_j的连锁不平衡系数，其中独立的连锁不平衡系数有9个；单倍型A_iB_j频率表示为p_ij；

其中，所述16个单倍型的频率包含相应的等位基因频率以及九个独立的连锁不平衡系数，其组成见表1；

表1二倍体配子的单倍型的基因频率的组成

其中，所述两个位点分别都具有35种基因型；

所述位点A有，4种A_iA_iA_iA_i（i=1，2，3，4），12种（i₁≠i₂=1，2，3，4），6种（i₁≠i₂=1，2，3，4），12种（i₁≠i₂≠i₃=1，2，3，4）和1种A₁A₂A₃A₄；

所述位点B有，4种B_jB_jB_jB_j，其中j=1，2，3，4；12种B_j1B_j1B_j1B_j2，其中j₁≠j₂=1，2，3，4；6种B_j1B_j1B_j2B_j2，其中j₁≠j₂=1，2，3，4；12种B_j1B_j1B_j2B_j3，其中j₁≠j₂≠j₃=1，2，3，4；1种B₁B₂B₃B₄；

由于单个位点，在遗传上能够识别的四倍体基因型只有15种，两个位点共产生15*15=225种（见下表A、表B），因此，两位点四倍体基因型只有225种，表示为A_u/B_v（u=1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234；v=1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234）。

表A：A位点的15种基因型

表B：B位点的15种基因型。

步骤2、假设从自然群体中随机抽取n个同源四倍体个体，用n_u/v表示不同基因型个体数量。

步骤3、基因型频率是由相应的配子频率计算而来，配子频率是由相应的单体型频率计算而来，其相应的基因型频率表达为p_u/v；如基因型为A₃/B₁₂的同源四倍体，其包含的基因型有A₃A₃A₃A₃B₁B₁B₁B₂，A₃A₃A₃A₃B₁B₁B₂B₂，A₃A₃A₃A₃B₁B₂B₂B₂，其基因型频率p_3/12＝2p³ ₃₁+3p₃₁ ²p₃₂ ²+2p₃₁p₃₂ ³。

步骤4、估算参数，未知单倍型频率Ω=（p_ij），可根据如下似然函数值算出来：

logL(Ω)＝Σn_ijlog(p_ij) （1）

其中n_ij是单倍型为A_iB_j的个体总数量；从式（1）可得到单倍型频率的最大似然估计，即

如果某种基因型频率不是单一的单倍型频率计算而来，而是多种（如A₃/B₁₂，其频率表示为p_3/12＝2_p ³ ₃₁+3p₃₁ ²p₃₂ ²+2p₃₁p₃₂ ³）；设n_u/v表示含有单倍型A_iB_j的基因型的个体数量，p_u/v表示含有单倍型A_iB_j的基因型频率，p_ij表示p_u/v中含有需要估算的单倍型A_iB_j的基因型频率；基因型中含有需要估算的单倍型并且只含有此单一的单倍型的同源四倍体的数量由m_u/v表示；

需要估算的单倍型为A_iB_j的个体数量为

其中

其中k是与u，v相关的自然数；根据公式（2）和（4）组成EM算法的一个循环；设置单倍型频率的初始值来计算在E步骤（式4）中φ_k的估算值，使用所计算出的φ_k的估算值带入M步骤（式2）来估算单倍型频率；重复进行以上E和M步骤，直至单体型频率估计收敛，得到单倍型频率的估算值；两个位点的等位基因的估算频率可通过单倍型频率得到，见表2:

表2A、B两个位点的等位基因频率的估算的推导公式

每两个位点间的连锁不平衡系数的估算也可通过单倍型频率得到，见表3:

表3A、B两位点间连锁不平衡系数的估算的推导公式

步骤5、假设检验是依据一定的概率原则，以较小的风险来判断估计分布与实际分布是否存在显著性差异，是否应当接受H0原假设选择的一种检验方法。

对于是否存在连锁不平衡系数，可以制定以下假设（H0是指同源四倍体各

位点不存在连锁不平衡系数）：

H₀：D₁₁=D₁₂=D₁₃=D₂₁=D₂₂=D₂₃=D₃₁=D₃₂=D₃₃=0；

H₁：Not H₀ （5）

通过最大似然估计，分别估算和波浪线对应的零假设的最大似然估计，抑扬符号对应的备择假设的最大似然估计；对数似然比的统计学检验可通过以下方程得到

LR是遵循χ²分布，自由度是9；

最终将得到的LR进行卡方检验，如果得到的LR值经卡方检验呈显著差异，我们则接受H₁假设，反之，若不呈显著差异，我们则接受H₀假设。

需要说明的是，虽然H₀与H₁是两个相互对立的假设，但是由于H₀假设所有参数都被限定，仅仅包含一种情况。在参数受到限制的时候，似然函数最大值必然会受到约束影响而降低。而H₁假设则自由度很高，没有任何的限制，因此似然函数最大值并不会受到约束。所以计算出来的L₁是恒定大于等于L₀的，因此得到的LR必定是一个非负数。最终计算出的LR是否大于零，也是检验模型的手段之一。

步骤6、上述进行了完全没有连锁不平衡系数和至少有一个系数的假设检验；在这个检验外，还可以9个连锁不平衡系数分别拿出来进行检验，H₀是指同源四倍体两个位点的其中一个连锁不平衡系数不存在的情况；例如：

H₀:D₁₁=0or D₁₂=0or D₁₃=0or D₂₁=0or D₂₂=0or D₂₃=0or D₃₁=0or D₃₂=0or D₃₃=0;

H₁：Not H₀

LR是遵循χ²分布，自由度是1；

通过上述检验，可以知道在该位点每一个基因型频率是否符合遗传平衡定律，通过这一数据了解群体的遗传结构；

除了上述假设检验，模拟实验还进行了在不同样本下的模型的假阳性分析和模型的功效分析，见表4:

表4模型的假阳性和功效

在假设检验中，估算得到的假设检验的假阳性很低，说明模型在估算出现假阳性的可能性很低，证明模型的可行性和准确性；估算得到的假设检验的功效整体较大，其中D₃₁=0的假设检验的功效明显较低，但并不影响整体的估算结果，连锁不平衡系数的估算是可以应用此模型进行的并可以得到较为准确的结果。

模拟实验

为了验证该同源四倍体两个位点分析模型的准确性与统计性质，使用计算机进行模拟实验；

假设A位点的A₁、A₂、A₃、A₄四种基因的基因频率为p₁、p₂、p₃、p₄分别等于0.0800、0.1000、0.1200、0.7000；同理，B位点的B₁、B₂、B₃、B₄四种基因的基因频率为q₁、q₂、q₃、q₄分别等于0.0900、0.0800、0.1300、0.7000。再假设所需估算的9个连锁不平衡系数D₁₁、D₁₂、D₁₃、D₂₁、D₂₂、D₂₃、D₃₁、D₃₂和D₃₃分别等于0.0100、0.0200、0.0100、0.0050、0.0100、0.0010、0.0030、0.0050和0.0070；这套参数假设需要先使用参数的估算里面所提到的基因频率与基因型频率的关系式，求出这套参数的粗略范围，再进行参数假设；表5给出的是此参数下1000次计算机模拟实验的参数估算结果与方差。

表51000次计算机模拟实验的参数估算平均值与标准差

根据估算出参数与真实数据的接近程度，可以认定模型本身是正确的；其中n代表的是样本数；表中是在不同样本数量下，估算出的参数值，括号内代表的是标准差；根据1000次模拟实验的平均值，可以看出，随着样本数量的增大，估算出来的参数越接近真实值，同时标准差越来越小，说明估算出的结果越来越准确；当然，如果模型本身没有错，样本的增大必然导致估算结果越来越准确，需要找出的是样本数量与精准度之间的一个均衡，即利用较少的样本估算出相对精准的结果；从上表可以看出，在样本数量为200的情况下，估算的平均值可以满足需求。

下面通过三组例子对本发明实施例的方案进行验证：

实验1：第一组参数下的假设检验结果和模拟实验结果。

假设A位点的A1、A2、A3、A4四种基因的基因频率为p1、p2、p3、p4分别等于0.0800、0.1000、0.1200、0.7000。同理，B位点的B1、B2、B3、B4四种基因的基因频率为q1、q2、q3、q4分别等于0.0900、0.0800、0.1300、0.7000。再假设所需估算的9个连锁不平衡系数D11、D12、D13、D21、D22、D23、D31、D32和D33分别等于0.0100、0.0200、0.0100、0.0050、0.0100、0.0010、0.0030、0.0050和0.0070。

假设检验结果：

所有D=0的假设检验结果

D₁₁=0的假设检验结果

D₁₂=0的假设检验结果

D₁₃=0的假设检验结果

D₂₁=0的假设检验结果

D₂₂=0的假设检验结果

D₂₃=0的假设检验结果

D₃₁=0的假设检验结果

D₃₂=0的假设检验结果

D₃₃=0的假设检验结果

模拟实验结果：

假阳性和功效结果：

实验2；第二组参数下的假设检验结果和模拟实验结果。

假设A位点的A₁、A₂、A₃、A₄四种基因的基因频率为p₁、p₂、p₃、p₄分别等于0.3000、0.2600、0.2000、0.2400。同理，B位点的B₁、B₂、B₃、B₄四种基因的基因频率为q₁、q₂、q₃、q₄分别等于0.2100、0.1800、0.4000、0.2100。再假设所需估算的9个连锁不平衡系数D₁₁、D₁₂、D₁₃、D₂₁、D₂₂、D₂₃、D₃₁、D₃₂和D₃₃分别等于0.0020、0.0300、-0.0100、0.0100、0.0030、0.0100、0.0050、-0.0090和0.0060。

假设检验结果：

所有D=0的假设检验结果

D₁₁=0的假设检验结果

D₁₂=0的假设检验结果

D₁₃=0的假设检验结果

D₂₁=0的假设检验结果

D₂₂=0的假设检验结果

D₂₃=0的假设检验结果

D₃₁=0的假设检验结果

D₃₂=0的假设检验结果

D₃₃=0的假设检验结果

模拟实验结果：

假阳性和功效结果：

实验3：第三组参数下的假设检验结果和模拟实验结果。

假设A位点的A₁、A₂、A₃、A₄四种基因的基因频率为p₁、p₂、p₃、p₄分别等于0.3600、0.4100、0.1200、0.1100。同理，B位点的B₁、B₂、B₃、B₄四种基因的基因频率为q₁、q₂、q₃、q₄分别等于0.0900、0.1400、0.3400、0.4300。再假设所需估算的9个连锁不平衡系数D₁₁、D₁₂、D₁₃、D₂₁、D₂₂、D₂₃、D₃₁、D₃₂和D₃₃分别等于0.0200、-0.0010、0.0040、-0.0070、0.0100、0.0300、-0.0080、-0.0020和-0.0090。

假设检验结果：

所有D=0的假设检验结果

D₁₁=0的假设检验结果

D₁₂=0的假设检验结果

D₁₃=0的假设检验结果

D₂₁=0的假设检验结果

D₂₂=0的假设检验结果

D₂₃=0的假设检验结果

D₃₁=0的假设检验结果

D₃₂=0的假设检验结果

D₃₃=0的假设检验结果

模拟实验结果：

假阳性和功效结果：

本发明实施例根据同源四倍体减数分裂的特点推导出计算分子标记之间连锁不平衡的统计模型，并指出检测各种连锁不平衡系数的计算程序。计算机模拟试验验证了模型有效性，指出了模型的适用条件；利用本研究发展的模型，我们可以对同源四倍体进行连锁不平衡检验，观测其连锁不平衡系数，了解其群体的结构。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种同源四倍体自然群体的连锁不平衡分析模型的构建方法，其特征在于，包括：

步骤1、对于来自自然群体中一组同源四倍体植物，设置一个位点A和一个位点B；

其中该位点A具有4个等位基因记为A₁，A₂，A₃，A₄；该位点B具有4个等位基因记为B₁，B₂，B₃，B₄；其中p₁，p₂，p₃，p₄为位点A的基因频率，且p₁+p₂+p₃+p₄＝1；其中q₁，q₂，q₃，q₄为位点B的基因频率，且q₁+q₂+q₃+q₄＝1；

两个位点产生16单倍型表示为A₁B₁，A₁B₂，A₁B₃，A₁B₄，A₂B₁，A₂B₂，A₂B₃，A₂B₄，A₃B₁，A₃B₂，A₃B₃，A₃B₄，A₄B₁，A₄B₂，A₄B₃，A₄B₄。基因连锁不平衡系数表示为D_ij，其代表A_i和B_j的连锁不平衡系数，i＝1，2，3，4；j＝1，2，3，4，下同；单倍型A_iB_j频率表示为p_ij；

步骤2、从所述自然群体中随机抽取n个同源四倍体个体，得到不同基因型个体数量n_u/v；

步骤3、根据单体型频率计算相应的配子频率，并根据配子频率计算计算相应的基因型频率为p_u/v；

如基因型为A₃/B₁₂的同源四倍体，其包含的基因型有A₃A₃A₃A₃B₁B₁B₁B₂，A₃A₃A₃A₃B₁B₁B₂B₂，A₃A₃A₃A₃B₁B₂B₂B₂，其基因型频率p_3/12＝2p³ ₃₁+3p₃₁ ²p₃₂ ²+2p₃₁p₃₂ ³；

步骤4、估算参数，未知单倍型频率Ω＝(p_ij)，根据如下似然函数值计算：

log L(Ω)＝∑n_ijlog(p_ij) (1)

由(1)获得单倍型频率的最大似然估计

${\hat{p}}_{ij}=\frac{1}{4n}\left(n_{ij}\right)---\left(2\right)$

让n_u/v表示含有单倍型A_iB_j的基因型的个体数量，p_u/v表示含有单倍型A_iB_j的基因型频率，p_ij表示p_u/v中含有需要估算的单倍型A_iB_j的基因型频率；基因型中含有需要估算的单倍型并且只含有此单一的单倍型的同源四倍体的数量由m_u/v表示；；

需要估算的单倍型为A_iB_j的个体数量为

$n_{ij}={\mathrm{\Σm}}_{u/v}+\underset{k}{\Σ}{\mathrm{\φ}}_k{n}_{u/v}---\left(3\right)$

其中

${\mathrm{\φ}}_k=\frac{p_{ij}}{p_{u/v}}---\left(4\right)$

其中k是与u，v相关的自然数；

根据公式(2)和公式(4)组成EM算法的一个循环；将单倍型频率的初始值来计算在E步骤的公式(4)中φ_k的估算值；然后使用所计算出的φ_k的估算值带入M步骤的公式(2)来估算单倍型频率；重复进行E和M步骤，直至单体型频率估计收敛，得到单倍型基因频率的估算值；

根据以下公式，通过单倍型基因频率获取两个位点的等位基因的估算频率；

$\begin{array}{cc}{\hat{p}}_1& {\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14}\\ {\hat{p}}_2& {\hat{p}}_{21}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24}\\ {\hat{p}}_3& {\hat{p}}_{31}+{\hat{p}}_{31}+{\hat{p}}_{33}+{\hat{p}}_{34}\\ {\hat{p}}_4& {\hat{p}}_{41}+{\hat{p}}_{42}+{\hat{p}}_{43}+{\hat{p}}_{44}\\ {\hat{q}}_1& {\hat{p}}_{11}+{\hat{p}}_{21}+{\hat{p}}_{31}+{\hat{p}}_{41}\\ {\hat{q}}_2& {\hat{p}}_{12}+{\hat{p}}_{22}+{\hat{p}}_{32}+{\hat{p}}_{42}\\ {\hat{q}}_3& {\hat{p}}_{13}+{\hat{p}}_{23}+{\hat{p}}_{33}+{\hat{p}}_{43}\\ {\hat{q}}_4& {\hat{p}}_{14}+{\hat{p}}_{24}+{\hat{p}}_{34}+{\hat{p}}_{44}\end{array}$

根据以下公式，通过单倍型基因频率计算每两个位点间的连锁不平衡系数；

$\begin{array}{cc}{\hat{D}}_{11}& {\hat{p}}_{11}-({\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14})({\hat{p}}_{11}+{\hat{p}}_{24}+{\hat{p}}_{34}+{\hat{p}}_{41})\\ {\hat{D}}_{12}& {\hat{p}}_{12}-({\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14})({\hat{p}}_{21}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})\\ {\hat{D}}_{13}& {\hat{p}}_{13}-({\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14})({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})\\ {\hat{D}}_{14}& -{\hat{D}}_{11}-{\hat{D}}_{12}-{\hat{D}}_{13}\\ {\hat{D}}_{21}& {\hat{p}}_{21}-({\hat{p}}_{21}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})({\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14})\\ {\hat{D}}_{22}& {\hat{p}}_{22}-({\hat{p}}_{21}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})({\hat{p}}_{24}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})\\ {\hat{D}}_{23}& {\hat{p}}_{23}-({\hat{p}}_{21}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})\end{array}$

$\begin{array}{cc}{\hat{D}}_{24}& -{\hat{D}}_{21}-{\hat{D}}_{22}-{\hat{D}}_{23}\\ {\hat{D}}_{31}& {\hat{p}}_{31}-({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})({\hat{p}}_{11}+{\hat{p}}_{12}+{\hat{p}}_{13}+{\hat{p}}_{14})\\ {\hat{D}}_{32}& {\hat{p}}_{32}-({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})({\hat{p}}_{24}+{\hat{p}}_{22}+{\hat{p}}_{23}+{\hat{p}}_{24})\\ {\hat{D}}_{33}& {\hat{p}}_{33}-({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})({\hat{p}}_{31}+{\hat{p}}_{32}+{\hat{p}}_{33}+{\hat{p}}_{34})\\ {\hat{D}}_{34}& -{\hat{D}}_{31}-{\hat{D}}_{32}-{\hat{D}}_{33}\\ {\hat{D}}_{41}& -{\hat{D}}_{11}-{\hat{D}}_{21}-{\hat{D}}_{31}\\ {\hat{D}}_{42}& -{\hat{D}}_{12}-{\hat{D}}_{22}-{\hat{D}}_{32}\\ {\hat{D}}_{43}& -{\hat{D}}_{13}-{\hat{D}}_{23}-{\hat{D}}_{33}\\ {\hat{D}}_{44}& {\hat{D}}_{11}+{\hat{D}}_{12}+{\hat{D}}_{13}+{\hat{D}}_{21}+{\hat{D}}_{22}+{\hat{D}}_{23}+{\hat{D}}_{31}+{\hat{D}}_{32}+{\hat{D}}_{33}\end{array}$

步骤5、根据以下公式验证估计分布与实际分布是否存在显著性差异；

H₀：D₁₁＝D₁₂＝D₁₃＝D₂₁＝D₂₂＝D₂₃＝D₃₁＝D₃₂＝D₃₃＝0；

H₁：Not H₀ (5)

其中H₀为同源四倍体各位点不存在连锁不平衡系数；

通过最大似然估计，分别估算和波浪线对应的零假设的最大似然估计，抑扬符号对应的备择假设的最大似然估计；

通过以下公式计算对数似然比的统计学检验

$LR=-2\[\ln L\left(\widetilde{\mathrm{\Ω}}\right)-\ln L\left(\widehat{\mathrm{\Ω}}\right)\]---\left(6\right)$

LR是遵循χ²分布，自由度是9；

将计算得到的LR进行卡方检验，判断得到的LR值经卡方检验是否呈显著差异，如果呈显著差异则最终结果为H₀，如果不呈显著差异则最终结果为H₁。

2.根据权利要求1所述的同源四倍体自然群体的连锁不平衡分析模型的构建方法，其特征在于，还包括：

步骤6、将9个连锁不平衡系数分别进行检验。

3.根据权利要求1所述的同源四倍体自然群体的连锁不平衡分析模型的构建方法，其特征在于，

所述步骤1中，所述16单倍型的频率包含以下等位基因频率以及九个独立的连锁不平衡系数：

$\begin{array}{cc}{p}_{11}& p_1{q}_1+D_{11}\\ p_{12}& p_1{q}_2+D_{12}\\ p_{13}& p_1{q}_3+D_{13}\\ p_{14}& p_1{q}_4-D_{11}-D_{12}-D_{13}\\ p_{21}& p_2{q}_1+D_{21}\end{array}$

$\begin{array}{cc}{p}_{22}& p_2{q}_2+D_{22}\\ p_{23}& p_2{q}_3+D_{23}\\ p_{24}& p_2{q}_4-D_{21}-D_{22}-D_{23}\\ p_{31}& p_3{q}_1+D_{31}\\ p_{32}& p_3{q}_2+D_{32}\\ p_{33}& p_3{q}_3+D_{33}\\ p_{34}& p_3{q}_4-D_{31}-D_{32}-D_{33}\\ p_{41}& p_4{q}_1-D_{11}-D_{21}-D_{31}\\ p_{42}& p_4{q}_2-D_{12}-D_{22}-D_{32}\\ p_{43}& p_4{q}_3-D_{13}-D_{23}-D_{33}\\ p_{44}& p_4{q}_4+D_{11}+D_{12}+D_{13}+D_{21}+D_{22}+D_{23}+D_{31}+D_{32}+D_{33}\end{array}.$

4.根据权利要求1所述的同源四倍体自然群体的连锁不平衡分析模型的构建方法，其特征在于，

所述步骤1中，所述位点A和位点B分别都具有35种基因型，在遗传上能够识别基因型共225种；

位点A有，4种A_iA_iA_iA_i，其中i＝1，2，3，4；12种其中i₁≠i₂＝1，2，3，4；6种其中i₁≠i₂＝1，2，3，4；12种其中i₁≠i₂≠i₃＝1，2，3，4；1种A₁A₂A₃A₄；

位点B有，4种B_jB_jB_jB_j，其中j＝1，2，3，4；12种其中j₁≠j₂＝1，2，3，4；6种其中j₁≠j₂＝1，2，3，4；12种其中j₁≠j₂≠j₃＝1，2，3，4；1种B₁B₂B₃B₄；

两个位点，在遗传上能识别的基因型有225种，表示为Au/Bv(u＝1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234；v＝1，2，3，4，12，13，14，23，24，34，123，124，134，234，1234)。