근교계수와 dii구하기 - Fortran 프로그램

 

근교계수와 dii 구하기

 

근교계수는 조상부터 구합니다.

 

조상의 근교계수가 있어야 dii를 구할 수 있습니다.

 

dii를 이용해서 근교계수를 구합니다.

 

Meuwissen and Luo(1992)의 알고리듬을 이용하여 근교계수 계산

 

Relationship Matrix(A)의 역행렬을 구성할 때 개체의 dii가 필요

: 개체의 근교계수 계산 -> 개체의 dii 계산 -> Relationship Matrix의 역행렬 구성순입니다.

 

여기서는 renumber된 혈통파일을 읽어 근교계수와 dii를 출력하는 프로그램

 

Pedigree

 

1 0 0
2 0 0
3 1 2
4 1 0
5 4 3
6 5 2

위 혈통을 ped.dat 로 저장

 

프로그램

 

PROGRAM dii

! program name : dii
! programmer : Sidong Kim
! usage : dii pedigree_file
! purpose : read pedigree file, calculate inbreeding coefficient and dii
! Pedigree file must be sorted by generation
! Animal ID must be a number
! date : 2009.4.4.
! update
! 2010.1.27. Park Byoungho

IMPLICIT NONE

! define max id-t-d vector length
INTEGER, PARAMETER :: MAXIDTD = 65536
! define max animal in the pedigree
INTEGER, PARAMETER :: MAXANIM = 100000

! define new data type for id-t-d vector
TYPE idtd_t
INTEGER ID
REAL T, D
END TYPE

! define new data type for pedigree
TYPE ped_t
INTEGER :: s, d
real :: f
END TYPE

! declare idtd vector with size = MAXIDTD
TYPE(idtd_t), DIMENSION(MAXIDTD) :: idtd
INTEGER :: no_of_idtd

! declare pedigree vector with size = MAXANIM
TYPE(ped_t), DIMENSION(MAXANIM) :: ped
INTEGER :: no_of_ped

!*****************************************
! main
!*****************************************

CHARACTER(256) :: pfn
INTEGER :: fostat
INTEGER :: ii, ss, dd, oss, odd
REAL :: of

! get pedigree filename - updated at 2010.1.27 by Park Byoungho
CALL get_command_argument(1, pfn) ! get pedigree data filename from the command line

