博文
全基因组选择介绍及实践-2：构建H矩阵
||
1, 编者自语
H矩阵作为一步法的入门技术, 是需要掌握的, 本文以一篇文献为例, 介绍如何从头构建H矩阵. 文章包括H矩阵推导过程和代码实现.
2, H矩阵定义
基因组选择中, GBLUP的一个挑战是, 在参考群构建时, 需要两步, 第一步根据系谱和表型数据, 计算出伪数据(pseudo-data)(比如, 根据系谱计算公牛的女儿产奶偏差作为表型值, 因为公牛没有产奶数据), 然后用基因组信息进行评估建模, 这就造成信息的损耗和偏离. 解决的方法是, 可以通过一种手段, 将系谱关系A矩阵和基因组信息构建的亲缘关系G矩阵合并为H矩阵, 这样就成了一步法(Single-setp).
As not all animals can be genotyped, a 2- or 3-step procedure has to be followed; first, a regular genetic evaluation is run; then, corrected phenotypes or pseudo-data are used in the second step, where the marker-assisted selection model is effectively applied (Guillaume et al., 2008; VanRaden et al., 2009). These phenotypes are daughter yield deviations (DYD) and yield deviations (YD) for dairy cattle.
3, H矩阵推导过程—三组进行推导
假设A矩阵, 包括三部分:
A is the numerator relationship matrix based on pedigree. Consider three types of animals in u: 1) ungenotyped ancestors with breeding values u1; 2) genotyped animals, with breeding values u2 (no ancestor is genotyped and phantom parents can be generated if necessary); and 3) ungenotyped animals with breeding values u3, which might descend from either one of the three types of animals.
编号1 是没有测序的个体, 祖先
编号2 是测序的个体(当代和后代)
编号3 是没有测序个体(当代和后代)
如果所有个体A矩阵按照世代排, 那么可以将其分为如下部分:
那么3组的方差, 3组和1组, 3组和2组的协方差为:
那么矩阵变为:
剖出A矩阵
4, H矩阵推导过程—二组进行推导
那么A矩阵和A逆矩阵可以剖分为:
编号1为非测序个体
编号2为测序个体
测序个体和非测序个体的方差以及协方差:
假定H矩阵为所有个体的矩阵:
那么H矩阵可以分为:
5, H逆矩阵推导过程
H矩阵还可以写为:
对其进行求逆:
6, 系谱数据构建A矩阵
R语言生成系谱
ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) ped_full
结果:
> ped_full ID Sire Dam 1 9 1 2 2 10 3 4 3 11 5 6 4 12 7 8 5 13 9 10 6 14 11 12 7 15 11 4 8 16 13 15 9 17 13 14
利用nadiv包计算A矩阵
因为nadiv中的makeA函数, 计算的A矩阵行名没有安装1~17编号, 需要生成一个id_r的编号对矩阵进行排序, 排序后的矩阵按照1 ~ 17编号.
ped = nadiv::prepPed(ped_full) A_mat = as.matrix(nadiv::makeA(ped)) id = row.names(A_mat) id_r = match(1:17,id) A_mat_sort = A_mat[id_r,id_r] Matrix::Matrix(A_mat_sort, sparse=TRUE) # 查看稀疏矩阵
结果:
> Matrix::Matrix(A_mat_sort, sparse=TRUE) 17 x 17 sparse Matrix of class "dsCMatrix" [[ suppressing 17 column names ‘1’, ‘2’, ‘3’ ... ]] 1 1.000 . . . . . . . 0.50 . . . 0.2500 . . 0.12500 0.12500 2 . 1.000 . . . . . . 0.50 . . . 0.2500 . . 0.12500 0.12500 3 . . 1.000 . . . . . . 0.500 . . 0.2500 . . 0.12500 0.12500 4 . . . 1.000 . . . . . 0.500 . . 0.2500 . 0.5000 0.37500 0.12500 5 . . . . 1.000 . . . . . 0.50 . . 0.250 0.2500 0.12500 0.12500 6 . . . . . 1.000 . . . . 0.50 . . 0.250 0.2500 0.12500 0.12500 7 . . . . . . 1.000 . . . . 0.50 . 0.250 . . 0.12500 8 . . . . . . . 1.000 . . . 0.50 . 0.250 . . 0.12500 9 0.500 0.500 . . . . . . 1.00 . . . 0.5000 . . 0.25000 0.25000 10 . . 0.500 0.500 . . . . . 1.000 . . 0.5000 . 0.2500 0.37500 0.25000 11 . . . . 0.500 0.500 . . . . 1.00 . . 0.500 0.5000 0.25000 0.25000 12 . . . . . . 0.500 0.500 . . . 1.00 . 0.500 . . 0.25000 13 0.250 0.250 0.250 0.250 . . . . 0.50 0.500 . . 1.0000 . 0.1250 0.56250 0.50000 14 . . . . 0.250 0.250 0.250 0.250 . . 0.50 0.50 . 1.000 0.2500 0.12500 0.50000 15 . . . 0.500 0.250 0.250 . . . 0.250 0.50 . 0.1250 0.250 1.0000 0.56250 0.18750 16 0.125 0.125 0.125 0.375 0.125 0.125 . . 0.25 0.375 0.25 . 0.5625 0.125 0.5625 1.06250 0.34375 17 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.25 0.250 0.25 0.25 0.5000 0.500 0.1875 0.34375 1.00000
和文章结果一致:
7, 基因组信息G矩阵
令个体9~12为测序个体, G矩阵为对角线为1, 非对角线为0.7的矩阵, 有行名和列名.
代码
G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 G
结果
> G 9 10 11 12 9 1.0 0.7 0.7 0.7 10 0.7 1.0 0.7 0.7 11 0.7 0.7 1.0 0.7 12 0.7 0.7 0.7 1.0
8, 根据公式计算H矩阵
# 提取子集 id_g = 9:12 id_ng = c(1:8,13:17) A22 = A[id_g,id_g] A12 = A[id_ng,id_g] A21 = A[id_g,id_ng] A22 = A[id_g,id_g] iA22 = solve(A22)
根据公式, 计算H11, H12, H21, H22
# 构建H矩阵 H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21 H12 = A12 %*% iA22 %*%G H21 = G %*% iA22 %*% A21 H22 = G H = cbind(rbind(H11,H21),rbind(H12,H22)) round(H,4)
结果
> round(H,2) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1 1.00 0.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 2 0.00 1.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 3 0.18 0.18 1.00 0.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.18 0.30 0.39 4 0.18 0.18 0.00 1.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.68 0.55 0.39 5 0.18 0.18 0.18 0.18 1.00 0.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 6 0.18 0.18 0.18 0.18 0.00 1.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 7 0.18 0.18 0.18 0.18 0.18 0.18 1.00 0.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 8 0.18 0.18 0.18 0.18 0.18 0.18 0.00 1.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 9 0.50 0.50 0.35 0.35 0.35 0.35 0.35 0.35 1.00 0.70 0.70 0.70 0.85 0.70 0.52 0.69 0.77 10 0.35 0.35 0.50 0.50 0.35 0.35 0.35 0.35 0.70 1.00 0.70 0.70 0.85 0.70 0.60 0.73 0.77 11 0.35 0.35 0.35 0.35 0.50 0.50 0.35 0.35 0.70 0.70 1.00 0.70 0.70 0.85 0.68 0.69 0.77 12 0.35 0.35 0.35 0.35 0.35 0.35 0.50 0.50 0.70 0.70 0.70 1.00 0.70 0.85 0.52 0.61 0.77 13 0.43 0.43 0.43 0.43 0.35 0.35 0.35 0.35 0.85 0.85 0.70 0.70 1.35 0.70 0.56 0.96 1.02 14 0.35 0.35 0.35 0.35 0.43 0.43 0.43 0.43 0.70 0.70 0.85 0.85 0.70 1.35 0.60 0.65 1.02 15 0.26 0.26 0.18 0.68 0.34 0.34 0.26 0.26 0.52 0.60 0.68 0.52 0.56 0.60 1.18 0.87 0.58 16 0.34 0.34 0.30 0.55 0.34 0.34 0.31 0.31 0.69 0.73 0.69 0.61 0.96 0.65 0.87 1.41 0.80 17 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.77 0.77 0.77 0.77 1.02 1.02 0.58 0.80 1.52
文章结果:
9, 代码汇总
# 构建系谱 ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) ped_full # 计算A矩阵 ped = nadiv::prepPed(ped_full) A_mat = as.matrix(nadiv::makeA(ped)) id = row.names(A_mat) id_r = match(1:17,id) A_mat_sort = A_mat[id_r,id_r] A = A_mat_sort Matrix::Matrix(A_mat_sort, sparse=TRUE) # 构建G矩阵 G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 G # 提取子集 id_g = 9:12 id_ng = c(1:8,13:17) A11 = A[id_ng,id_ng] A22 = A[id_g,id_g] A12 = A[id_ng,id_g] A21 = A[id_g,id_ng] A22 = A[id_g,id_g] iA22 = solve(A22) # 构建H矩阵 H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21 H12 = A12 %*% iA22 %*%G H21 = G %*% iA22 %*% A21 H22 = G H = cbind(rbind(H11,H21),rbind(H12,H22)) id = rownames(H) id_r = match(1:17,id) H = H[id_r,id_r] round(H,2)
10, 编写一个函数
这里编写一个H_matrix函数, 参数有:
G矩阵, 包括行名和列名
ped, 系谱文件
id_g, 是基因型id
id_full, 是所有个体的id
H_matrix = function(G = G, ped = ped,id_g, id_full){
ped = nadiv::prepPed(ped_full)
A_mat = as.matrix(nadiv::makeA(ped))
A = A_mat[id_full,id_full]
id_ng = setdiff(id_full,id_g)
A11 = A[id_ng,id_ng]
A22 = A[id_g,id_g]
A12 = A[id_ng,id_g]
A21 = A[id_g,id_ng]
A22 = A[id_g,id_g]
iA22 = solve(A22)
# 构建H矩阵
H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21
H12 = A12 %*% iA22 %*%G
H21 = G %*% iA22 %*% A21
H22 = G
H = cbind(rbind(H11,H21),rbind(H12,H22))
id = rownames(H)
id_r = match(1:17,id)
H = H[id_r,id_r]
return(H)
}测试
# pedigree ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) # genotype G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 # id_g id_g = 9:12 # id_full id_full = 1:17 H = H_matrix(G,ped_full,id_g,id_full) round(H,2)
结果
> round(H,2) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1 1.00 0.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 2 0.00 1.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 3 0.18 0.18 1.00 0.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.18 0.30 0.39 4 0.18 0.18 0.00 1.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.68 0.55 0.39 5 0.18 0.18 0.18 0.18 1.00 0.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 6 0.18 0.18 0.18 0.18 0.00 1.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 7 0.18 0.18 0.18 0.18 0.18 0.18 1.00 0.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 8 0.18 0.18 0.18 0.18 0.18 0.18 0.00 1.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 9 0.50 0.50 0.35 0.35 0.35 0.35 0.35 0.35 1.00 0.70 0.70 0.70 0.85 0.70 0.52 0.69 0.77 10 0.35 0.35 0.50 0.50 0.35 0.35 0.35 0.35 0.70 1.00 0.70 0.70 0.85 0.70 0.60 0.73 0.77 11 0.35 0.35 0.35 0.35 0.50 0.50 0.35 0.35 0.70 0.70 1.00 0.70 0.70 0.85 0.68 0.69 0.77 12 0.35 0.35 0.35 0.35 0.35 0.35 0.50 0.50 0.70 0.70 0.70 1.00 0.70 0.85 0.52 0.61 0.77 13 0.43 0.43 0.43 0.43 0.35 0.35 0.35 0.35 0.85 0.85 0.70 0.70 1.35 0.70 0.56 0.96 1.02 14 0.35 0.35 0.35 0.35 0.43 0.43 0.43 0.43 0.70 0.70 0.85 0.85 0.70 1.35 0.60 0.65 1.02 15 0.26 0.26 0.18 0.68 0.34 0.34 0.26 0.26 0.52 0.60 0.68 0.52 0.56 0.60 1.18 0.87 0.58 16 0.34 0.34 0.30 0.55 0.34 0.34 0.31 0.31 0.69 0.73 0.69 0.61 0.96 0.65 0.87 1.41 0.80 17 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.77 0.77 0.77 0.77 1.02 1.02 0.58 0.80 1.52
搞定.
公众号ID：R-breeding
11, 数据及结果来源
paper:Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information.[J]. Journal of Dairy Science, 2009, 92(9):4656-63.
https://blog.sciencenet.cn/blog-2577109-1165174.html
上一篇：全基因组选择介绍及实践-1
下一篇：育种4.0世代的到来个人应该准备什么
