단형질 임의회귀모형의 육종가 분할하기

 

단형질 임의회귀모형의 육종가 분할하기

 

조건 : 단형질 임의회귀모형 육종가 추정프로그램의 결과 다음과 같은 파일이 생긴다.

 

sol_animal.dat - 개체 육종가의 해 sol_pe.dat - 개체 영구환경효과의 해 sol_htd.dat - herd-test-day 효과의 해 sol_freg.dat- fixed regression(days-in-milk)효과의 해

 

결과파일로 ptbv_strrm.rst라는 파일이 생긴다.

 

역행렬을 구하기 위하여 inverse.f95라는 파일이 필요

 

sol_animal.dat-5.83171366463021701E-002 5.51949593578108383E-002 -4.41815610161045771E-002
-7.27854215878254718E-002 -3.04895988807174119E-002 -2.43928678312897995E-002
0.13108791818763491 -2.47079050896106614E-002 6.85809683543737447E-002
0.34455565510303948 6.28153836460515393E-003 -0.31640906759476300
-0.45374052814589877 -5.20171296475242101E-002 0.27982278588533233
-0.54855497787185925 7.30057170955604601E-002 0.19457594145119805
0.85180185558070165 -9.50264187408815866E-003 -0.31306198881963443
0.22084487785983828 1.26955535468971717E-002 -1.74370132185816401E-002
sol_pe.dat 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
-0.64862564652799926 -0.36004453207614845 -1.4718223028472115
-0.77616922037832869 0.13701234773550813 0.96882278786722387
-1.9926586979646272 0.98509553381424153 -6.93086269042940128E-002
3.5187699375965571 -1.0509720917831973 -0.40475540593421461
-0.10135823945903098 0.28891043192047283 0.97707271902248471
sol_htd.dat21.851215151168674
19.355854125686076
20.325140095758933
20.008023326152944
18.081103566222183
14.775098104209940
14.873509605460567
14.936790948941656
12.269442512214830
11.765061264213903
sol_freg.dat

-2.02663490753281781E-012
-0.62535641942759423
-0.13458176861200177
0.34790193427196797
-0.42176752253868177

프로그램PROGRAM PTBV_STRRM ! program name
! - partitioning animal solutions from the single trait random regression model
! prgrammer
! - Park Byoungho
! update
! - 2011.5.26 : 1st edition USE gi

IMPLICIT NONE

! data dictionary
INTEGER, PARAMETER :: NO_OF_ANIMAL = 8 ! define new data type for pedigree
TYPE ped_t
INTEGER :: sire, dam
REAL(KIND = 8) :: alpha_par
REAL(KIND = 8) :: alpha_prog
REAL(KIND = 8) :: alpha_anim
REAL(KIND = 8), DIMENSION(3,3) :: qrq ! Q'R-1Q
REAL(KIND = 8), DIMENSION(3,3) :: w1 ! w1 of w1 * pa
REAL(KIND = 8), DIMENSION(3) :: pa ! pa of w1 * pa
REAL(KIND = 8), DIMENSION(3,3) :: w2 ! w2 of w2 * yd
REAL(KIND = 8), DIMENSION(3) :: yield_devi ! yd of w2 * yd
REAL(KIND = 8), DIMENSION(3,3) :: w3 ! w3 of w3 * pc
REAL(KIND = 8), DIMENSION(3) :: n_pc ! nemerator of progeny contribution
REAL(KIND = 8), DIMENSION(3) :: pc ! pc of w3 * pc
END TYPE REAL(KIND = 8), PARAMETER :: animal_vr(3,3) = & ! animal genetic variance and covariance
RESHAPE((/ 3.297, 0.594, -1.381, &
0.594, 0.921, -0.289, &
-1.381, -0.289, 1.005/),(/3,3/))
REAL(KIND = 8), DIMENSION(3,3) :: animal_vr_i ! (inverse of animal and permanent environment variance-covariance) * residual variance REAL(KIND = 8), PARAMETER :: residual_vr = 3.710 ! residual variance

INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file
TYPE(ped_t), DIMENSION(NO_OF_ANIMAL) :: pedi ! pedigree array

