一种异源六倍体遗传连锁分析方法

摘要

本发明公开一种异源六倍体遗传连锁分析方法，其推导是基于异源多倍体减数分裂过程中染色体的优先配对现象；优先配对包括三种类型；分别获取：六条染色体两两优先配对的1种情况，六条染色体中的两条染色体优先配对且其他四条染色体自由配对的6种情况，六条染色体全部随机配对的8种情况；生成基于异源六倍体优先配对因子的分子标记共分离基础模型，并建立共分离概率矩阵；根据概率矩阵中的每一个元素，构建过渡矩阵；统计实际数据中每一种配子的真实数量，并根据可能的配子情况统计出M矩阵；根据总概率矩阵、过渡矩阵、M矩阵计算异源六倍体减数分裂过程中优先配对因子的估算值；通过不断迭代重新计算重组率r；以此构建异源六倍体遗传连锁图。

说明书

技术领域

本发明涉及生物建模技术领域，特别是指一种异源六倍体遗传连锁分析方 法。

背景技术

多倍体在被子植物中是很常见的，大约有70%的被子植物在他们各自的进 化历程中发生过多倍体化(Masterson,1994)，多倍体化在植物进化过程中起到 了至关重要的作用。植物多倍体化被认为是植物进化的主要动力，因此多倍体 研究对植物进化研究具有指导意义。

多倍体按照起源可以分为同源多倍体和异源多倍体。同源多倍体的多个染 色体组是由同一物种经染色体加倍形成的；异源多倍体的染色体组由不同物种 得来。根据同源染色体倾向于优先配对的原则，多倍体的多条染色体在减数分 裂过程中，异源多倍体的染色体会产生优先配对现象。根据对多倍体染色体优 先配对因子的估算，有助于探究多倍体的物种起源，并且对多倍体的进化研究 有重要意义。

2001年国际上首次把优先配对现象融入多倍体统计模型之中。但是在近 10年，对多倍体优先配对现象及优先配对因子的估算，还仅仅停留在对四倍 体的分析上，而对六倍体优先配对因子的估算尚属空白。其主要原因在于：现 有估算四倍体优先配对因子的模型仅仅涉及一个参数，估算系统相对简单；而 六倍体由于具有六条染色体，若想估算其优先配对因子至少需要同时估算两个 参数，这就为估算增大了难度。

发明内容

本发明要解决的技术问题是提供一种异源六倍体遗传连锁分析方法，能够 对多倍体优先配对因子进行估算。

为了解决上述问题，本发明实施例提出了一种异源六倍体遗传连锁分析方 法，包括：

步骤1、根据异源六倍体的六条染色体两两优先配对、两条优先配对四条 随机配对、六条全部随机配对的三种配对方式，分别获取：六条染色体两两优 先配对的1种情况，六条染色体中的两条染色体优先配对且其他四条染色体自 由配对的6种情况，六条染色体全部随机配对的8种情况；

步骤2、根据步骤1的三种配对方式，确定每种配对方式的单个亲本的三 倍体配子；

步骤3、针对每一种配对方式中的每一种情况分别构建每一种配子产生的 概率矩阵作为概率表；针对每一种情况的每一种配子产生的概率表，计算每一 种配对方式中的全部配子的概率总表，作为六倍体优先配对因子估算的基础模 型；根据每种配对方式的配子的概率总表，计算每一优先配对情况下的配子产 生概率，并针对所有配子生成总概率矩阵；

步骤4、根据总概率矩阵中的每一元素，构建过渡矩阵；统计实际数据中 每一种配子的真实数量，并根据可能的配子情况统计出M矩阵；通过迭代重新 计算重组率r；

步骤5、根据总概率矩阵、过渡矩阵、M矩阵计算异源六倍体减数分裂过 程中优先配对因子的估算值。

作为上述技术方案的优选，所述步骤2具体为：

步骤21、在异源六条染色体两两配对的15种配对方式下，每一种配对情 况产生8种三倍体的配子，经过合并共得到20种三倍体的配子：

