基于GBLUP/ssGBLUP的育种值计算 SNP 效应。

估计SNP效应

一步法中的 \(\mathbf{H}^{-1}\) 公式为 \[ \mathbf{H}^{-1}=\mathbf{A}^{-1}+\left[\begin{array}{lc} 0 & 0 \\ 0 & \mathbf{G}^{-1}-\mathbf{A}_{22}^{-1} \end{array}\right] \] 其中 \(\mathbf{G}\) 是校正到与 \(\mathbf{A}_{22}\) 矩阵规模相同后的矩阵。

我们将个体育种值分为基因型个体 (\(\mathbf{a_{g}}\)) 和 非基因型个体 (\(\mathbf{a_{n}}\)) ， 基因型个体的育种值可以视为所有 snp 效应之和，即： \[ \mathbf{a_{g}} = \mathbf{Z}\mathbf{u} \] 这里 \(\mathbf{Z}\) 是 SNP 效应的系数矩阵（中心化的系数矩阵），\(\mathbf{u}\) 是 SNP 效应向量。

因此，个体加性效应的方差为 \[ \operatorname{var}(\mathbf{a_{g}}) = \operatorname{var}(\mathbf{Z}\mathbf{u}) = \mathbf{Z} \mathbf{D} \mathbf{Z}^{\prime} \sigma_{u}^{2} = \mathbf{G}^{*}\sigma_{a}^{2} \] 这里 \(\mathbf{D}\) 是位点的权重的对角矩阵（GBLUP 中为 \(\mathbf{I}\) ），\(\sigma_{u}^{2}\) 是位点的加性方差，\(\mathbf{G}^{*}\) 是有权重的基因组亲缘关系矩阵。

基因型个体的育种值 \(\mathbf{a_{g}}\) 和 SNP 效应 \(\mathbf{u}\) 的联合协方差矩阵为 \[ \operatorname{var}\left[\begin{array}{lc} \mathbf{a_{g}} \\ \mathbf{u} \end{array}\right] = \left[\begin{array}{lc} \mathbf{ZDZ^{\prime}} & \mathbf{ZD^{\prime}} \\ \mathbf{DZ^{\prime}} & \mathbf{D} \end{array}\right] \sigma_{u}^{2} \] 因此 \[ \mathbf{G}^{*} = \frac{\operatorname{var}(\mathbf{a_{g}})} {\sigma_{a}^{2}} = \frac{\operatorname{var}(\mathbf{Zu})} {\sigma_{a}^{2}} = \mathbf{ZDZ^{\prime}} \frac{\sigma_{u}^{2}}{\sigma_{a}^{2}} = \mathbf{ZDZ^{\prime}}\lambda \] 其中 \(\lambda\) 是一个方差比值，基于 VanRaden 2009 ，存在公式 \[ \lambda = \frac{\sigma_{u}^{2}}{\sigma_{a}^{2}} = \frac{1}{\sum_{i=1}^{M}2p_{i}(1-p_{i})} \] 这里 \(M\) 是位点数目，\(p_{i}\) 是等位基因频率。

基于 Stranden & Garrick (2009) ，可以推导出公式 \[ \hat{\mathbf{u}} = \lambda \mathbf{DZ^{\prime}G^{*-1}\mathbf{\hat{a}_{g}}} = \mathbf{DZ^{\prime}(ZDZ^{\prime})^{-1}\mathbf{\hat{a}_{g}}} \] 我们可以检查一下， 没有问题 \[ \mathbf{\hat{a}_{g}} = \mathbf{Zu} =\mathbf{ZDZ^{\prime}(ZDZ^{\prime})^{-1}\mathbf{\hat{a}_{g}}} = \mathbf{\hat{a}_{g}} \]

估计SNP效应的方差

每一个 SNP 效应的方差估计值为 （Zhang, 2010） \[ \hat{\sigma}_{u,i}^{2} = \hat{u}_{i}^{2}2p_{i}(1-p_{i}) \] 根据上面的公式可以构建一个估计一步法的 \(\mathbf{D}\) 矩阵的算法，定义 \(t\) 是迭代次数，\(i\) 第 \(i\) 个 SNP ，算法流程如下（这也是 blupf90 进行 ssGWAS 的流程）：

1. \(t=0\) ，\(\mathbf{D}_{(t)} = \mathbf{I}\) ， \(\mathbf{G}^{*}_{(t)} = \mathbf{Z} \mathbf{D}_{(t)} \mathbf{Z}^{\prime} \lambda\)
2. 使用 ssGBLUP 计算 \(\mathbf{\hat{a}_{g}}\)
3. 计算 \(\hat{\mathbf{u}_{(t)}} = \lambda \mathbf{D_{(t)}Z^{\prime}G^{*-1}_{(t)}\mathbf{\hat{a}_{g}}}\)
4. 对每一个位点，计算 \(d^{*}_{i(t+1)} = \hat{u}_{i(t)}^{2}2p_{i}(1-p_{i})\) ，得到 \(\mathbf{D}^{*}_{(t+1)}\) 矩阵
5. 标准化 \(\mathbf{D}_{(t+1)} = \frac{\operatorname{tr} (\mathbf{D}_{(0)})}{\operatorname{tr} (\mathbf{D}^{*}_{(t+1)})} \mathbf{D}^{*}_{(t+1)}\)
6. 计算 \(\mathbf{G}^{*}_{(t+1)} = \mathbf{ZD_{(t+1)}Z^{\prime}} \lambda\)
7. \(t = t+1\)
8. 退出程序，或循环跳到第二步或第三步。

如果循环跳到第三步，则每一次计算均使用相同的 \(\mathbf{\hat{a}_{g}}\) 。

具体是跳到第二步或第三步 ，还有循环次数均需要根据具体情况确定。

参考文献

1. Wang H, Misztal I, Aguilar I, et al. Genome-wide association mapping including phenotypes from relatives without genotypes[J]. Genetics research, 2012, 94(2): 73-83.
2. VanRaden P M, Van Tassell C P, Wiggans G R, et al. Invited review: Reliability of genomic predictions for North American Holstein bulls[J]. Journal of dairy science, 2009, 92(1): 16-24.
3. Strandén I, Garrick D J. Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[J]. Journal of dairy science, 2009, 92(6): 2971-2975.
4. Zhang Z, Liu J, Ding X, et al. Best linear unbiased prediction of genomic breeding values using a trait-specific marker-derived relationship matrix[J]. PloS one, 2010, 5(9): e12648.