INTEGER :: htd ! each herd-test-day
INTEGER :: days ! days in milk
INTEGER :: pre_animal = 0 ! 이전 개체 저장
REAL(KIND = 8) :: yield ! 형질의 값 REAL(KIND = 8), DIMENSION(10) :: mat_obs = 0.0 ! matrix for observation temporarily
REAL(KIND = 8), DIMENSION(10) :: mat_x1b1 = 0.0 ! matrix for herd - test - day X solution
REAL(KIND = 8), DIMENSION(10,5) :: mat_x2 = 0.0 ! design matrix for days-in-milk(fixed regression)
REAL(KIND = 8), DIMENSION(10,3) :: mat_q = 0.0 ! design matrix for animal
REAL(KIND = 8), DIMENSION(10,3) :: mat_z = 0.0 ! design matrix for pe
REAL(KIND = 8), DIMENSION(3,3) :: qrq = 0.0 ! Q'R-1Q
REAL(KIND = 8), DIMENSION(3,3) :: qrq_i = 0.0 ! inverse of Q'R-1Q
REAL(KIND = 8), DIMENSION(3,3) :: diag = 0.0 ! Q'R-1Q +G-1 alpha
REAL(KIND = 8), DIMENSION(3,3) :: diag_i = 0.0 ! inverse of Q'R-1Q +G-1 alpha
!REAL(KIND = 8), DIMENSION(3,3) :: w1 = 0.0 ! w1 of w1 * pa ! Legendre polynomial
REAL(KIND = 8), PARAMETER :: lambda(5,5) = &
RESHAPE((/SQRT(.5) , 0.0 , 0.0 , 0.0 , 0.0, &
0.0 , SQRT(1.5) , 0.0 , 0.0 , 0.0, &
-SQRT(2.5)*.5 , 0.0 , SQRT(2.5)*1.5 , 0.0 , 0.0, &
0.0 , -SQRT(3.5)*1.5, 0.0 , SQRT(3.5)*2.5, 0.0, &
SQRT(4.5)*3./8., 0.0 , -SQRT(4.5)*30./8., 0.0 , SQRT(4.5)*35./8./),(/5,5/))

INTEGER :: dim_min, dim_max, dim_range ! minimun, maximum and range of days in milk
REAL(KIND = 8) :: sut ! standardized unit of time
REAL(KIND = 8), DIMENSION(5) :: sut_poly ! polynomial of standardized unit of time

REAL(KIND = 8), DIMENSION(500,5) :: dim_lp ! Legendre polynomial of days in milk, for x2 design matrix

REAL(KIND = 8), DIMENSION(10) :: sol_htd = 0 ! herd - test - day solution
REAL(KIND = 8), DIMENSION(5) :: sol_dim = 0 ! days-in-milk solution
REAL(KIND = 8), DIMENSION(NO_OF_ANIMAL,3) :: sol_ani = 0 ! animal solution
REAL(KIND = 8), DIMENSION(NO_OF_ANIMAL,3) :: sol_pe = 0 ! permanent environment solution

INTEGER :: i,j ! for loop
INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message CHARACTER(LEN=256), PARAMETER :: pedi_filename = 'strrm.ped' ! pedigree file name
CHARACTER(LEN=256), PARAMETER :: data_filename = 'strrm.dat' ! data file name
CHARACTER(LEN=256), PARAMETER :: htd_filename = 'sol_htd.dat' ! herd-test-day solution file name
CHARACTER(LEN=256), PARAMETER :: freg_filename = 'sol_freg.dat' ! fixed regression solution file name
CHARACTER(LEN=256), PARAMETER :: ani_filename = 'sol_animal.dat' ! animal solution file name
CHARACTER(LEN=256), PARAMETER :: pe_filename = 'sol_pe.dat' ! permanent environment solution file name !pedi 초기화
!pedi%sire = 0
!pedi%dam = 0
!pedi%alpha_par = 0.0
!pedi%alpha_prog = 0.0
!pedi%alpha_anim = 0.0
!pedi%w1 = 0.0
!pedi%yield_devi = 0.0 ! inverse of animal genetic variance and covariance
CALL geninv(animal_vr, animal_vr_i)