作为上述技术方案的优选，所述步骤3具体为：

步骤31、在异源六条染色体两两优先配对的配对方式下，通过一条染色 体上两个标记位点之间的重组率推导出20种配子的概率；这些概率是由重组 率以及优先配对因子决定的；

步骤32、计算15种配对情况中的每一种配对情况的配子的概率矩阵：

${\mathrm{Pf}}_1=\left(\begin{array}{cccccccc}{\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3\\ \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)\\ \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3\end{array}\right);$

其中Pf₁至Pf₁₅在数值上相等;

步骤33、将15个概率矩阵合并生成20种配子的概率总表以作为异源六 倍体优先配对因子估算的基础模型；

作为上述技术方案的优选，所述步骤4具体为：

步骤41、在六条染色体两两配对的配对方式下，根据所述概率总表计算 20种配子的总概率矩阵P；

$P=\underset{k=1}{\overset{15}{\Sigma }}{A}_k$

其中A_k为第k种优先配对情况中各种配子的概率矩阵；

$A_k=\left(\begin{array}{cc}{\mathrm{Pf}}_k\times (\frac{1}{15}+{\theta }_1),& k=1;\\ {\mathrm{Pf}}_k\times (\frac{1}{15}+\frac{1}{6}\times {\theta }_2),& k=2,···,7;\\ {\mathrm{Pf}}_k\times (\frac{1}{15}-\frac{1}{8}\times {\theta }_1-\frac{1}{8}\times {\theta }_2),& k=8,···,15\end{array}\right)$

其中Pf_k（k=1，…，15）为所述15种情况下的第k种优先配对情况下配 子产生的概率矩阵。

步骤42、构建过渡矩阵Q；

总概率矩阵P当中的每一个元素都是由若干个加号连接的关于重组率r 以及优先配对因子θ的多项式；把加号连接的每一个多项式作为独立的个体， 把每一个独立多项式中r的幂指数乘以该多项式，形成新的矩阵作为过渡矩阵 Q：

步骤43、构建过渡矩阵Ψ；

${\Psi }_{\mathrm{ij}}=\frac{P_{\mathrm{ij}}}{Q_{\mathrm{ij}}};$

步骤44、统计实际数据中每一种配子的真实数量，按照20种理论上可能 产生的配子情况，统计得出M矩阵；

步骤45、根据EM算法估算重组率r：

其中N为染色体的数量；

作为上述技术方案的优选，所述步骤5具体为：

步骤51、在六条染色体两两配对的配对方式下，根据P矩阵、Q矩阵、M 矩阵、Ψ矩阵，计算所述配对方式下的优先配对因子θ₁，θ₂；

${\hat{\theta }}_1=\Sigma \frac{M\times \Psi }{N};$

${\hat{\theta }}_2=\Sigma \frac{M\times \Psi }{N};$

步骤52、根据θ₁，θ₂计算θ₃；(θ₁+θ₂+θ₃)=1。

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

本发明提供了一种异源六倍体遗传连锁分析方法，通过实验可以看出，采 用本发明的方法获得的参数估算的结果与其真值非常接近，均值相差很小，方 差也很小，说明本发明的估算方法具有很好的准确度和精确度。

附图说明

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

具体实施方式

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

本发明提供了一种异源六倍体遗传连锁分析方法，其流程如图1所示的， 包括：

步骤A、构建估算异源六倍体杂交群体优先配对因子的统计模型；构建步 骤为：

步骤1、定义优先配对因子：

把异源六倍体的六条染色体分别标号为：1、2、3、4、5、6。这六条染色 体配对会有三种配对方式：①六条染色体两两优先配对、②两条优先配对四条 随机配对、③六条全部随机配对的三种配对方式。其中：

①六条染色体两两优先配对会产生1种配对情况：（1）1和2配对，3和 4配对，则5和6自然配对，把这种情况命名为全部亲和。把这种情况的优先 配对因子定义为θ₁。

