一种基于混合线性模型的多性状关联分析方法

摘要

本发明公开了一种基于混合线性模型的多性状关联分析方法，该方法包括：构建统计遗传模型；确定效应显著的单位点SNP标记；确定互作效应显著的二互作上位性SNP标记；估算遗传效应。本发明基于多变量混合线性模型的多性状全基因组关联分析方法综合利用了多个遗传相关性状的变异信息，与单性状分析方法相比，具有较高的分析功效和较低的假阳性，QTS位置估算更为准确，效应估计更为稳健。

说明书

技术领域

本发明涉及多性状联合定位技术领域，尤其涉及一种基于混合线性模型的多性状 关联分析方法。

背景技术

全基因组关联分析已经变成一种标准的探索复杂性状遗传结构和解释数量性状 遗传变异基础的有效方法。关联分析方法中的主要问题在于解释数据的依赖性，包括个体 之间的依赖和位点之间的依赖性。混合线性模型同时包含固定效应和随机效应，可以有效 的解释大数据中存在的群体结构，包括群体分层和亲缘关系。

然而，大多数的关联分析方法都只针对单个数量性状，不能综合考虑多个遗传相 关的性状进行联合分析，无法剖析多个性状遗传相关的分子机理，无法分析遗传位点的多 效性等，例如：

申请公布号CN103632067A的发明专利申请文献公开了一种基于混合线性模型的 种子数量性状位点定位方法，该方法包括：(1)统计遗传模型的建立；(2)全基因组扫描显著 的标记区间：基于模型在全基因组范围内，通过表型和每个标记区间做基于HendersonIII 的F检验，搜索得到所有可能存在QTL的候选标记区间；(3)在全基因组范围内搜索显著的数 量性状位点，将步骤(2)中得到的候选标记区间作为协变量，然后以1cm为步长，在全基因组 范围内做基于HendersonIII的F检验，搜索显著的QTL位点；(4)在全基因组范围内搜索显著 的二互作标记区间：将步骤(2)中得到的候选标记区间作为协变量，然后在全基因组范围内 做基于HendersonIII的F检验，搜索得到显著的二互作标记区间；(5)搜索显著的二互作上 位性的位点：以步骤(3)得到的QTL以及步骤(4)得到的显著的互作标记区间作为协变量，然 后在所述显著的互作标记区间中做基于HendersonIII的F检验，搜索得到显著的二互作上 位性的位点；(6)遗传参数的估算：通过步骤(3)得到的显著的QTL位点和步骤(5)得到的显 著的二互作上位性位点，获取得到模型中各种效应系数，然后通过模型，计算得到这些位点 的效应以及估算每个位点的遗传率。

目前，已有一些基于多变量混合线性模型的多性状关联分析方法，如MTMM(Krote etal.，2012，NatureGenetics，44(9)：1066-1071)，GEMMA(Zhouetal.,2014.Nature Methods,11(4),407-409)等。但是，这些方法或者计算量大，或者无法分析三个或三个以上 的形状，或者无法分析基因间互作、基因与环境间互作效应。

因此，有必要再探究一种更为高效的能够进行多性状关联分析的定位方法，已解 决上述问题。

发明内容

本发明提供了一种基于混合线性模型的多性状联合定位方法，该方法分析功效较 高、假阳性较低，QTS位置估算更准确，效应估计更稳健。

一种基于混合线性模型的多性状关联分析方法，包括：

(1)构建统计遗传模型：

假设一个自然群体由n个个体组成，在p个不同环境中进行田间试验，m个相关性状 的遗传变异受s个QTS位点和t对二互作上位性的调控；环境k中针对第i个性状的第j个株系 的表型观测值y_ijk表示为：

$\left(\begin{array}{c}{y}_{ijk}={\mu }_i+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+\underset{l}{\overset{s}{\Sigma }}{a}_{il}{x}_{lj}+\underset{l,h\in (1,\mathrm{2...},s),l<h}{\overset{t}{\Sigma }}{\mathrm{aa}}_{ilh}{x}_{lj}{x}_{hj}\\ +e_{ik}+\underset{l}{\overset{s}{\Sigma }}{\mathrm{ae}}_{ilk}{x}_{ljk}+\underset{l,h\in (1,\mathrm{2...},s),l<h}{\overset{t}{\Sigma }}{\mathrm{aae}}_{ilhk}{x}_{ljk}{x}_{hjk}+{ϵ}_{ijk}\end{array}\right)$

