R을 이용하여 단형질 개체 모형(single trait animal model)에 대한 blupf90 분석 결과 정리하기

blupf90을 이용하여 단형질 개체 모형의 해를 구하는 방법을 아래에서 설명한 바 있다.

 

2014/01/28 - [Animal Breeding/BLUPF90] - blupf90으로 단형질 개체 모형(Single Trait Animal Model)의 육종가(Breeding Value) 구하기

 

여기서는 blupf90이 구한 결과를 R로 정리하는 것을 설명한다.

 

blupf90은 육종가(breeding value)와 육종가의 SE(SEP : Standard Error of Prediction)를 제공하는데 SE를 이용하여 정확도(accuracy)를 구해야 한다. 또한 renumber된 개체 ID이기 때문에 원래 ID를 찾아 주어야 한다. 여기서는 육종가 분할 및 DYD 구하는 것은 생략한다. 육종가 분할 및 DYD는 다음을 참고한다.

 

2020/05/18 - [Animal Breeding/R for Genetic Evaluation] - R을 이용한 단형질 개체 모형(Single Trait Animal Model)의 육종가 구하기 프로그래밍

 

blupf90이 구한 solutions는 다음과 같다.

 

trait/effect level  solution          s.e.
   1   1         1          4.35850233          4.88082357
   1   1         2          3.40443010          5.66554023
   1   2         1         -0.24945855          4.20407502
   1   2         2         -0.00866312          4.13608581
   1   2         3          0.17687209          4.20610397
   1   2         4          0.18261469          4.11029997
   1   2         5         -0.18573210          4.13814812
   1   2         6          0.09844458          4.34094096
   1   2         7         -0.01877010          4.43664612
   1   2         8         -0.04108420          4.27297922

 

renumber된 혈통과 원래 ID는 renadd02.ped에 나와 있다.

 

1 2 5 1 0 2 1 0 0 7
7 0 0 3 0 0 0 0 2 2
2 6 0 2 0 1 1 1 0 4
3 6 7 1 0 2 1 0 1 6
6 0 0 3 0 0 0 2 0 1
4 8 3 1 0 2 1 0 0 8
8 0 0 3 0 0 0 2 0 3
5 8 7 1 0 2 1 0 1 5

 

1, 2, 3열이 리넘버된 혈통이고 10열이 원래 ID이다.

 

단형질 개체 모형(single trait animal model)에 대한 blupf90의 결과로부터 정확도를 계산하고 원래 ID를 붙이는 R 프로그램이다.

 

# blupf90의 결과를 R로 정리하기 - Date :2020-06-03

# blupf90의 육종가와 SE(standard error of prediction)읽어 

#정확도를 계산하고, 원래 ID를 붙임

# set working directory 
setwd("D:/temp/stam_blupf90") 

# print working directory 
getwd()

# sigma_a2
sigma_a = 20

sigma_a

# read solutions
solutions = read.table("solutions", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)

solutions

# read renumbered pedigree
pedi = read.table("renadd02.ped", sep = "", stringsAsFactors = FALSE)

pedi

# read renumbered data
data = read.table("renf90.dat", sep = "", stringsAsFactors = FALSE)

data


# read fixed solutions from solutions
sol_fixed = solutions[solutions$effect == 1, 4]

sol_fixed

# read animal solutions from solutions
sol_animal = solutions[solutions$effect == 2, 3:5]

rownames(sol_animal) = sol_animal[,1]

sol_animal

# reliability
rel = 1 - (sol_animal$se^2 / sigma_a)

rel

# accuracy
acc = sqrt(rel)

acc

# 육종가 분할하기, DYD 구하기 생략

# result
stam_result = data.frame(animal = sol_animal$level, animal_bv = sol_animal$solution, rel = rel, acc = acc, sep = sol_animal$se)

stam_result

# join original animal id
aid_renum = pedi[,c(1,10)]

colnames(aid_renum) = c("animal", "orig_aid")

aid_renum

# 조인하기 위한 패키지 실행 
#install.packages("plyr") 
library(plyr)

stam_result_final = join(x = stam_result, y = aid_renum, by = "animal", type = "left")

stam_result_final = stam_result_final[order(stam_result_final$orig_aid),]

stam_result_final

 

다음은 실행 결과이다.

 

> # blupf90의 결과를 R로 정리하기 - Date :2020-06-03
> 
> # blupf90의 육종가와 SE(standard error of prediction)읽어 
> 
> #정확도를 계산하고, 원래 ID를 붙임
> 
> # set working directory 
> setwd("D:/temp/stam_blupf90") 
> 
> # print working directory 
> getwd()
[1] "D:/temp/stam_blupf90"
> 
> # sigma_a2
> sigma_a = 20
> 
> sigma_a
[1] 20
> 
> # read solutions
> solutions = read.table("solutions", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)
> 
> solutions
   trait effect level    solution       se