②两条优先配对四条随机配对，其余四条染色体自由组合：这样的配对情 况有6种，把这些配对情况的优先配对因子定义为θ₂。

③异源六倍体的六条染色体在减数分裂过程中不出现亲和配对，全部随机 配对：这样的配对情况有8种，把这些配对情况的优先配对因子定义为θ₃。

步骤2、估算异源六倍体亲本减数分裂产生的配子类型。

在第一种配对方式下，如果异源六倍体的六条染色体中1和2配对，3和 4配对，则5和6自然配对。即会产生如下1种配对情况⑴1|2| 3|4| 5|6|； 同时这个六倍体亲本将产生8种3倍体的配子：1|3|5|，1|3|6|，2|3|5|， 2|3|6|，1|4|5|，1|4|6|，2|4|5|，2|4|6|；其中1|代表1号染色体，以此 类推。如果考虑两位点的情况，位点A会产生8种配子，位点B会产生8中配 子，那么综合考虑两个位点的信息则会产生64种配子类型。这种情况下配子 的概率如下表所示：

表1.情况⑴对应产生两位点64种配子概率

在第二种配对方式下，如果六倍体的六条染色体中某两条染色体优先配 对，其余四条染色体随机组合配对。则会产生如下6种配对情况；

⑵1|2| 3|5| 4|6|，⑶1|2| 3|6| 4|5|，（1号2号染色体优先配对）

⑷1|5|3|4|2|6|，⑸1|6|3|4|2|5|，（3号4号染色体优先配对）

⑹1|3|2|4|5|6|，⑺1|4|2|3|5|6|。（5号6号染色体优先配对）

同时每一种配对情况六倍体亲本将对应产生8种3倍体的配子：

情况⑵对应的8种三倍体配子为：1|3|4|，1|3|6|，2|3|4|，2|3|6|， 1|4|5|，1|5|6|，2|4|5|，2|5|6|；情况⑶对应的8种三倍体配子为：1|3|4|， 1|3|5|，2|3|4|，2|3|5|，1|6|4|，1|6|5|，2|6|4|，2|6|5|；情况⑷对应的 8种三倍体配子为：1|2|3|，1|3|6|，2|3|5|，3|5|6|，1|2|4|，1|4|6|，2|4|5|， 4|5|6|；情况⑸对应的8种三倍体配子为：1|3|2|，1|3|5|，1|4|2|，1|4|5|， 6|3|2|，3|5|6|，2|4|6|，4|5|6|；情况⑹对应的8种三倍体配子为：1|2|5|， 1|2|6|，1|4|5|，1|4|6|，2|3|5|，2|3|6|，3|4|5|，3|4|6|；情况⑺对应的 8种三倍体配子为：1|2|5|，1|2|6|，1|3|5|，1|3|6|，2|4|5|，2|4|6|，3|4|5|， 3|4|6|；

在第三种配对方式下，六倍体的六条染色体全部随机组合配对。则会产生 如下8种配对情况；

⑼1|3|2|5|4|6|，⑼1|3|2|6|4|5|，⑽1|4|2|5|3|6|，

⑾1|4|2|6|3|5|，⑿1|5|2|3|4|6|，⒀1|5|3|6|2|4|，

⒁1|6|2|3|4|5|，⒂1|6|3|5|2|4|。

同时每一种配对情况六倍体亲本将对应产生8种3倍体的配子：