式中，μ_i是性状i的群体均值；b_ic是性状i第c个协变量的效应，系数为x_jkc；a_il是性 状i第l个QTS的加性效应，系数为x_lj；aa_ilh是性状i第l个QTS与第h个QTS之间的加加上位性， 系数为x_ljx_hj；e_ik是性状i在环境k下的效应；ae_ilk是a_il与第k个环境的互作效应，系数为x_ljk； aae_ilhk是aa_ilh与第k个环境的互作效应，系数为x_ljkx_hjk；ε_ijk是性状i株系j在环境k下的随机 残差；

(2)确定效应显著的单位点SNP标记：

在全基因组范围内，逐一检测各SNP标记，通过Lambda统计量和置换检验方法，获 取效应显著的单位点SNP标记；

$y_{ijk}={\mu }_{ik}+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+a_{ilk}{x}_{lj}+{ϵ}_{ijk}$

式中，y_ijk是环境k中第i个性状的第j个株系的表型观测值；μ_ik是环境k中性状i的 群体均值；b_ic是性状i第c个协变量的效应，系数为x_jkc；a_ilk是环境k中性状i位点l的加性效 应，系数为x_lj；ε_ijk是株系j在环境k下第i个性状的随机残差；

(3)确定互作效应显著的二互作上位性SNP标记：

将步骤(2)中获得的单位点SNP标记作为协变量，通过Lambda统计量和置换检验方 法，获取互作效应显著的二互作上位性SNP标记；

$y_{ijk}={\mu }_{ik}+{\mathrm{aa}}_{ilhk}{x}_{hj}{x}_{lj}+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+\underset{r=1}{\overset{s}{\Sigma }}{a}_{irk}{x}_{rj}+{ϵ}_{ijk}$

式中，y_ijk是环境k中针对第i个性状的第j个株系的表型观测值；μ_ik是环境k中性状 i的群体均值；aa_ilhk是性状i在环境k下第l个QTS与第h个QTS之间的加加上位性效应，系数 为x_ljx_hj；b_ic是性状i第c个协变量的效应，系数为x_jkc；a_irk是在步骤(2)中获得的所述单位点 SNP标记的加性效应，系数为x_rj；ε_ijk是株系j在环境k下第i个性状的随机残差；

(4)估算遗传效应：

针对步骤(2)得到的单位点SNP标记和步骤(3)得到的二互作上位性SNP标记，采用 向前选择法剔除假阳性的单位点SNP标记和二互作上位性SNP标记，得到显著的单位点QTS 和二互作上位性QTS的效应系数，将所述效应系数代入步骤(1)中构建全模型，计算得到所 述单位点QTS和二互作上位性QTS的遗传效应。

具体地，所述的遗传效应为单位点QTS的加性效应，二互作上位性QTS的加加上位 性效应，单位点QTS与环境的互作效应以及二互作上位性QTS与环境的互作效应。

步骤(2)中，所述Lambda统计量的计算方法为：采用矩阵向量的方式表示表型观测 值Y，

Y＝W_QB_Q+W_MB_M+E_ε

式中，Y是一个n×m阶表型观测值矩阵；B_Q＝(a_1la_2l…a_ml)是p×m阶矩阵，由检 测位点l各环境下的加性效应组成，其中a_il＝(a_i1la_i2l…a_ipl)^T；Β_Μ是协变量矩阵，包含 了各个环境下的群体均值和各协变量的效应；Β_Q和Β_Μ的系数矩阵分别是W_Q和W_M；E_ε是m个 性状的随机残差矩阵；

其中，W＝[W_Q,W_M],Q＝Y^TY-Y^TW(W^TW)⁺W^TY,Q₁＝Y^TY-Y^TW_M(W_M^TW_M)⁺W_M^TY，构建Wilks’ Lambda统计量λ＝|Q|/|Q₁|；在原假设H₀:B_Q＝0的情况下，统计量服从Lambda分布。

