# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.3 p18
앞에서는 SNP effect를 구한 후, 다시 DGV를 구하였는데 여기서는 GRM(Genomic Relationship Matrix)를 이용하여 DGV를 바로 구하는 방법에 관하여 설명을 한다.
SNP genotype을 이용하여 NRM(Numerator Relationship Matrix)와 유사한 행렬을 구할 수 있다. 이 행렬에 약 0.90 또는 0.95 또는 0.99의 가중치를 주고 나머지를 NRM에 주어 improved non-singular matrix Gwt를 구할 수 있다.
Data
animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
13 0 0 1 558 9 2 0 1 1 0 0 0 2 1 2
14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0
Pedigree
animal sire dam
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 13 4
16 15 2
17 15 5
18 14 6
19 14 9
20 14 9
21 1 3
22 14 8
23 14 11
24 14 10
25 14 7
26 14 12
R Program
# Mixed Linear Model(GBLUP Model) for Computing SNP Effect
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.2 p183
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function to make numerator relationship matrix
makenrm = function(pedi) {
n = nrow(pedi)
nrm = matrix(c(0), nrow = n, ncol = n)
for (i in 1:n) {
animal = pedi[i, 1]
sire = pedi[i, 2]
dam = pedi[i, 3]
if (sire == 0 & dam == 0) {
# 부모를 모를 경우
nrm[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# 부만 알고 있을 경우
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, sire])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1
} else if (sire == 0 & dam != 0) {
# 모만 알고 있을 경우
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, dam])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1
} else {
# 부와 모를 알고 있을 경우
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
}
}
return(nrm)
}
# set working directory
setwd(".")
# print working directory
getwd()
# read data
data = read.table("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data
# read pedigree : animal sire dam
pedi = read.table("snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
pedi
# variance component and ratio
sigma_a = 35.241
sigma_e = 245
sigma_a
sigma_e
# SNP
M_all = as.matrix(data[, 7:16])
M_all
# mean of each SNP(allele frequency)
snp_mean = colMeans(M_all) / 2
snp_mean
# sum of 2pq(heterozygote ratio)
sum2pq = 0
for (i in 1:10) {
sum2pq = sum2pq + 2 * snp_mean[i] * (1 - snp_mean[i])
}
sum2pq
# alpha
alpha = sigma_e / sigma_a
alpha
# observation
y = data[1:8, 6]
y
# design matrix of fixed effect(general mean)
X = matrix(rep(1, 8), 8, 1)
X
# design matrix for DGV
W = cbind(diag(8), matrix(c(0), 8, 6))
W
# genomic relationship matrix, G
Z = sweep(M_all, 2, snp_mean * 2)
Z
G = Z %*% t(Z) / sum2pq
round(G, 3)
# Numerater Relationship Matrix
A = makenrm(pedi)
A
A22 = A[13:26, 13:26]
A22
# weighted G
wt = 0.99
Gw = wt * G + (1 - wt) * A22
round(Gw, 3)
Gwi = ginv(Gw)
round(Gwi, 3)
# Left Hand Side
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% W),
cbind(t(W) %*% X, t(W) %*% W + Gwi * alpha))
round(LHS, 3)
# Right Hand Side
RHS = rbind(t(X) %*% y,
t(W) %*% y)
round(RHS, 3)
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)
# solution
sol = gi_LHS %*% RHS
round(sol, 3)
RStuio에서 실행한 화면
프로그램 실행 결과
> # Mixed Linear Model(GBLUP Model) for Computing SNP Effect
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
>
> # Raphael Mrode
>
> # Example 11.2 p183
>
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # function to make numerator relationship matrix
> makenrm = function(pedi) {
+ n = nrow(pedi)
+ nrm = matrix(c(0), nrow = n, ncol = n)
+
+ for (i in 1:n) {
+ animal = pedi[i, 1]
+ sire = pedi[i, 2]
+ dam = pedi[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 부모를 모를 경우
+ nrm[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # 부만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, sire])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1
+ } else if (sire == 0 & dam != 0) {
+ # 모만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, dam])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1
+ } else {
+ # 부와 모를 알고 있을 경우
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
+ }
+ }
+ return(nrm)
+ }
>
> # set working directory
> setwd(".")
>
> # print working directory
> getwd()
[1] "D:/users/bhpark/2020/job/공부하기/07_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/19_mixed linear model GBLUP for computing SNP effects"
>
> # read data
> data = read.table("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> data
animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
1 13 0 0 1 558 9.0 2 0 1 1 0 0 0 2 1 2
2 14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
3 15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
4 16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
5 17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
6 18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
7 19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
8 20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
9 21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
10 22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
11 23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
12 24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
13 25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
14 26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0
>
> # read pedigree : animal sire dam
> pedi = read.table("snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> pedi
animal sire dam
1 1 0 0
2 2 0 0
3 3 0 0
4 4 0 0
5 5 0 0
6 6 0 0
7 7 0 0
8 8 0 0
9 9 0 0
10 10 0 0
11 11 0 0
12 12 0 0
13 13 0 0
14 14 0 0
15 15 13 4
16 16 15 2
17 17 15 5
18 18 14 6
19 19 14 9
20 20 14 9
21 21 1 3
22 22 14 8
23 23 14 11
24 24 14 10
25 25 14 7
26 26 14 12
>
> # variance component and ratio
> sigma_a = 35.241
> sigma_e = 245
> sigma_a
[1] 35.241
> sigma_e
[1] 245
>
> # SNP
> M_all = as.matrix(data[, 7:16])
> M_all
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
[1,] 2 0 1 1 0 0 0 2 1 2
[2,] 1 0 0 0 0 2 0 2 1 0
[3,] 1 1 2 1 1 0 0 2 1 2
[4,] 0 0 2 1 0 1 0 2 2 1
[5,] 0 1 1 2 0 0 0 2 1 2
[6,] 1 1 0 1 0 2 0 2 2 1
[7,] 0 0 1 1 0 2 0 2 2 0
[8,] 0 1 1 0 0 1 0 2 2 0
[9,] 2 0 0 0 0 1 2 2 1 2
[10,] 0 0 0 1 1 2 0 2 0 0
[11,] 0 1 1 0 0 1 0 2 2 1
[12,] 1 0 0 0 1 1 0 2 0 0
[13,] 0 0 0 1 1 2 0 2 1 0
[14,] 1 0 1 1 0 2 0 1 0 0
>
> # mean of each SNP(allele frequency)
> snp_mean = colMeans(M_all) / 2
> snp_mean
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
0.32142857 0.17857143 0.35714286 0.35714286 0.14285714 0.60714286 0.07142857 0.96428571 0.57142857 0.39285714
>
> # sum of 2pq(heterozygote ratio)
> sum2pq = 0
> for (i in 1:10) {
+ sum2pq = sum2pq + 2 * snp_mean[i] * (1 - snp_mean[i])
+ }
> sum2pq
snp1
3.538265
>
> # alpha
> alpha = sigma_e / sigma_a
> alpha
[1] 6.95213
>
> # observation
> y = data[1:8, 6]
> y
[1] 9.0 13.4 12.7 15.4 5.9 7.7 10.2 4.8
>
> # design matrix of fixed effect(general mean)
> X = matrix(rep(1, 8), 8, 1)
> X
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
[7,] 1
[8,] 1
>
> # design matrix for DGV
> W = cbind(diag(8), matrix(c(0), 8, 6))
> W
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0
>
> # genomic relationship matrix, G
> Z = sweep(M_all, 2, snp_mean * 2)
> Z
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
[1,] 1.3571429 -0.3571429 0.2857143 0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
[2,] 0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143
[3,] 0.3571429 0.6428571 1.2857143 0.2857143 0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
[4,] -0.6428571 -0.3571429 1.2857143 0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857
[5,] -0.6428571 0.6428571 0.2857143 1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
[6,] 0.3571429 0.6428571 -0.7142857 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 0.2142857
[7,] -0.6428571 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 -0.7857143
[8,] -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 -0.7857143
[9,] 1.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 -0.2142857 1.8571429 0.07142857 -0.1428571 1.2142857
[10,] -0.6428571 -0.3571429 -0.7142857 0.2857143 0.7142857 0.7857143 -0.1428571 0.07142857 -1.1428571 -0.7857143
[11,] -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857
[12,] 0.3571429 -0.3571429 -0.7142857 -0.7142857 0.7142857 -0.2142857 -0.1428571 0.07142857 -1.1428571 -0.7857143
[13,] -0.6428571 -0.3571429 -0.7142857 0.2857143 0.7142857 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143
[14,] 0.3571429 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 -0.92857143 -1.1428571 -0.7857143
> G = Z %*% t(Z) / sum2pq
> round(G, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.472 -0.446 0.988 0.059 0.685 -0.163 -0.708 -0.547 0.887 -0.789 -0.203 -0.143 -0.829 -0.264
[2,] -0.446 0.745 -0.930 -0.446 -0.950 0.180 0.200 0.079 0.099 0.402 -0.143 0.483 0.362 0.362
[3,] 0.988 -0.930 1.634 0.422 1.048 -0.365 -0.627 -0.183 0.120 -0.708 0.160 -0.345 -0.748 -0.466
[4,] 0.059 -0.446 0.422 0.907 0.402 -0.163 0.422 0.301 -0.526 -0.506 0.362 -0.708 -0.264 -0.264
[5,] 0.685 -0.950 1.048 0.402 1.593 -0.102 -0.365 -0.203 -0.183 -0.446 0.140 -0.647 -0.486 -0.486
[6,] -0.163 0.180 -0.365 -0.163 -0.102 0.745 0.200 0.079 0.099 -0.163 0.140 -0.365 0.079 -0.203
[7,] -0.708 0.200 -0.627 0.422 -0.365 0.200 0.786 0.382 -0.728 0.140 0.160 -0.345 0.382 0.099
[8,] -0.547 0.079 -0.183 0.301 -0.203 0.079 0.382 0.826 -0.567 -0.264 0.604 -0.183 -0.022 -0.304
[9,] 0.887 0.099 0.120 -0.526 -0.183 0.099 -0.728 -0.567 2.280 -0.526 -0.224 0.120 -0.567 -0.284
[10,] -0.789 0.402 -0.708 -0.506 -0.446 -0.163 0.140 -0.264 -0.526 1.190 -0.486 0.705 0.867 0.584
[11,] -0.203 -0.143 0.160 0.362 0.140 0.140 0.160 0.604 -0.224 -0.486 0.665 -0.405 -0.244 -0.526
[12,] -0.143 0.483 -0.345 -0.708 -0.647 -0.365 -0.345 -0.183 0.120 0.705 -0.405 1.068 0.382 0.382
[13,] -0.829 0.362 -0.748 -0.264 -0.486 0.079 0.382 -0.022 -0.567 0.867 -0.244 0.382 0.826 0.261
[14,] -0.264 0.362 -0.466 -0.264 -0.486 -0.203 0.099 -0.304 -0.284 0.584 -0.526 0.382 0.261 1.109
>
> # Numerater Relationship Matrix
> A = makenrm(pedi)
> A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 1.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00
[2,] 0.0 1.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.50 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[3,] 0.0 0.0 1.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00
[4,] 0.0 0.0 0.0 1.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[5,] 0.0 0.0 0.0 0.00 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.50 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[6,] 0.0 0.0 0.0 0.00 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.50 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[7,] 0.0 0.0 0.0 0.00 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.50 0.00
[8,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.50 0.00 0.00 0.00 0.00
[9,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.50 0.50 0.0 0.00 0.00 0.00 0.00 0.00
[10,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.50 0.00 0.00
[11,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.50 0.00 0.00 0.00
[12,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.50
[13,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[14,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 1.0 0.0 0.00 0.00 0.50 0.50 0.50 0.0 0.50 0.50 0.50 0.50 0.50
[15,] 0.0 0.0 0.0 0.50 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.50 0.0 1.0 0.50 0.50 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[16,] 0.0 0.5 0.0 0.25 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.5 1.00 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[17,] 0.0 0.0 0.0 0.25 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.5 0.25 1.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[18,] 0.0 0.0 0.0 0.00 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 1.00 0.25 0.25 0.0 0.25 0.25 0.25 0.25 0.25
[19,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 1.00 0.50 0.0 0.25 0.25 0.25 0.25 0.25
[20,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.50 1.00 0.0 0.25 0.25 0.25 0.25 0.25
[21,] 0.5 0.0 0.5 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 1.0 0.00 0.00 0.00 0.00 0.00
[22,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 1.00 0.25 0.25 0.25 0.25
[23,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 1.00 0.25 0.25 0.25
[24,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 1.00 0.25 0.25
[25,] 0.0 0.0 0.0 0.00 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 0.25 1.00 0.25
[26,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 0.25 0.25 1.00
> A22 = A[13:26, 13:26]
> A22
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 1.0 0.0 0.00 0.00 0.50 0.50 0.50 0 0.50 0.50 0.50 0.50 0.50
[3,] 0.50 0.0 1.0 0.50 0.50 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[4,] 0.25 0.0 0.5 1.00 0.25 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[5,] 0.25 0.0 0.5 0.25 1.00 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.5 0.0 0.00 0.00 1.00 0.25 0.25 0 0.25 0.25 0.25 0.25 0.25
[7,] 0.00 0.5 0.0 0.00 0.00 0.25 1.00 0.50 0 0.25 0.25 0.25 0.25 0.25
[8,] 0.00 0.5 0.0 0.00 0.00 0.25 0.50 1.00 0 0.25 0.25 0.25 0.25 0.25
[9,] 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 1 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 1.00 0.25 0.25 0.25 0.25
[11,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 1.00 0.25 0.25 0.25
[12,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 1.00 0.25 0.25
[13,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 0.25 1.00 0.25
[14,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 0.25 0.25 1.00
>
> # weighted G
> wt = 0.99
> Gw = wt * G + (1 - wt) * A22
> round(Gw, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.468 -0.441 0.983 0.061 0.681 -0.161 -0.701 -0.541 0.878 -0.781 -0.201 -0.141 -0.821 -0.261
[2,] -0.441 0.748 -0.921 -0.441 -0.941 0.183 0.203 0.084 0.099 0.403 -0.136 0.483 0.363 0.363
[3,] 0.983 -0.921 1.627 0.423 1.043 -0.361 -0.621 -0.181 0.118 -0.701 0.158 -0.341 -0.741 -0.461
[4,] 0.061 -0.441 0.423 0.908 0.401 -0.161 0.418 0.298 -0.521 -0.501 0.358 -0.701 -0.261 -0.261
[5,] 0.681 -0.941 1.043 0.401 1.587 -0.101 -0.361 -0.201 -0.181 -0.441 0.138 -0.641 -0.481 -0.481
[6,] -0.161 0.183 -0.361 -0.161 -0.101 0.748 0.201 0.081 0.099 -0.159 0.141 -0.359 0.081 -0.199
[7,] -0.701 0.203 -0.621 0.418 -0.361 0.201 0.788 0.383 -0.721 0.141 0.161 -0.339 0.381 0.101
[8,] -0.541 0.084 -0.181 0.298 -0.201 0.081 0.383 0.828 -0.561 -0.259 0.601 -0.179 -0.019 -0.299
[9,] 0.878 0.099 0.118 -0.521 -0.181 0.099 -0.721 -0.561 2.267 -0.521 -0.221 0.118 -0.561 -0.281
[10,] -0.781 0.403 -0.701 -0.501 -0.441 -0.159 0.141 -0.259 -0.521 1.188 -0.479 0.701 0.860 0.581
[11,] -0.201 -0.136 0.158 0.358 0.138 0.141 0.161 0.601 -0.221 -0.479 0.668 -0.399 -0.239 -0.519
[12,] -0.141 0.483 -0.341 -0.701 -0.641 -0.359 -0.339 -0.179 0.118 0.701 -0.399 1.068 0.381 0.381
[13,] -0.821 0.363 -0.741 -0.261 -0.481 0.081 0.381 -0.019 -0.561 0.860 -0.239 0.381 0.828 0.261
[14,] -0.261 0.363 -0.461 -0.261 -0.481 -0.199 0.101 -0.299 -0.281 0.581 -0.519 0.381 0.261 1.108
> Gwi = ginv(Gw)
> round(Gwi, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 36.374 -27.981 -8.256 -7.273 -16.994 1.865 1.034 23.580 2.958 39.443 15.760 -17.655 0.710 5.177
[2,] -27.981 88.271 31.300 -7.935 26.742 -30.838 -16.075 -3.254 -1.422 -44.402 -27.862 -31.454 38.437 -6.093
[3,] -8.256 31.300 17.521 -13.840 10.826 -13.961 10.967 -4.374 1.150 -15.361 -3.670 -9.610 12.371 -1.151
[4,] -7.273 -7.935 -13.840 59.030 1.183 35.036 -44.042 13.642 4.069 8.965 -16.849 24.984 -3.933 4.327
[5,] -16.994 26.742 10.826 1.183 16.560 -3.741 6.055 -10.356 1.661 -23.795 -5.410 8.979 8.462 0.465
[6,] 1.865 -30.838 -13.961 35.036 -3.741 33.909 -9.403 4.611 4.196 17.932 1.618 33.086 -19.287 5.158
[7,] 1.034 -16.075 10.967 -44.042 6.055 -9.403 76.599 -30.644 2.333 -4.681 31.963 24.525 -23.251 0.443
[8,] 23.580 -3.254 -4.374 13.642 -10.356 4.611 -30.644 38.819 3.866 31.071 -11.692 -26.201 11.746 3.079
[9,] 2.958 -1.422 1.150 4.069 1.661 4.196 2.333 3.866 3.572 4.523 2.296 3.933 1.097 2.907
[10,] 39.443 -44.402 -15.361 8.965 -23.795 17.932 -4.681 31.071 4.523 65.241 12.420 -9.861 -22.528 2.734
[11,] 15.760 -27.862 -3.670 -16.849 -5.410 1.618 31.963 -11.692 2.296 12.420 33.290 12.724 -8.689 6.556
[12,] -17.655 -31.454 -9.610 24.984 8.979 33.086 24.525 -26.201 3.933 -9.861 12.724 65.530 -29.485 5.120
[13,] 0.710 38.437 12.371 -3.933 8.462 -19.287 -23.251 11.746 1.097 -22.528 -8.689 -29.485 46.431 4.505
[14,] 5.177 -6.093 -1.151 4.327 0.465 5.158 0.443 3.079 2.907 2.734 6.556 5.120 4.505 6.133
>
> # Left Hand Side
> LHS = rbind(
+ cbind(t(X) %*% X, t(X) %*% W),
+ cbind(t(W) %*% X, t(W) %*% W + Gwi * alpha))
> round(LHS, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 8 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000
[2,] 1 253.878 -194.527 -57.394 -50.564 -118.147 12.966 7.191 163.929 20.563 274.211 109.565 -122.742 4.935 35.992
[3,] 1 -194.527 614.671 217.604 -55.166 185.914 -214.391 -111.754 -22.623 -9.886 -308.686 -193.698 -218.671 267.222 -42.358
[4,] 1 -57.394 217.604 122.805 -96.215 75.261 -97.060 76.244 -30.406 7.998 -106.789 -25.516 -66.811 86.008 -8.002
[5,] 1 -50.564 -55.166 -96.215 411.383 8.227 243.574 -306.187 94.843 28.289 62.324 -117.135 173.693 -27.340 30.084
[6,] 1 -118.147 185.914 75.261 8.227 116.124 -26.006 42.094 -71.996 11.545 -165.427 -37.608 62.421 58.832 3.233
[7,] 1 12.966 -214.391 -97.060 243.574 -26.006 236.736 -65.371 32.055 29.169 124.667 11.249 230.015 -134.084 35.857
[8,] 1 7.191 -111.754 76.244 -306.187 42.094 -65.371 533.526 -213.039 16.218 -32.543 222.212 170.498 -161.643 3.080
[9,] 1 163.929 -22.623 -30.406 94.843 -71.996 32.055 -213.039 270.876 26.875 216.008 -81.283 -182.155 81.657 21.406
[10,] 0 20.563 -9.886 7.998 28.289 11.545 29.169 16.218 26.875 24.832 31.445 15.965 27.341 7.624 20.209
[11,] 0 274.211 -308.686 -106.789 62.324 -165.427 124.667 -32.543 216.008 31.445 453.565 86.347 -68.557 -156.619 19.006
[12,] 0 109.565 -193.698 -25.516 -117.135 -37.608 11.249 222.212 -81.283 15.965 86.347 231.435 88.456 -60.406 45.576
[13,] 0 -122.742 -218.671 -66.811 173.693 62.421 230.015 170.498 -182.155 27.341 -68.557 88.456 455.571 -204.985 35.597
[14,] 0 4.935 267.222 86.008 -27.340 58.832 -134.084 -161.643 81.657 7.624 -156.619 -60.406 -204.985 322.794 31.319
[15,] 0 35.992 -42.358 -8.002 30.084 3.233 35.857 3.080 21.406 20.209 19.006 45.576 35.597 31.319 42.639
>
> # Right Hand Side
> RHS = rbind(t(X) %*% y,
+ t(W) %*% y)
> round(RHS, 3)
[,1]
[1,] 79.1
[2,] 9.0
[3,] 13.4
[4,] 12.7
[5,] 15.4
[6,] 5.9
[7,] 7.7
[8,] 10.2
[9,] 4.8
[10,] 0.0
[11,] 0.0
[12,] 0.0
[13,] 0.0
[14,] 0.0
[15,] 0.0
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
> round(gi_LHS, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0.139 -0.018 0.019 -0.025 -0.029 -0.027 -0.010 -0.008 -0.015 0.013 0.036 -0.020 0.035 0.022 0.023
[2,] -0.018 0.152 -0.033 0.085 0.007 0.054 -0.009 -0.064 -0.051 0.095 -0.088 -0.019 -0.016 -0.086 -0.026
[3,] 0.019 -0.033 0.074 -0.085 -0.048 -0.091 0.014 0.012 0.003 0.016 0.039 -0.017 0.054 0.031 0.038
[4,] -0.025 0.085 -0.085 0.158 0.044 0.088 -0.031 -0.054 -0.007 0.002 -0.070 0.023 -0.033 -0.070 -0.046
[5,] -0.029 0.007 -0.048 0.044 0.111 0.041 -0.015 0.055 0.037 -0.061 -0.063 0.044 -0.087 -0.032 -0.031
[6,] -0.027 0.054 -0.091 0.088 0.041 0.164 0.000 -0.025 -0.013 -0.032 -0.039 0.019 -0.071 -0.040 -0.050
[7,] -0.010 -0.009 0.014 -0.031 -0.015 0.000 0.095 0.021 0.008 0.015 -0.030 0.020 -0.053 0.002 -0.032
[8,] -0.008 -0.064 0.012 -0.054 0.055 -0.025 0.021 0.085 0.037 -0.081 0.007 0.016 -0.047 0.035 0.007
[9,] -0.015 -0.051 0.003 -0.007 0.037 -0.013 0.008 0.037 0.100 -0.062 -0.043 0.076 -0.024 -0.014 -0.044
[10,] 0.013 0.095 0.016 0.002 -0.061 -0.032 0.015 -0.081 -0.062 0.296 -0.067 -0.020 0.008 -0.068 -0.040
[11,] 0.036 -0.088 0.039 -0.070 -0.063 -0.039 -0.030 0.007 -0.043 -0.067 0.154 -0.067 0.092 0.107 0.073
[12,] -0.020 -0.019 -0.017 0.023 0.044 0.019 0.020 0.016 0.076 -0.020 -0.067 0.088 -0.052 -0.036 -0.071
[13,] 0.035 -0.016 0.054 -0.033 -0.087 -0.071 -0.053 -0.047 -0.024 0.008 0.092 -0.052 0.138 0.049 0.047
[14,] 0.022 -0.086 0.031 -0.070 -0.032 -0.040 0.002 0.035 -0.014 -0.068 0.107 -0.036 0.049 0.099 0.028
[15,] 0.023 -0.026 0.038 -0.046 -0.031 -0.050 -0.032 0.007 -0.044 -0.040 0.073 -0.071 0.047 0.028 0.151
>
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
[,1]
[1,] 9.943
[2,] 0.070
[3,] 0.111
[4,] 0.048
[5,] 0.258
[6,] -0.493
[7,] -0.355
[8,] 0.143
[9,] -0.227
[10,] 0.027
[11,] 0.114
[12,] -0.238
[13,] 0.142
[14,] 0.054
[15,] 0.350
>
책의 결과와 값이 다르다. 책에는 A(NRM)을 섞었다는 설명이 없지만 G에 0.99, A에 0.01 가중치를 주었을 때 비슷한 값을 나타내었다. 그래도 같은 결과값이 나오지 않은 이유는 책에 나오는 5세대까지 추적하여 얻은 A22에 0.01의 가중치를 주어 구한 값으로 생각된다. 하지만 우리는 혈통이 없으므로 책과 같은 A22를 구할 수 없고 그래서 책과 같은 결과값을 구하지 못한 것으로 생각된다. 책에 나오는 A22를 입력하여 사용하면 같은 값이 나올까?
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| GBLUP Model with residual polygenic effects (0) | 2020.12.29 |
|---|---|
| SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0) | 2020.12.24 |
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis (0) | 2020.12.19 |
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis (0) | 2020.12.19 |
| Fixed Effect Model for SNP Effects, weighted analysis (0) | 2020.12.16 |