# GBLUP Model with residual polygenic effects

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.

# Raphael Mrode

# Example 11.2 p189

 

GBLUP Model은 유전체 자료를 가지고 개체들의 DGV를 구하는 모형이다. 여기서 유전체 자료가 모든 상가적 유전 효과를 설명할 수 있다고 보기 어려워 polygenic effect를 추가한다. 간단한 설명은 다음 포스팅을 참고한다.

2020/12/24 - [Animal Breeding/R for Genetic Evaluation] - SNP-BLUP Model with Polygenic Effects(unweighted analysis)

 

SNP-BLUP Model with Polygenic Effects(unweighted analysis)

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.5 p189 SNP effect로 모든 상가적 유전 분산을 설명할 수 없을 수도 있다. 그래서 설명..

bhpark.tistory.com

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

 

program

# GBLUP Model with residual polygenic effects

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.

# Raphael Mrode

# Example 11.2 p189

# MASS packages
if (!("MASS" %in% installed.packages())) {
  install.packages("MASS", dependencies = TRUE)
}
library(MASS)

# function to make inverse of numerator relationship matrix
ainv = function(pedi) {
  n = nrow(pedi)
  Ainv = 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) {
      # both parents unknown
      alpha = 1
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
    } else if (sire != 0 & dam == 0) {
      # sire known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
    } else if (sire == 0 & dam != 0) {
      # dam known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    } else {
      # both parents known
      alpha = 2
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
      Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
      Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    }
  }
  return(Ainv)
}

# 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

# ratio of polygenic effect
rpe = 0.1

# variance of residual polygenic effect
sigma_u = sigma_a * rpe
sigma_u

# variance ratio, alpha
alpha1 = sigma_e / sigma_u
alpha1

alpha2 = sigma_e / (sigma_a - sigma_u)
alpha2

# 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 = sum(2 * snp_mean * (1 - snp_mean))
sum2pq

# 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 polygenic effects
W1 = cbind(matrix(c(0), 8, 12), diag(8), matrix(c(0), 8, 6))
W1
# inverse matrix of NRM 
Ai = ainv(pedi) 
Ai

# design matrix for DGV
W2 = cbind(diag(8), matrix(c(0), 8, 6))
W2

# genomic relationship matrix, G
Z = sweep(M_all, 2, snp_mean * 2)
round(Z, 3)
G = Z %*% t(Z) / sum2pq
round(G, 3)

# Numerator 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)  %*% W1,               t(X)  %*% W2), 
  cbind(t(W1) %*% X, t(W1) %*% W1 + Ai * alpha1, t(W1) %*% W2),
  cbind(t(W2) %*% X, t(W2) %*% W1,               t(W2) %*% W2 + Gwi * alpha2))
round(LHS, 3)

# Right Hand Side
RHS = rbind(t(X)  %*% y, 
            t(W1) %*% y,
            t(W2) %*% 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)

 

RStudio에서 실행한 화면

 

프로그램 실행 결과