进一步地，步骤(2)中，在全基因组范围内进行搜索，通过所述Lambda统计量的计 算方法计算每一个SNP标记的Lambda统计量，并通过置换检验方法确定效应显著的阈值；然 后，判断Lambda统计量与阈值的大小关系；当某一单位点SNP标记的Lambda统计量大于或等 于阈值时，该SNP标记为效应显著的单位点SNP标记。

步骤(3)中，所述Lambda统计量的计算方法为：采用矩阵向量的方式表示表型观测 值Y，

Y＝W_QB_Q+W_MB_M+E_ε

式中，Y是一个n×m阶表型观测值矩阵；B_Q＝(aa_1lhaa_2lh…aa_mlh)是p×m阶矩阵， 由检测位点l和位点h各环境下的加加上位性效应组成，其中aa_ilh＝(aa_i1lhaa_i2lh… aa_iplh)^T；Β_Μ是协变量矩阵，包含了各个环境下的群体均值和步骤(2)中获得的显著的单位 点SNP标记；Β_Q和Β_Μ的系数矩阵分别是W_Q和W_M；E_ε是m个性状的随机残差矩阵；

其中，W＝[W_Q,W_M],Q＝Y^TY-Y^TW(W^TW)⁺W^TY,Q₁＝Y^TY-Y^TW_M(W_M^TW_M)⁺W_M^TY，构建Wilks’ Lambda统计量λ＝|Q|/|Q₁|；在原假设H₀:B_QQ＝0的情况下，统计量服从Lambda分布。

进一步地，步骤(3)中，在步骤(2)获得的单位点SNP标记的基础上进行搜索，通过 所述Lambda统计量的计算方法计算每一对二互作上位性SNP标记的Lambda统计量；并通过 置换检验方法确定效应显著的阈值；然后，判断Lambda统计量与阈值的大小关系；当某一对 二互作上位性SNP标记的Lambda统计量大于或等于阈值时，该对SNP标记为互作效应显著的 二互作上位性SNP标记。

步骤(1)中，所述表型观测值采用矩阵向量进行表示；

$\left(\begin{array}{c}Y=1\mu +X_C{B}_C+X_A{B}_A+X_{AA}{B}_{AA}+U_E{E}_E+\underset{l=1}{\overset{s}{\Sigma }}{U}_{A_{l}E}{E}_{A_{l}E}+\underset{l<h}{\overset{t}{\underset{l,h\in (1,2,...,s),}{\Sigma }}}{U}_{{\mathrm{AA}}_{lh}E}{E}_{{\mathrm{AA}}_{lh}E}+{\mathrm{IE}}_{ϵ}\\ =XB+\underset{u=1}{\overset{r+1}{\Sigma }}{U}_u{E}_u\end{array}\right)$

