补充理解（2020-12-03）

补充一些理论知识，即为什么要采取本文所述的方式构建G矩阵。主要来源于这篇ssGBLUP综述文章，写的非常好。

如何理解，Z矩阵是对M矩阵进行了中心化。中心化即减掉均值，使总体均值为0。

先回顾一下，对M矩阵中心化的方法： \[Z = M-P\] M矩阵包括了每个样本、每个位点次等位基因的含量。其中每个元素都减掉了1，即M=M*-1，由0、1、2变为-1、0和1形式，便于计算。 \[P = 2(p_{i}-0.5)\] P矩阵包括了每个位点次等位基因的频率。这里减掉0.5，是因为M*矩阵每一个SNP位点等位基因含量减掉了1，所以M*每一个SNP位点等位基因含量均值即2倍的次等位基因频率$p_{i}$需要减掉1，由于每个位点的的等位基因含量为2个等位基因之和，也就是减掉2*0.5。

为什么这样做，就对M矩阵进行了中心化？

利用正文中的SNP位点1来做说明。3个样本包括2种基因型：AA、AA、TT。以T作为MAF等位基因，那么转换为M矩阵，则是0、0、2。

假定T等位基因的效应是a，发生频率是p，那么TT的频率是$p^{2}$，带有该基因型的个体的育种值u=2a；

AA基因型的频率为$q^{2}$，q=1-p，带有该基因型的个体的育种值u=0a=0；

TA基因型的频率为2pq，带有该基因型的个体的育种值u=1a；

该位点可解释的方差为$Var(u)=E(u^{2})-E(u)^{2}$；

在群体水平上，根据TT、TA和AA三种基因型的频率计算育种值u的平均值：$E(u)=2ap^{2}+0q^{2}+a2pq\\=2ap^{2}+2apq\\=2ap^{2}+2ap(1-p)=2pa$；

因此，$E(u)^{2}=(2pa)^{2}=4p^{2}a^{2}$；

根据基因型TT、TA、AA的频率，$E(u^{2})=p^{2}(2a)^{2}+2pq(1a)^{2}+q^{2}(0a)^{2}=4p^{2}a^{2}+2pqa^{2}$；

代入$E(u^{2})$和$E(u)^{2}$推导结果，方差$Var(u)=(4p^{2}a^{2}+2pqa^{2})-4p^{2}a^{2}=2pqa^{2}$。上述两个该公式，参考Leggara的文章中的Table 4。请注意，表中第3个基因型应为aa，PDF文件中错写为Aa。

两个体i、j的育种值可以表述为$m_{i}a、m_{j}a$，其中m表示MAF基因的含量；

减掉群体均值$2pa$，对其进行中心化：$u_{i}=m_{i}a-2pa=(m_{i}-2p)a=z_{i}a \\ u_{j}=m_{j}a-2pa=(m_{j}-2p)a=z_{j}a$；

这样基本就推导出了Z为什么是M-P这个公式；

两个体育种值之间的协方差：$Cov(u_{i},u_{j})=(z_{i}a-2pa)(z_{j}a-2pa)=(z_{i}-2p)(z_{j}-2p)\sigma_{a}^{2}$；

前述假定等位基因含量为0、1、2形式，进一步减掉1，则变为-1、0、 1，p=0；

因此，$Cov(u_{i},u_{j})=z_{i}z_{j}\sigma_{a}^{2}$；

因此，$u_{i}$和$u_{j}$之间的相关性表述为$r(u_{i},u_{j})=\frac{Cov(u_{i},u_{j})}{2p_{i}q_{j}\sigma_{a}^{2}}=\frac{z_{i}z_{j}\sigma_{a}^{2}}{2p_{i}q_{j}\sigma_{a}^{2}}=\frac{z_{i}z_{j}}{2p_{i}q_{j}}$。

从上边两个公式来看，G矩阵中元素表示的是个体间育种值的相关系数。

接下来推导G矩阵的计算公式

从1个SNP位点推广到多个SNP位点，假设所有位点的效应值都是一样的。那么$u=Za$，其中Z为中心化后M矩阵，a为等位基因的效应值。

育种值u方差为$Var(u)=Var(Za)=ZVar(a)Z'=ZZ'\sigma_{a}^{2}$；

