Let's Study Animal Breeding

전체 글

Mixed Linear Model(GBLUP model) for Computing DGV 2020.12.21
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기 2020.12.19
blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기 2020.12.19
Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis 2020.12.19
Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis 2020.12.19

Mixed Linear Model(GBLUP model) for Computing DGV

투정이 2020. 12. 21. 20:41

2020. 12. 21. 20:41

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

# Raphael Mrode

# Example 11.3 p18

앞에서는 SNP effect를 구한 후, 다시 DGV를 구하였는데 여기서는 GRM(Genomic Relationship Matrix)를 이용하여 DGV를 바로 구하는 방법에 관하여 설명을 한다.

SNP genotype을 이용하여 NRM(Numerator Relationship Matrix)와 유사한 행렬을 구할 수 있다. 이 행렬에 약 0.90 또는 0.95 또는 0.99의 가중치를 주고 나머지를 NRM에 주어 improved non-singular matrix Gwt를 구할 수 있다.

Data

animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
13 0 0 1 558 9 2 0 1 1 0 0 0 2 1 2
14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0

Pedigree

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

R Program

# Mixed Linear Model(GBLUP Model) for Computing SNP Effect

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

# Raphael Mrode

# Example 11.2 p183

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

# function to make numerator relationship matrix
makenrm = function(pedi) {
  n = nrow(pedi)
  nrm = matrix(c(0), nrow = n, ncol = n)
  
  for (i in 1:n) {
    animal = pedi[i, 1]
    sire = pedi[i, 2]
    dam = pedi[i, 3]
    
    if (sire == 0 & dam == 0) {
      # 부모를 모를 경우
      nrm[animal, animal] = 1
    } else if (sire != 0 & dam == 0) {
      # 부만 알고 있을 경우
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, sire])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1
    } else if (sire == 0 & dam != 0) {
      # 모만 알고 있을 경우
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, dam])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1
    } else {
      # 부와 모를 알고 있을 경우
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
    }
  }
  return(nrm)
}

# set working directory 
setwd(".") 

# print working directory 
getwd()

