# Social Interaction Model
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 8.1 p125
# Introduction to BLUPF90 suite programs
# Standard Edition
# Yutaka Masuda
# University of Georgia
# September 2019
참고
R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기
# Social Interaction Model # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 8.1 p125 하나의 펜(pen) 또는 케이지에 여러 마리를 검정할..
bhpark.tistory.com
자료
7 1 4 1 1 5.5 8 9 1
8 1 4 1 2 9.8 7 9 1
9 2 5 1 2 4.9 7 8 2
10 1 4 2 1 8.23 11 12 1
11 2 5 2 2 7.5 10 12 2
12 3 6 2 2 10 10 11 3
13 2 5 3 1 4.5 14 15 2
14 3 6 3 2 8.4 13 15 3
15 3 6 3 1 6.4 13 14 3
1열 : animal
2열 : sire
3열 : dam
4열 : pen
5열 : sex
6열 : 성장율
7 - 8열 : 같은 펜에 있는 동기축 개체
9열 : fullsib 효과
주어지는 자료는 6열까지이다. fullsib 효과야 2열과 3열을 붙이면 되는 거니까 어렵지 않다. 문제는 동기축을 자료로 붙이는 일이다. 그리고 리넘된 자료로 해야 하므로 renumf90을 적당히 돌리고, 같은 펜에 있는 동기축의 renumber된 번호를 붙여야 한다. 여기서는 일단 수동으로 자료를 준비하여 분석한다.
혈통
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 1 4
8 1 4
9 2 5
10 1 4
11 2 5
12 3 6
13 2 5
14 3 6
15 3 6
blupf90 파라미터 파일
DATAFILE
social2_data.txt
NUMBER_OF_TRAITS
1
NUMBER_OF_EFFECTS
6
OBSERVATION(S)
6
WEIGHT(S)
EFFECTS:
5 2 cross # sex
1 15 cross # direct genetic
7 0 cross # associative genetic 1
8 15 cross # associative genetic 2
9 3 cross # common environmental
4 3 cross # random group (pen)
RANDOM_RESIDUAL VALUES
48.48
RANDOM_GROUP
2 3
RANDOM_TYPE
add_animal
FILE
social_pedi.txt
(CO)VARIANCES
25.7 2.25
2.25 3.60
RANDOM_GROUP
5
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.5
RANDOM_GROUP
6
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.12
OPTION solv_method FSPAK
보통의 유전평가는 renumf90 프로그램으로 renumber를 하고, renumf90이 만들어낸 파라미터 파일을 약간 수정하여 분석하지만 여기서는 renumer되었으므로 blupf90용 파라미터 파일을 이용한다.
파라미터 파일이 약간 tricky하다.
EFFECTS:
5 2 cross # sex
1 15 cross # direct genetic
7 0 cross # associative genetic 1
8 15 cross # associative genetic 2
sex 효과, direct genetic 효과는 일반적이다. 그러나 associative genetic 효과의 design matrix는 일반적인 design matrix와는 다르다. 일반적인 design matrix는 한 행에 1이 하나만 들어가나, associative genetic 효과의 design matrix는 한 행에 같은 pen(group)에 있는 동기축들을 1로 해야 한다. 그래서 7열은 레벨을 만들지 말고, 8열로만 레벨을 만들어 한 효과로 처리하고 동기축을 1로 한다는 의미이다.
실행화면
blupf90 blupf90_social.par | tee blupf90_social_01.log
로그
blupf90_social.par
BLUPF90 ver. 1.68
Parameter file: blupf90_social.par
Data file: social2_data.txt
Number of Traits 1
Number of Effects 6
Position of Observations 6
Position of Weight (1) 0
Value of Missing Trait/Observation 0
EFFECTS
# type position (2) levels [positions for nested]
1 cross-classified 5 2
2 cross-classified 1 15
3 cross-classified 7 0
4 cross-classified 8 15
5 cross-classified 9 3
6 cross-classified 4 3
Residual (co)variance Matrix
48.480
correlated random effects 2 3
Type of Random Effect: additive animal
Pedigree File: social_pedi.txt
trait effect (CO)VARIANCES
1 2 25.70 2.250
1 3 2.250 3.600
Random Effect(s) 5
Type of Random Effect: diagonal
trait effect (CO)VARIANCES
1 5 12.50
Random Effect(s) 6
Type of Random Effect: diagonal
trait effect (CO)VARIANCES
1 6 12.12
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 = 9
# equations = 38
G
25.700 2.2500
2.2500 3.6000
G
12.500
G
12.120
read 9 records in 0.1406250 s, 139
nonzeroes
read 15 additive pedigrees
finished peds in 0.1406250 s, 250 nonzeroes
solutions stored in file: "solutions"
solutions
trait/effect level solution
1 1 1 6.00414149
1 1 2 8.24268672
1 2 1 0.29558614
1 2 2 -0.48338296
1 2 3 0.18779681
1 2 4 0.29558614
1 2 5 -0.48338296
1 2 6 0.18779681
1 2 7 0.12478232
1 2 8 0.52180202
1 2 9 -0.87358211
1 2 10 0.53576023
1 2 11 -0.48764109
1 2 12 0.39919156
1 2 13 -0.57230863
1 2 14 0.15302685
1 2 15 0.19896884
1 4 1 -0.04446498
1 4 2 0.02782802
1 4 3 0.01663697
1 4 4 -0.04446498
1 4 5 0.02782802
1 4 6 0.01663697
1 4 7 -0.07597626
1 4 8 -0.09883240
1 4 9 0.00894718
1 4 10 -0.00305127
1 4 11 0.08331352
1 4 12 0.05970749
1 4 13 0.01905137
1 4 14 0.00474261
1 4 15 0.00209776
1 5 1 0.33277773
1 5 2 -0.51533364
1 5 3 0.18255591
1 6 1 -0.26871653
1 6 2 0.35903349
1 6 3 -0.09031696
그러나 리넘버된 개체의 동기축(same pen or same group)을 옆에 추가한다는 것이 쉽지는 않은 일이다. 그래서 다음과 같은 자료에서 renumf90 과 R을 이용해서 자료를 준비해 본다.
자료
a7 a1 a4 1 1 5.5
a8 a1 a4 1 2 9.8
a9 a2 a5 1 2 4.9
a10 a1 a4 2 1 8.23
a11 a2 a5 2 2 7.5
a12 a3 a6 2 2 10
a13 a2 a5 3 1 4.5
a14 a3 a6 3 2 8.4
a15 a3 a6 3 1 6.4
animal, sire, dam, sex, growth
혈통
a1 0 0
a2 0 0
a3 0 0
a4 0 0
a5 0 0
a6 0 0
a7 a1 a4
a8 a1 a4
a9 a2 a5
a10 a1 a4
a11 a2 a5
a12 a3 a6
a13 a2 a5
a14 a3 a6
a15 a3 a6
renumf90 파라미터 파일
여기서는 associative effect 는 없고, sex, animal, pen, common environment effect(full-sib) 효과만 있는 것으로 한다.
# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
COMBINE 7 2 3
DATAFILE
social1_data.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
RESIDUAL_VARIANCE
48.48
EFFECT
5 cross numer
EFFECT
1 cross alpha
RANDOM
animal
FILE
social1_pedi.txt
FILE_POS
1 2 3
PED_DEPTH
0
(CO)VARIANCES
25.7
EFFECT
4 cross numer
RANDOM
diagonal
(CO)VARIANCES
12.12
EFFECT
7 cross alpha
RANDOM
diagonal
(CO)VARIANCES
12.5
OPTION sol se
공통 환경 효과는 2열과 3열을 combine해서 만들었다.
실행
renumf90 renumf90_social1.par
실행결과
renadd02.ped
1 12 15 1 0 2 1 0 0 a14
2 10 13 1 0 2 1 0 0 a8
12 0 0 3 0 0 0 3 0 a3
3 11 14 1 0 2 1 0 0 a9
13 0 0 3 0 0 0 0 3 a4
14 0 0 3 0 0 0 0 3 a5
4 12 15 1 0 2 1 0 0 a15
15 0 0 3 0 0 0 0 3 a6
5 10 13 1 0 2 1 0 0 a10
10 0 0 3 0 0 0 3 0 a1
6 11 14 1 0 2 1 0 0 a11
7 12 15 1 0 2 1 0 0 a12
8 10 13 1 0 2 1 0 0 a7
11 0 0 3 0 0 0 3 0 a2
9 11 14 1 0 2 1 0 0 a13
renf90.dat
5.5 1 8 1 1
9.8 2 2 1 1
4.9 2 3 1 2
8.23 1 5 2 1
7.5 2 6 2 2
10 2 7 2 3
4.5 1 9 3 2
8.4 2 1 3 3
6.4 1 4 3 3
renf90.tables
Effect group 1 of column 1 with 2 levels, effect # 1
Value # consecutive number
1 4 1
2 5 2
Effect group 3 of column 1 with 3 levels, effect # 3
Value # consecutive number
1 3 1
2 3 2
3 3 3
Effect group 4 of column 1 with 3 levels, effect # 4
Value # consecutive number
a1a4 3 1
a2a5 3 2
a3a6 3 3
renf90.par
# BLUPF90 parameter file created by RENUMF90
DATAFILE
renf90.dat
NUMBER_OF_TRAITS
1
NUMBER_OF_EFFECTS
4
OBSERVATION(S)
1
WEIGHT(S)
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
2 2 cross
3 15 cross
4 3 cross
5 3 cross
RANDOM_RESIDUAL VALUES
48.480
RANDOM_GROUP
2
RANDOM_TYPE
add_animal
FILE
renadd02.ped
(CO)VARIANCES
25.700
RANDOM_GROUP
3
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.120
RANDOM_GROUP
4
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.500
OPTION sol se
위의 renf90.dat와 renf90.par는 바로 사용할 수가 없다. renf90.dat를 읽어서 same pen( or same group)의 동기축을 추가하는 프로그램이 필요하다.
# Analyse the social interaction model using blupf90
# Add competitors' renum ID in the group after running renumf90
# set working directory
setwd(".")
# print working directory
getwd()
# read renf90.dat after running renumf90
data_01 = read.table("renf90.dat", header = FALSE, sep = "", col.names = c("grow", "sex", "raid", "pen", "ce"))
data_01
# raid by pen
data_02 = data_01[order(data_01$pen),]
data_02
rabp_r = length(unique(data_02$pen))
rabp_r
rabp_c = max(table(data_02$pen)) + 1
rabp_c
rabp = matrix(0, rabp_r, rabp_c)
rabp
j = 1
k = 1
rabp[j, k] = rabp[j, k] = data_02$pen[1]
k = k + 1
rabp[j, k] = rabp[j, k] = data_02$raid[1]
for (i in 2 : nrow(data_02)) {
# i = 2
if(data_02$pen[i - 1] == data_02$pen[i]) {
k = k + 1
rabp[j, k] = data_02$raid[i]
}
else {
j = j + 1
k = 1
rabp[j, k] = rabp[j, k] = data_02$pen[i]
k = k + 1
rabp[j, k] = rabp[j, k] = data_02$raid[i]
}
}
rabp
# data_02에 competitor 열 추가
data_03 = data_02
for (i in 1 : (rabp_c - 2)) {
print(i)
data_03 = cbind(data_03, 0)
}
data_03
# add competitors
for (i in 1 : nrow(data_03)) {
# pen number of data
pn = data_03$pen[i]
# row number of pen number in rabp
rn = which(rabp[,1] == pn)
# number of columns of data_01
nc = ncol(data_01)
for (j in 2 : ncol(rabp)) {
# if raid is different from competitor
if (data_03$raid[i] != rabp[rn, j]) {
nc = nc + 1
data_03[i, nc] = rabp[rn, j]
}
}
}
data_03
write.table(data_03, "renf90_02.dat", sep=" ", row.names = FALSE, col.names = FALSE, quote = FALSE)
결과로 생긴 renf90_02.dat는 다음과 같다.
5.5 1 8 1 1 2 3
9.8 2 2 1 1 8 3
4.9 2 3 1 2 8 2
8.23 1 5 2 1 6 7
7.5 2 6 2 2 5 7
10 2 7 2 3 5 6
4.5 1 9 3 2 1 4
8.4 2 1 3 3 9 4
6.4 1 4 3 3 9 1
blupf90을 실행하기 위하여 renf90.par를 수정한다.
# BLUPF90 parameter file created by RENUMF90
DATAFILE
renf90_02.dat
NUMBER_OF_TRAITS
1
NUMBER_OF_EFFECTS
6
OBSERVATION(S)
1
WEIGHT(S)
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
2 2 cross
3 15 cross
6 0 cross # associative genetic 1
7 15 cross # associative genetic 2
4 3 cross
5 3 cross
RANDOM_RESIDUAL VALUES
48.480
RANDOM_GROUP
2 3
RANDOM_TYPE
add_animal
FILE
renadd02.ped
(CO)VARIANCES
25.7 2.25
2.25 3.60
RANDOM_GROUP
5
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.120
RANDOM_GROUP
6
RANDOM_TYPE
diagonal
FILE
(CO)VARIANCES
12.500
OPTION sol se
실행 결과 solutions
trait/effect level solution s.e.
1 1 1 6.00414149 5.87722318
1 1 2 8.24268672 5.56132372
1 2 1 0.15302685 4.57296859
1 2 2 0.52180202 4.52043322
1 2 3 -0.87358211 4.57571923
1 2 4 0.19896884 4.53216874
1 2 5 0.53576023 4.60038547
1 2 6 -0.48764109 4.57506398
1 2 7 0.39919156 4.55704542
1 2 8 0.12478232 4.60557709
1 2 9 -0.57230863 4.56574172
1 2 10 0.29558614 4.86839186
1 2 11 -0.48338296 4.86816815
1 2 12 0.18779681 4.85815710
1 2 13 0.29558614 4.86839186
1 2 14 -0.48338296 4.86816815
1 2 15 0.18779681 4.85815710
1 4 1 0.00474261 1.84114013
1 4 2 -0.09883240 1.83170549
1 4 3 0.00894718 1.88441482
1 4 4 0.00209776 1.83375166
1 4 5 -0.00305127 1.86104141
1 4 6 0.08331352 1.88132183
1 4 7 0.05970749 1.85378900
1 4 8 -0.07597626 1.83835336
1 4 9 0.01905137 1.88519935
1 4 10 -0.04446498 1.87204137
1 4 11 0.02782802 1.89561822
1 4 12 0.01663697 1.87131689
1 4 13 -0.04446498 1.87204137
1 4 14 0.02782802 1.89561822
1 4 15 0.01663697 1.87131689
1 5 1 -0.26871653 3.16958408
1 5 2 0.35903349 3.07874277
1 5 3 -0.09031696 3.19664443
1 6 1 0.33277773 3.25808128
1 6 2 -0.51533364 3.19778089
1 6 3 0.18255591 3.23629084
renadd02.ped를 보고 original ID를 찾아볼 수 있다. 순서가 바뀌었을 뿐 결과가 이전과 같다.
'Animal Breeding > BLUPF90' 카테고리의 다른 글
| inbupgf90을 이용한 근교계수 계산과 교배계획 세우기 (0) | 2023.05.25 |
|---|---|
| blupf90으로 single-step GBLUP Model 풀기 (0) | 2021.01.05 |
| blupf90으로 GBLUP Model with residual polygenic effects 풀기 (0) | 2020.12.29 |
| blupf90으로 gblup model 풀기 (0) | 2020.12.29 |
| blupf90으로 SNP-BLUP Model(with polygenic effects unweighted analysis) 풀기 (0) | 2020.12.24 |