근교계수(inbreeding coefficient)와 Dii 구하기

혈통자료

 

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

 

위 자료를 ped.dat로 저장
자료 설명:animal, sire, dam

 

프로그램

 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! author : Kim Sidong
! date : 2009.3
! purpose : calculate inbreeding coefficient and Dii
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PROGRAM dii 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
REAL :: D
END TYPE

! define new data type for pedigree
TYPE ped_t
INTEGER :: s ! sire
INTEGER :: d ! dam
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(LEN = 256) :: pfn, ofn ! pedigree and output file name
INTEGER :: status ! status - open file and read data
INTEGER :: ii, ss, dd, oss, odd ! animal, sire, dam, old sire, old dam
REAL :: of ! old F

WRITE(*,*) 'Pedigree fine name : '
READ(*,*) pfn
WRITE(*,*) 'Output fine name : '
READ(*,*) ofn

OPEN(UNIT = 12, FILE = pfn, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status)
IF(status /= 0) STOP 'Pedigree File open error '
OPEN(UNIT = 13, FILE = ofn, STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status)
IF(status /= 0) STOP 'Output File open error'

! read pedigree and store it to ped array
no_of_ped = 1
DO
READ(UNIT = 12, FMT = *, IOSTAT = status), ii, ss, dd ! animal, sire, dam
IF(status /= 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
ped(ii)%s = ss ! store sire to ii th of ped array, animal number is index number
ped(ii)%d = dd ! store dam to ii th of ped array, animal number is index number
ped(ii)%f = 0.0 ! initialize
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
oss = ped(ii)%s
odd = ped(ii)%d
of = ped(ii)%f
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 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 ! out of range or process all animal
s = ped(idtd(offset)%id)%s ! sire of animal
d = ped(idtd(offset)%id)%d ! dam of animal

IF(s > 0) THEN ! if sire is known
idtd(no_of_idtd)%id = s ! add sire to idtd array
idtd(no_of_idtd)%t = idtd(offset)%t / 2.0 ! calculate t of sire
idtd(no_of_idtd)%d = getDii(s) ! calculate d of dam
no_of_idtd = no_of_idtd + 1 ! move next
IF (no_of_idtd == MAXIDTD) THEN ! out of range
WRITE(*,*) 'Increase MAXIDTD'
EXIT
END IF
END IF IF(d > 0) THEN ! if dam is known
idtd(no_of_idtd)%id = d ! add dam to idtd array
idtd(no_of_idtd)%t = idtd(offset)%t / 2.0 ! calculate t of dam
idtd(no_of_idtd)%d = getDii(d) ! calculate d of dam
no_of_idtd = no_of_idtd + 1 ! move next
IF (no_of_idtd == MAXIDTD) THEN ! out of range
WRITE(*,*) 'Increase MAXIDTD'
EXIT
END IF
END IF

offset = offset + 1 ! process next animal of idtd
END DO
END SUBROUTiNE trace

SUBROUTINE getF(id, f)
INTEGER, INTENT(IN) :: id
REAL, INTENT(OUT) :: f
INTEGER :: i, oid
REAL :: t, od

CALL trace(ID)

CALL gnomesort()

f = -1.0 ! inbreeding coefficient
i = 1
t = 0.0
oid = idtd(i)%id ! old id
od = idtd(i)%d ! old 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

FUNCTION getDii(id)
INTEGER, INTENT(IN) :: id ! animal
REAL :: getDii
INTEGER :: s, d, p ! sire, dam, sire + dam

s = ped(id)%s
d = ped(id)%d
p = s + d

IF (s == 0 .AND. d == 0) THEN ! both parents are 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 ! both 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

 

프로그램의 손쉬운 실행을 위하여ped.dat
ped.out
위를 dii.par로 저장