情况⑻对应的8种三倍体配子为：1|2|4|，1|2|6|，1|4|5|，1|5|6|， 2|3|4|，2|3|6|，3|4|5|，3|5|6|；情况⑼对应的8种三倍体配子为：1|2|4|， 1|2|5|，1|4|6|，1|5|6|，2|3|4|，2|3|5|，3|4|6|，3|4|5|；情况⑽对应的 8种三倍体配子为：1|2|3|，1|2|6|，1|3|5|，1|5|6|，2|3|4|，2|4|6|，3|4|5|， 4|5|6|；情况⑾对应的8种三倍体配子为：1|2|3|，1|2|5|，1|3|6|，1|5|6|， 2|3|4|，2|4|5|，3|4|6|，4|5|6|；情况⑿对应的8种三倍体配子为：1|2|4|， 1|2|6|，1|3|4|，1|3|6|，2|4|5|，2|5|6|，3|4|5|，3|5|6|；情况⒀对应的 8种三倍体配子为：1|2|3|，1|3|4|，1|2|6|，1|4|6|，2|3|5|，3|4|5|，2|5|6|， 4|5|6|；情况⒁对应的8种三倍体配子为：1|2|4|，1|2|5|，1|3|4|，1|3|5|， 2|4|6|，2|5|6|，3|4|6|，3|5|6|；情况⒂对应的8种三倍体配子为：1|2|3|， 1|3|4|，1|2|5|，1|4|5|，2|3|6|，3|4|6|，2|5|6|，4|5|6|；

因此一共会产生如下所示的20种三倍体的配子：

步骤3、构建估算六倍体杂交群体优先配对因子的统计模型。

步骤31、首先需要构建每一种配子产生的概率矩阵：

根据一条染色体上两个标记位点之间的重组率能够推导出这20种配子的 概率；这些概率是由重组率以及优先配对因子决定的。

根据步骤1可知六倍体减数分裂理论上会产生3倍体的配子，每一种优先 配对情况都对应产生8种配子。例如，全部优先配对情况下（1|2|3|4| 5|6|）会产生8种配子：1|3|5|，1|3|6|，1|4|5|，1|4|6|，2|3|5|，2|3|6|， 2|4|5|，2|4|6|。如果考虑两个标记位点，那么配对情况⑴产生的配子概率情 况如下表：

表1.情况⑴对应产生两位点64种配子概率

情况⑵-⑻对应的64种配子概率分别由表2-表8给出：

表2.情况⑵对应产生两位点64种配子概率

表3.情况⑶对应产生两位点64种配子概率

表4.情况⑷对应产生两位点64种配子概率

表5.情况⑸对应产生两位点64种配子概率

表6.情况⑹对应产生两位点64种配子概率

表7.情况⑺对应产生两位点64种配子概率

情况⑻-⒂对应的64种配子概率分别由表8-表15给出：

表8.情况⑴对应产生两位点64种配子概率

表9.情况⑼对应产生两位点64种配子概率

表10.情况⑽对应产生两位点64种配子概率

表11.情况⑾对应产生两位点64种配子概率

表12.情况⑿对应产生两位点64种配子概率

表13.情况⒀对应产生两位点64种配子概率

表14.情况⒁对应产生两位点64种配子概率

表15.情况⒂对应产生两位点64种配子概率

步骤32、根据步骤1可知，异源六倍体减数分裂染色体配对方式共有15 种配对情况，对每一种情况按照步骤2计算对应配子的概率表，这样会将得到 15个8*8的概率矩阵。将这15个概率矩阵合并，得到全部20种配子的概率 总表；这个概率总表就是异源六倍体优先配对因子估算的基础模型;

${\mathrm{Pf}}_1=\left(\begin{array}{cccccccc}{\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3\\ \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)\\ \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{1-r}{2}\right)}^3& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2\\ {\left(\frac{r}{2}\right)}^3& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{r}{2}\right)}^2\times \left(\frac{1-r}{2}\right)& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& \frac{r}{2}\times {\left(\frac{1-r}{2}\right)}^2& {\left(\frac{1-r}{2}\right)}^3\end{array}\right)$

Pf₂…Pf₁₅在数值上与Pf₁相等。

位点A共可能产生20种配子，位点B共可能产生20种配子，则考虑两个 位点的情况，六倍体共可能产生400种配子类型；该400种配子类型为一20*20 的矩阵。

步骤4、估算重组率r；