# read data
data = read.table("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data

# read pedigree : animal sire dam
pedi = read.table("snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
pedi

# variance component and ratio
sigma_a = 35.241
sigma_e = 245
sigma_a
sigma_e

# SNP 
M_all = as.matrix(data[, 7:16])
M_all

#  mean of each SNP(allele frequency)
snp_mean = colMeans(M_all) / 2
snp_mean

#  sum of 2pq(heterozygote ratio)
sum2pq = 0
for (i in 1:10) {
  sum2pq = sum2pq + 2 * snp_mean[i] * (1 - snp_mean[i])
}
sum2pq

# alpha
alpha = sigma_e  / sigma_a
alpha

# observation
y = data[1:8, 6]
y

# design matrix of fixed effect(general mean)
X = matrix(rep(1, 8), 8, 1)
X

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

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

# Numerater Relationship Matrix
A = makenrm(pedi)
A
A22 = A[13:26, 13:26]
A22

# weighted G
wt = 0.99
Gw = wt * G + (1 - wt) * A22
round(Gw, 3)
Gwi = ginv(Gw)
round(Gwi, 3)

# Left Hand Side
LHS = rbind(
  cbind(t(X) %*% X, t(X) %*% W), 
  cbind(t(W) %*% X, t(W) %*% W + Gwi * alpha))
round(LHS, 3)

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

# Solutions

# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)

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

RStuio에서 실행한 화면

프로그램 실행 결과

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

책의 결과와 값이 다르다. 책에는 A(NRM)을 섞었다는 설명이 없지만 G에 0.99, A에 0.01 가중치를 주었을 때 비슷한 값을 나타내었다. 그래도 같은 결과값이 나오지 않은 이유는 책에 나오는 5세대까지 추적하여 얻은 A22에 0.01의 가중치를 주어 구한 값으로 생각된다. 하지만 우리는 혈통이 없으므로 책과 같은 A22를 구할 수 없고 그래서 책과 같은 결과값을 구하지 못한 것으로 생각된다. 책에 나오는 A22를 입력하여 사용하면 같은 값이 나올까?

저작자표시

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

GBLUP Model with residual polygenic effects (0)	2020.12.29
SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0)	2020.12.24
Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis (0)	2020.12.19
Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis (0)	2020.12.19
Fixed Effect Model for SNP Effects, weighted analysis (0)	2020.12.16

blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기

투정이 2020. 12. 19. 12:47

2020. 12. 19. 12:47

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

# Raphael Mrode

# Example 11.2 p183

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

2020/12/19 - [Animal Breeding/R for Genetic Evaluation] - Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis

Mixed Linear Model for Computing SNP Effects(SNP-BLUP Model), unweighted analysis

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.2 p183 SNP effect를 임의 효과로 다루어 SNP effect를 추정한다. 각 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

1 ~ 3 : animal, sire, dam

4 : general mean

5 : EDC(using weight)

6 : Fat DYD

7 : EDC 역수

8 - 17 : SNP1 ~ SNP10의 coding하고 평균을 0으로 scaling한 값

(7 - 17 컬럼은 원래의 자료에서 계산을 하여 입력하여야 한다.)

* 계산 방법은 위 포스팅을 참고

Renumf90 Parameter File

# 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)
7
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
4 cross alpha
RANDOM
diagonal
RANDOM_REGRESSION
data
RR_POSITION
8 9 10 11 12 13 14 15 16 17
(CO)VARIANCES
9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96
OPTION solv_method FSPAK

가중치로 7열을 주었다.

나머지는 설명은 다음 포스팅 참조

2020/12/19 - [Animal Breeding/BLUPF90] - blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기

blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.2 p183 간단한 설명은 다음 포스팅을 참고한다. 2020/12/19 - [Animal Breeding/R f..

bhpark.tistory.com

Renumf90 실행 화면

renumf90 실행 로그

 RENUMF90 version 1.145
 renumf90_snpblup_snp_w.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     
 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
 table with            1  elements sorted
 added count
 Effect group            2  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:diag      

 Wrote parameter file "renf90.par"
 Wrote renumbered data "renf90.dat" 8 records

renumf90 실행 결과로 생성된 파일

renf90.dat

 9 0.00179211 1.3571429 -0.3571429 0.2857143 0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1
 13.4 0.00138504 0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143 1 1
 12.7 0.00333333 0.3571429 0.6428571 1.2857143 0.2857143 0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1
 15.4 0.01369863 -0.6428571 -0.3571429 1.2857143 0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857 1 1
 5.9 0.01923077 -0.6428571 0.6428571 0.2857143 1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1
 7.7 0.01149425 0.3571429 0.6428571 -0.7142857 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 0.2142857 1 1
 10.2 0.01562500 -0.6428571 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 -0.7857143 1 1
 4.8 0.00970874 -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 -0.7857143 1 1

renf90.par

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

가중치 열이 2열로 재배치.

위에서 생성된 renf90.dat와 renf90.par를 이용하여 blupf90 실행

blupf90 실행 화면

blupf90 실행 로그

renf90.par
     BLUPF90 ver. 1.68

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

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

 Residual (co)variance Matrix
  245.00    

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

 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 =           11
 G
  9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  9.9600    
 read            8  records in   0.1718750      s,                      66 
  nonzeroes
 finished peds in   0.1718750      s,                      66  nonzeroes
 left hand side
      0.0003     -0.0001      0.0001      0.0001      0.0001     -0.0001     -0.0000     -0.0000      0.0000      0.0002      0.0001
     -0.0001      0.1005     -0.0000     -0.0001     -0.0001      0.0000      0.0000      0.0000     -0.0000     -0.0001      0.0000
      0.0001     -0.0000      0.1005     -0.0000      0.0000     -0.0000     -0.0001     -0.0000      0.0000      0.0000      0.0001
      0.0001     -0.0001     -0.0000      0.1006      0.0000     -0.0000     -0.0001     -0.0000      0.0000      0.0001      0.0000
      0.0001     -0.0001      0.0000      0.0000      0.1006     -0.0000     -0.0001     -0.0000      0.0000      0.0000      0.0001
     -0.0001      0.0000     -0.0000     -0.0000     -0.0000      0.1004     -0.0000      0.0000     -0.0000     -0.0000      0.0000
     -0.0000      0.0000     -0.0001     -0.0001     -0.0001     -0.0000      0.1006      0.0000     -0.0000      0.0001     -0.0002
     -0.0000      0.0000     -0.0000     -0.0000     -0.0000      0.0000      0.0000      0.1004     -0.0000     -0.0000     -0.0000
      0.0000     -0.0000      0.0000      0.0000      0.0000     -0.0000     -0.0000     -0.0000      0.1004      0.0000      0.0000
      0.0002     -0.0001      0.0000      0.0001      0.0000     -0.0000      0.0001     -0.0000      0.0000      0.1006     -0.0001
      0.0001      0.0000      0.0001      0.0000      0.0001      0.0000     -0.0002     -0.0000      0.0000     -0.0001      0.1006
 right hand side:
    0.00   -0.00    0.00    0.00    0.00   -0.00   -0.00   -0.00    0.00    0.00
    0.00
 solution:
    9.12    0.00   -0.00    0.00   -0.00    0.00    0.00    0.00   -0.00    0.00
   -0.00
 solutions stored in file: "solutions"

blupf90 실행 결과 : solutions

trait/effect level  solution
   1   1         1          9.12441059
   1   2         1          0.00004354
   1   3         1         -0.00440133
   1   4         1          0.00439877
   1   5         1         -0.00104827
   1   6         1          0.00048476
   1   7         1          0.00229456
   1   8         1          0.00000000
   1   9         1         -0.00000000
   1  10         1          0.00179833
   1  11         1         -0.00125139

실행 결과가 책과 다른데 이것은 책에서 EDC의 역수가 EDC를 가중치로 주어 방정식을 푼 것으로 보인다. 이 결과를 이용한 DGV 계산은 생략한다.

저작자표시

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

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(unweighted analysis) 풀기 (0)	2020.12.19
blupf90으로 Fixed Effect Model for SNP Effects(weighted analysis) 풀기 (0)	2020.12.17
blupf90으로 Fixed Effect Model for SNP Effects(unweighted analysis) 풀기 (0)	2020.12.17

blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기

투정이 2020. 12. 19. 12:21

2020. 12. 19. 12:21

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

# Raphael Mrode

# Example 11.2 p183

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

2020/12/19 - [Animal Breeding/R for Genetic Evaluation] - Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis

Mixed Linear Model for Computing SNP Effects(SNP-BLUP Model), unweighted analysis

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.2 p183 SNP effect를 임의 효과로 다루어 SNP effect를 추정한다. 각 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

1 ~ 3 : animal, sire, dam

4 : general mean

5 : EDC(using weight)

6 : Fat DYD

7 : EDC 역수

8 - 17 : SNP1 ~ SNP10의 coding하고 평균을 0으로 scaling한 값

(7 - 17 컬럼은 원래의 자료에서 계산을 하여 입력하여야 한다.)

* 계산 방법은 위 포스팅을 참고

Renumf90 Parameter File

# 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
4 cross alpha
RANDOM
diagonal
RANDOM_REGRESSION
data
RR_POSITION
8 9 10 11 12 13 14 15 16 17
(CO)VARIANCES
9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96
OPTION solv_method FSPAK

첫째

EFFECT
4 cross alpha

general mean의 추정을 위하여 넣은 effect

두번째

EFFECT
4 cross alpha

SNP 효과들을 회귀 변수로 다루기 위한 effect. 공변이를 임의 효과로 다루므로 일종의 random regression model이다. 이 효과에 대해서 nested 하는데 결국 모두 1이므로 전체에 대해서 회귀 계수(SNP effect)를 추정하게 된다.

RANDOM
diagonal

대각 성분에만 분산을 더하여 효과 추정

RANDOM_REGRESSION
data

공변이들이 자료에 있다는 의미

RR_POSITION
8 9 10 11 12 13 14 15 16 17

자료에 있는 공변이들의 위치

(CO)VARIANCES
9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.96

10개 random SNP effect의 분산 행렬. 여기서 9.96은 상가적 유전분산(35.242)을 각 SNP의 2p(1-p)의 합(3.5382)으로 나누어준 값이다.

Renumf90 실행 화면

renumf90 실행 로드

 RENUMF90 version 1.145
 renumf90_mlm_snp_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     
 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
 table with            1  elements sorted
 added count
 Effect group            2  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:diag      

 Wrote parameter file "renf90.par"
 Wrote renumbered data "renf90.dat" 8 records

renumf90으로 생성된 파일

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
 13.4 0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143 0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143 1 1
 12.7 0.3571429 0.6428571 1.2857143 0.2857143 0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1
 15.4 -0.6428571 -0.3571429 1.2857143 0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 0.2142857 1 1
 5.9 -0.6428571 0.6428571 0.2857143 1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571 1.2142857 1 1
 7.7 0.3571429 0.6428571 -0.7142857 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 0.2142857 1 1
 10.2 -0.6428571 -0.3571429 0.2857143 0.2857143 -0.2857143 0.7857143 -0.1428571 0.07142857 0.8571429 -0.7857143 1 1
 4.8 -0.6428571 0.6428571 0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857 0.8571429 -0.7857143 1 1

renf90.par

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

위에서 생성된 renf90.dat와 renf90.par를 이용하여 blupf90 실행

blupf90 실행 화면

blupf90 실행 로그

renf90.par
     BLUPF90 ver. 1.68

 Parameter file:             renf90.par
 Data file:                  renf90.dat
 Number of Traits             1
 Number of Effects           11
 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  covariable             2         1    13
    3  covariable             3         1    13
    4  covariable             4         1    13
    5  covariable             5         1    13
    6  covariable             6         1    13
    7  covariable             7         1    13
    8  covariable             8         1    13
    9  covariable             9         1    13
   10  covariable            10         1    13
   11  covariable            11         1    13

 Residual (co)variance Matrix
  245.00    

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

 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 =           13
 # equations =           11
 G
  9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600      0.0000    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      9.9600    
  0.0000    
  0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000      0.0000    
  9.9600    
 read            8  records in   0.2343750      s,                      66 
  nonzeroes
 finished peds in   0.2343750      s,                      66  nonzeroes
 left hand side
      0.0327     -0.0006      0.0047      0.0093      0.0052     -0.0052     -0.0070     -0.0047      0.0023      0.0117      0.0070
     -0.0006      0.1162     -0.0021     -0.0042     -0.0016      0.0016     -0.0040      0.0001     -0.0000     -0.0063      0.0080
      0.0047     -0.0021      0.1092      0.0013      0.0028      0.0013     -0.0051     -0.0007      0.0003      0.0017      0.0051
      0.0093     -0.0042      0.0013      0.1194      0.0056      0.0026     -0.0142     -0.0013      0.0007      0.0033      0.0102
      0.0052     -0.0016      0.0028      0.0056      0.1130     -0.0003     -0.0093     -0.0007      0.0004     -0.0002      0.0134
     -0.0052      0.0016      0.0013      0.0026     -0.0003      0.1048     -0.0030      0.0007     -0.0004     -0.0039      0.0030
     -0.0070     -0.0040     -0.0051     -0.0142     -0.0093     -0.0030      0.1264      0.0010     -0.0005      0.0057     -0.0219
     -0.0047      0.0001     -0.0007     -0.0013     -0.0007      0.0007      0.0010      0.1011     -0.0003     -0.0017     -0.0010
      0.0023     -0.0000      0.0003      0.0007      0.0004     -0.0004     -0.0005     -0.0003      0.1006      0.0008      0.0005
      0.0117     -0.0063      0.0017      0.0033     -0.0002     -0.0039      0.0057     -0.0017      0.0008      0.1127     -0.0057
      0.0070      0.0080      0.0051      0.0102      0.0134      0.0030     -0.0219     -0.0010      0.0005     -0.0057      0.1264
 right hand side:
    0.32    0.00    0.01    0.12    0.04   -0.04   -0.05   -0.05    0.02    0.11
    0.07
 solution:
    9.94    0.09   -0.31    0.26   -0.08    0.11    0.14   -0.00    0.00   -0.06
   -0.02
 solutions stored in file: "solutions"

blupf90 실행 결과 : solutions

trait/effect level  solution
   1   1         1          9.94394207
   1   2         1          0.08702092
   1   3         1         -0.31079215
   1   4         1          0.26246002
   1   5         1         -0.08027711
   1   6         1          0.11020813
   1   7         1          0.13908022
   1   8         1         -0.00000000
   1   9         1          0.00000000
   1  10         1         -0.06069044
   1  11         1         -0.01580233

general mean effect(1개) 추정하고, 10개의 SNP effect를 추정

저작자표시

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

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으로 Fixed Effect Model for SNP Effects(weighted analysis) 풀기 (0)	2020.12.17
blupf90으로 Fixed Effect Model for SNP Effects(unweighted analysis) 풀기 (0)	2020.12.17
blupf90으로 multiple trait random regression model의 육종가 구하기 (0)	2014.04.01

Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis

투정이 2020. 12. 19. 11:55

2020. 12. 19. 11:55

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

# Raphael Mrode

# Example 11.2 p183

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

2020/12/19 - [Animal Breeding/R for Genetic Evaluation] - Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis

Mixed Linear Model for Computing SNP Effects(SNP-BLUP Model), unweighted analysis

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.2 p183 SNP effect를 임의 효과로 다루어 SNP effect를 추정한다. 각 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

프로그램

# Mixed Linear Model for Computing SNP Effect, weighted analysis

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

# Raphael Mrode

# Example 11.2 p183

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

# set working directory 
setwd(".") 

# print working directory 
getwd()

# read data
data = read.table("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data

# 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])
}
s2pq

# alpha
alpha = sigma_e * s2pq / sigma_a
alpha

# observation
y = data[1:8, 6]
y

# design matrix of fixed effect(general mean)
X = matrix(rep(1, 8), 8, 1)
X

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

# Weight : R from EDC
R = diag(data[1:8, 5])
R
Rinv = ginv(R)
Rinv

# Left Hand Side
LHS = rbind(
  cbind(t(X) %*% Rinv %*% X, t(X) %*% Rinv %*% Z), 
  cbind(t(Z) %*% Rinv %*% X, t(Z) %*% Rinv %*% Z + diag(10) * alpha))
LHS

# Right Hand Side
RHS = rbind(t(X) %*% Rinv %*% y, 
            t(Z) %*% Rinv %*% 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[2:11]
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에서 실행하는 화면

실행 결과

> # Mixed Linear Model for Computing SNP Effect, weighted analysis
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
> 
> # Raphael Mrode
> 
> # Example 11.2 p183
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+   install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # 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/16_mixed linear model SNPBLUP for computing SNP effects_weighted"
> 
> # 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
> 
> # 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])
+ }
> s2pq
    snp1 
3.538265 
> 
> # alpha
> alpha = sigma_e * s2pq / sigma_a
> alpha
    snp1 
24.59848 
> 
> # 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 SNP effect
> M = as.matrix(data[1:8, 7:16])
> Z = sweep(M, 2, snp_mean * 2)
> Z
        snp1       snp2       snp3       snp4       snp5       snp6       snp7       snp8       snp9      snp10
1  1.3571429 -0.3571429  0.2857143  0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
2  0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143  0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143
3  0.3571429  0.6428571  1.2857143  0.2857143  0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
4 -0.6428571 -0.3571429  1.2857143  0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857  0.8571429  0.2142857
5 -0.6428571  0.6428571  0.2857143  1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
6  0.3571429  0.6428571 -0.7142857  0.2857143 -0.2857143  0.7857143 -0.1428571 0.07142857  0.8571429  0.2142857
7 -0.6428571 -0.3571429  0.2857143  0.2857143 -0.2857143  0.7857143 -0.1428571 0.07142857  0.8571429 -0.7857143
8 -0.6428571  0.6428571  0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857  0.8571429 -0.7857143
> 
> # Weight : R from EDC
> R = diag(data[1:8, 5])
> R
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]  558    0    0    0    0    0    0    0
[2,]    0  722    0    0    0    0    0    0
[3,]    0    0  300    0    0    0    0    0
[4,]    0    0    0   73    0    0    0    0
[5,]    0    0    0    0   52    0    0    0
[6,]    0    0    0    0    0   87    0    0
[7,]    0    0    0    0    0    0   64    0
[8,]    0    0    0    0    0    0    0  103
> Rinv = ginv(R)
> Rinv
            [,1]        [,2]        [,3]       [,4]       [,5]       [,6]     [,7]        [,8]
[1,] 0.001792115 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[2,] 0.000000000 0.001385042 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[3,] 0.000000000 0.000000000 0.003333333 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[4,] 0.000000000 0.000000000 0.000000000 0.01369863 0.00000000 0.00000000 0.000000 0.000000000
[5,] 0.000000000 0.000000000 0.000000000 0.00000000 0.01923077 0.00000000 0.000000 0.000000000
[6,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.01149425 0.000000 0.000000000
[7,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.015625 0.000000000
[8,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.009708738
> 
> # Left Hand Side
> LHS = rbind(
+   cbind(t(X) %*% Rinv %*% X, t(X) %*% Rinv %*% Z), 
+   cbind(t(Z) %*% Rinv %*% X, t(Z) %*% Rinv %*% Z + diag(10) * alpha))
> LHS
                            snp1          snp2         snp3          snp4          snp5          snp6          snp7          snp8          snp9         snp10
       0.076267880 -0.0292324941  0.0165285648  0.025943492  0.0299278126 -1.845749e-02 -0.0121950399 -0.0108954114  0.0054477057  0.0396312095  1.398055e-02
snp1  -0.029232494 24.6279258974 -0.0028682259 -0.020567675 -0.0149681047  9.542617e-03  0.0095592165  0.0041760706 -0.0020880353 -0.0168110757  7.872053e-04
snp2   0.016528565 -0.0028682259 24.6207119740 -0.004921569  0.0113384107 -2.579590e-03 -0.0160931411 -0.0023612235  0.0011806118  0.0007965454  1.724097e-02
snp3   0.025943492 -0.0205676752 -0.0049215685 24.636989205  0.0111223220 -3.126712e-03 -0.0205867825 -0.0037062131  0.0018531066  0.0129343415  9.602677e-03
snp4   0.029927813 -0.0149681047  0.0113384107  0.011122322 24.6396792508 -7.598423e-03 -0.0258437486 -0.0042754018  0.0021377009  0.0004520377  3.606263e-02
snp5  -0.018457489  0.0095426174 -0.0025795900 -0.003126712 -0.0075984226  2.460613e+01 -0.0005633219  0.0026367842 -0.0013183921 -0.0117993932  5.317478e-05
snp6  -0.012195040  0.0095592165 -0.0160931411 -0.020586783 -0.0258437486 -5.633219e-04 24.6530639051  0.0017421486 -0.0008710743  0.0180342684 -4.347322e-02
snp7  -0.010895411  0.0041760706 -0.0023612235 -0.003706213 -0.0042754018  2.636784e-03  0.0017421486 24.6000355316 -0.0007782437 -0.0056616014 -1.997222e-03
snp8   0.005447706 -0.0020880353  0.0011806118  0.001853107  0.0021377009 -1.318392e-03 -0.0008710743 -0.0007782437 24.5988681661  0.0028308007  9.986111e-04
snp9   0.039631209 -0.0168110757  0.0007965454  0.012934342  0.0004520377 -1.179939e-02  0.0180342684 -0.0056616014  0.0028308007 24.6361259751 -1.650383e-02
snp10  0.013980555  0.0007872053  0.0172409722  0.009602677  0.0360626336  5.317478e-05 -0.0434732183 -0.0019972221  0.0009986111 -0.0165038270  2.465204e+01
> 
> # Right Hand Side
> RHS = rbind(t(X) %*% Rinv %*% y, 
+             t(Z) %*% Rinv %*% y)
> RHS
             [,1]
       0.69592505
snp1  -0.26572369
snp2   0.04235790
snp3   0.34506266
snp4   0.24713577
snp5  -0.15650240
snp6  -0.05461040
snp7  -0.09941786
snp8   0.04970893
snp9   0.40602373
snp10  0.09651420
> 
> # Solutions
> 
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
> gi_LHS
              [,1]          [,2]          [,3]          [,4]          [,5]          [,6]          [,7]          [,8]          [,9]         [,10]         [,11]
 [1,] 13.147293052  1.556141e-02 -8.807918e-03 -1.380434e-02 -1.592688e-02  9.837634e-03  6.465596e-03  5.807560e-03 -2.903780e-03 -2.113444e-02 -7.423874e-03
 [2,]  0.015561405  4.062281e-02 -5.712222e-06  1.751802e-05  5.776172e-06 -4.076705e-06 -8.056763e-06  8.744768e-18  3.433872e-18  2.674969e-06 -1.014599e-05
 [3,] -0.008807918 -5.712222e-06  4.062216e-02  1.740921e-05 -7.949681e-06 -2.339667e-06  2.211834e-05 -8.721898e-18  5.350707e-18  1.280970e-05 -2.336232e-05
 [4,] -0.013804342  1.751802e-05  1.740921e-05  4.060396e-02 -1.520532e-06 -5.199698e-06  2.706639e-05  4.538403e-18 -2.510182e-18  8.726196e-07 -7.949932e-06
 [5,] -0.015926881  5.776172e-06 -7.949681e-06 -1.520532e-06  4.060441e-02  6.012267e-07  3.457150e-05  1.866691e-17 -2.210756e-19  2.482223e-05 -5.028287e-05
 [6,]  0.009837634 -4.076705e-06 -2.339667e-06 -5.199698e-06  6.012267e-07  4.064766e-02  5.778862e-06  1.549308e-17 -7.883759e-18  3.634628e-06 -5.651215e-06
 [7,]  0.006465596 -8.056763e-06  2.211834e-05  2.706639e-05  3.457150e-05  5.778862e-06  4.056633e-02 -5.221111e-18 -1.382570e-18 -4.006944e-05  6.776774e-05
 [8,]  0.005807560  2.030423e-17 -1.004242e-17  8.459318e-18  1.655272e-17 -1.081492e-17 -8.314475e-18  4.065292e-02  2.593318e-17 -1.023724e-17 -6.425168e-18
 [9,] -0.002903780 -6.247715e-18  3.719745e-18  3.142069e-18 -1.326496e-17 -8.469482e-18  9.578672e-18  9.534626e-18  4.065292e-02 -2.255649e-18  5.814246e-18
[10,] -0.021134439  2.674969e-06  1.280970e-05  8.726196e-07  2.482223e-05  3.634628e-06 -4.006944e-05  5.082198e-20  2.004080e-18  4.062485e-02  3.906646e-05
[11,] -0.007423874 -1.014599e-05 -2.336232e-05 -7.949932e-06 -5.028287e-05 -5.651215e-06  6.776774e-05 -1.295537e-18 -7.271990e-18  3.906646e-05  4.056904e-02
> 
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
        [,1]
 [1,]  9.124
 [2,]  0.000
 [3,] -0.004
 [4,]  0.004
 [5,] -0.001
 [6,]  0.000
 [7,]  0.002
 [8,]  0.000
 [9,]  0.000
[10,]  0.002
[11,] -0.001
> 
> # mean effect of solutions
> mean_eff = sol[1]
> round(mean_eff, 3)
[1] 9.124
> 
> # snp effect of solutions
> snp_eff = sol[2:11]
> round(snp_eff, 3)
 [1]  0.000 -0.004  0.004 -0.001  0.000  0.002  0.000  0.000  0.002 -0.001
> 
> # 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.002
 [2,]  0.002
 [3,] -0.002
 [4,]  0.008
 [5,] -0.008
 [6,] -0.003
 [7,]  0.007
 [8,]  0.001
 [9,] -0.003
[10,] -0.001
[11,]  0.000
[12,] -0.002
[13,]  0.001
[14,]  0.003
>

책의 값과 다른 결과 값이 나온다. EDC의 역수를 가중치로 주어야 하는데 EDC 값을 그대로 주어서 나온 것으로 보인다. EDC의 역수를 가중치로 주자.

저작자표시

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

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 Effects, unweighted analysis (0)	2020.12.19
Fixed Effect Model for SNP Effects, weighted analysis (0)	2020.12.16
Fixed Effect Model for SNP Effects, unweighted analysis (0)	2020.12.16

Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis

투정이 2020. 12. 19. 11:41

2020. 12. 19. 11:41

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

# Raphael Mrode

# Example 11.2 p183

SNP effect를 임의 효과로 다루어 SNP effect를 추정한다. 각 SNP effect의 분산은 1) 상가적 유전 분산을 SNP의 수로 나누어 사용할 수도 있고 2) 상가적 유전 분산을 각 SNP의 2p(1-p)를 합한 값(s2pq)으로 나누어 사용하기도 한다. 두번째 경우가 대립 유전자 빈도의 차이를 고려한다는 점에서 많이 쓰이고 있다. 두번째를 사용할 경우 alpha = sigma_e * s2pq / sigma_a 로 계산하여 사용한다. 여기서는 혈통이 필요없다.

추정한 SNP effect를 이용하여 SNP genotype을 가지고 있는 개체의 DGV를 계산한다.

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

Program

# Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effect, unweighted analysis

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

# Raphael Mrode

# Example 11.2 p183

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

# set working directory 
setwd(".") 

# print working directory 
getwd()

# read data
data = read.table("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data

# 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])
}
s2pq

# alpha
alpha = sigma_e * s2pq / sigma_a
alpha

# observation
y = data[1:8, 6]
y

# design matrix of fixed effect(general mean)
X = matrix(rep(1, 8), 8, 1)
X

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

# Left Hand Side
LHS = rbind(
  cbind(t(X) %*% X, t(X) %*% Z), 
  cbind(t(Z) %*% X, t(Z) %*% Z + diag(10) * alpha))
LHS

# Right Hand Side
RHS = rbind(t(X) %*% 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[2:11]
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에서 실행하는 화면

실행 결과

> # Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effect, unweighted analysis
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
> 
> # Raphael Mrode
> 
> # Example 11.2 p183
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+   install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # 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/15_mixed linear model SNPBLUP for computing SNP effects_unweighted"
> 
> # 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
> 
> # 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])
+ }
> s2pq
    snp1 
3.538265 
> 
> # alpha
> alpha = sigma_e * s2pq / sigma_a
> alpha
    snp1 
24.59848 
> 
> # 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 SNP effect
> M = as.matrix(data[1:8, 7:16])
> Z = sweep(M, 2, snp_mean * 2)
> Z
        snp1       snp2       snp3       snp4       snp5       snp6       snp7       snp8       snp9      snp10
1  1.3571429 -0.3571429  0.2857143  0.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
2  0.3571429 -0.3571429 -0.7142857 -0.7142857 -0.2857143  0.7857143 -0.1428571 0.07142857 -0.1428571 -0.7857143
3  0.3571429  0.6428571  1.2857143  0.2857143  0.7142857 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
4 -0.6428571 -0.3571429  1.2857143  0.2857143 -0.2857143 -0.2142857 -0.1428571 0.07142857  0.8571429  0.2142857
5 -0.6428571  0.6428571  0.2857143  1.2857143 -0.2857143 -1.2142857 -0.1428571 0.07142857 -0.1428571  1.2142857
6  0.3571429  0.6428571 -0.7142857  0.2857143 -0.2857143  0.7857143 -0.1428571 0.07142857  0.8571429  0.2142857
7 -0.6428571 -0.3571429  0.2857143  0.2857143 -0.2857143  0.7857143 -0.1428571 0.07142857  0.8571429 -0.7857143
8 -0.6428571  0.6428571  0.2857143 -0.7142857 -0.2857143 -0.2142857 -0.1428571 0.07142857  0.8571429 -0.7857143
> 
> # Left Hand Side
> LHS = rbind(
+   cbind(t(X) %*% X, t(X) %*% Z), 
+   cbind(t(Z) %*% X, t(Z) %*% Z + diag(10) * alpha))
> LHS
                        snp1        snp2       snp3        snp4        snp5       snp6        snp7        snp8        snp9      snp10
       8.0000000 -0.14285714  1.14285714  2.2857143  1.28571429 -1.28571429 -1.7142857 -1.14285714  0.57142857  2.85714286  1.7142857
snp1  -0.1428571 28.47603006 -0.52040816 -1.0408163 -0.39795918  0.39795918 -0.9693878  0.02040816 -0.01020408 -1.55102041  1.9693878
snp2   1.1428571 -0.52040816 26.76174435  0.3265306  0.68367347  0.31632653 -1.2448980 -0.16326531  0.08163265  0.40816327  1.2448980
snp3   2.2857143 -1.04081633  0.32653061 29.2515403  1.36734694  0.63265306 -3.4897959 -0.32653061  0.16326531  0.81632653  2.4897959
snp4   1.2857143 -0.39795918  0.68367347  1.3673469 27.68011170 -0.08163265 -2.2755102 -0.18367347  0.09183673 -0.04081633  3.2755102
snp5  -1.2857143  0.39795918  0.31632653  0.6326531 -0.08163265 25.68011170 -0.7244898  0.18367347 -0.09183673 -0.95918367  0.7244898
snp6  -1.7142857 -0.96938776 -1.24489796 -3.4897959 -2.27551020 -0.72448980 30.9658260  0.24489796 -0.12244898  1.38775510 -5.3673469
snp7  -1.1428571  0.02040816 -0.16326531 -0.3265306 -0.18367347  0.18367347  0.2448980 24.76174435 -0.08163265 -0.40816327 -0.2448980
snp8   0.5714286 -0.01020408  0.08163265  0.1632653  0.09183673 -0.09183673 -0.1224490 -0.08163265 24.63929537  0.20408163  0.1224490
snp9   2.8571429 -1.55102041  0.40816327  0.8163265 -0.04081633 -0.95918367  1.3877551 -0.40816327  0.20408163 27.61888721 -1.3877551
snp10  1.7142857  1.96938776  1.24489796  2.4897959  3.27551020  0.72448980 -5.3673469 -0.24489796  0.12244898 -1.38775510 30.9658260
> 
> # Right Hand Side
> RHS = rbind(t(X) %*% y, 
+             t(Z) %*% y)
> RHS
        [,1]
       79.10
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)
> gi_LHS
               [,1]          [,2]          [,3]          [,4]          [,5]          [,6]          [,7]          [,8]          [,9]         [,10]
 [1,]  1.388294e-01 -5.663085e-05 -4.970969e-03 -9.042826e-03 -4.643044e-03  6.969389e-03  5.829371e-03  5.807560e-03 -2.903780e-03 -1.427606e-02
 [2,] -5.663085e-05  3.548365e-02  7.785237e-04  1.479041e-03  7.596326e-04 -4.419903e-04  8.832420e-04  2.243621e-17 -1.284610e-17  1.772102e-03
 [3,] -4.970969e-03  7.785237e-04  3.771277e-02  2.541537e-04 -4.906823e-04 -6.784168e-04  1.063631e-03  2.168404e-18 -4.824700e-18 -1.381252e-04
 [4,] -9.042826e-03  1.479041e-03  2.541537e-04  3.555059e-02 -8.504422e-04 -1.229133e-03  3.172508e-03  1.517883e-18 -4.770490e-18 -3.301998e-04
 [5,] -4.643044e-03  7.596326e-04 -4.906823e-04 -8.504422e-04  3.694391e-02  5.655121e-05  1.781453e-03 -1.338990e-17  6.478108e-18  3.572256e-04
 [6,]  6.969389e-03 -4.419903e-04 -6.784168e-04 -1.229133e-03  5.655121e-05  3.941053e-02  9.383621e-04  2.276825e-18 -4.391019e-18  5.733858e-04
 [7,]  5.829371e-03  8.832420e-04  1.063631e-03  3.172508e-03  1.781453e-03  9.383621e-04  3.414696e-02  7.155734e-18 -3.713392e-18 -2.095373e-03
 [8,]  5.807560e-03  8.548257e-18  2.059984e-18 -1.702197e-17 -6.559423e-18 -1.095044e-17 -1.084202e-18  4.065292e-02  2.574980e-18  7.155734e-18
 [9,] -2.903780e-03  2.554651e-18 -2.141299e-18  7.806256e-18  1.208885e-17  1.040834e-17 -6.803369e-18 -9.215718e-19  4.065292e-02  1.192622e-18
[10,] -1.427606e-02  1.772102e-03 -1.381252e-04 -3.301998e-04  3.572256e-04  5.733858e-04 -2.095373e-03 -2.168404e-18  4.878910e-18  3.802145e-02
[11,] -5.999035e-03 -2.241301e-03 -1.064947e-03 -1.808286e-03 -3.287558e-03 -9.713208e-04  4.937702e-03 -2.547875e-18  5.068645e-18  1.999292e-03
              [,11]
 [1,] -5.999035e-03
 [2,] -2.241301e-03
 [3,] -1.064947e-03
 [4,] -1.808286e-03
 [5,] -3.287558e-03
 [6,] -9.713208e-04
 [7,]  4.937702e-03
 [8,] -2.005774e-18
 [9,] -1.095044e-17
[10,]  1.999292e-03
[11,]  3.427246e-02
> 
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
        [,1]
 [1,]  9.944
 [2,]  0.087
 [3,] -0.311
 [4,]  0.262
 [5,] -0.080
 [6,]  0.110
 [7,]  0.139
 [8,]  0.000
 [9,]  0.000
[10,] -0.061
[11,] -0.016
> 
> # mean effect of solutions
> mean_eff = sol[1]
> round(mean_eff, 3)
[1] 9.944
> 
> # snp effect of solutions
> snp_eff = sol[2:11]
> round(snp_eff, 3)
 [1]  0.087 -0.311  0.262 -0.080  0.110  0.139  0.000  0.000 -0.061 -0.016
> 
> # 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.070
 [2,]  0.111
 [3,]  0.045
 [4,]  0.253
 [5,] -0.495
 [6,] -0.357
 [7,]  0.145
 [8,] -0.224
 [9,]  0.027
[10,]  0.114
[11,] -0.240
[12,]  0.143
[13,]  0.054
[14,]  0.354
>

저작자표시

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

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
Fixed Effect Model for SNP Effects, weighted analysis (0)	2020.12.16
Fixed Effect Model for SNP Effects, unweighted analysis (0)	2020.12.16
임의 회귀 모형의 해를 이용한 개체 육종가 그래프 및 305일령 육종가 계산 (0)	2020.11.26

PREV 이전 1 2 3 4 5 6 7 ···50 NEXT 다음

Let's Study Animal Breeding

전체 글

Mixed Linear Model(GBLUP model) for Computing DGV

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

blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(weighted analysis) 풀기

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

blupf90으로 Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects(unweighted analysis) 풀기

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

Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis

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

Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis

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

+ Recent posts

티스토리툴바