! open pedigree filename
OPEN(UNIT = 12, FILE = pfn, STATUS = 'OLD', ACTION = 'READ', IOSTAT = fostat)
IF(fostat /= 0) STOP 'File open error (input)'
!open output filename
OPEN(UNIT = 13, FILE = TRIM(pfn) // '.ico', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = fostat)
IF(fostat /= 0) STOP 'File open error (output)'
WRITE(*,*) 'Output file : ', TRIM(pfn), '.ico'

! read pedigree file and store it into ped array
no_of_ped = 1
DO
! read animal, sire, dam from pedigree file
READ(UNIT = 12, FMT = *, IOSTAT = fostat), ii, ss, dd
IF(fostat /= 0) EXIT

IF (ii .NE. no_of_ped) THEN
WRITE(*,*) 'Pedigree should be numbered and sorted in chronological order'
STOP 'Pedigree file error'
END IF

! store sire and dam into ped array
ped(ii)%s = ss
ped(ii)%d = dd
ped(ii)%f = 0.0

no_of_ped = no_of_ped + 1
END DO

CLOSE(12)

ii = 1
oss = -1
odd = -1
of = 0.0

DO
IF(ii >= no_of_ped) EXIT
! check animal with the same sire and the sam dam
IF (oss == ped(ii)%s .AND. odd == ped(ii)%d) THEN
! just give F value of previous animal
ped(ii)%f = of
ELSE
!for new animal, F should be calculated
CALL getF(ii, ped(ii)%f)
! store values for an animal with same sire and dam of animal, ii
of = ped(ii)%f
oss = ped(ii)%s
odd = ped(ii)%d
END IF

WRITE(UNIT = 13, FMT = '(1X, 3(I6, 1X), F12.10, 1X, F12.10)'), ii, ped(ii)%s, ped(ii)%d, 1.0/getDii(ii), ped(ii)%f

ii = ii + 1
END DO

CLOSE(13)

!************************************
! subroutine
!************************************

CONTAINS

SUBROUTINE getF(id, f)
! purpose : calculate inbreeding coefficient
INTEGER, INTENT(IN) :: id
REAL, INTENT(OUT) :: f
INTEGER :: i, oid
REAL :: t, od

CALL trace(ID)

CALL gnomesort()

f = -1.0
i = 1
t = 0.0
oid = idtd(i)%id
od = idtd(i)%d

DO
IF(i >= no_of_idtd) EXIT

IF(oid .NE. idtd(i)%id) THEN
f = f + t**2 * od
t = idtd(i)%t
oid = idtd(i)%id
od = idtd(i)%d
ELSE
t = t + idtd(i)%t
END IF

i = i + 1
END DO

! process last line
f = f + t**2 * od
END SUBROUTiNE getF

SUBROUTINE trace(id)
INTEGER, INTENT(IN) :: id
INTEGER :: offset, s, d

no_of_idtd = 1
idtd(no_of_idtd)%id = id
idtd(no_of_idtd)%t = 1.0
idtd(no_of_idtd)%d = getDii(id)
no_of_idtd = no_of_idtd + 1

offset = 1
DO
IF (offset >= MAXIDTD .OR. offset >= no_of_idtd) EXIT

! find sire and dam
s = ped(idtd(offset)%id)%s
d = ped(idtd(offset)%id)%d

! add sire to idtd array
IF(s > 0) THEN
idtd(no_of_idtd)%id = s
idtd(no_of_idtd)%t = idtd(offset)%t / 2.0
idtd(no_of_idtd)%d = getDii(s)
no_of_idtd = no_of_idtd + 1
IF (no_of_idtd == MAXIDTD) THEN
WRITE(*,*) 'Increase MAXIDTD'
EXIT
END IF
END IF

! add dam to idtd array
IF(d > 0) THEN
idtd(no_of_idtd)%id = d
idtd(no_of_idtd)%t = idtd(offset)%t / 2.0
idtd(no_of_idtd)%d = getDii(d)
no_of_idtd = no_of_idtd + 1
IF (no_of_idtd == MAXIDTD) THEN
WRITE(*,*) 'Increase MAXIDTD'
EXIT
END IF
END IF

offset = offset + 1
END DO
END SUBROUTiNE trace

FUNCTION getDii(id)
INTEGER, INTENT(IN) :: id
REAL :: getDii
INTEGER :: p, s, d

! find sire and dam
s = ped(id)%s
d = ped(id)%d
p = s + d
IF (s == 0 .AND. d == 0) THEN ! no parents is unknown
getDii = 1.0
ELSE IF( s == 0 .OR. d == 0) THEN ! one parent is known
getDii = 0.75 - 0.25 * ped(p)%f
ELSE ! parents are known
getDii = 0.5 - 0.25 * (ped(s)%f + ped(d)%f)
END IF
END FUNCTION getDii

! see wikipedia for gnome sort
SUBROUTINE gnomesort()
INTEGER :: i, j, max_idtd
TYPE(idtd_t) :: temp
i = 2
j = 3
max_idtd = no_of_idtd - 1
DO
IF (i > max_idtd) EXIT
IF (idtd(i-1)%id >= idtd(i)%id) THEN
i = j
j = j + 1
ELSE
temp = idtd(i-1)
idtd(i-1) = idtd(i)
idtd(i) = temp
i = i - 1
if (i == 1) i = 2
END IF
END DO
END SUBROUTINE gnomesort

END PROGRAM dii

위 프로그램을 dii.f95로 저장

 

컴파일 방법 : gfortran dii.f95 -o dii

 

결과

 

1 0 0 1.0000000000 0.0000000000
2 0 0 1.0000000000 0.0000000000
3 1 2 2.0000000000 0.0000000000
4 1 0 1.3333333731 0.0000000000
5 4 3 2.0000000000 0.1250000000
6 5 2 2.1333334446 0.1250000000

4th col : dii

5th col : inbreeding coefficient

 

프로그램 컴파일 및 실행화면

 

 

 

파일

 

 

1264649670_dii.zip

다운로드

1264649670_dii.zip