式中，Y是一个n×m阶性状观测值矩阵；μ＝(μ₁μ₂…μ_m)是一个1×m阶行向量； $B_C=\left(\begin{array}{cccc}{b}_{C_1}& b_{C_2}& ...& b_{C_q}\end{array}\right)$是一个q×m阶向量，表示m个性状q个协变量的效应；X_C是Β_C的系 数矩阵；$B_A=\left(\begin{array}{cccc}{b}_{A_1}& b_{A_2}& ...& b_{A_m}\end{array}\right)$是一个s×m阶矩阵，表示m个性状s个QTS的加性效应矩 阵，其中$b_{A_k}={\left(\begin{array}{cccc}{a}_{k1}& a_{k2}& ...& a_{ks}\end{array}\right)}^T,$表示性状k显著QTS的加性效应向量；X_A是Β_A的系数矩阵； $B_{AA}=\left(\begin{array}{cccc}{b}_{{\mathrm{AA}}_1}& b_{{\mathrm{AA}}_2}& ...& b_{{\mathrm{AA}}_m}\end{array}\right)$是一个t×m阶矩阵，表示m个性状t对QTS的加加上位性效应矩阵， 其中$b_{{\mathrm{AA}}_k}={\left(\begin{array}{cccc}{\mathrm{aa}}_{k1}& {\mathrm{aa}}_{k2}& ...& {\mathrm{aa}}_{kt}\end{array}\right)}^T;$X_AA是B_AA的系数矩阵；环境效应矩阵E_E＝(e₁e₂…e_m)，其 中e_k＝(e_k1e_k2…e_kp)^T是性状k的环境效应向量，随机效应，服从多元正态分布 是性状k的环境效应方差；U_E是E_E的系数；$E_{A_{l}E}=\left(\begin{array}{cccc}{e}_{A_{1l}E}& e_{A_{2l}E}& ...& e_{A_{ml}E}\end{array}\right)$是m个性状第l个QTS的加性与环境互作效应矩阵，是随机效应，其中$e_{A_{kl}E}={\left(\begin{array}{cccc}{\mathrm{ae}}_{kl1}& {\mathrm{ae}}_{kl2}& ...& {\mathrm{ae}}_{klp}\end{array}\right)}^T;$是的系数；$E_{{\mathrm{AA}}_{lh}E}=\left(\begin{array}{cccc}{e}_{{\mathrm{AA}}_{1lh}E}& e_{{\mathrm{AA}}_{2lh}E}& ...& e_{{\mathrm{AA}}_{mlh}E}\end{array}\right)$是m个性状的加加上位性与环境互作 效应矩阵，随机效应，其中$e_{{\mathrm{AA}}_{klh}E}={\left(\begin{array}{cccc}{\mathrm{aae}}_{klh1}& {\mathrm{aae}}_{klh2}& ...& {\mathrm{aae}}_{klhp}\end{array}\right)}^T;$是的系数； $E_{ϵ}=\left(\begin{array}{cccc}{e}_{{ϵ}_1}& e_{{ϵ}_2}& ...& e_{{ϵ}_m}\end{array}\right)$是m个性状的随机残差效应矩阵；I是单位矩阵。

进一步地，步骤(4)中，所述效应系数代入步骤(1)中构建全模型进行计算时，采用 蒙特卡罗马尔可夫链方法进行吉布斯抽样，获得遗传效应参数的分布样本，根据所述分布 样本估算单位点QTS和二互作上位性QTS的遗传效应值，并进行显著性检验。

与现有技术相比，本发明具有以下有益效果：

(1)本发明基于多变量混合线性模型的多性状全基因组关联分析方法综合利用了 多个遗传相关性状的变异信息，与单性状分析方法相比，具有较高的分析功效和较低的假 阳性，QTS位置估算更为准确，效应估计更为稳健。

(2)本发明采用包括上位性及位点与环境互作效应的模型，更符合复杂性状的遗 传特征，全模型方法能更准确地估算QTS各项效应分量，也易于解释各项效应分量的生物学 及育种含义。

附图说明

图1为本发明基于混合线性模型的多性状联合定位方法的流程图。

具体实施方式

下面结合具体实施方式对本发明作进一步的阐释。

一种基于混合线性模型的多性状联合定位方法，包括以下步骤：

(一)统计遗传模型的建立：

假设一自然群体由n个个体组成，在p个环境中进行了田间试验，调查每个个体的 性状表型值，对各遗传材料进行基因型检测。

假设s个QTS位点控制了m个相关性状的遗传变异，其中部分QTS参与t对二互作上 位性。

对于第i个表型性状(i＝1,2,…,m)，第l个QTS的加性效应记作a_il，第l个QTS与第h 个QTS之间的加加上位性记作aa_ilh，都作固定效应处理；ae_ilk和aae_ilhk分别是a_il、aa_ilh与第k 个环境的互作效应，都作随机效应。

那么，对于m个相互关联的表型性状，株系j在环境k中的观察值可以用如下的混合 线性模型表示：

$\begin{array}{c}{y}_{ijk}={\mu }_i+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+\underset{l}{\overset{s}{\Sigma }}{a}_{il}{x}_{lj}+\underset{l,h\in (1,\mathrm{2...},s),l<h}{\overset{t}{\Sigma }}{\mathrm{aa}}_{ilh}{x}_{lj}{x}_{hj}\\ +e_{ik}+\underset{l}{\overset{s}{\Sigma }}{\mathrm{ae}}_{ilk}{x}_{ljk}+\underset{l,h\in (1,\mathrm{2...},s),l<h}{\overset{t}{\Sigma }}{\mathrm{aae}}_{ilhk}{x}_{ljk}{x}_{hjk}+{ϵ}_{ijk}\end{array}---\left(1\right)$