步骤41、根据前述的概率总表，设这15种情况下的第k种优先配对情况 下，配子产生的概率为Pf_k（k=1，…，15）；设A_k为第k种优先配对情况中各 种配子的概率矩阵；

$A_k=\left(\begin{array}{cc}{\mathrm{Pf}}_k\times (\frac{1}{15}+{\theta }_1),& k=1;\\ {\mathrm{Pf}}_k\times (\frac{1}{15}+\frac{1}{6}\times {\theta }_2),& k=2,···,7;\\ {\mathrm{Pf}}_k\times (\frac{1}{15}-\frac{1}{8}\times {\theta }_1-\frac{1}{8}\times {\theta }_2),& k=8,···,15\end{array}\right)$

则全部20种配子的总概率矩阵P为

$P=\underset{k=1}{\overset{15}{\Sigma }}{A}_k.$

步骤42、构建过渡矩阵Q；Q矩阵是根据计算系统临时定义的矩阵，没有 实际生物学含义。根据估算参数不同，Q矩阵形式也不同，分为Qr矩阵、Qθ ₁矩阵和Qθ₂矩阵。

总概率矩阵P当中的每一个元素都是由若干个加号连接的关于重组率r 以及优先配对因子θ的多项式。把加号连接的每一个多项式看成独立的个体， 要估算重组率r，就要把每一个独立多项式中r的幂指数乘以该多项式，得到

估算Qθ₁矩阵和Qθ₂矩阵的方法与估算Qr矩阵类似，这里不再赘述。 并统称过渡矩阵为Q矩阵。

步骤43、构建过渡矩阵Ψ。Ψ矩阵是根据计算系统临时定义的矩阵，没 有实际生物学含义。

${\Psi }_{\mathrm{ij}}=\frac{P_{\mathrm{ij}}}{Q_{\mathrm{ij}}}$

步骤44、统计实际数据中每一种配子的真实数量，按照20种理论上可能 产生的配子情况，统计得出M矩阵。

步骤45、根据EM算法估算重组率r：

$\hat{r}=\Sigma \frac{M\times \Psi }{3\times N}$

其中N为染色体的数量；

步骤5、由步骤4可以得出P矩阵、Q矩阵、M矩阵、Ψ矩阵，由此计算。 则

${\hat{\theta }}_1=\Sigma \frac{M\times \Psi }{N};{\hat{\theta }}_2=\Sigma \frac{M\times \Psi }{N}.$

其中，步骤5具体为：

步骤51、在异源六条染色体两两配对的配对方式下，根据P矩阵、Q矩阵、 M矩阵、Ψ矩阵，计算所述配对方式下的优先配对因子θ₁，θ₂；

${\hat{\theta }}_1=\Sigma \frac{M\times \Psi }{N};$

${\hat{\theta }}_2=\Sigma \frac{M\times \Psi }{N};$

步骤52、根据θ₁，θ₂计算θ₃；(θ1+θ₂+θ₃)=1。

根据步骤1-5得出异源六倍体减数分裂过程中优先配对因子的估算值。

以下通过计算机模拟实验验证算法准确性及精度。

假设优先配对因子θ_i（i=1,2,3）满足以下两个条件：（1）(0﹤θi＜ 1,i=1,2,3).（2）(θ₁+θ₂+θ₃)=1。

假设一个亲本分别为P₁（1|2|3|4|5|6|）和P₂（1|1|1|1|1|1|）。假设他 们子代样本量分别为n=200,400,800。

假设两个标记位点A和B之间的重组率有四种，分别是两个高度连锁（r= 0.05,0.15）和两个低连锁程度（r=0.25,0.35）。

针对两位点全信息量连锁相两位点部分信息连锁相 以及显性标记连锁相这三种有代表性的连锁相，其 中每一种连锁相都对应上述参数做1000次计算机模拟实验，充分验证估算六 倍体染色体优先配对因子模型的准确度与精确度。

其中，在样本量n=200时，各参数的估算结果如下表所示的:

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