根据方差公式，如果位点x和y间是相互独立的，那么Var(x+y)=Var(x)+Var(y)；

由于假设位点间是互不连锁的，因此n个独立SNP位点的方差为$Var(u)=\sigma_{u}^{2}=\sum_{i=1}^{n}2pq\sigma_{a}^{2}$；

因此$\sigma_{a}^{2}=\frac{\sigma_{u}^{2}}{\sum_{i=1}^{n}2pq}$；

$Var(u)=ZZ'\sigma_{a}^{2}=ZZ'\frac{\sigma_{u}^{2}}{\sum_{i=1}^{n}2pq}=\frac{ZZ'}{\sum_{i=1}^{n}2pq}\sigma_{u}^{2}$；

定义$G=\frac{ZZ'}{\sum_{i=1}^{n}2pq}$，那么$Var(u)=G\sigma_{u}^{2}$。

原文（2018-12-09）

自己的选择育种知识，又到了大规模系统学习和更新的时候了。

自从2005年开始学习水产动物选择育种知识以来，采取囫囵吞枣模式，虽然对于REML算法和BLUP算法的细节知识仍然一知半解，好歹相关的软件如ASReml、R等可以以较为熟练的使用，对于选择育种体系的理解也逐渐深入，常规的遗传评估工作算是可以胜任。

全基因组选择方法建立以来，尽管迅速在奶牛、猪等大动物上展现了优越的性能，然而对于虾等个体小、世代周期短的动物，如何应用，仍然是一个难题，也是我正在努力研究的方向。

然而，基于高通量SNP标记信息，获得的亲缘关系更加准确，这一结论已经毋容置疑。

因此，自这个blog开始，将逐渐学习并记录基于SNP标记的基因组选择评估方法，逐渐更新自己的知识体系。随着年龄的增大，知识体系的更新会慢很多，且行且记吧。

计算G矩阵的方法有很多种，目前引用较多的是Vanraden于2008年发表在Journal of Dairy Science上的文章“Efficient methods to compute genomic predictions”中所建立的方法，目前该论文已经被引用了超过1500次。

一些背景知识：

IBD (Identical by descent)，血缘同源。

IBS (Identical by state)，状态同源。

通常我们会讲，如果一段DNA在两个或者多个个体中均有一致的核苷酸序列，那么可以定义该DNA片段属性为状态同源，即IBS。

在两个或者多个个体中的一个IBS片段，如果遗传自一个共同的祖先个体，且没有发生重组，那么该片段也是血缘同源的，即IBD。

一段DNA，如果是IBD，那么肯定也是IBS；如果片段不是IBD，也有可能是IBS，因为不同个体中发生了相同的突变，或者重组没有改变片段等原因所致。

关于IBD和IBS的定义，参见WIKIPEDIA。

个体间的亲缘关系(relationship)，可以通过以下两种方法估计：

通过系谱进行计算，利用处于IBD状态的等位基因估计个体间的亲缘关系

通过分子标记，利用处于IBS状态的基因估计个体间共享染色体片段的百分比，估计个体间的相似性

下边的SNP例子，主要来自Fikret Isik的一个Presentation： “Genomic Relationships and GBLUP”，但幻灯片中P、Z矩阵的计算结果是错误的。

1. 构建M矩阵

M矩阵是一个n×m矩阵，n表示个体数（the number of individuals），m表示位点数（the number of loci）。

譬如有3个二倍体个体，4个SNP位点。对于每个位点，小写等位基因表示其出现频率最低。

rm(list=ls()) # 加载synbreed包 is_installed <- require(xtable)

## Loading required package: xtable

if (!is_installed) { install.packages("xtable") library(xtable) } genotypes <- data.frame( snp1 = c("AA", "AA", "tt"), snp2 = c("Ct", "Ct", "CC"), snp3 = c("GG", "Ga", "GG"), snp4 = c("Ag", "AA", "AA"), row.names = c("Ind1", "Ind2", "Ind3"), stringsAsFactors = FALSE ) print(xtable(genotypes),type="html")

snp1 snp2 snp3 snp4

Ind1 AA Ct GG Ag

Ind2 AA Ct Ga AA

Ind3 tt CC GG AA

把基因型矩阵转换为基因含量矩阵，具体参数用基因型中次等位基因的个数表示，即0、1、2。因此我们称这个矩阵为MAF矩阵。

转换用到了synbreed包中的codeGeno函数。

# 加载synbreed包 is_installed <- require(synbreed)

## Loading required package: synbreed

if (!is_installed) { install.packages("synbreed") library(synbreed) } # 把小写字母转换为大写字母 genotypes <- apply(genotypes, 2, function(x) toupper(x)) # 创建codeGeno函数需要的数据类型 genotypes_gp <- create.gpData(geno = genotypes) # codeGeno默认计算MAF等位基因的个数 M <- codeGeno(genotypes_gp,reference.allele = "minor")$geno M

## snp1 snp2 snp3 snp4 ## Ind1 0 1 0 1 ## Ind2 0 1 1 0 ## Ind3 2 0 0 0

为了便于计算，通常会将M矩阵减1，每个位点等位基因含量变为-1、0、1三种形式。减1的目的也是为了将每列的均值标准化为0。如果三种等位基因型出现的频率是相同的，那么-1后，每列的均值就是0。

M <- M - 1 M

## snp1 snp2 snp3 snp4 ## Ind1 -1 0 -1 0 ## Ind2 -1 0 0 -1 ## Ind3 1 -1 -1 -1

2. MM’矩阵

M %*% t(M)

## Ind1 Ind2 Ind3 ## Ind1 2 1 0 ## Ind2 1 2 0 ## Ind3 0 0 4

MM’对角线数字为2、2、4，分别表示三个个体基因型为纯合的位点数。譬如Ind1，四个位点中有两个位点的基因型为纯合类型（AA和GG），而Ind3个体四个位点基因型全部为纯合类型。

非对角线元素，表示亲缘个体间共享的等位基因数。原英文：off-diagonals measure the number of alleles shared by relatives，不是很理解。

3. M’M矩阵

t(M) %*% M

## snp1 snp2 snp3 snp4 ## snp1 3 -1 0 0 ## snp2 -1 1 1 1 ## snp3 0 1 2 1 ## snp4 0 1 1 2

M’M对角线数字为3、1、2、2，分别表示四个位点纯合基因型的个体数。譬如位点snp1，有3个个体基因型全部为纯合；位点snp2，有1个个体基因型为纯合。

非对角线元素表示什么意思，还不是很清楚。原英文：off-diagonals measure the number of times alleles at different loci were inherited by the same individual.

4. 构建P矩阵

根据下述公式构建P矩阵的每一列：

\[P_{i} = 2(p_{i}-0.5)\]

其中i表示位点编号，$p_{i}$表示每个位点次等位基因频率。P矩阵中反应的是每个snp位点（每列）两个等位基因频率偏离均值0.5的程度。以snp1，次等位基因为t，频率为p=2/6=0.333333；那么数字化后=c(0, 0, 2), 再减1后为c(-1, -1, 1)，相加为-1而不是0，如何标准化为0？因此这里需要用到偏离矩阵，P=2(0.333333-0.5)= -0.333334。在下一部分构建Z矩阵时, M-P = -1 - 3*(-0.333334) ≈ 0。

#计算每个位点的MAF. p_lower <- (apply(M,2,sum)+nrow(M))/(nrow(M)*2) p_upper <- 2*(p_lower-0.5) p_matrix <- matrix(p_upper,byrow=T,nrow=nrow(M),ncol=ncol(M)) P <- p_matrix P

## [,1] [,2] [,3] [,4] ## [1,] -0.3333333 -0.3333333 -0.6666667 -0.6666667 ## [2,] -0.3333333 -0.3333333 -0.6666667 -0.6666667 ## [3,] -0.3333333 -0.3333333 -0.6666667 -0.6666667

每一列，即每个位点的次等位基因频率都是相同的数字。

5. 构建Z矩阵

Z矩阵实际上是中心化后的M矩阵。

原文中M-P的作用：

Sets means values of the allele effects to 0. 从下边Z矩阵可以看出，每一列SNP位点，等位基因效应的和为0。

Gives more credit to rare alleles than to common alleles when calculating genomic relationships. 意思是以MAF基因频率去进行中心化？

Allele frequencies in P should be from the unselected base population rather than those that occur after selection or inbreeding. An earlier or later base population can lead to greater or fewer relationships and to more or less inbreeding.

The genomic inbreeding coefficient is greater if the individual is homozygous for rare alleles than if homozygous for common alleles.

Z <- M - P Z

## snp1 snp2 snp3 snp4 ## Ind1 -0.6666667 0.3333333 -0.3333333 0.6666667 ## Ind2 -0.6666667 0.3333333 0.6666667 -0.3333333 ## Ind3 1.3333333 -0.6666667 -0.3333333 -0.3333333

6. 构建G矩阵

构建方法，在文章中提到了3种：

6.1 第一种方法：

\[G = \frac{ZZ'}{2\sum{p_{i}(1-p_{i})}}\]

Z为SNP标记的设计矩阵；公式中，$p_{i}$为每个位点的次等位基因频率。公式中分母部分，是为了使G矩阵与A矩阵的尺度相一致。$\sum{p_{i}(1-p_{i})}$ scales G to be analogous to the numerator relationship matrix A.

d <- 2*sum(p_lower*(1 - p_lower)) G <- Z %*% t(Z) / d print(xtable(G),type="html")

Ind1 Ind2 Ind3

Ind1 0.77 0.08 -0.85

Ind2 0.08 0.77 -0.85

Ind3 -0.85 -0.85 1.69

邓飞在他的科学网博客，利用sommer包对上述结果进行了验证。在此，同样利用sommer包的A.mat()函数进行验证。可以看出，计算结果是一致的。sommer 包现在核心函数已经用C++重写，性能可期待。

suppressPackageStartupMessages(is_installed <- require(sommer)) if (!is_installed) { install.packages("sommer") library(sommer) } print(xtable(A.mat(M)),type="html")

Ind1 Ind2 Ind3

Ind1 0.77 0.08 -0.85

Ind2 0.08 0.77 -0.85

Ind3 -0.85 -0.85 1.69

个体j的基因组近交系数等于$G_{jj}-1$。

diag(G) - 1

## Ind1 Ind2 Ind3 ## -0.2307692 -0.2307692 0.6923077

Ind1和Ind2个体的近交系数为负值，这表明基因组近交系数跟传统的基于系谱的近交系数并不在一个尺度上。

个体间的基因组亲缘关系，譬如j、k个体间的亲缘关系，= $G_{jk}/\sqrt{G_{jj}G_{kk}}$。

现在计算Ind1和Ind2个体间的亲缘关系，= $G_{12}/\sqrt{G_{11}G_{22}}$ = 0.13。这个系数，跟经典的Relationship coefficient of Wright （1922）非常相似。

由于对角线元素一般都在1左右，因此j与k的基因组亲缘关系（Genomic relationship）也就基本约等于G矩阵中相对应的值。

如果把G阵与传统A阵（Additive genetic relationship matrix 或 numerator relationship matrix）类比，那么j与k的共亲系数（coancestry coefficient）同样是基因组亲缘系数的一半。

6.2 第二种方法

\[G = ZDZ'\]

D是一个对角矩阵，其中$D_{ii}$等于 \[D_{ii} = \frac{1}{m[2p_{i}(1-p_{i})]}\]

第二种方法是第一种方法的一个变体，利用每个位点期望方差的倒数，对标记效应值进行加权。

D_ii <- 1/(ncol(M)*(2*p_lower*(1 - p_lower))) # 构建对角矩阵 D <- diag(D_ii,ncol(M)) G <- Z %*% D %*% t(Z) print(xtable(G),type="html")

Ind1 Ind2 Ind3

Ind1 0.81 -0.09 -0.72

Ind2 -0.09 0.81 -0.72

Ind3 -0.72 -0.72 1.45

参考文献：

Vanraden P M. Efficient methods to compute genomic predictions.[J]. Journal of Dairy Science, 2008, 91(11):4414-4423.

Wright, S. 1922. Coefficients of inbreeding and relationship. Am.Nat. 56:330–338.

如何构建G矩阵-基因组亲缘关系矩阵（Genomic relationships matrix）

补充理解（2020-12-03）

原文（2018-12-09）

1. 构建M矩阵

2. MM’矩阵

3. M’M矩阵

4. 构建P矩阵

5. 构建Z矩阵

6. 构建G矩阵

6.1 第一种方法：

6.2 第二种方法