式(1)中，y_ijk是环境k中第i个性状的第j个株系的表型观测值；μ_i是性状i的群体 均值；b_ic是性状i第c个协变量的效应，这些协变量可以是离散或连续的，比如年份，地点，主 成分分析的一个或多个特征向量等，系数为x_jkc；a_il是性状i第l个QTS的加性效应,系数为 x_lj；aa_ilh是性状i第l个QTS与第h个QTS之间的加加上位性，系数为x_ljx_hj；e_ik是第i个性状在 环境k下的效应，是随机效应；ae_ilk是a_il与第k个环境的互作效应，系数为x_ljk；aae_ilhk是aa_ilh与第k个环境的互作效应，系数为x_ljkx_hjk；ε_ijk是性状i株系j在环境k下的随机残差。

上述模型(1)可改写成矩阵向量的表示形式：

$\begin{array}{c}Y=1\mu +X_C{B}_C+X_A{B}_A+X_{AA}{B}_{AA}+U_E{E}_E+\underset{l=1}{\overset{s}{\Sigma }}{U}_{A_{l}E}{E}_{A_{l}E}+\underset{l<h}{\overset{t}{\underset{l,h\in (1,2,...,s),}{\Sigma }}}{U}_{{\mathrm{AA}}_{lh}E}{E}_{{\mathrm{AA}}_{lh}E}+{\mathrm{IE}}_{ϵ}\\ =XB+\underset{u=1}{\overset{r+1}{\Sigma }}{U}_u{E}_u\end{array}---\left(2\right)$

式(2)中，Y是一个n×m阶性状观测值矩阵；μ＝(μ₁μ₂…μ_m)是一个1×m阶行向 量；$B_C=\left(\begin{array}{cccc}{b}_{C_1}& b_{C_2}& ...& b_{C_q}\end{array}\right)$是一个q×m阶向量，表示m个性状q个协变量的效应；X_C是Β_C的 系数矩阵；$B_A=\left(\begin{array}{cccc}{b}_{A_1}& b_{A_2}& ...& b_{A_m}\end{array}\right)$是一个s×m阶矩阵，表示m个性状s个QTS的加性效应矩 阵，其中$b_{A_k}={\left(\begin{array}{cccc}{a}_{k1}& a_{k2}& ...& a_{ks}\end{array}\right)}^T,$表示性状k显著QTS的加性效应向量；X_A是Β_A的系数矩阵； $B_{AA}=\left(\begin{array}{cccc}{b}_{{\mathrm{AA}}_1}& b_{{\mathrm{AA}}_2}& ...& b_{{\mathrm{AA}}_m}\end{array}\right)$是一个t×m阶矩阵，表示m个性状t对QTS的加加上位性效应矩 阵，其中$b_{{\mathrm{AA}}_k}={\left(\begin{array}{cccc}{\mathrm{aa}}_{k1}& {\mathrm{aa}}_{k2}& ...& {\mathrm{aa}}_{kt}\end{array}\right)}^T;$X_AA是B_AA的系数矩阵；环境效应矩阵E_E＝(e₁e₂… e_m)，其中e_k＝(e_k1e_k2…e_kp)^T是性状k的环境效应向量，随机效应，服从多元正态分布 是性状k的环境效应方差；U_E是E_E的系数；$E_{A_{l}E}=\left(\begin{array}{cccc}{e}_{A_{1l}E}& e_{A_{2l}E}& ...& e_{A_{ml}E}\end{array}\right)$是m 个性状第l个QTS的加性与环境互作效应矩阵，是随机效应，其中$e_{A_{kl}E}={\left(\begin{array}{cccc}{\mathrm{ae}}_{kl1}& {\mathrm{ae}}_{kl2}& ...& {\mathrm{ae}}_{klp}\end{array}\right)}^T;$是的系数；$E_{{\mathrm{AA}}_{lh}E}=\left(\begin{array}{cccc}{e}_{{\mathrm{AA}}_{1lh}E}& e_{{\mathrm{AA}}_{2lh}E}& ...& e_{{\mathrm{AA}}_{mlh}E}\end{array}\right)$是m个性状的加加上位性与环境互作 效应矩阵，随机效应，其中$e_{{\mathrm{AA}}_{klh}E}={\left(\begin{array}{cccc}{\mathrm{aae}}_{klh1}& {\mathrm{aae}}_{klh2}& ...& {\mathrm{aae}}_{klhp}\end{array}\right)}^T;$是的系数； $E_{ϵ}=\left(\begin{array}{cccc}{e}_{{ϵ}_1}& e_{{ϵ}_2}& ...& e_{{ϵ}_m}\end{array}\right)$是m个性状的随机残差效应矩阵；I是单位矩阵。