! open pedigree file
OPEN(UNIT = 20, FILE = pedi_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF ! 자료 파일 열기
OPEN(UNIT = 30, FILE = data_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF

! herd-test-day solution 파일 열기
OPEN(UNIT = 40, FILE = htd_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF

! fixed regression solution 파일 열기
OPEN(UNIT = 50, FILE = freg_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF

! animal solution 파일 열기
OPEN(UNIT = 60, FILE = ani_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF

! permanent environment 파일 열기
OPEN(UNIT = 70, FILE = pe_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF ! htd 해 읽기
i = 1
DO
READ(UNIT = 40, FMT = *, IOSTAT = status) sol_htd(i)
IF (status /= 0) EXIT ! end of file
i = i + 1
END DO

! WRITE(*,*) sol_htd
CLOSE(40)

! days in milk 해 읽기
i = 1
DO
READ(UNIT = 50, FMT = *, IOSTAT = status) sol_dim(i)
IF (status /= 0) EXIT ! end of file
i = i + 1
END DO

! WRITE(*,*) sol_dim
CLOSE(50) ! animal 해 읽기
i = 1
DO
READ(UNIT = 60, FMT = *, IOSTAT = status) sol_ani(i,:)
IF (status /= 0) EXIT ! end of file
i = i + 1
END DO

! DO i = 1, 8
! WRITE(*,*) sol_ani(i,:)
! END DO
CLOSE(60)

! 영구 환경 효과 해 읽기
i = 1
DO
READ(UNIT = 70, FMT = *, IOSTAT = status) sol_pe(i,:)
IF (status /= 0) EXIT ! end of file
i = i + 1
END DO

! DO i = 1, 8
! WRITE(*,*) sol_pe(i,:)
! END DO
CLOSE(70) !===================================================
!===================================================
! 혈통 파일을 읽어서 계산하기 시작
!===================================================
!=================================================== ! alpha_par - 양쪽 부모를 알면 1, 한 부모만 알면 2/3, 양 부모를 모르면 0.5
! alpha_prog - 개체의 배우자가 있으면 1, 없으면 2/3
! parent average 계산
! nemerator of progeny contribution 계산
DO
READ(UNIT = 20, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file
pedi(animal)%sire = sire
pedi(animal)%dam = dam

! 양 부모 모르면
IF (sire == 0 .AND. dam == 0) THEN
pedi(animal)%alpha_par = 0.5
pedi(animal)%pa = 0
! 부만 알면
ELSE IF (sire /= 0 .AND. dam == 0) THEN
pedi(animal)%alpha_par = 2./3.
pedi(sire)%alpha_prog = pedi(sire)%alpha_prog + 2./3.
pedi(animal)%pa = sol_ani(sire,:) / 2.
pedi(sire)%n_pc = pedi(sire)%n_pc + 2./3. * 2. * sol_ani(animal,:)
! 모만 알면
ELSE IF (sire == 0 .AND. dam /= 0) THEN
pedi(animal)%alpha_par = 2./3.
pedi(dam)%alpha_prog = pedi(dam)%alpha_prog + 2./3.
pedi(animal)%pa = sol_ani(dam,:) / 2.
pedi(dam)%n_pc = pedi(dam)%n_pc + 2./3. * 2. * sol_ani(animal,:)
! 양 부모 알면
ELSE IF (sire /= 0 .AND. dam /= 0) THEN
pedi(animal)%alpha_par = 1
pedi(sire)%alpha_prog = pedi(sire)%alpha_prog + 1.
pedi(dam)%alpha_prog = pedi(dam)%alpha_prog + 1.
pedi(animal)%pa = (sol_ani(sire, :) + sol_ani(dam, :))/ 2.
pedi(sire)%n_pc = pedi(sire)%n_pc + 2. * sol_ani(animal,:) - sol_ani(dam, :)
pedi(dam)%n_pc = pedi(dam)%n_pc + 2. * sol_ani(animal,:) - sol_ani(sire, :)
END IF
END DO
CLOSE(20)

! alpah_par와 alpha_prog를 이용하여 alpha_anim계산, n_pc와 alpha_prog를 이용하여 PC계산
! alpha_anim = 2 * alpha_par + 0.5 * alpha_prog
DO i = 1, NO_OF_ANIMAL
pedi(i)%alpha_anim = 2* pedi(i)%alpha_par + 0.5 * pedi(i)%alpha_prog
IF (pedi(i)%alpha_prog /= 0.0) THEN
! pedi(i)%pc = 0.0
!ELSE
pedi(i)%pc = pedi(i)%n_pc / pedi(i)%alpha_prog
END IF
END DO !===================================================
!===================================================
! 혈통 파일을 읽어서 계산하기 끝
!===================================================
!===================================================

 

!===================================================
! 자료 파일 처리할 때 필요한 르장드르 폴리노미얼 만들기 시작
!=================================================== ! prepare Legendre polynomial of standardized unit of time(days in milk)
! open setup result file
OPEN(UNIT = 25, FILE = 'strrm_setup.rst', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg) IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Setup result file open failed -- error message : ', error_msg
STOP
END IF

READ(UNIT = 25, FMT ='(/ I12 / I12 / I12)', IOSTAT = status) dim_min, dim_max, dim_range
CLOSE(25)

!initialize
dim_lp = 0.0

! prepare Legendre polynomial of standardized unit of time(days in milk)
DO i = dim_min, dim_max
sut = -1. + 2. * (i - dim_min) / dim_range ! standardized unit of time
sut_poly(1) = 1.0 ! polynomial of standardized unit of time
sut_poly(2) = sut
sut_poly(3) = sut_poly(2) * sut
sut_poly(4) = sut_poly(3) * sut
sut_poly(5) = sut_poly(4) * sut
dim_lp(i,:) = MATMUL(sut_poly, lambda)
END DO
!===================================================
! 자료 파일 처리할 때 필요한 르장드르 폴리노미얼 만들기 끝
!=================================================== !===================================================
!===================================================
! 자료 파일을 읽어 개체의 yield deviation 만들기 시작
!===================================================
!===================================================
DO
READ(UNIT = 30, FMT = *, IOSTAT = status) htd, days, animal, yield
IF (status /= 0) EXIT ! end of file

!이전 개체와 같은 계속 읽으면서 design matrix를 만든다
IF (pre_animal == animal .OR. pre_animal == 0 ) THEN
mat_obs(htd) = yield
mat_x1b1(htd) = sol_htd(htd)
mat_x2(htd,:) = dim_lp(days,:)

! 이전 개체와 다른 개체와 나왔으므로 이전 개체의 yield deviation 계산
ELSE
mat_q = mat_x2(:, 1:3)
mat_z = mat_x2(:, 1:3) ! Q'R-1Q 계산
qrq = MATMUL(TRANSPOSE(mat_q),mat_q)/residual_vr
! Q'R-1Q 저장
pedi(pre_animal)%qrq = qrq !diag, diag_i 계산
diag = qrq + animal_vr_i * pedi(pre_animal)%alpha_anim
CALL geninv(diag, diag_i) ! w2 of w2 * yd
pedi(pre_animal)%w2 = MATMUL(diag_i, qrq) ! inverse of Q'R-1Q
CALL geninv(qrq, qrq_i) ! yield deviation 저장
pedi(pre_animal)%yield_devi = MATMUL(MATMUL(qrq_i, TRANSPOSE(mat_q)), &
mat_obs - mat_x1b1 - MATMUL(mat_x2, sol_dim) - MATMUL(mat_z, sol_pe(pre_animal,:))) / residual_vr mat_obs = 0.0
mat_obs(htd) = yield

mat_x1b1 = 0.0
mat_x1b1(htd) = sol_htd(htd)

mat_x2 = 0.0
mat_x2(htd,:) = dim_lp(days,:)

END IF
pre_animal = animal
END DO
CLOSE(30)

!==========================
!== 마지막 개체
!==========================
mat_q = mat_x2(:, 1:3)
mat_z = mat_x2(:, 1:3) !Q'R-1Q 계산
qrq = MATMUL(TRANSPOSE(mat_q),mat_q)/residual_vr
! Q'R-1Q 저장
pedi(pre_animal)%qrq = qrq !diag, diag_i 계산
diag = qrq + animal_vr_i * pedi(pre_animal)%alpha_anim
CALL geninv(diag, diag_i)

! w2 of w2 * yd
pedi(pre_animal)%w2 = MATMUL(diag_i, qrq) ! inverse of Q'R-1Q
CALL geninv(qrq, qrq_i) ! yield deviation을 저장
pedi(pre_animal)%yield_devi = MATMUL(MATMUL(qrq_i, TRANSPOSE(mat_q)), &
mat_obs - mat_x1b1 - MATMUL(mat_x2, sol_dim) - MATMUL(mat_z, sol_pe(pre_animal,:))) / residual_vr

!===================================================
!===================================================
! 자료 파일을 읽어 개체의 yield deviation 만들기 끝
!===================================================
!=================================================== !===================================================
!===================================================
! pedi를 읽으면서 다시 w1과 w3 계산 시작
!===================================================
!===================================================
DO i = 1, NO_OF_ANIMAL !diag, diag_i 계산
diag = pedi(i)%qrq + animal_vr_i * pedi(i)%alpha_anim
CALL geninv(diag, diag_i) ! w1 계산
pedi(i)%w1 = MATMUL(diag_i, 2 * animal_vr_i * pedi(i)%alpha_par)

! w3 계산
pedi(i)%w3 = MATMUL(diag_i, 0.5 * animal_vr_i * pedi(i)%alpha_prog)

END DO
!===================================================
!===================================================
! pedi를 읽으면서 다시 w1과 w3 계산 끝
!===================================================
!===================================================

 

OPEN(UNIT = 80, FILE = 'ptbv_strrm.rst', STATUS = 'REPLACE', ACTION = 'WRITE')
DO i = 1, NO_OF_ANIMAL
!WRITE(80,*) i , pedi(i)%w3
!WRITE(80,*) i, pedi(i)%pc
WRITE(80,*) i, MATMUL(pedi(i)%w1, pedi(i)%pa) + MATMUL(pedi(i)%w2, pedi(i)%yield_devi) + MATMUL(pedi(i)%w3, pedi(i)%pc), &
MATMUL(pedi(i)%w1, pedi(i)%pa) , MATMUL(pedi(i)%w2, pedi(i)%yield_devi) , MATMUL(pedi(i)%w3, pedi(i)%pc)
END DO
CLOSE(80)

END PROGRAM

 

컴파일 및 실행

 

 

실행결과 파일

ptbv_strrm.rst

 

1열 : 개체 번호

3,4,5 :분할한 육종가의 합 - 검증용

6,7,8 : 부모에서 받은 능력

9,10,11:자기 자신의 능력

12,13,14 : 자손에서 받은 능력

(아래를 에디터에 카피하여 볼 것)

 

1 -5.83169591329880596E-002 5.51949976832910860E-002 -4.41816415917402755E-002 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 -5.83169591329880596E-002 5.51949976832910860E-002 -4.41816415917402755E-002
2 -7.27851319067628344E-002 -3.04895592085095771E-002 -2.43929926892826471E-002 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 -7.27851319067628344E-002 -3.04895592085095771E-002 -2.43929926892826471E-002
3 0.13108810533859794 -2.47078706317781262E-002 6.85808823893621911E-002 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.13108810533859794 -2.47078706317781262E-002 6.85808823893621911E-002
4 0.34455516670178793 6.28155459564510264E-003 -0.31640890902907554 -2.84852419927414564E-002 8.57884745632065022E-003 -2.64305963876728181E-002 0.26020689846632500 2.11876409039978633E-002 -0.24430365959618697 0.11283351022820441 -2.34849337646734109E-002 -4.56746530452157434E-002
5 -0.45374108313145534 -5.20171058672236689E-002 0.27982296843841797 1.51343394000211205E-002 -1.49646744960539538E-002 1.47284797031553838E-002 -0.39182897983379933 -6.45193632822226443E-002 0.23759842992643523 -7.70464426976771449E-002 2.74669319110529275E-002 2.74960588088273318E-002
6 -0.54855454789824498 7.30055573642011818E-002 0.19457584078331497 -0.12210634257138972 5.65594236527098594E-003 5.57255574072271290E-002 -0.42644820532685523 6.73496149989301907E-002 0.13885028337608785 0.0000000000000000 0.0000000000000000 0.0000000000000000
7 0.85180209170182319 -9.50267113537227662E-003 -0.31306206354249089 6.85799265470027342E-002 -2.09585824630550387E-002 -4.13130150283222211E-002 0.72921007533089621 1.18797999783033879E-002 -0.28916562221182918 5.40120898239242545E-002 -4.23888650620626091E-004 1.74165736976605208E-002
8 0.22084430320370835 1.26955792689948041E-002 -1.74368298374452946E-002 0.12102323201761608 -1.68087093228478750E-002 -5.04367831183003060E-002 9.98210711860922778E-002 2.95042885918426791E-002 3.29999532808550114E-002 0.0000000000000000 0.0000000000000000 0.0000000000000000

 

 

관련 파일

1306412035_육종가분할_단형질임의회귀모형.zip

다운로드

파일