> # GBLUP Model with residual polygenic effects
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
> 
> # Raphael Mrode
> 
> # Example 11.2 p189
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+   install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # function to make inverse of numerator relationship matrix
> ainv = function(pedi) {
+   n = nrow(pedi)
+   Ainv = 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) {
+       # both parents unknown
+       alpha = 1
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+     } else if (sire != 0 & dam == 0) {
+       # sire known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+     } else if (sire == 0 & dam != 0) {
+       # dam known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     } else {
+       # both parents known
+       alpha = 2
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+       Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+       Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     }
+   }
+   return(Ainv)
+ }
> 
> # 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/23_GBLUP with polygenic effect"
> 
> # 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
> 
> # ratio of polygenic effect
> rpe = 0.1
> 
> # variance of residual polygenic effect
> sigma_u = sigma_a * rpe
> sigma_u
[1] 3.5241
> 
> # variance ratio, alpha
> alpha1 = sigma_e / sigma_u
> alpha1
[1] 69.5213
> 
> alpha2 = sigma_e / (sigma_a - sigma_u)
> alpha2
[1] 7.724588
> 
> # 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 = sum(2 * snp_mean * (1 - snp_mean))
> sum2pq
[1] 3.538265
> 
> # 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 polygenic effects
> W1 = cbind(matrix(c(0), 8, 12), diag(8), matrix(c(0), 8, 6))
> W1
     [,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,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0
[2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
[3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0
[4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0
[5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0
[6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0
[7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0
[8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0
> # inverse matrix of NRM 
> Ai = ainv(pedi) 
> Ai
      [,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.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [2,]  0.0  1.5  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5    -1     0     0     0     0     0     0     0     0     0     0
 [3,]  0.5  0.0  1.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [4,]  0.0  0.0  0.0  1.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
 [5,]  0.0  0.0  0.0  0.0  1.5  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5     0    -1     0     0     0     0     0     0     0     0     0
 [6,]  0.0  0.0  0.0  0.0  0.0  1.5  0.0  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0    -1     0     0     0     0     0     0     0     0
 [7,]  0.0  0.0  0.0  0.0  0.0  0.0  1.5  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0    -1     0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.5    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0    -1     0     0     0     0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    2   0.0   0.0   0.0   0.0   1.0   0.0     0     0     0    -1    -1     0     0     0     0     0     0
[10,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   1.5   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0    -1     0     0
[11,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   1.5   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0    -1     0     0     0
[12,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   1.5   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0     0    -1
[13,]  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   1.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
[14,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5    1   0.5   0.5   0.5   0.0   5.0   0.0     0     0    -1    -1    -1     0    -1    -1    -1    -1    -1
[15,]  0.0  0.5  0.0 -1.0  0.5  0.0  0.0  0.0    0   0.0   0.0   0.0  -1.0   0.0   3.0    -1    -1     0     0     0     0     0     0     0     0     0
[16,]  0.0 -1.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     2     0     0     0     0     0     0     0     0     0     0
[17,]  0.0  0.0  0.0  0.0 -1.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     0     2     0     0     0     0     0     0     0     0     0
[18,]  0.0  0.0  0.0  0.0  0.0 -1.0  0.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     2     0     0     0     0     0     0     0     0
[19,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     2     0     0     0     0     0     0     0
[20,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     2     0     0     0     0     0     0
[21,] -1.0  0.0 -1.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0     2     0     0     0     0     0
[22,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     2     0     0     0     0
[23,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0  -1.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     2     0     0     0
[24,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0  -1.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     2     0     0
[25,]  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     2     0
[26,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0  -1.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     0     2
> 
> # design matrix for DGV
> W2 = cbind(diag(8), matrix(c(0), 8, 6))
> W2
     [,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)
> round(Z, 3)
        snp1   snp2   snp3   snp4   snp5   snp6   snp7   snp8   snp9  snp10
 [1,]  1.357 -0.357  0.286  0.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
 [2,]  0.357 -0.357 -0.714 -0.714 -0.286  0.786 -0.143  0.071 -0.143 -0.786
 [3,]  0.357  0.643  1.286  0.286  0.714 -1.214 -0.143  0.071 -0.143  1.214
 [4,] -0.643 -0.357  1.286  0.286 -0.286 -0.214 -0.143  0.071  0.857  0.214
 [5,] -0.643  0.643  0.286  1.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
 [6,]  0.357  0.643 -0.714  0.286 -0.286  0.786 -0.143  0.071  0.857  0.214
 [7,] -0.643 -0.357  0.286  0.286 -0.286  0.786 -0.143  0.071  0.857 -0.786
 [8,] -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143  0.071  0.857 -0.786
 [9,]  1.357 -0.357 -0.714 -0.714 -0.286 -0.214  1.857  0.071 -0.143  1.214
[10,] -0.643 -0.357 -0.714  0.286  0.714  0.786 -0.143  0.071 -1.143 -0.786
[11,] -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143  0.071  0.857  0.214
[12,]  0.357 -0.357 -0.714 -0.714  0.714 -0.214 -0.143  0.071 -1.143 -0.786
[13,] -0.643 -0.357 -0.714  0.286  0.714  0.786 -0.143  0.071 -0.143 -0.786
[14,]  0.357 -0.357  0.286  0.286 -0.286  0.786 -0.143 -0.929 -1.143 -0.786
> 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
> 
> # Numerator 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)  %*% W1,               t(X)  %*% W2), 
+   cbind(t(W1) %*% X, t(W1) %*% W1 + Ai * alpha1, t(W1) %*% W2),
+   cbind(t(W2) %*% X, t(W2) %*% W1,               t(W2) %*% W2 + Gwi * alpha2))
> round(LHS, 3)
      [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]   [,10]   [,11]   [,12]   [,13]   [,14]   [,15]   [,16]   [,17]   [,18]   [,19]   [,20]   [,21]   [,22]   [,23]
 [1,]    8   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   0.000   0.000
 [2,]    0 104.282   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000
 [3,]    0   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761 -69.521   0.000   0.000   0.000   0.000   0.000   0.000
 [4,]    0  34.761   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000
 [5,]    0   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 [6,]    0   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000 -69.521   0.000   0.000   0.000   0.000   0.000
 [7,]    0   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000
 [8,]    0   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 [9,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521
[10,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043   0.000   0.000   0.000   0.000  69.521   0.000   0.000   0.000   0.000 -69.521 -69.521   0.000   0.000
[11,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
[12,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
[13,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
[14,]    1   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 105.282   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000
[15,]    1   0.000   0.000   0.000   0.000   0.000  34.761  34.761  34.761  69.521  34.761  34.761  34.761   0.000 348.606   0.000   0.000   0.000 -69.521 -69.521 -69.521   0.000 -69.521
[16,]    1   0.000  34.761   0.000 -69.521  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000 209.564 -69.521 -69.521   0.000   0.000   0.000   0.000   0.000
[17,]    1   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521 140.043   0.000   0.000   0.000   0.000   0.000   0.000
[18,]    1   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000 140.043   0.000   0.000   0.000   0.000   0.000
[19,]    1   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000 140.043   0.000   0.000   0.000   0.000
[20,]    1   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 140.043   0.000   0.000   0.000
[21,]    1   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000 140.043   0.000   0.000
[22,]    0 -69.521   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043   0.000
[23,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043
[24,]    0   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
        [,24]   [,25]   [,26]   [,27]    [,28]    [,29]    [,30]    [,31]    [,32]    [,33]    [,34]    [,35]   [,36]    [,37]    [,38]    [,39]    [,40]   [,41]
 [1,]   0.000   0.000   0.000   0.000    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,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [3,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [4,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [5,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [6,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [7,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [8,]   0.000   0.000 -69.521   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [9,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[10,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[11,]   0.000 -69.521   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[12,] -69.521   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[13,]   0.000   0.000   0.000 -69.521    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[14,]   0.000   0.000   0.000   0.000    1.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[15,] -69.521 -69.521 -69.521 -69.521    0.000    1.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[16,]   0.000   0.000   0.000   0.000    0.000    0.000    1.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[17,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    1.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[18,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    1.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[19,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    1.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[20,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    1.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[21,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    1.000   0.000    0.000    0.000    0.000    0.000   0.000
[22,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[23,]   0.000   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
[24,] 139.043   0.000   0.000   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000   0.000    0.000    0.000    0.000    0.000   0.000
 [ getOption("max.print") 에 도달했습니다 -- 17 행들을 생략합니다 ]
> 
> # Right Hand Side
> RHS = rbind(t(X)  %*% y, 
+             t(W1) %*% y,
+             t(W2) %*% y)
> round(RHS, 3)
      [,1]
 [1,] 79.1
 [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,]  9.0
[15,] 13.4
[16,] 12.7
[17,] 15.4
[18,]  5.9
[19,]  7.7
[20,] 10.2
[21,]  4.8
[22,]  0.0
[23,]  0.0
[24,]  0.0
[25,]  0.0
[26,]  0.0
[27,]  0.0
[28,]  9.0
[29,] 13.4
[30,] 12.7
[31,] 15.4
[32,]  5.9
[33,]  7.7
[34,] 10.2
[35,]  4.8
[36,]  0.0
[37,]  0.0
[38,]  0.0
[39,]  0.0
[40,]  0.0
[41,]  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]  [,16]  [,17]  [,18]  [,19]  [,20]  [,21] [,22]  [,23]  [,24]  [,25]  [,26]  [,27]
 [1,]  0.142 0.000 -0.001 0.000 -0.002 -0.001 -0.001 0.000 0.000 -0.002 0.000 0.000 0.000 -0.003 -0.005 -0.004 -0.003 -0.003 -0.004 -0.004 -0.004 0.000 -0.002 -0.002 -0.002 -0.002 -0.002
 [2,]  0.000 0.014  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.007  0.000  0.000  0.000  0.000  0.000
 [3,] -0.001 0.000  0.014 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.007  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [4,]  0.000 0.000  0.000 0.014  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.007  0.000  0.000  0.000  0.000  0.000
 [5,] -0.002 0.000  0.000 0.000  0.014  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.007  0.004  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [6,] -0.001 0.000  0.000 0.000  0.000  0.014  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.007  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [7,] -0.001 0.000  0.000 0.000  0.000  0.000  0.014 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.007  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [8,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.014 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.007  0.000
 [9,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.014  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.007  0.000  0.000  0.000  0.000
[10,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.014 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.007  0.007 0.000  0.000  0.000  0.000  0.000  0.000
[11,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.014 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.007  0.000  0.000
[12,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.014 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.007  0.000  0.000  0.000
[13,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.014  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.007
[14,] -0.003 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.014  0.000  0.007  0.004  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[15,] -0.005 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.014  0.000  0.000  0.000  0.007  0.007  0.007 0.000  0.007  0.007  0.007  0.007  0.007
[16,] -0.004 0.000  0.000 0.000  0.007  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.007  0.000  0.014  0.007  0.007  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[17,] -0.003 0.000  0.007 0.000  0.004  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.004  0.000  0.007  0.014  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[18,] -0.003 0.000  0.000 0.000  0.004  0.007  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.004  0.000  0.007  0.004  0.014  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[19,] -0.004 0.000  0.000 0.000  0.000  0.000  0.007 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.014  0.004  0.004 0.000  0.004  0.004  0.004  0.004  0.004
[20,] -0.004 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.007 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.014  0.007 0.000  0.004  0.004  0.004  0.004  0.004
[21,] -0.004 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.007 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.007  0.014 0.000  0.004  0.004  0.004  0.004  0.004
[22,]  0.000 0.007  0.000 0.007  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.014  0.000  0.000  0.000  0.000  0.000
[23,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.007  0.000 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.014  0.004  0.004  0.004  0.004
[24,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.007 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.004  0.014  0.004  0.004  0.004
       [,28]  [,29]  [,30]  [,31]  [,32]  [,33]  [,34]  [,35]  [,36]  [,37]  [,38]  [,39]  [,40]  [,41]
 [1,] -0.017  0.018 -0.023 -0.027 -0.025 -0.009 -0.007 -0.013  0.012  0.033 -0.018  0.032  0.020  0.021
 [2,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [3,]  0.000  0.000  0.000 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [4,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [5,] -0.001  0.001 -0.001  0.000 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [6,]  0.000  0.000  0.000  0.000 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [7,]  0.000  0.000  0.000  0.000  0.000 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [8,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [9,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[10,]  0.001  0.000  0.001  0.000  0.001  0.000 -0.001 -0.001  0.001  0.000  0.000  0.000  0.000  0.000
[11,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[12,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[13,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[14,] -0.002  0.001 -0.002  0.000 -0.001  0.001  0.001  0.001 -0.001  0.001  0.000  0.000  0.001  0.000
[15,]  0.002 -0.002  0.003  0.001  0.002 -0.001 -0.001 -0.001  0.000 -0.001  0.000 -0.001 -0.001 -0.001
[16,] -0.002  0.002 -0.003 -0.001 -0.002  0.001  0.001  0.001  0.000  0.001  0.000  0.000  0.001  0.001
[17,] -0.001  0.001 -0.002 -0.001 -0.001  0.001  0.000  0.000  0.000  0.001  0.000  0.001  0.001  0.000
[18,] -0.001  0.002 -0.002  0.000 -0.002  0.001  0.001  0.001  0.000  0.001  0.000  0.001  0.001  0.001
[19,]  0.001 -0.001  0.002  0.001  0.001 -0.001 -0.001  0.000  0.000 -0.001  0.000  0.000 -0.001  0.000
[20,]  0.002 -0.001  0.002  0.000  0.002  0.000 -0.001 -0.001  0.001 -0.001  0.000  0.000 -0.001 -0.001
[21,]  0.002 -0.001  0.002  0.000  0.002  0.000 -0.001 -0.001  0.001 -0.001  0.000 -0.001 -0.001  0.000
[22,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[23,]  0.001 -0.001  0.001  0.001  0.001  0.000  0.000  0.000  0.000 -0.001  0.000 -0.001 -0.001  0.000
[24,]  0.001 -0.001  0.001  0.001  0.001  0.000  0.000  0.000  0.000 -0.001  0.000 -0.001 -0.001  0.000
 [ getOption("max.print") 에 도달했습니다 -- 17 행들을 생략합니다 ]
> 
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
        [,1]
 [1,]  9.936
 [2,]  0.000
 [3,]  0.037
 [4,]  0.000
 [5,]  0.025
 [6,] -0.026
 [7,] -0.014
 [8,]  0.000
 [9,]  0.000
[10,] -0.034
[11,]  0.000
[12,]  0.000
[13,]  0.000
[14,]  0.011
[15,]  0.001
[16,]  0.043
[17,]  0.077
[18,] -0.017
[19,] -0.020
[20,] -0.016
[21,] -0.052
[22,]  0.000
[23,]  0.000
[24,]  0.000
[25,]  0.000
[26,]  0.000
[27,]  0.000
[28,]  0.055
[29,]  0.107
[30,]  0.032
[31,]  0.228
[32,] -0.454
[33,] -0.317
[34,]  0.132
[35,] -0.201
[36,]  0.023
[37,]  0.107
[38,] -0.214
[39,]  0.132
[40,]  0.053
[41,]  0.318
>

 

저작자표시

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기  (0) 2021.04.15
Single-step GBLUP(ssGBLUP)  (0) 2021.01.05
SNP-BLUP Model with Polygenic Effects(unweighted analysis)  (0) 2020.12.24
Mixed Linear Model(GBLUP model) for Computing DGV  (0) 2020.12.21
Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis  (0) 2020.12.19

blupf90으로 gblup model을 푸는 것은 다음 포스팅을 참고한다.

2020/12/29 - [Animal Breeding/BLUPF90] - blupf90으로 gblup model 풀기

 

blupf90으로 gblup model 풀기

blupf90은 ssGBLUP을 위하여 만든 프로그램이다. 그래서 gblup model을 다루기 위해서는 약간의 트릭이 필요하다. 그래서 renumf90으로 renumbering한 후 pregsf90을 실행하고 blupf90을 실행할 수 없다. pregsf90..

bhpark.tistory.com

 

위 포스팅에서 pregsf90 단계까지는 같은 자료를 가지고 똑같이 진행한다.

 

blupf90 실행을 위한 파라미터 파일

# BLUPF90 parameter file created by RENUMF90
DATAFILE
gblup_data.txt
NUMBER_OF_TRAITS
1
NUMBER_OF_EFFECTS
3
OBSERVATION(S)
6
WEIGHT(S)
 
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
4  1 cross 
1 26 cross # residual polygenic effect
8 14 cross # new ID (renumbered only for genotyped animals)
RANDOM_RESIDUAL VALUES
245.00    
RANDOM_GROUP
2
RANDOM_TYPE
add_animal
FILE
gblup_pedi.txt                                                            
(CO)VARIANCES
3.5241
RANDOM_GROUP
3
RANDOM_TYPE
user_file   
FILE
Gi                                                                    
(CO)VARIANCES
31.717
OPTION solv_method FSPAK

 

polygenic effect를 추가하여 NUMBER_OF_EFFECTS가 3이 되었다. 3개의 effects는 다음과 같다.

4  1 cross  
1 26 cross # residual polygenic effect 
8 14 cross # new ID (renumbered only for genotyped animals)

 

둘째 효과가 polygenic effect이다. 추가된 effect에 대한 설명은 다음과 같다.

RANDOM_GROUP
2
RANDOM_TYPE
add_animal
FILE
gblup_pedi.txt                                                            
(CO)VARIANCES
3.5241

나머지는 이전 분석과 같다.

 

blupf90을 실행한 화면

 

blupf90 실행 로그

blupf90_gblup_polygenic.par
     BLUPF90 ver. 1.68

 Parameter file:             blupf90_gblup_polygenic.par
 Data file:                  gblup_data.txt
 Number of Traits             1
 Number of Effects            3
 Position of Observations      6
 Position of Weight (1)        0
 Value of Missing Trait/Observation           0

EFFECTS
 #  type                position (2)        levels   [positions for nested]
    1  cross-classified       4         1
    2  cross-classified       1        26
    3  cross-classified       8        14

 Residual (co)variance Matrix
  245.00    

 Random Effect(s)    2
 Type of Random Effect:      additive animal
 Pedigree File:              gblup_pedi.txt                                                                                                                                                                                                                                            
 trait   effect    (CO)VARIANCES
  1       2     3.524    

 Random Effect(s)    3
 Type of Random Effect:      user defined from file
 User File:                  Gi                                                                                                                                                                                                                                                        
 trait   effect    (CO)VARIANCES
  1       3     31.72    

 REMARKS
  (1) Weight position 0 means no weights utilized
  (2) Effect positions of 0 for some effects and traits means that such
      effects are missing for specified traits
 

 * The limited number of OpenMP threads = 4

 * solving method (default=PCG):FSPAK               
 
 Data record length =            8
 # equations =           41
 G
  3.5241    
 G
  31.717    
 read            8  records in   0.1562500      s,                      41 
  nonzeroes
  read           26  additive pedigrees
 g_usr_inv: read          105  elements
 largest row, column, diagonal:          14          14          14
 user defined matrix
     36.3741    -27.9809     -8.2555     -7.2732    -16.9944      1.8651      1.0343     23.5797      2.9578     39.4427     15.7599    -17.6553      0.7099      5.1771
    -27.9809     88.2709     31.3003     -7.9351     26.7420    -30.8381    -16.0748     -3.2541     -1.4220    -44.4016    -27.8617    -31.4539     38.4374     -6.0928
     -8.2555     31.3003     17.5206    -13.8396     10.8256    -13.9611     10.9670     -4.3736      1.1504    -15.3606     -3.6702     -9.6102     12.3714     -1.1511
     -7.2732     -7.9351    -13.8396     59.0298      1.1834     35.0358    -44.0422     13.6423      4.0692      8.9647    -16.8488     24.9842     -3.9326      4.3273
    -16.9944     26.7420     10.8256      1.1834     16.5595     -3.7407      6.0548    -10.3559      1.6606    -23.7952     -5.4096      8.9787      8.4624      0.4650
      1.8651    -30.8381    -13.9611     35.0358     -3.7407     33.9085     -9.4030      4.6108      4.1957     17.9323      1.6181     33.0855    -19.2868      5.1577
      1.0343    -16.0748     10.9670    -44.0422      6.0548     -9.4030     76.5989    -30.6437      2.3328     -4.6811     31.9632     24.5245    -23.2509      0.4430
     23.5797     -3.2541     -4.3736     13.6423    -10.3559      4.6108    -30.6437     38.8191      3.8658     31.0707    -11.6919    -26.2014     11.7456      3.0790
      2.9578     -1.4220      1.1504      4.0692      1.6606      4.1957      2.3328      3.8658      3.5719      4.5230      2.2964      3.9327      1.0966      2.9069
     39.4427    -44.4016    -15.3606      8.9647    -23.7952     17.9323     -4.6811     31.0707      4.5230     65.2412     12.4202     -9.8613    -22.5282      2.7339
     15.7599    -27.8617     -3.6702    -16.8488     -5.4096      1.6181     31.9632    -11.6919      2.2964     12.4202     33.2898     12.7236     -8.6889      6.5557
    -17.6553    -31.4539     -9.6102     24.9842      8.9787     33.0855     24.5245    -26.2014      3.9327     -9.8613     12.7236     65.5298    -29.4853      5.1203
      0.7099     38.4374     12.3714     -3.9326      8.4624    -19.2868    -23.2509     11.7456      1.0966    -22.5282     -8.6889    -29.4853     46.4310      4.5049
      5.1771     -6.0928     -1.1511      4.3273      0.4650      5.1577      0.4430      3.0790      2.9069      2.7339      6.5557      5.1203      4.5049      6.1332
 finished peds in   0.1562500      s,                     191  nonzeroes
 solutions stored in file: "solutions"

 

blupf90 실행 결과 : solutions

trait/effect level  solution
   1   1         1          9.93629996
   1   2         1          0.00000000
   1   2         2          0.03710016
   1   2         3          0.00000000
   1   2         4          0.02506699
   1   2         5         -0.02564345
   1   2         6         -0.01365606
   1   2         7         -0.00000000
   1   2         8         -0.00000000
   1   2         9         -0.03406764
   1   2        10         -0.00000000
   1   2        11         -0.00000000
   1   2        12         -0.00000000
   1   2        13          0.01065311
   1   2        14          0.00054688
   1   2        15          0.04292705
   1   2        16          0.07711376
   1   2        17         -0.01700165
   1   2        18         -0.02021065
   1   2        19         -0.01570340
   1   2        20         -0.05188500
   1   2        21          0.00000000
   1   2        22          0.00027344
   1   2        23          0.00027344
   1   2        24          0.00027344
   1   2        25          0.00027344
   1   2        26          0.00027344
   1   3         1          0.05511872
   1   3         2          0.10731930
   1   3         3          0.03187816
   1   3         4          0.22808352
   1   3         5         -0.45376632
   1   3         6         -0.31731567
   1   3         7          0.13243836
   1   3         8         -0.20059690
   1   3         9          0.02264487
   1   3        10          0.10681153
   1   3        11         -0.21392987
   1   3        12          0.13204457
   1   3        13          0.05330732
   1   3        14          0.31846507

 

GBLUP with polygenic effect_blupf90.zip
0.00MB

저작자표시

'Animal Breeding > BLUPF90' 카테고리의 다른 글

blupf90으로 social interaction model 풀기  (1) 2021.04.21
blupf90으로 single-step GBLUP Model 풀기  (0) 2021.01.05
blupf90으로 gblup model 풀기  (0) 2020.12.29
blupf90으로 SNP-BLUP Model(with polygenic effects unweighted analysis) 풀기  (0) 2020.12.24
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기  (0) 2020.12.19

blupf90은 ssGBLUP을 위하여 만든 프로그램이다. 그래서 gblup model을 다루기 위해서는 약간의 트릭이 필요하다. 그래서 renumf90으로 renumbering한 후 renumf90이 만들어낸 파라미터 파일을 이용하여 pregsf90을 실행하고 blupf90을 실행할 수 없다. pregsf90을 실행할 파라미터 파일을 직접 만들어 주어야 한다. 

 

GBLUP과 ssGBLUP의 차이는 유전체 자료를 가지고 있지 않은 개체가 평가 모형에 포함되느냐 아니냐의 차이다. 유전체 자료(SNP data)를 이용해서 NRM과 비슷한 행렬을 만든다. 그것을 GRM(genomic relationship matrix)라고 한다. 이 GRM의 효율을 높이기 위하여 NRM을 약간 섞는데 NRM 1%, 5%, 10%정도 섞을 수 있다. 이상은 Mrode 책에 나오는 설명인데 예제에서는 정확하게 NRM을 몇 % 섞었는지 설명이 없다. 아래 예제에서는 NRM을 1% 섞는 것으로 한다. 또한 책에서는 주어진 혈통 그 이상을 알고 있어 계산한 NRM이 주어지는데 본 예제에서는 주어진 혈통 이외의 혈통을 알 수가 없으므로 책에 주어진 NRM과 같은 NRM을 구할 수가 없다. 결론적으로 책과 같은 해를 구할 수가 없다. 자세한 사항은 다음 글을 참고한다.

masuday.github.io/blupf90_tutorial/mrode_c11ex113_gblup.html

 

Numerical examples from Mrode (2014)

Numerical examples from Mrode (2014) Yutaka Masuda September 2019 Back to index.html. GBLUP Model The SNP marker model is equivalent to a model, so-called GBLUP, under some assumptions. The model includes a genomic relationship matrix, \(\mathbf{G}\). This

masuday.github.io

 

gblup_snpgenotype.txt

13 20110002120000000000000000000000000000000000000000
14 10000202100000000000000000000000000000000000000000
15 11211002120000000000000000000000000000000000000000
16 00210102210000000000000000000000000000000000000000
17 01120002120000000000000000000000000000000000000000
18 11010202210000000000000000000000000000000000000000
19 00110202200000000000000000000000000000000000000000
20 01100102200000000000000000000000000000000000000000
21 20000122120000000000000000000000000000000000000000
22 00011202000000000000000000000000000000000000000000
23 01100102210000000000000000000000000000000000000000
24 10001102000000000000000000000000000000000000000000
25 00011202100000000000000000000000000000000000000000
26 10110201000000000000000000000000000000000000000000

 

10개의 SNP 이외의 genotype이 있는데 pregsf90이 50개 미만의 SNP자료를 읽지를 못하여 부가적으로 붙인 것인다. 결과에 영향을 미치지 않는다.

gblup_snpgenotype_XrefID.txt

13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26

 

gblup_data.txt

13 0 0 1 558 9 0.001792115 1
14 0 0 1 722 13.4 0.001385042 2
15 13 4 1 300 12.7 0.003333333 3
16 15 2 1 73 15.4 0.01369863 4
17 15 5 1 52 5.9 0.019230769 5
18 14 6 1 87 7.7 0.011494253 6
19 14 9 1 64 10.2 0.015625 7
20 14 9 1 103 4.8 0.009708738 8

관측치는 6열이고 8열은 1열(개체 ID)의 새로운 ID이다.

 

gblup_pedi.txt

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

 

pregsf90_gblup.par

# BLUPF90 parameter file created by RENUMF90
DATAFILE
gblup_data.txt
NUMBER_OF_TRAITS
1
NUMBER_OF_EFFECTS
2
OBSERVATION(S)
6
WEIGHT(S)
 
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
4         1 cross 
1        26 cross 
RANDOM_RESIDUAL VALUES
245.00    
RANDOM_GROUP
2
RANDOM_TYPE
add_animal   
FILE
gblup_pedi.txt
(CO)VARIANCES
35.250    
OPTION SNP_file gblup_snpgenotype.txt gblup_snpgenotype_XrefID.txt
OPTION no_quality_control
OPTION AlphaBeta 0.99 0.01
OPTION tunedG 0
OPTION saveAscii
OPTION saveG
OPTION saveGInverse
OPTION createGimA22i 0

 

pregsf90 실행 화면

 

pregsf90 실행 로그

pregsf90_gblup.par
     preGSf90 ver. 1.19

 Parameter file:             pregsf90_gblup.par
 Data file:                  gblup_data.txt
 Number of Traits             1
 Number of Effects            2
 Position of Observations      6
 Position of Weight (1)        0
 Value of Missing Trait/Observation           0

EFFECTS
 #  type                position (2)        levels   [positions for nested]
    1  cross-classified       4         1
    2  cross-classified       1        26

 Residual (co)variance Matrix
  245.00    

 Random Effect(s)    2
 Type of Random Effect:      additive animal
 Pedigree File:              gblup_pedi.txt                                                                                                                                                                                                                                            
 trait   effect    (CO)VARIANCES
  1       2     35.25    

 REMARKS
  (1) Weight position 0 means no weights utilized
  (2) Effect positions of 0 for some effects and traits means that such
      effects are missing for specified traits
 
 

 Options read from parameter file for genomic
 * SNP format: BLUPF90 standard (text)
 * SNP file: gblup_snpgenotype.txt
 * SNP Xref file: gblup_snpgenotype_XrefID.txt
 * NOT CreateGimA22i (default=.false.)
 * Save G Matrix (default=.false.)
 * Save G Inverse matrix (default=.false.)
 * No Quality Control Checks !!!!! (default .false.):  T
 * Create a tuned G (default = 2): 0
 * AlphaBeta defaults alpha=0.95, beta=0.05) :  0.99  0.01
 * Matrix in Ascii format(default=binary)
 
 *--------------------------------------------------------------*
 *                 Genomic Library: Dist Version 1.281          *
 *                                                              *
 *             Optimized OpenMP Version -  6 threads            *
 *                                                              *
 *  Modified relationship matrix (H) created for effect:   2    *
 *--------------------------------------------------------------*
 
 Read 26 animals from pedigree file: "D:\users\bhpark\2020\job\공부하기\07_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode\20_gblup for computing SNP effects_blupf90\gblup_pedi.txt"
 Number of Genotyped Animals: 14

 Creating A22 
    Extracting subset of: 26 pedigrees from: 26 ...elapsed time:     0.0000
    Calculating A22 Matrix by Colleau OpenMP...elapsed time: .0004
    Numbers of threads=1 6

 Reading SNP file
    Column position in file for the first marker: 4
    Format to read SNP file: (3x,400000i1)                                     
    Number of SNPs: 50
    Format: integer genotypes (0 to 5) to double-precision array
    Number of Genotyped animals: 14
    Reading SNP file elapsed time: .00
 
 Statistics of alleles frequencies in the current population
    N:             50
    Mean:       0.079
    Min:        0.000
    Max:        0.964
    Var:        0.038
 
 
 Quality Control - Monomorphic SNPs Exist - NOT REMOVED: 40

 Genotypes missings (%):  0.000

 Calculating G Matrix 
    Dgemm MKL #threads=     1    6 Elapsed omp_get_time:     0.0005
 
 Scale by Sum(2pq). Average:   3.53826530612245     

 Detecting samples with similar genotypes
    elapsed time=     0.0
 
 Blend G as alpha*G + beta*A22: (alpha,beta)     0.990     0.010

 Frequency - Diagonal of G
    N:          14
    Mean:        1.131
    Min:         0.668
    Max:         2.267
    Range:       0.080
    Class:     20
 
  #Class       Class   Count
       1  0.6681           1
       2  0.7480           5
       3  0.8280           1
       4  0.9079           0
       5  0.9879           0
       6   1.068           2
       7   1.148           1
       8   1.228           0
       9   1.308           0
      10   1.388           0
      11   1.468           1
      12   1.547           1
      13   1.627           1
      14   1.707           0
      15   1.787           0
      16   1.867           0
      17   1.947           0
      18   2.027           0
      19   2.107           0
      20   2.187           1
      21   2.267           0
 

 Check for diagonal of genomic relationship matrix

 ** High Diagonal of genotype 3  1.63 Not Removed
 ** High Diagonal of genotype 9  2.27 Not Removed
 ** Low Diagonal of genotype 11  0.67 Not Removed

 Check for diagonal of genomic relationship matrix, genotypes not removed: 3

 ------------------------------
  Final Pedigree-Based Matrix 
 ------------------------------
 
 Statistic of Rel. Matrix A22
                            N      Mean       Min       Max       Var
     Diagonal              14     1.000     1.000     1.000     0.000
     Off-diagonal         182     0.148     0.000     0.500     0.032
 

 ----------------------
  Final Genomic Matrix 
 ----------------------
 
 Statistic of Genomic Matrix
                            N      Mean       Min       Max       Var
     Diagonal              14     1.131     0.668     2.267     0.210
     Off-diagonal         182    -0.085    -0.941     1.043     0.203
 

 Correlation of Genomic Inbreeding and Pedigree Inbreeding
     Variance of Y is Zero !!
     Correlation:     0.0000
 
 Diagonal elements
    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =        NaN       NaN

    Correlation diagonal elements G & A       NaN
 
 All elements - Diagonal / Off-Diagonal
    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =     -0.300     1.443

    Correlation all elements G & A     0.733
 
 Off-Diagonal
    Using 84 elements from A22 >= .02000

    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =     -0.215     1.226

    Correlation Off-Diagonal elements G & A     0.361
 
 ***********************************************************************
 * CORRELATION FOR OFF-DIAGONALS G & A22 IS LOW THAN  0.50  !!!!!  *
 * MISIDENTIFIED GENOMIC SAMPLES OR POOR QUALITY GENOMIC DATA *
 ***********************************************************************
 Saving G in file: "G"
       elapsed time=     0.0

 Creating A22-inverse 
    Inverse LAPACK MKL dpotrf/i  #threads=    1    6 Elapsed omp_get_time:     0.0019

 ----------------------
  Final A22 Inv Matrix 
 ----------------------
 
 Statistic of Inv. Rel. Matrix A22
                            N      Mean       Min       Max       Var
     Diagonal              14     1.536     1.000     3.500     0.364
     Off-diagonal         182    -0.082    -0.667     0.000     0.046
 
 
 Creating G-inverse 
    Inverse LAPACK MKL dpotrf/i  #threads=    1    6 Elapsed omp_get_time:     0.0000

 --------------------------
  Final Genomic Inv Matrix 
 --------------------------
 
 Statistic of Inv. Genomic Matrix
                            N      Mean       Min       Max       Var
     Diagonal              14    41.948     3.572    88.271   687.247
     Off-diagonal         182    -0.425   -44.402    39.443   305.004
 
 
 Saving G-inverse in file: "Gi"
       elapsed time=     0.0

*--------------------------------------------------*
* Setup Genomic Done !!!, elapsed time:      0.043 *
*--------------------------------------------------*

 

자세한 설명은 생략. 실제 자료가 아니라 문제가 막 튀어 나오는데 QC를 하지 않으므로 그냥 PASS

pregsf90 실행 결과로 생긴 파일

freqdata.count : allele frequency

         1  0.321429
         2  0.178571
         3  0.357143
         4  0.357143
         5  0.142857
         6  0.607143
         7  0.071429
         8  0.964286
         9  0.571429
        10  0.392857
        11  0.000000
        12  0.000000
        13  0.000000
        14  0.000000
        15  0.000000
        16  0.000000
        17  0.000000
        18  0.000000
        19  0.000000
        20  0.000000
        21  0.000000
        22  0.000000
        23  0.000000
        24  0.000000
        25  0.000000
        26  0.000000
        27  0.000000
        28  0.000000
        29  0.000000
        30  0.000000
        31  0.000000
        32  0.000000
        33  0.000000
        34  0.000000
        35  0.000000
        36  0.000000
        37  0.000000
        38  0.000000
        39  0.000000
        40  0.000000
        41  0.000000
        42  0.000000
        43  0.000000
        44  0.000000
        45  0.000000
        46  0.000000
        47  0.000000
        48  0.000000
        49  0.000000
        50  0.000000

 

sum2pq

   3.53826530612245

 

NRM을 1% 섞은 Genomic Relationship Matrix

1 1 1.467519840782
1 2 -.441110314271
2 2 .748038939835
1 3 .982865906929
2 3 -.920764248235
3 3 1.627404485436
1 4 .061029200220
2 4 -.441110314271
3 4 .423269650637
4 4 .907923584490
1 5 .680582198257
2 5 -.940749828817
3 5 1.042822648674
4 5 .400784070111
5 5 1.587433324273
1 6 -.161312186125
2 6 .183442683655
3 6 -.361167991943
4 6 -.161312186125
5 6 -.101355444379
6 6 .748038939835
1 7 -.700922861835
2 7 .203428264237
3 7 -.620980539508
4 7 .418269650749
5 7 -.361167991943
6 7 .200928264293
7 7 .788010100999
1 8 -.541038217180
2 8 .083514780746
3 8 -.181297766707
4 8 .298356167258
5 8 -.201283347289
6 8 .081014780801
7 8 .383298489473
8 8 .827981262162
1 9 .877938004131
2 9 .098500361439
3 9 .118485942021
4 9 -.521052636598
5 9 -.181297766707
6 9 .098500361439
7 9 -.720908442417
8 9 -.561023797762
9 9 2.266943064056
1 10 -.780865184162
2 10 .403284070055
3 10 -.700922861835
4 10 -.501067056016
5 10 -.441110314271
6 10 -.158812186181
7 10 .140971522547
8 10 -.258740089090
9 10 -.521052636598
10 10 1.187721712636
1 11 -.201283347289
2 11 -.136326605655
3 11 .158457103185
4 11 .358312909003
5 11 .138471522603
6 11 .140971522547
7 11 .160957103129
8 11 .600639875930
9 11 -.221268927871
10 11 -.478581475490
11 11 .668096617507
1 12 -.141326605543
2 12 .483226392383
3 12 -.341182411362
4 12 -.700922861835
5 12 -.640966120089
6 12 -.358667991999
7 12 -.338682411417
8 12 -.178797766763
9 12 .118485942021
10 12 .700567778839
11 12 -.398639153163
12 12 1.067808229145
1 13 -.820836345326
2 13 .363312908891
3 13 -.740894022999
4 13 -.261240089034
5 13 -.481081475435
6 13 .081014780801
7 13 .380798489529
8 13 -.018913122108
9 13 -.561023797762
10 13 .860452423494
11 13 -.238754508508
12 13 .380798489529
13 13 .827981262162
1 14 -.261240089034
2 14 .363312908891
3 14 -.461095894853
4 14 -.261240089034
5 14 -.481081475435
6 14 -.198783347345
7 14 .101000361383
8 14 -.298711250254
9 14 -.281225669616
10 14 .580654295348
11 14 -.518552636654
12 14 .380798489529
13 14 .260885006038
14 14 1.107779390308

NRM을 구성할 때 혈통을 이용

 

NRM을 1% 섞은 GRM의 inverse

1 1 36.374112147756
1 2 -27.980861774302
2 2 88.270926666259
1 3 -8.255548353247
2 3 31.300275470378
3 3 17.520566628931
1 4 -7.273154470752
2 4 -7.935058588596
3 4 -13.839579561710
4 4 59.029844228979
1 5 -16.994409615350
2 5 26.742040202839
3 5 10.825563830254
4 5 1.183435492773
5 5 16.559517365607
1 6 1.865090988759
2 6 -30.838117848755
3 6 -13.961143487346
4 6 35.035830728961
5 6 -3.740681212886
6 6 33.908530891321
1 7 1.034309816800
2 7 -16.074769262009
3 7 10.966959887452
4 7 -44.042226922475
5 7 6.054790042008
6 7 -9.402982983010
7 7 76.598931919425
1 8 23.579731230504
2 8 -3.254062300818
3 8 -4.373575563548
4 8 13.642319923702
5 8 -10.355941361690
6 8 4.610816099161
7 8 -30.643656293298
8 8 38.819147341585
1 9 2.957758933413
2 9 -1.421961970816
3 9 1.150380418263
4 9 4.069157034358
5 9 1.660646618910
6 9 4.195695275841
7 9 2.332803703845
8 9 3.865774059947
9 9 3.571863763864
1 10 39.442736583671
2 10 -44.401612414027
3 10 -15.360569585413
4 10 8.964674694062
5 10 -23.795163924540
6 10 17.932250568059
7 10 -4.681081006798
8 10 31.070735197126
9 10 4.523044730378
10 10 65.241156282393
1 11 15.759887763419
2 11 -27.861655984766
3 11 -3.670240835263
4 11 -16.848814124536
5 11 -5.409628976497
6 11 1.618086379708
7 11 31.963192006762
8 11 -11.691869904844
9 11 2.296427605333
10 11 12.420172214465
11 11 33.289780498461
1 12 -17.655277777373
2 12 -31.453869451043
3 12 -9.610207612292
4 12 24.984206447650
5 12 8.978713862467
6 12 33.085527157687
7 12 24.524531681789
8 12 -26.201393052725
9 12 3.932687423384
10 12 -9.861323926268
11 12 12.723588395898
12 12 65.529756384831
1 13 .709886622910
2 13 38.437384536377
3 13 12.371399723371
4 13 -3.932589146850
5 13 8.462389615927
6 13 -19.286778706349
7 13 -23.250880956380
8 13 11.745624274467
9 13 1.096598339716
10 13 -22.528242543689
11 13 -8.688864312408
12 13 -29.485256595986
13 13 46.430966020324
1 14 5.177116824003
2 14 -6.092801764317
3 14 -1.151050225076
4 14 4.327310182797
5 14 .465043915303
6 14 5.157724634274
7 14 .443047275678
8 14 3.079022874246
9 14 2.906885066300
10 14 2.733882533018
11 14 6.555679888597
12 14 5.120296531692
13 14 4.504947248223
14 14 6.133181220876

 

blupf90 실행 화면

 

blupf90 실행 로그

blupf90_gblup.par
     BLUPF90 ver. 1.68

 Parameter file:             blupf90_gblup.par
 Data file:                  gblup_data.txt
 Number of Traits             1
 Number of Effects            2
 Position of Observations      6
 Position of Weight (1)        0
 Value of Missing Trait/Observation           0

EFFECTS
 #  type                position (2)        levels   [positions for nested]
    1  cross-classified       4         1
    2  cross-classified       8        14

 Residual (co)variance Matrix
  245.00    

 Random Effect(s)    2
 Type of Random Effect:      user defined from file
 User File:                  Gi                                                                                                                                                                                                                                                        
 trait   effect    (CO)VARIANCES
  1       2     35.25    

 REMARKS
  (1) Weight position 0 means no weights utilized
  (2) Effect positions of 0 for some effects and traits means that such
      effects are missing for specified traits
 

 * The limited number of OpenMP threads = 4

 * solving method (default=PCG):FSPAK               
 
 Data record length =            8
 # equations =           15
 G
  35.250    
 read            8  records in   0.1562500      s,                      17 
  nonzeroes
 g_usr_inv: read          105  elements
 largest row, column, diagonal:          14          14          14
 user defined matrix
     36.3741    -27.9809     -8.2555     -7.2732    -16.9944      1.8651      1.0343     23.5797      2.9578     39.4427     15.7599    -17.6553      0.7099      5.1771
    -27.9809     88.2709     31.3003     -7.9351     26.7420    -30.8381    -16.0748     -3.2541     -1.4220    -44.4016    -27.8617    -31.4539     38.4374     -6.0928
     -8.2555     31.3003     17.5206    -13.8396     10.8256    -13.9611     10.9670     -4.3736      1.1504    -15.3606     -3.6702     -9.6102     12.3714     -1.1511
     -7.2732     -7.9351    -13.8396     59.0298      1.1834     35.0358    -44.0422     13.6423      4.0692      8.9647    -16.8488     24.9842     -3.9326      4.3273
    -16.9944     26.7420     10.8256      1.1834     16.5595     -3.7407      6.0548    -10.3559      1.6606    -23.7952     -5.4096      8.9787      8.4624      0.4650
      1.8651    -30.8381    -13.9611     35.0358     -3.7407     33.9085     -9.4030      4.6108      4.1957     17.9323      1.6181     33.0855    -19.2868      5.1577
      1.0343    -16.0748     10.9670    -44.0422      6.0548     -9.4030     76.5989    -30.6437      2.3328     -4.6811     31.9632     24.5245    -23.2509      0.4430
     23.5797     -3.2541     -4.3736     13.6423    -10.3559      4.6108    -30.6437     38.8191      3.8658     31.0707    -11.6919    -26.2014     11.7456      3.0790
      2.9578     -1.4220      1.1504      4.0692      1.6606      4.1957      2.3328      3.8658      3.5719      4.5230      2.2964      3.9327      1.0966      2.9069
     39.4427    -44.4016    -15.3606      8.9647    -23.7952     17.9323     -4.6811     31.0707      4.5230     65.2412     12.4202     -9.8613    -22.5282      2.7339
     15.7599    -27.8617     -3.6702    -16.8488     -5.4096      1.6181     31.9632    -11.6919      2.2964     12.4202     33.2898     12.7236     -8.6889      6.5557
    -17.6553    -31.4539     -9.6102     24.9842      8.9787     33.0855     24.5245    -26.2014      3.9327     -9.8613     12.7236     65.5298    -29.4853      5.1203
      0.7099     38.4374     12.3714     -3.9326      8.4624    -19.2868    -23.2509     11.7456      1.0966    -22.5282     -8.6889    -29.4853     46.4310      4.5049
      5.1771     -6.0928     -1.1511      4.3273      0.4650      5.1577      0.4430      3.0790      2.9069      2.7339      6.5557      5.1203      4.5049      6.1332
 finished peds in   0.1562500      s,                     114  nonzeroes
 left hand side
      0.0327      0.0041      0.0041      0.0041      0.0041      0.0041      0.0041      0.0041      0.0041      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000
      0.0041      1.0360     -0.7938     -0.2342     -0.2063     -0.4821      0.0529      0.0293      0.6689      0.0839      1.1189      0.4471     -0.5009      0.0201      0.1469
      0.0041     -0.7938      2.5082      0.8880     -0.2251      0.7586     -0.8748     -0.4560     -0.0923     -0.0403     -1.2596     -0.7904     -0.8923      1.0904     -0.1728
      0.0041     -0.2342      0.8880      0.5011     -0.3926      0.3071     -0.3961      0.3111     -0.1241      0.0326     -0.4358     -0.1041     -0.2726      0.3510     -0.0327
      0.0041     -0.2063     -0.2251     -0.3926      1.6787      0.0336      0.9939     -1.2494      0.3870      0.1154      0.2543     -0.4780      0.7088     -0.1116      0.1228
      0.0041     -0.4821      0.7586      0.3071      0.0336      0.4739     -0.1061      0.1718     -0.2938      0.0471     -0.6750     -0.1535      0.2547      0.2401      0.0132
      0.0041      0.0529     -0.8748     -0.3961      0.9939     -0.1061      0.9660     -0.2668      0.1308      0.1190      0.5087      0.0459      0.9386     -0.5471      0.1463
      0.0041      0.0293     -0.4560      0.3111     -1.2494      0.1718     -0.2668      2.1771     -0.8693      0.0662     -0.1328      0.9068      0.6957     -0.6596      0.0126
      0.0041      0.6689     -0.0923     -0.1241      0.3870     -0.2938      0.1308     -0.8693      1.1053      0.1097      0.8814     -0.3317     -0.7433      0.3332      0.0873
      0.0000      0.0839     -0.0403      0.0326      0.1154      0.0471      0.1190      0.0662      0.1097      0.1013      0.1283      0.0651      0.1116      0.0311      0.0825
      0.0000      1.1189     -1.2596     -0.4358      0.2543     -0.6750      0.5087     -0.1328      0.8814      0.1283      1.8508      0.3523     -0.2798     -0.6391      0.0776
      0.0000      0.4471     -0.7904     -0.1041     -0.4780     -0.1535      0.0459      0.9068     -0.3317      0.0651      0.3523      0.9444      0.3610     -0.2465      0.1860
      0.0000     -0.5009     -0.8923     -0.2726      0.7088      0.2547      0.9386      0.6957     -0.7433      0.1116     -0.2798      0.3610      1.8590     -0.8365      0.1453
      0.0000      0.0201      1.0904      0.3510     -0.1116      0.2401     -0.5471     -0.6596      0.3332      0.0311     -0.6391     -0.2465     -0.8365      1.3172      0.1278
      0.0000      0.1469     -0.1728     -0.0327      0.1228      0.0132      0.1463      0.0126      0.0873      0.0825      0.0776      0.1860      0.1453      0.1278      0.1740
 right hand side:
    0.32    0.04    0.05    0.05    0.06    0.02    0.03    0.04    0.02    0.00
    0.00    0.00    0.00    0.00    0.00
 solution:
    9.94    0.07    0.11    0.05    0.26   -0.49   -0.36    0.14   -0.23    0.03
    0.11   -0.24    0.14    0.05    0.35
 solutions stored in file: "solutions"

 

blupf90 실행 결과 : solutions

trait/effect level  solution
   1   1         1          9.94327493
   1   2         1          0.06981142
   1   2         2          0.11059485
   1   2         3          0.04780743
   1   2         4          0.25776306
   1   2         5         -0.49280367
   1   2         6         -0.35505536
   1   2         7          0.14266450
   1   2         8         -0.22698269
   1   2         9          0.02673761
   1   2        10          0.11379233
   1   2        11         -0.23776246
   1   2        12          0.14174468
   1   2        13          0.05374991
   1   2        14          0.35057024

 

gblup for computing SNP effects_blupf90.zip
0.00MB

저작자표시

'Animal Breeding > BLUPF90' 카테고리의 다른 글

blupf90으로 single-step GBLUP Model 풀기  (0) 2021.01.05
blupf90으로 GBLUP Model with residual polygenic effects 풀기  (0) 2020.12.29
blupf90으로 SNP-BLUP Model(with polygenic effects unweighted analysis) 풀기  (0) 2020.12.24
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기  (0) 2020.12.19
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기  (0) 2020.12.19

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.

# Raphael Mrode

# Example 11.5 p189

간단한 모델 설명은 다음 포스팅을 참고한다.

2020/12/24 - [Animal Breeding/R for Genetic Evaluation] - SNP-BLUP Model with Polygenic Effects(unweighted analysis)

 

SNP-BLUP Model with Polygenic Effects(unweighted analysis)

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.5 p189 SNP effect로 모든 상가적 유전 분산을 설명할 수 없을 수도 있다. 그래서 설명..

bhpark.tistory.com

 

Data

13 0 0 1 558 9 0.00179211 1.3571429 -0.3571429 0.2857143 0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
14 0 0 1 722 13.4 0.00138504 0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143
15 13 4 1 300 12.7 0.00333333 0.3571429 0.6428571 1.2857143 0.2857143 0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
16 15 2 1 73 15.4 0.01369863 -0.6428571 -0.3571429 1.2857143 0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857
17 15 5 1 52 5.9 0.01923077 -0.6428571 0.6428571 0.2857143 1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857
18 14 6 1 87 7.7 0.01149425 0.3571429 0.6428571 -0.7142857 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 0.2142857
19 14 9 1 64 10.2 0.01562500 -0.6428571 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 -0.7857143
20 14 9 1 103 4.8 0.00970874 -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 -0.7857143

 

Pedigree

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

 

Renumf90을 위한 파라미터 파일

# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
snp_data2.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
 
WEIGHT(S)
 
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
1 cross alpha
RANDOM
animal
FILE
snp_pedi2.txt
FILE_POS
1 2 3
PED_DEPTH
0
(CO)VARIANCES
3.5241
EFFECT
4 cross alpha
RANDOM
diagonal
RANDOM_REGRESSION
data
RR_POSITION
8 9 10 11 12 13 14 15 16 17
(CO)VARIANCES
8.963969 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 8.963969 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 8.963969 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 8.963969 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 8.963969 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 8.963969 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 8.963969 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 8.963969 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 8.963969 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 8.963969
OPTION solv_method FSPAK

 

Renumf90 실행 화면

 

renumf90 실행 로그

 RENUMF90 version 1.150
 renumf90_snpblup_polygenic_uw.par
 datafile:snp_data2.txt
 traits:           6
 R
   245.0    

 Processing effect  1 of type cross     
 item_kind=alpha     

 Processing effect  2 of type cross     
 item_kind=alpha     
 pedigree file name  "snp_pedi2.txt"
 positions of animal, sire, dam, alternate dam, yob, and group     1     2     3     0     0     0     0
 all pedigrees to be included
 Reading (CO)VARIANCES:           1 x           1

 Processing effect  3 of type cross     
 item_kind=alpha     
 Reading (CO)VARIANCES:          10 x          10

 Maximum size of character fields: 20

 Maximum size of record (max_string_readline): 800

 Maximum number of fields for input file (max_field_readline): 100

 Pedigree search method (ped_search): convention

 Order of pedigree animals (animal_order): default

 Order of UPG (upg_order): default

 Missing observation code (missing): 0

 hash tables for effects set up
 first 3 lines of the data file (up to 70 characters)
    13 0 0 1 558 9 0.00179211 1.3571429 -0.3571429 0.2857143 0.2857143 -0.
    14 0 0 1 722 13.4 0.00138504 0.3571429 -0.3571429 -0.7142857 -0.714285
    15 13 4 1 300 12.7 0.00333333 0.3571429 0.6428571 1.2857143 0.2857143
 read            8  records
 table with            1  elements sorted
 added count
 Effect group            1  of column            1  with            1  levels
 table expanded from        10000  to        10000  records
 added count
 Effect group            2  of column            1  with            8  levels
 table with            1  elements sorted
 added count
 Effect group            3  of column            1  with            1  levels
 table expanded from        10000  to        10000  records
 wrote statistics in file "renf90.tables"

 Basic statistics for input data  (missing value code is '0')
 Pos  Min         Max         Mean        SD                 N
   6    4.8000      15.400      9.8875      3.7434           8

 random effect   2
 type:animal    
 opened output pedigree file "renadd02.ped"
 read           26  pedigree records
 loaded           18  parent(s) in round            0

 Pedigree checks
 
 Number of animals with records                  =            8
 Number of parents without records               =           18
 Total number of animals                         =           26

 random effect   3
 type:diag      

 Wrote parameter file "renf90.par"
 Wrote renumbered data "renf90.dat" 8 records
 Wrote field information "renf90.fields" for 14 fields in data

 

renumf90 실행으로 생긴 파일

renadd02.ped

26 3 20 1 0 2 0 0 0 26
1 0 0 3 0 0 1 1 0 13
21 9 11 1 0 2 0 0 0 21
13 0 0 3 0 0 0 0 1 5
2 3 17 1 0 2 1 0 0 19
3 0 0 3 0 0 1 8 0 14
22 3 16 1 0 2 0 0 0 22
11 0 0 3 0 0 0 0 1 3
16 0 0 3 0 0 0 0 1 8
4 1 12 1 0 2 1 2 0 15
23 3 19 1 0 2 0 0 0 23
18 0 0 3 0 0 0 0 1 10
9 0 0 3 0 0 0 1 0 1
14 0 0 3 0 0 0 0 1 6
5 4 10 1 0 2 1 0 0 16
24 3 18 1 0 2 0 0 0 24
19 0 0 3 0 0 0 0 1 11
12 0 0 3 0 0 0 0 1 4
6 4 13 1 0 2 1 0 0 17
17 0 0 3 0 0 0 0 2 9
25 3 15 1 0 2 0 0 0 25
20 0 0 3 0 0 0 0 1 12
7 3 17 1 0 2 1 0 0 20
10 0 0 3 0 0 0 0 1 2
8 3 14 1 0 2 1 0 0 18
15 0 0 3 0 0 0 0 1 7

renf90.dat

 9 1.3571429 -0.3571429 0.2857143 0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1 1
 13.4 0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143 1 3 1
 12.7 0.3571429 0.6428571 1.2857143 0.2857143 0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 4 1
 15.4 -0.6428571 -0.3571429 1.2857143 0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857 1 5 1
 5.9 -0.6428571 0.6428571 0.2857143 1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 6 1
 7.7 0.3571429 0.6428571 -0.7142857 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 0.2142857 1 8 1
 10.2 -0.6428571 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 -0.7857143 1 2 1
 4.8 -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 -0.7857143 1 7 1

 

renf90.par

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
          12
OBSERVATION(S)
    1
WEIGHT(S)
 
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 12         1 cross 
 13        26 cross 
 2          1 cov 14
 3          1 cov 14
 4          1 cov 14
 5          1 cov 14
 6          1 cov 14
 7          1 cov 14
 8          1 cov 14
 9          1 cov 14
 10          1 cov 14
 11          1 cov 14
RANDOM_RESIDUAL VALUES
   245.00    
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_animal   
 FILE
renadd02.ped                                                                    
(CO)VARIANCES
   3.5241    
 RANDOM_GROUP
     3     4     5     6     7     8     9    10    11    12
 RANDOM_TYPE
 diagonal     
 FILE
                                                                                
(CO)VARIANCES
   8.9640       0.0000       0.0000       0.0000       0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       8.9640       0.0000       0.0000       0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       8.9640       0.0000       0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000       8.9640       0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000       0.0000       8.9640       0.0000       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000       0.0000       0.0000       8.9640       0.0000    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000       0.0000       0.0000       0.0000       8.9640    
   0.0000       0.0000       0.0000    
   0.0000       0.0000       0.0000       0.0000       0.0000       0.0000       0.0000    
   8.9640       0.0000       0.0000    
   0.0000       0.0000       0.0000       0.0000       0.0000       0.0000       0.0000    
   0.0000       8.9640       0.0000    
   0.0000       0.0000       0.0000       0.0000       0.0000       0.0000       0.0000    
   0.0000       0.0000       8.9640    
OPTION solv_method FSPAK

 

renf90.fields

field     variable  origfield     group    column random    effect     file
    1        trait          6         0         0 *         cov        *
    2    covariate          8         0         0 *         cov        *
    3    covariate          9         0         0 *         cov        *
    4    covariate         10         0         0 *         cov        *
    5    covariate         11         0         0 *         cov        *
    6    covariate         12         0         0 *         cov        *
    7    covariate         13         0         0 *         cov        *
    8    covariate         14         0         0 *         cov        *
    9    covariate         15         0         0 *         cov        *
   10    covariate         16         0         0 *         cov        *
   11    covariate         17         0         0 *         cov        *
   12   renumbered          4         1         1 *         cross      *                                                                               
   13   renumbered          1         2         1 animal    cross      renadd02.ped                                                                    
   14   renumbered          4         3         1 diagonal  cross      *

 

blupf90 실행 화면

 

blupf90 실행 로그

renf90.par
     BLUPF90 ver. 1.70

 Parameter file:             renf90.par
 Data file:                  renf90.dat
 Number of Traits             1
 Number of Effects           12
 Position of Observations      1
 Position of Weight (1)        0
 Value of Missing Trait/Observation           0

EFFECTS
 #  type                position (2)        levels   [positions for nested]
    1  cross-classified      12         1
    2  cross-classified      13        26
    3  covariable             2         1    14
    4  covariable             3         1    14
    5  covariable             4         1    14
    6  covariable             5         1    14
    7  covariable             6         1    14
    8  covariable             7         1    14
    9  covariable             8         1    14
   10  covariable             9         1    14
   11  covariable            10         1    14
   12  covariable            11         1    14

 Residual (co)variance Matrix
  245.00    

 Random Effect(s)    2
 Type of Random Effect:      additive animal
 Pedigree File:              renadd02.ped                                                                                                                                                                                                                                              
 trait   effect    (CO)VARIANCES
  1       2     3.524    

 correlated random effects     3  4  5  6  7  8  9 10 11 12
 Type of Random Effect:      diagonal
 trait   effect    (CO)VARIANCES
  1       3     8.964       0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000    
  1       4     0.000       8.964       0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000    
  1       5     0.000       0.000       8.964       0.000       0.000       0.000       0.000       0.000       0.000       0.000    
  1       6     0.000       0.000       0.000       8.964       0.000       0.000       0.000       0.000       0.000       0.000    
  1       7     0.000       0.000       0.000       0.000       8.964       0.000       0.000       0.000       0.000       0.000    
  1       8     0.000       0.000       0.000       0.000       0.000       8.964       0.000       0.000       0.000       0.000    
  1       9     0.000       0.000       0.000       0.000       0.000       0.000       8.964       0.000       0.000       0.000    
  1      10     0.000       0.000       0.000       0.000       0.000       0.000       0.000       8.964       0.000       0.000    
  1      11     0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000       8.964       0.000    
  1      12     0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000       0.000       8.964    

 REMARKS
  (1) Weight position 0 means no weights utilized
  (2) Effect positions of 0 for some effects and traits means that such
      effects are missing for specified traits
 

 * The limited number of OpenMP threads = 4

 * solving method (default=PCG):FSPAK               
 
 Data record length =           14
 # equations =           37
 G
  3.5241    
 G
  8.9640      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      8.9640      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      8.9640      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      8.9640      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      8.9640      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      8.9640      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      8.9640      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      8.9640      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      8.9640    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  8.9640    
 read            8  records in   0.2187500      s,                     162 
  nonzeroes
  read           26  additive pedigrees
 finished peds in   0.2187500      s,                     215  nonzeroes
 solutions stored in file: "solutions"

 

blupf90 실행 결과 : solutions

trait/effect level  solution
   1   1         1          9.93695036
   1   2         1          0.01067048
   1   2         2         -0.01577057
   1   2         3          0.00049787
   1   2         4          0.04298274
   1   2         5          0.07718088
   1   2         6         -0.01695885
   1   2         7         -0.05195194
   1   2         8         -0.02022509
   1   2         9          0.00000000
   1   2        10          0.03712634
   1   2        11          0.00000000
   1   2        12          0.02509833
   1   2        13         -0.02563348
   1   2        14         -0.01364935
   1   2        15         -0.00000000
   1   2        16         -0.00000000
   1   2        17         -0.03411019
   1   2        18         -0.00000000
   1   2        19         -0.00000000
   1   2        20         -0.00000000
   1   2        21          0.00000000
   1   2        22          0.00024894
   1   2        23          0.00024894
   1   2        24          0.00024894
   1   2        25          0.00024894
   1   2        26          0.00024894
   1   3         1          0.07836105
   1   4         1         -0.28018933
   1   5         1          0.23400811
   1   6         1         -0.07435738
   1   7         1          0.09844808
   1   8         1          0.12723456
   1   9         1          0.00000000
   1  10         1          0.00000000
   1  11         1         -0.05409385
   1  12         1         -0.01787671

 

첫째 효과는 general mean, 둘째 효과는 개체의 residual polygenic effects, 3 ~ 12까지의 효과는 10개의 snp effect이다. 

snpblup_polygenic_uw_blupf90.zip
0.00MB

 

저작자표시

'Animal Breeding > BLUPF90' 카테고리의 다른 글

blupf90으로 GBLUP Model with residual polygenic effects 풀기  (0) 2020.12.29
blupf90으로 gblup model 풀기  (0) 2020.12.29
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기  (0) 2020.12.19
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기  (0) 2020.12.19
blupf90으로 Fixed Effect Model for SNP Effects(weighted analysis) 풀기  (0) 2020.12.17

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.

# Raphael Mrode

# Example 11.5 p189

SNP effect로 모든 상가적 유전 분산을 설명할 수 없을 수도 있다. 그래서 설명하지 못하는 부분을 위하여 polygenic effect를 넣는다. 여기서는 polygenic effect의 분산을 상가적 유전분산의 10%인 3.5241로 하였다. 그래서 MME에 들어가는 polygenic effect의 분산 비율 alpha1 = 245 / 3.5241 = 69.521이다. MME에 들어가는 SNP effect의 분산 비율 alpha2를 계산하기 위하여 먼저 SNP effect가 설명하는 분산을 계산한다. 35.241 - 3.5241 = 31.7169이다. 개개의 SNP의 분산은 sum of 2pq로 나눈다. 즉 31.7169 / 3.5383 = 8.964이다. 마지막으로 alpha2 = 245 / 8.964 = 27.332이다.

 

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

 

Progmam

# SNP-BLUP Model with residual polygenic effect, unweighted analysis

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.

# Raphael Mrode

# Example 11.5 p189

# MASS packages
if (!("MASS" %in% installed.packages())) {
  install.packages("MASS", dependencies = TRUE)
}
library(MASS)

# function to make inverse of numerator relationship matrix
ainv = function(pedi) {
  n = nrow(pedi)
  Ainv = 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) {
      # both parents unknown
      alpha = 1
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
    } else if (sire != 0 & dam == 0) {
      # sire known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
    } else if (sire == 0 & dam != 0) {
      # dam known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    } else {
      # both parents known
      alpha = 2
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
      Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
      Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    }
  }
  return(Ainv)
}


# 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)
#s2pq = 0
#for (i in 1:10) {
#  s2pq = s2pq + 2 * snp_mean[i] * (1 - snp_mean[i])
#}
sum2pq = sum(2 * snp_mean * (1 - snp_mean))
sum2pq

# ratio of polygenic effect
rpe = 0.1

# variance of residual polygenic effect
sigma_u = sigma_a * rpe
sigma_u

# variance ratio, alpha
alpha1 = sigma_e / sigma_u
alpha1

alpha2 =  sum2pq * sigma_e / (sigma_a - sigma_u)
alpha2

# 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 polygenic effect
W = cbind(matrix(rep(0, 8 * 12), 8, 12), diag(8), matrix(rep(0, 8 * 6), 8, 6))
W

# inverse matrix of NRM
ai = ainv(pedi)
ai

# design matrix for SNP effect
M = as.matrix(data[1:8, 7:16])
Z = sweep(M, 2, snp_mean * 2)
round(Z, 3)

# Left Hand Side
LHS = rbind(
  cbind(t(X) %*% X, t(X) %*% W,               t(X) %*% Z), 
  cbind(t(W) %*% X, t(W) %*% W + ai * alpha1, t(W) %*% Z),
  cbind(t(Z) %*% X, t(Z) %*% W,               t(Z) %*% Z + diag(10) * alpha2))
round(LHS, 3)

# Right Hand Side
RHS = rbind(t(X) %*% y, 
            t(W) %*% y,
            t(Z) %*% y)
RHS

# Solutions

# generalized inverse of LHS
gi_LHS = ginv(LHS)
gi_LHS

# solution
sol = gi_LHS %*% RHS
round(sol, 3)

# mean effect of solutions
mean_eff = sol[1]
round(mean_eff, 3)

# snp effect of solutions
snp_eff = sol[28:37]
round(snp_eff, 3)

# calculate DGV
# genotype of animals
# design matrix of animals
P2 = sweep(M_all, 2, snp_mean * 2)
round(P2, 3)

# DGV
dgv = P2 %*% snp_eff
round(dgv, 3)

 

프로그램을 RStudio에서 실행한 화면

 

프로그램 실행 결과

> # SNP-BLUP Model with residual polygenic effect, unweighted analysis
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
> 
> # Raphael Mrode
> 
> # Example 11.5 p189
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+   install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # function to make inverse of numerator relationship matrix
> ainv = function(pedi) {
+   n = nrow(pedi)
+   Ainv = 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) {
+       # both parents unknown
+       alpha = 1
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+     } else if (sire != 0 & dam == 0) {
+       # sire known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+     } else if (sire == 0 & dam != 0) {
+       # dam known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     } else {
+       # both parents known
+       alpha = 2
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+       Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+       Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     }
+   }
+   return(Ainv)
+ }
> 
> 
> # 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)
> #s2pq = 0
> #for (i in 1:10) {
> #  s2pq = s2pq + 2 * snp_mean[i] * (1 - snp_mean[i])
> #}
> sum2pq = sum(2 * snp_mean * (1 - snp_mean))
> sum2pq
[1] 3.538265
> 
> # ratio of polygenic effect
> rpe = 0.1
> 
> # variance of residual polygenic effect
> sigma_u = sigma_a * rpe
> sigma_u
[1] 3.5241
> 
> # variance ratio, alpha
> alpha1 = sigma_e / sigma_u
> alpha1
[1] 69.5213
> 
> alpha2 =  sum2pq * sigma_e / (sigma_a - sigma_u)
> alpha2
[1] 27.33164
> 
> # 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 polygenic effect
> W = cbind(matrix(rep(0, 8 * 12), 8, 12), diag(8), matrix(rep(0, 8 * 6), 8, 6))
> W
     [,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,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0
[2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
[3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0
[4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0
[5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0
[6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0
[7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0
[8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0
> 
> # inverse matrix of NRM
> ai = ainv(pedi)
> ai
      [,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.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [2,]  0.0  1.5  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5    -1     0     0     0     0     0     0     0     0     0     0
 [3,]  0.5  0.0  1.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [4,]  0.0  0.0  0.0  1.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
 [5,]  0.0  0.0  0.0  0.0  1.5  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5     0    -1     0     0     0     0     0     0     0     0     0
 [6,]  0.0  0.0  0.0  0.0  0.0  1.5  0.0  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0    -1     0     0     0     0     0     0     0     0
 [7,]  0.0  0.0  0.0  0.0  0.0  0.0  1.5  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0    -1     0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.5    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0    -1     0     0     0     0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    2   0.0   0.0   0.0   0.0   1.0   0.0     0     0     0    -1    -1     0     0     0     0     0     0
[10,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   1.5   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0    -1     0     0
[11,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   1.5   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0    -1     0     0     0
[12,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   1.5   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0     0    -1
[13,]  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   1.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
[14,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5    1   0.5   0.5   0.5   0.0   5.0   0.0     0     0    -1    -1    -1     0    -1    -1    -1    -1    -1
[15,]  0.0  0.5  0.0 -1.0  0.5  0.0  0.0  0.0    0   0.0   0.0   0.0  -1.0   0.0   3.0    -1    -1     0     0     0     0     0     0     0     0     0
[16,]  0.0 -1.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     2     0     0     0     0     0     0     0     0     0     0
[17,]  0.0  0.0  0.0  0.0 -1.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     0     2     0     0     0     0     0     0     0     0     0
[18,]  0.0  0.0  0.0  0.0  0.0 -1.0  0.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     2     0     0     0     0     0     0     0     0
[19,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     2     0     0     0     0     0     0     0
[20,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     2     0     0     0     0     0     0
[21,] -1.0  0.0 -1.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0     2     0     0     0     0     0
[22,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     2     0     0     0     0
[23,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0  -1.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     2     0     0     0
[24,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0  -1.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     2     0     0
[25,]  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     2     0
[26,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0  -1.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     0     2
> 
> # design matrix for SNP effect
> M = as.matrix(data[1:8, 7:16])
> Z = sweep(M, 2, snp_mean * 2)
> round(Z, 3)
    snp1   snp2   snp3   snp4   snp5   snp6   snp7  snp8   snp9  snp10
1  1.357 -0.357  0.286  0.286 -0.286 -1.214 -0.143 0.071 -0.143  1.214
2  0.357 -0.357 -0.714 -0.714 -0.286  0.786 -0.143 0.071 -0.143 -0.786
3  0.357  0.643  1.286  0.286  0.714 -1.214 -0.143 0.071 -0.143  1.214
4 -0.643 -0.357  1.286  0.286 -0.286 -0.214 -0.143 0.071  0.857  0.214
5 -0.643  0.643  0.286  1.286 -0.286 -1.214 -0.143 0.071 -0.143  1.214
6  0.357  0.643 -0.714  0.286 -0.286  0.786 -0.143 0.071  0.857  0.214
7 -0.643 -0.357  0.286  0.286 -0.286  0.786 -0.143 0.071  0.857 -0.786
8 -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143 0.071  0.857 -0.786
> 
> # Left Hand Side
> LHS = rbind(
+   cbind(t(X) %*% X, t(X) %*% W,               t(X) %*% Z), 
+   cbind(t(W) %*% X, t(W) %*% W + ai * alpha1, t(W) %*% Z),
+   cbind(t(Z) %*% X, t(Z) %*% W,               t(Z) %*% Z + diag(10) * alpha2))
> round(LHS, 3)
                                                                                                                                                                                            
       8.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   0.000   0.000
       0.000 104.282   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000
       0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761 -69.521   0.000   0.000   0.000   0.000   0.000   0.000
       0.000  34.761   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000
       0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000 -69.521   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043   0.000   0.000   0.000   0.000  69.521   0.000   0.000   0.000   0.000 -69.521 -69.521   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 104.282   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       1.000   0.000   0.000   0.000  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 105.282   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       1.000   0.000   0.000   0.000   0.000   0.000  34.761  34.761  34.761  69.521  34.761  34.761  34.761   0.000 348.606   0.000   0.000   0.000 -69.521 -69.521 -69.521   0.000 -69.521
       1.000   0.000  34.761   0.000 -69.521  34.761   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000 209.564 -69.521 -69.521   0.000   0.000   0.000   0.000   0.000
       1.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521 140.043   0.000   0.000   0.000   0.000   0.000   0.000
       1.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000 140.043   0.000   0.000   0.000   0.000   0.000
       1.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000 140.043   0.000   0.000   0.000   0.000
       1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 140.043   0.000   0.000   0.000
       1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000 140.043   0.000   0.000
       0.000 -69.521   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000 139.043
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
       0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 -69.521   0.000 -69.521   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
                                        snp1   snp2   snp3   snp4   snp5   snp6   snp7   snp8   snp9  snp10
        0.000   0.000   0.000   0.000 -0.143  1.143  2.286  1.286 -1.286 -1.714 -1.143  0.571  2.857  1.714
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000 -69.521   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000 -69.521   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
      -69.521   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000 -69.521  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  1.357 -0.357  0.286  0.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
      -69.521 -69.521 -69.521 -69.521  0.357 -0.357 -0.714 -0.714 -0.286  0.786 -0.143  0.071 -0.143 -0.786
        0.000   0.000   0.000   0.000  0.357  0.643  1.286  0.286  0.714 -1.214 -0.143  0.071 -0.143  1.214
        0.000   0.000   0.000   0.000 -0.643 -0.357  1.286  0.286 -0.286 -0.214 -0.143  0.071  0.857  0.214
        0.000   0.000   0.000   0.000 -0.643  0.643  0.286  1.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
        0.000   0.000   0.000   0.000  0.357  0.643 -0.714  0.286 -0.286  0.786 -0.143  0.071  0.857  0.214
        0.000   0.000   0.000   0.000 -0.643 -0.357  0.286  0.286 -0.286  0.786 -0.143  0.071  0.857 -0.786
        0.000   0.000   0.000   0.000 -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143  0.071  0.857 -0.786
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
      139.043   0.000   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000 139.043   0.000   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000 139.043   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
        0.000   0.000   0.000 139.043  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [ getOption("max.print") 에 도달했습니다 -- 10 행들을 생략합니다 ]
> 
> # Right Hand Side
> RHS = rbind(t(X) %*% y, 
+             t(W) %*% y,
+             t(Z) %*% y)
> RHS
        [,1]
       79.10
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
        9.00
       13.40
       12.70
       15.40
        5.90
        7.70
       10.20
        4.80
        0.00
        0.00
        0.00
        0.00
        0.00
        0.00
snp1    0.95
snp2    2.85
snp3   29.60
snp4   10.30
snp5   -9.90
snp6  -13.25
snp7  -11.30
snp8    5.65
snp9   26.80
snp10  16.15
> 
> # 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]  [,16]  [,17]  [,18]  [,19]  [,20]  [,21] [,22]  [,23]  [,24]  [,25]  [,26]  [,27]
 [1,]  0.142 0.000 -0.001 0.000 -0.002 -0.001 -0.001 0.000 0.000 -0.002 0.000 0.000 0.000 -0.003 -0.005 -0.004 -0.003 -0.003 -0.004 -0.004 -0.004 0.000 -0.002 -0.002 -0.002 -0.002 -0.002
 [2,]  0.000 0.014  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.007  0.000  0.000  0.000  0.000  0.000
 [3,] -0.001 0.000  0.014 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.007  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [4,]  0.000 0.000  0.000 0.014  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.007  0.000  0.000  0.000  0.000  0.000
 [5,] -0.002 0.000  0.000 0.000  0.014  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.007  0.004  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [6,] -0.001 0.000  0.000 0.000  0.000  0.014  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.007  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [7,] -0.001 0.000  0.000 0.000  0.000  0.000  0.014 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.007  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
 [8,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.014 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.007  0.000
 [9,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.014  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.007  0.000  0.000  0.000  0.000
[10,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.014 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.007  0.007 0.000  0.000  0.000  0.000  0.000  0.000
[11,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.014 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.007  0.000  0.000
[12,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.014 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.007  0.000  0.000  0.000
[13,]  0.000 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.014  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.007
[14,] -0.003 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.014  0.000  0.007  0.004  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[15,] -0.005 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.014  0.000  0.000  0.000  0.007  0.007  0.007 0.000  0.007  0.007  0.007  0.007  0.007
[16,] -0.004 0.000  0.000 0.000  0.007  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.007  0.000  0.014  0.007  0.007  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[17,] -0.003 0.000  0.007 0.000  0.004  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.004  0.000  0.007  0.014  0.004  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[18,] -0.003 0.000  0.000 0.000  0.004  0.007  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.004  0.000  0.007  0.004  0.014  0.000  0.000  0.000 0.000  0.000  0.000  0.000  0.000  0.000
[19,] -0.004 0.000  0.000 0.000  0.000  0.000  0.007 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.014  0.004  0.004 0.000  0.004  0.004  0.004  0.004  0.004
[20,] -0.004 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.007 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.014  0.007 0.000  0.004  0.004  0.004  0.004  0.004
[21,] -0.004 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.007 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.007  0.014 0.000  0.004  0.004  0.004  0.004  0.004
[22,]  0.000 0.007  0.000 0.007  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.014  0.000  0.000  0.000  0.000  0.000
[23,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.007  0.000 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.014  0.004  0.004  0.004  0.004
[24,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.007 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.004  0.014  0.004  0.004  0.004
[25,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.007 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.004  0.004  0.014  0.004  0.004
[26,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.007 0.000  0.000 0.000 0.000 0.000  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.004  0.004  0.004  0.014  0.004
[27,] -0.002 0.000  0.000 0.000  0.000  0.000  0.000 0.000 0.000  0.000 0.000 0.000 0.007  0.000  0.007  0.000  0.000  0.000  0.004  0.004  0.004 0.000  0.004  0.004  0.004  0.004  0.014
       [,28]  [,29]  [,30]  [,31]  [,32]  [,33] [,34]  [,35]  [,36]  [,37]
 [1,]  0.000 -0.005 -0.008 -0.004  0.006  0.005 0.005 -0.003 -0.013 -0.006
 [2,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [3,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [4,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [5,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [6,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [7,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [8,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [9,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[10,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[11,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[12,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[13,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[14,] -0.001  0.000  0.000  0.000  0.000  0.001 0.000  0.000  0.000 -0.001
[15,]  0.000  0.000  0.001  0.000  0.000 -0.001 0.000  0.000  0.000  0.001
[16,]  0.000  0.000 -0.001  0.000  0.000  0.001 0.000  0.000  0.000 -0.001
[17,]  0.000  0.000 -0.001  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[18,]  0.000  0.000  0.000  0.000  0.000  0.001 0.000  0.000  0.000 -0.001
[19,]  0.000  0.000  0.001  0.000  0.000 -0.001 0.000  0.000  0.000  0.000
[20,]  0.000  0.000  0.000  0.000  0.000 -0.001 0.000  0.000  0.000  0.001
[21,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.001
[22,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[23,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[24,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[25,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[26,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
[27,]  0.000  0.000  0.000  0.000  0.000  0.000 0.000  0.000  0.000  0.000
 [ getOption("max.print") 에 도달했습니다 -- 10 행들을 생략합니다 ]
> 
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
        [,1]
 [1,]  9.937
 [2,]  0.000
 [3,]  0.037
 [4,]  0.000
 [5,]  0.025
 [6,] -0.026
 [7,] -0.014
 [8,]  0.000
 [9,]  0.000
[10,] -0.034
[11,]  0.000
[12,]  0.000
[13,]  0.000
[14,]  0.011
[15,]  0.000
[16,]  0.043
[17,]  0.077
[18,] -0.017
[19,] -0.020
[20,] -0.016
[21,] -0.052
[22,]  0.000
[23,]  0.000
[24,]  0.000
[25,]  0.000
[26,]  0.000
[27,]  0.000
[28,]  0.078
[29,] -0.280
[30,]  0.234
[31,] -0.074
[32,]  0.098
[33,]  0.127
[34,]  0.000
[35,]  0.000
[36,] -0.054
[37,] -0.018
> 
> # mean effect of solutions
> mean_eff = sol[1]
> round(mean_eff, 3)
[1] 9.937
> 
> # snp effect of solutions
> snp_eff = sol[28:37]
> round(snp_eff, 3)
 [1]  0.078 -0.280  0.234 -0.074  0.098  0.127  0.000  0.000 -0.054 -0.018
> 
> # calculate DGV
> # genotype of animals
> # design matrix of animals
> P2 = sweep(M_all, 2, snp_mean * 2)
> round(P2, 3)
        snp1   snp2   snp3   snp4   snp5   snp6   snp7   snp8   snp9  snp10
 [1,]  1.357 -0.357  0.286  0.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
 [2,]  0.357 -0.357 -0.714 -0.714 -0.286  0.786 -0.143  0.071 -0.143 -0.786
 [3,]  0.357  0.643  1.286  0.286  0.714 -1.214 -0.143  0.071 -0.143  1.214
 [4,] -0.643 -0.357  1.286  0.286 -0.286 -0.214 -0.143  0.071  0.857  0.214
 [5,] -0.643  0.643  0.286  1.286 -0.286 -1.214 -0.143  0.071 -0.143  1.214
 [6,]  0.357  0.643 -0.714  0.286 -0.286  0.786 -0.143  0.071  0.857  0.214
 [7,] -0.643 -0.357  0.286  0.286 -0.286  0.786 -0.143  0.071  0.857 -0.786
 [8,] -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143  0.071  0.857 -0.786
 [9,]  1.357 -0.357 -0.714 -0.714 -0.286 -0.214  1.857  0.071 -0.143  1.214
[10,] -0.643 -0.357 -0.714  0.286  0.714  0.786 -0.143  0.071 -1.143 -0.786
[11,] -0.643  0.643  0.286 -0.714 -0.286 -0.214 -0.143  0.071  0.857  0.214
[12,]  0.357 -0.357 -0.714 -0.714  0.714 -0.214 -0.143  0.071 -1.143 -0.786
[13,] -0.643 -0.357 -0.714  0.286  0.714  0.786 -0.143  0.071 -0.143 -0.786
[14,]  0.357 -0.357  0.286  0.286 -0.286  0.786 -0.143 -0.929 -1.143 -0.786
> 
> # DGV
> dgv = P2 %*% snp_eff
> round(dgv, 3)
        [,1]
 [1,]  0.055
 [2,]  0.108
 [3,]  0.029
 [4,]  0.224
 [5,] -0.456
 [6,] -0.319
 [7,]  0.135
 [8,] -0.198
 [9,]  0.023
[10,]  0.107
[11,] -0.216
[12,]  0.133
[13,]  0.053
[14,]  0.321
>

 

책의 결과 값과 SNP effect 값은 동일하나 polygenic effect가 다르게 나왔다. blupf90으로도 확인을 했는데 역시 다르다. 저자의 오타인가 보다. DGV도 다르다.

 

자료, 혈통, 프로그램 파일 : 

21_SNPBLUP with polygenic effect_unweighted.zip
0.00MB

 

저작자표시

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

Single-step GBLUP(ssGBLUP)  (0) 2021.01.05
GBLUP Model with residual polygenic effects  (0) 2020.12.29
Mixed Linear Model(GBLUP model) for Computing DGV  (0) 2020.12.21
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

+ Recent posts

티스토리툴바