IOD를 이용한 육종가 구하기 프로그램을 아래 이미 작성한 바 있다. 지금 소개하는 방법은 1986년에 Schaeffer와 Kennedy가 소개한 방법으로 여러 나라의 국가 단위 유전평가 프로그램 작성시 사용하는 방법이다.
이 방법에 대한 자세한 설명을 위 논문이나 R. A. Mrode의 Linear Models for the prediction of animal breeding values 2nd Edition의 13장을 참고하기 바란다
본 프로그램을 실행하려면 UNIX sort 명령어가 실행될 수 있어야 한다. 리눅스 계열이라면 문제없겠지만 윈도우즈 명령창에서 실행하려면 원도우즈 명령창에서 UNIX sort명령어가 실행될 수 있도록 미리 프로그램을 설치하여야 한다.
참조(http://bhpark.tistory.com/58#)
혈통 파일 준비
1 0 0
2 0 0
3 0 0
4 1 0
5 3 2
6 1 2
7 4 5
8 3 6
위 파일을 stam_pedi.dat로 저장
혈통 파일은 renumbering되어 있어야 하며, 부모를 모를 경우 0으로 코딩한다.
첫째 : 개체
둘째 : 아비
셋째 : 어미
검정 파일 준비
4 1 4.5
5 2 2.9
6 2 3.9
7 1 3.5
8 1 5.0
위 검정 파일을 stam_perf.dat로 저장
첫째 열 : 개체
둘째 열 : 고정효과
셋째 열 : 관측치
고정효과가 1개인 모형만 고려하며 고정효과가 더 있을 경우 소스를 수정하길 바란다.
프로그램 소스
!================================================
!Single Trait Animal Model의 Solution 찾기
!
! Programmer : Park Byoungho
! Date : 2012.8.22.
!================================================
PROGRAM stam_sol
IMPLICIT NONE
TYPE perf
INTEGER :: ani, f1
REAL :: trt
END TYPE
REAL, PARAMETER :: epsilon = 1.e-8 ! convergence
INTEGER, PARAMETER :: round_max = 1000 ! maximum iteration round
REAL :: alpha ! variance ratio
INTEGER :: round ! iteration round
INTEGER :: i, j, k, l ! loop
INTEGER :: ind, sire, dam ! when reading pedigree
INTEGER :: no_of_pedi ! number of pedigree file
INTEGER :: no_of_perf ! number of performance file
INTEGER :: no_of_f1 ! number of equations for fixed 1
INTEGER :: no_of_ani ! number of equations for animal
REAL, ALLOCATABLE :: fixed1(:), animal(:) ! for solutions of fixed 1 and animal
REAL, ALLOCATABLE :: fixed1_old(:), animal_old(:) ! for solutions of fixed 1 and animal
REAL :: diag, arhs ! diagonal, adjusted right hand side
INTEGER :: f1_lvl_cur, f1_lvl_pre ! previous and current level of fixed 1
INTEGER, ALLOCATABLE :: pedi(:,:) ! for pedigree data
TYPE(perf), ALLOCATABLE :: perf_ani(:) ! for performance file sorted by animal
TYPE(perf), ALLOCATABLE :: perf_f1(:) ! for performance file sorted by fixed 1
REAL :: criteria
INTEGER :: status ! error number
CHARACTER(LEN = 256) :: pedi_filename ! pedigree data file name
CHARACTER(LEN = 256) :: perf_filename ! performance data file name
CHARACTER(LEN = 256) :: variance_ratio ! variance ratio
CHARACTER(LEN = 40) :: error_msg ! error message
WRITE(*,*) "=============================================================="
WRITE(*,*) " Single Trait Animal Model Solver"
WRITE(*,*) " Model y = fixed + animal + e"
WRITE(*,*) ""
WRITE(*,*) " Usage : stam_sol pedigree_file performance_file variance_ratio"
WRITE(*,*) "==============================================================="
IF (COMMAND_ARGUMENT_COUNT() /= 3) THEN
WRITE(*,*) "3 arguments are needed"
STOP
END IF
! 아규먼트에서 파일 네임 얻어 오기
CALL GET_COMMAND_ARGUMENT(1, pedi_filename) ! get pedigree data file name from the command line
CALL GET_COMMAND_ARGUMENT(2, perf_filename) ! get performance data file name from the command line
CALL GET_COMMAND_ARGUMENT(3, variance_ratio) ! get variance ratio from the command line
! 문자열을 숫자로 변환
READ(variance_ratio,*) alpha
WRITE(*,*) "Pedigree file name : ", TRIM(pedi_filename)
WRITE(*,*) "Performance file name : ", TRIM(perf_filename)
WRITE(*,*) "Variance ratio : ", alpha
! 혈통 파일 열기
OPEN(UNIT = 20, FILE = pedi_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF
! pedifgree Output 파일 열기
OPEN(UNIT = 30, FILE = 'stam_01.txt', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'pedigree output file open failed -- error message : ', error_msg
STOP
END IF
! 혈통 파일을 읽어 분석에 필요한 혈통 파일 만들기
no_of_pedi = 0
DO
READ(UNIT = 20, FMT = *, IOSTAT = status) ind, sire, dam
IF (status /= 0) EXIT ! reach end of file
WRITE(30,'(4i8)') ind, 1, sire, dam
no_of_pedi = no_of_pedi + 1
IF (sire /= 0) THEN
WRITE(30,'(4i8)') sire, 2, ind, dam
no_of_pedi = no_of_pedi + 1
END IF
IF (dam /= 0) THEN
WRITE(30,'(4i8)') dam, 2, ind, sire
no_of_pedi = no_of_pedi + 1
END IF
END DO
CLOSE(20)
CLOSE(30)
!WRITE(*,*) no_of_pedi
! 외부 명령어(sort)를 사용하여 혈통 파일 정렬하기
CALL SYSTEM('sort -n stam_01.txt -o stam_02.txt')
! 배열 할당
ALLOCATE(pedi(no_of_pedi,4))
pedi = 0
! 혈통 파일 열기
OPEN(UNIT = 40, FILE = 'stam_02.txt', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'stam_02.txt file open failed -- error message : ', error_msg
STOP
END IF
! 혈통 자료를 다시 읽어, 배열에 저장
DO i = 1, no_of_pedi
READ(UNIT = 40, FMT = *, IOSTAT = status) pedi(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(40)
!DO i = 1, no_of_pedi
! WRITE(*,*) pedi(i,:)
!END DO
! 검정 자료 정렬 하기
CALL SYSTEM('sort -n ' // perf_filename // '-o stam_03.txt')
CALL SYSTEM('sort -nk2 ' // perf_filename // '-o stam_04.txt')
! performance data file 열기
OPEN(UNIT = 50, FILE = 'stam_03.txt', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'performance data file (sorted by animal) open failed -- error message : ', error_msg
STOP
END IF
! data file 열기
OPEN(UNIT = 60, FILE = 'stam_04.txt', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'data file (sorted by fixed 1) open failed -- error message : ', error_msg
STOP
END IF
! Solution Output 파일 열기
OPEN(UNIT = 70, FILE = 'stam_05.txt', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'solution output file open failed -- error message : ', error_msg
STOP
END IF
! Solution Output 파일 열기
OPEN(UNIT = 80, FILE = 'stam_06.txt', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
! file open failed
IF (status /= 0) THEN
WRITE (*,'(1X, A, A)') 'running information file open failed -- error message : ', error_msg
STOP
END IF
WRITE(80,*) "Pedigree file name : ", TRIM(pedi_filename)
WRITE(80,*) "Performance file name : ", TRIM(perf_filename)
WRITE(80,*) "Variance ratio : ", alpha
! 능력검정 파일 읽어 저장하기
no_of_perf = 0
DO
READ(UNIT = 50, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_perf = no_of_perf + 1
END DO
! 배열 할당
ALLOCATE(perf_ani(no_of_perf))
! 자료 파일을 다시 읽어, 배열에 저장(sorted by animal)
REWIND(50)
DO i = 1, no_of_perf
READ(UNIT = 50, FMT = *, IOSTAT = status) perf_ani(i)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(50)
!DO i = 1, no_of_perf
! WRITE(*,*) perf_ani(i)
!END DO
! 배열 할당
ALLOCATE(perf_f1(no_of_perf))
! 자료 파일을 다시 읽어, 배열에 저장(sorted by fixed 1)
DO i = 1, no_of_perf
READ(UNIT = 60, FMT = *, IOSTAT = status) perf_f1(i)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(60)
!DO i = 1, no_of_perf
! WRITE(*,*) perf_f1(i)
!END DO
! number of equations for fixed 1
no_of_f1 = 1
DO i = 2, no_of_perf
IF (perf_f1(i-1)%f1 /= perf_f1(i)%f1) THEN
no_of_f1 = no_of_f1 + 1
END IF
END DO
WRITE(*,*) "level of fixed effect : ", no_of_f1
WRITE(80,*) "level of fixed effect : ", no_of_f1
! 배열 할당
ALLOCATE(fixed1(no_of_f1))
fixed1 = 0.0
ALLOCATE(fixed1_old(no_of_f1))
fixed1_old = 0.0
! number of equations for animal
no_of_ani = 1
DO i = 2, no_of_pedi
IF (pedi(i-1,1) /= pedi(i,1)) THEN
no_of_ani = no_of_ani + 1
END IF
END DO
WRITE(*,*) "number of animals : ", no_of_ani
WRITE(80,*) "number of animals : ", no_of_ani
! 배열 할당
ALLOCATE(animal(no_of_ani))
animal = 0.0
ALLOCATE(animal_old(no_of_ani))
animal_old = 0.0
WRITE(*,*) ' Round ', 'convergence criteria'
WRITE(80,*) ' Round ', 'convergence criteria'
! fixed 1 초기값 계산
j = 1
DO i = 1, no_of_f1
diag = 0.0
arhs = 0.0
DO WHILE ( i == perf_f1(j)%f1 )
diag = diag + 1
arhs = arhs + perf_f1(j)%trt
j = j + 1
END DO
fixed1(i) = arhs/diag
END DO
fixed1_old = fixed1
!WRITE(*,*) fixed1
! animal 초기값 계산
DO i = 1, no_of_perf
animal(perf_ani(i)%ani) = perf_ani(i)%trt - fixed1(perf_ani(i)%f1)
END DO
animal_old = animal_old
!WRITE(*,*) animal
! Iteration on data에 의한 solution 구하기
DO round = 1, round_max
! fixed effect
j = 1
DO i = 1, no_of_f1
diag = 0.0
arhs = 0.0
DO WHILE ( i == perf_f1(j)%f1 )
diag = diag + 1.
arhs = arhs + perf_f1(j)%trt - animal(perf_f1(j)%ani)
j = j + 1
END DO
fixed1(i) = arhs/diag
END DO
!WRITE(*,*) fixed1
! animal effect
j = 1
k = 1
! 개체의 해를 하나씩 구함
DO i = 1, no_of_ani
diag = 0.
arhs = 0.
! 개체순으로 정렬된 파일 처리
DO WHILE ( i == perf_ani(j)%ani)
diag = diag + 1.
arhs = arhs + perf_ani(j)%trt - fixed1(perf_ani(j)%f1)
j = j + 1
END DO
! 혈통 파일 처리
DO WHILE ( i == pedi(k,1) )
! 코드가 1일 때 - 부 모
IF ( pedi(k,2) == 1 ) THEN
! 양쪽 부모를 모를 때
IF ( pedi(k,3) == 0 .AND. pedi(k,4) == 0 ) THEN
diag = diag + alpha
! 아비만 알 때
ELSE IF ( pedi(k,3) /= 0 .AND. pedi(k,4) == 0 ) THEN
diag = diag + 4./3. * alpha
arhs = arhs + 2./3. * alpha * animal(pedi(k,3))
! 어미만 알 때
ELSE IF ( pedi(k,3) == 0 .AND. pedi(k,4) /= 0 ) THEN
diag = diag + 4./3. * alpha
arhs = arhs + 2./3. * alpha * animal(pedi(k,4))
! 양쪽 모두 알 때
ELSE IF ( pedi(k,3) /= 0 .AND. pedi(k,4) /= 0 ) THEN
diag = diag + 2. * alpha
arhs = arhs + alpha * (animal(pedi(k,3)) + animal(pedi(k,4)))
END IF
! 코드가 2일 때 - 배우자 자손
ELSE IF ( pedi(k,2) == 2 ) THEN
! 배우자를 모를 때
IF (pedi(k,4) == 0) THEN
diag = diag + 1./3. * alpha
arhs = arhs + 2./3. * alpha * animal(pedi(k,3))
! 배우자를 알 때
ELSE IF (pedi(k,4) /= 0) THEN
diag = diag + 1./2. * alpha
arhs = arhs + alpha * (animal(pedi(k,3)) - animal(pedi(k,4))/2.)
END IF
END IF
k = k + 1
END DO
animal(i) = arhs / diag
END DO
!WRITE(*,*) animal
criteria = (SUM((fixed1 - fixed1_old)**2) + SUM((animal - animal_old)**2)) / &
&(SUM(fixed1**2) + SUM(animal**2))
WRITE(*,*) round, criteria
WRITE(80,*) round, criteria
fixed1_old = fixed1
animal_old = animal
! convergen check
IF ( criteria < epsilon) THEN
EXIT
END IF
END DO
! 해 출력
WRITE(70,*) "Single Trait Animal Solver"
! writing fixed1 effect
DO i = 1, no_of_f1
WRITE(70,'(a8,i8,f12.4)') 'f1', i, fixed1(i)
END DO
! writing animal effect
DO i = 1, no_of_ani
WRITE(70,'(a8,i8,f12.4)') 'animal', i, animal(i)
END DO
CLOSE(70)
! result information
WRITE(*,*) "stam_01.txt : pedigree file (unsorted)"
WRITE(*,*) "stam_02.txt : pedigree file (sorted)"
WRITE(*,*) "stam_03.txt : performance file (sorted by animal)"
WRITE(*,*) "stam_04.txt : performance file (sorted by fixed 1 effect)"
WRITE(*,*) "stam_05.txt : solution file"
WRITE(*,*) "stam_06.txt : running information file"
WRITE(80,*) "stam_01.txt : pedigree file (unsorted)"
WRITE(80,*) "stam_02.txt : pedigree file (sorted)"
WRITE(80,*) "stam_03.txt : performance file (sorted by animal)"
WRITE(80,*) "stam_04.txt : performance file (sorted by fixed 1 effect)"
WRITE(80,*) "stam_05.txt : solution file"
WRITE(80,*) "stam_06.txt : running information file"
CLOSE(80)
END PROGRAM stam_sol
위 파일을 stam_sol.f95로 저장
프로그램 컴파일 및 실행
우분투에서 컴파일일하고 실행한 화면
결과 파일
stam_01.txt - 분석에 사용하기 위하면 만든 혈통 파일
1 1 0 0
2 1 0 0
3 1 0 0
4 1 1 0
1 2 4 0
5 1 3 2
3 2 5 2
2 2 5 3
6 1 1 2
1 2 6 2
2 2 6 1
7 1 4 5
4 2 7 5
5 2 7 4
8 1 3 6
3 2 8 6
6 2 8 3
stam_02.txt - stam_01을 정렬한 혈통 파일
1 1 0 0
1 2 4 0
1 2 6 2
2 1 0 0
2 2 5 3
2 2 6 1
3 1 0 0
3 2 5 2
3 2 8 6
4 1 1 0
4 2 7 5
5 1 3 2
5 2 7 4
6 1 1 2
6 2 8 3
7 1 4 5
8 1 3 6
stam_03.txt - 개체 순으로 정렬된 능력 검정 파일
4 1 4.5
5 2 2.9
6 2 3.9
7 1 3.5
8 1 5.0
stam_04.txt - 고정효과 순으로 정렬된 능력 검정 파일
4 1 4.5
7 1 3.5
8 1 5.0
5 2 2.9
6 2 3.9
stam_05.txt - 고정효과와 개체의 해
Single Trait Animal Solver
f1 1 4.3591
f1 2 3.4052
animal 1 0.0980
animal 2 -0.0193
animal 3 -0.0415
animal 4 -0.0091
animal 5 -0.1863
animal 6 0.1762
animal 7 -0.2500
animal 8 0.1821
stam_06.txt - 실행 화면
Pedigree file name : stam_pedi.dat
Performance file name : stam_perf.dat
Variance ratio : 2.0000000
level of fixed effect : 2
number of animals : 8
Round convergence criteria
1 2.06746776E-02
2 3.40968557E-03
3 3.08418588E-04
4 1.03507882E-04
5 2.29280631E-05
6 4.30842601E-06
7 7.19477498E-07
8 1.25288381E-07
9 3.24400240E-08
10 1.45340202E-08
11 8.40427372E-09
stam_01.txt : pedigree file (unsorted)
stam_02.txt : pedigree file (sorted)
stam_03.txt : performance file (sorted by animal)
stam_04.txt : performance file (sorted by fixed 1 effect)
stam_05.txt : solution file
stam_06.txt : running information file
관련 파일
혹시라도 프로그램을 개선하면 알려주세요....
'Animal Breeding > Fortran program' 카테고리의 다른 글
| 혈통 grouping하기 (0) | 2012.03.04 |
|---|---|
| 혈통 파일 트레이스(trace) 하기 (0) | 2012.01.25 |
| 단형질 임의회귀모형의 육종가 분할하기 (0) | 2011.05.26 |
| Single Trait Random Regression Model - 305 육종가 계산 (0) | 2010.02.23 |
| Single Trait Random Regression Model (0) | 2010.02.11 |