1      1      1     1  4.35850233 4.880824
2      1      1     2  3.40443010 5.665540
3      1      2     1 -0.24945855 4.204075
4      1      2     2 -0.00866312 4.136086
5      1      2     3  0.17687209 4.206104
6      1      2     4  0.18261469 4.110300
7      1      2     5 -0.18573210 4.138148
8      1      2     6  0.09844458 4.340941
9      1      2     7 -0.01877010 4.436646
10     1      2     8 -0.04108420 4.272979
> 
> # read renumbered pedigree
> pedi = read.table("renadd02.ped", sep = "", stringsAsFactors = FALSE)
> 
> pedi
  V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1  1  2  5  1  0  2  1  0  0   7
2  7  0  0  3  0  0  0  0  2   2
3  2  6  0  2  0  1  1  1  0   4
4  3  6  7  1  0  2  1  0  1   6
5  6  0  0  3  0  0  0  2  0   1
6  4  8  3  1  0  2  1  0  0   8
7  8  0  0  3  0  0  0  2  0   3
8  5  8  7  1  0  2  1  0  1   5
> 
> # read renumbered data
> data = read.table("renf90.dat", sep = "", stringsAsFactors = FALSE)
> 
> data
   V1 V2 V3
1 4.5  1  2
2 2.9  2  5
3 3.9  2  3
4 3.5  1  1
5 5.0  1  4
> 
> 
> # read fixed solutions from solutions
> sol_fixed = solutions[solutions$effect == 1, 4]
> 
> sol_fixed
[1] 4.358502 3.404430
> 
> # read animal solutions from solutions
> sol_animal = solutions[solutions$effect == 2, 3:5]
> 
> rownames(sol_animal) = sol_animal[,1]
> 
> sol_animal
  level    solution       se
1     1 -0.24945855 4.204075
2     2 -0.00866312 4.136086
3     3  0.17687209 4.206104
4     4  0.18261469 4.110300
5     5 -0.18573210 4.138148
6     6  0.09844458 4.340941
7     7 -0.01877010 4.436646
8     8 -0.04108420 4.272979
> 
> # reliability
> rel = 1 - (sol_animal$se^2 / sigma_a)
> 
> rel
[1] 0.11628766 0.14463971 0.11543447 0.15527171 0.14378651 0.05781158 0.01580856 0.08708243
> 
> # accuracy
> acc = sqrt(rel)
> 
> acc
[1] 0.3410098 0.3803153 0.3397565 0.3940453 0.3791919 0.2404404 0.1257321 0.2950973
> 
> # 육종가 분할하기, DYD 구하기 생략
> 
> # result
> stam_result = data.frame(animal = sol_animal$level, animal_bv = sol_animal$solution, rel = rel, acc = acc, sep = sol_animal$se)
> 
> stam_result
  animal   animal_bv        rel       acc      sep
1      1 -0.24945855 0.11628766 0.3410098 4.204075
2      2 -0.00866312 0.14463971 0.3803153 4.136086
3      3  0.17687209 0.11543447 0.3397565 4.206104
4      4  0.18261469 0.15527171 0.3940453 4.110300
5      5 -0.18573210 0.14378651 0.3791919 4.138148
6      6  0.09844458 0.05781158 0.2404404 4.340941
7      7 -0.01877010 0.01580856 0.1257321 4.436646
8      8 -0.04108420 0.08708243 0.2950973 4.272979
> 
> # join original animal id
> aid_renum = pedi[,c(1,10)]
> 
> colnames(aid_renum) = c("animal", "orig_aid")
> 
> aid_renum
  animal orig_aid
1      1        7
2      7        2
3      2        4
4      3        6
5      6        1
6      4        8
7      8        3
8      5        5
> 
> # 조인하기 위한 패키지 실행 
> #install.packages("plyr") 
> library(plyr)
> 
> stam_result_final = join(x = stam_result, y = aid_renum, by = "animal", type = "left")
> 
> stam_result_final = stam_result_final[order(stam_result_final$orig_aid),]
> 
> stam_result_final
  animal   animal_bv        rel       acc      sep orig_aid
6      6  0.09844458 0.05781158 0.2404404 4.340941        1
7      7 -0.01877010 0.01580856 0.1257321 4.436646        2
8      8 -0.04108420 0.08708243 0.2950973 4.272979        3
2      2 -0.00866312 0.14463971 0.3803153 4.136086        4
5      5 -0.18573210 0.14378651 0.3791919 4.138148        5
3      3  0.17687209 0.11543447 0.3397565 4.206104        6
1      1 -0.24945855 0.11628766 0.3410098 4.204075        7
4      4  0.18261469 0.15527171 0.3940453 4.110300        8
>