由于，上述模型中各QTS位置未知；所以，需要采用全基因组一维和二维扫描，根据 效应显著性检验的Lambda统计量分布特征，鉴别各显著的QTS，从而确定QTS的数目和位置， 在此基础上，构建上述全模型估算QTS的各项效应分量。

(二)全基因组一维扫描搜索单位点效应显著的QTS

首先，全基因组一维检测每个SNP标记的效应。当对位点l检验显著性时，可以用如 下的混合线性模型表示：

$y_{ijk}={\mu }_{ik}+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+a_{ilk}{x}_{lj}+{ϵ}_{ijk}---\left(3\right)$

式(3)中，y_ijk是环境k中第i个性状的第j个株系的表型观测值；μ_ik是环境k中性状i 的群体均值；b_ic是性状i第c个协变量的效应，系数为x_jkc；a_ilk是环境k中性状i位点l的加性 效应，系数为x_lj；ε_ijk是株系j在环境k下第i个性状的随机残差；

或者，采用矩阵向量进行表示：

Y＝W_QB_Q+W_MB_M+E_ε(4)

式(4)中，Y是一个n×m阶表型观测值矩阵；B_Q＝(a_1la_2l…a_ml)是p×m阶矩阵，由 检测位点l各环境下的加性效应组成，其中a_il＝(a_i1la_i2l…a_ipl)^T；Β_Μ是协变量矩阵，包 含了各个环境下的群体均值和各协变量的效应。Β_Q和Β_Μ的系数矩阵分别是W_Q和W_M；E_ε是m 个性状的随机残差矩阵。

其中，W＝[W_Q,W_M]，Q＝Y^TY-Y^TW(W^TW)⁺W^TY，Q₁＝Y^TY-Y^TW_M(W_M^TW_M)⁺W_M^TY，构建Wilks’ Lambda统计量λ＝|Q|/|Q₁|，在原假设H₀:B_Q＝0的情况下，统计量服从Lambda分布，即广义的 F分布。

对全基因组进行一维搜索后，计算每一个SNP位点的Lambda统计量，当某一SNP的 Lambda统计量值超过设定的阈值时，选取该SNP为显著的SNP位点，可采用置换检验方法确 定显著性的阈值。

(三)全基因组二维扫描搜索互作效应显著的成对QTS

由于，存在互作的QTS位点不一定存在单位点效应，所以，需要采用全基因组二维 扫描策略搜索所有可能的二互作上位性。在扫描中把单位点效应显著的QTS作为协变量放 入模型中，控制背景遗传效应，逐个扫描成对的上位性位点。

假设在一维扫描中搜索得到s个单位点效应显著的QTS，那么检验基因组中位点l 和h之间上位性的显著性可以用以下混合线性模型：

$y_{ijk}={\mu }_{ik}+{\mathrm{aa}}_{ilhk}{x}_{lj}{x}_{hj}+\underset{c}{\overset{q}{\Sigma }}{b}_{ic}{x}_{jkc}+\underset{r=1}{\overset{s}{\Sigma }}{a}_{irk}{x}_{rj}+{ϵ}_{ijk}---\left(5\right)$

式(5)中，y_ijk是环境k中第i个性状的第j个株系的表型观测值；μ_ik是环境k中性状i 的群体均值；aa_ilhk是性状i环境k下第l个QTS与第h个QTS之间的加加上位性效应；系数为 x_ljx_hj；b_ic是性状i第c个协变量的效应，系数为x_jkc；a_irk是在步骤(2)中检测到的显著的SNP 标记加性效应；系数为x_rj；ε_ijk是株系j在环境k下第i个性状的随机残差。

或者，采用矩阵向量进行表示：

Y＝W_QB_Q+W_MB_M+E_ε(6)

式(6)中，Y是一个n×m阶表型观测值矩阵；B_QQ＝(aa_1laa_2l…aa_ml)是p×m阶矩 阵，由检测位点l和位点h在各环境下的加加上位性效应组成，其中aa_ilh＝(aa_i1lhaa_i2lh… aa_iplh)^T；Β_Μ是协变量矩阵，包含了各个环境下的群体均值和步骤(2)中检测得到的显著的 SNP标记位点。Β_Q和Β_Μ的系数矩阵分别是W_Q和W_M；E_ε是m个性状的随机残差矩阵。

其中，W＝[W_Q,W_M]，Q＝Y^TY-Y^TW(W^TW)⁺W^TY，Q₁＝Y^TY-Y^TW_M(W_M^TW_M)⁺W_M^TY，构建Wilks’ Lambda统计量λ＝|Q|/|Q₁|，在原假设H₀:B_QQ＝0的情况下，统计量服从Lambda分布，即广义 的F分布。

对全基因组进行二维搜索后，计算检验位点h与位点l间的Lambda统计量，从而可 以获得Lambda统计量的分布曲线，采用置换检验法获得阈值，对每一超过设定阈值的统计 量点作为显著的二互作位点。

(四)遗传效应估算

针对所有鉴别出的效应显著的单位点SNP标记和二互作上位性效应显著的成对 SNP标记，采用向前选择法剔除假阳性的单位点SNP标记和二互作上位性SNP标记，得到单位 点QTS和二互作上位性QTS的应系数，将所述效应系数代入步骤(1)中构建全模型，获得最终 的多性状QTS分析全模型。

基于上述最终的全模型，采用蒙特卡罗马尔可夫链方法进行吉布斯抽样(SmithA FMandRobertsGO，1993，JournaloftheRoyalStatisticalSociety.SeriesB (Methodological)，p3-23)获取各项效应参数和方差参数的分布样本，用抽取的样本估算 QTS各项效应和方差，并进行显著性检验。

应用例1

(1)区域试验与DNA重测序

以一个衍生于超级稻协优9308的重组自交系(RIL)群体为研究对象。

试验采用随机区组试验设计，群体包含138个株系，2009年在海南省陵水县和浙江 省杭州市两个不同的地点进行试验，共276个观察样本，调查了三个品质性状。

利用BWA软件对DNA重测序数据进行序列连配，利用SAMtools进行SNP的挖掘，共有 701867个SNP用来后续的多性状全基因组关联分析。具体的材料和方法参照(Xuetal., 2015，BioMedResearchInternational)在QTS定位分析中，共对3个性状进行关联分析，直 链淀粉含量(AC)、糊化温度(GT)和凝胶稠度(GC)。这些性状都对水稻品质有十分重要的影 响，且相互间存在相关性。

(2)QTS作图软件及定位结果分析

本方法实施所采用的QTS作图软件是基于上述实施例所述方法所编写的 QTXNetwork软件；即：采用混合线性模型进行主效QTS以及上位性的定位，并估算主效QTS， 上位性以及环境互作的遗传参数。利用SNP标记、混合线性模型以及表型性状进行QTS定位 分析。

结果发现(如表1所示)，在第4、第6、第7和第8条染色体上分别检测到1个显著的 QTS。其中，Q6-44对三个性状同时起作用，说明具有一因多效，但是效应的方向和大小不一 致，而且对GT和GC存在基因与环境互作效应；Q7-73对GT和GC同时起作用，但是此位点对GC 的作用非常小。此外，还检测到一对上位性Q10-21/Q10-81，同时控制AC，GT和GC，且位点在 不同的环境对GT的效应大小不一样。

表1水稻数据多性状联合QTS定位的分析结果(性状(1)AC,(2)GT,(3)GC)

附：1、显著性水平，“^*”-p<0.05，“^**”-p<0.01；2、在结果中，SNPi/SNPj表示一个 二互作上位性结果。

最后，还需要特别注意的是，以上所举例子仅是本发明的具体实施例子。显然，本 发明不仅仅限于以上实施例子，还可以有许多变通的情况。本领域的技术人员从本发明公 开的内容直接推导出或联想到的所有变通情况，均认为是本发明的保护范围。