혈통 grouping하기

혈통 grouping하기

 

혈통이 혈연관계에 따라 몇 개의 그룹으로 나누어 지기도 한다. 한 개체를 기준으로 부모, 자식을 반복적으로 끝없이 연결할 때 연결되는 개체들이 하나의 그룹이 될 것이다.

 

예를 들어 다음과 같은 혈통을 보자

 

8 5 7
9 0 0
10 0 0
11 0 0
12 9 10
13 12 11
1 0 0
2 0 0
3 0 0
4 0 0
5 1 2
6 3 4
7 5 6
14 0 0

그림으로 그려보면 다음과 같다.

 

 

세 개의 그룹으로 나누어 진 것을 알 수 있다.

 

만일 1과 9의 부가 15이면 2 개의 그룹으로 나누어 진다. 단 15는 혈통의 개체에 등장하지 않는다.

 

혈통과 그림은 다음과 같다.

 

8 5 7
9 15 0
10 0 0
11 0 0
12 9 10
13 12 11
1 15 0
2 0 0
3 0 0
4 0 0
5 1 2
6 3 4
7 5 6
14 0 0

 

비록 15가 혈통의 개체에는 나타나지 않지만 그림에는 나타나고 15때문에 두 그룹이 연결된 것을 볼 수 있다.

 

프로그램

 

!===============================================================================
! 혈통 그루핑하기
! - 혈연 관계가 있는 개체들끼리 그루핑
! 프로그래머 : 박병호
! 사용 로직 : 그놈소트, 이진검색 이용
! 날짜 : 2012.3.4.
!===============================================================================PROGRAM grouping IMPLICIT NONE ! define new data type for pedigree
TYPE pedi_t
CHARACTER(16) :: a, s, d ! animal, sire, dam
INTEGER :: gn, f ! group number, flag
END TYPE

TYPE po_range_t
CHARACTER(16) :: p ! parent
INTEGER :: s, e, f ! start, end, flag
END TYPE

CHARACTER(256) :: pfn, ofn ! pedigree file name, output file name
INTEGER :: no_of_pedi ! number of animals in pedigree
INTEGER :: no_of_ubp ! number of sire and dam which don't appear in animals
INTEGER :: no_of_sd ! number of sire and dam
INTEGER :: no_of_usd ! number of unique sire and dam
TYPE(pedi_t), ALLOCATABLE :: pedi_asd(:) ! pedigree : animal - sire - dam - group
CHARACTER(16), ALLOCATABLE :: pedi_po(:,:) ! pedigree : parent - offsping
TYPE(po_range_t), ALLOCATABLE :: po_range(:) ! parent offsprng range
CHARACTER(16), ALLOCATABLE :: grp_pedi(:) ! temporary animal array
INTEGER :: cur_pos, ist_pos ! current position and insert position position in grp_pedi
INTEGER :: grp_number ! gropu number
INTEGER :: pos_a, pos_s, pos_d, pos_p, pos_o ! position of animal, sire, dam, parent and offspring
INTEGER :: rs, re ! start and end of po_range
INTEGER :: status ! io status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER :: i, j ! loop

IF (command_argument_count() /= 2) THEN
WRITE(*,*) "================================================================="
WRITE(*,*) " Usage : grouping pedigfree_file_name"
WRITE(*,*) "================================================================="
STOP
END IF !=============================================================================
! Reading pedigree file ...
!=============================================================================

! get pedigree filename
CALL get_command_argument(1, pfn) ! get pedigree data filename from the command line ! 자료 파일 열기
OPEN(UNIT = 20, FILE = pfn, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

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

! pedi 파일 자료 개수 알아내기
no_of_pedi = 0
DO
READ(UNIT = 20, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_pedi = no_of_pedi + 1
END DO ! 배열 할당
ALLOCATE(pedi_asd(no_of_pedi))
ALLOCATE(grp_pedi(no_of_pedi)) ! 혈통 자료를 다시 읽어, 배열에 저장
REWIND(20)
no_of_sd = 0
DO i = 1, no_of_pedi
READ(UNIT = 20, FMT = *, IOSTAT = status) pedi_asd(i)%a, pedi_asd(i)%s, pedi_asd(i)%d
IF (status /= 0) EXIT ! reach end of file
pedi_asd(i)%gn = 0
pedi_asd(i)%f = 0

! 자식을 가진 부모 개수 구하기(중복)
IF (pedi_asd(i)%s /= "0") THEN
no_of_sd = no_of_sd + 1
END IF ! 자식을 가진 부모 개수 구하기(중복)
IF (pedi_asd(i)%d /= "0") THEN
no_of_sd = no_of_sd + 1
END IF END DO !=============================================================================
! 부모 - 자식 배열 만들기
! Making parent - offspring array ...
!=============================================================================

! 배열 할당
ALLOCATE(pedi_po(no_of_sd, 2))
pedi_po = "" ! 혈통 배열을 다시 읽으며 부모 - 자손 배열 저장
j = 1
DO i = 1, no_of_pedi
IF (pedi_asd(i)%s /= "0") THEN
pedi_po(j,1) = pedi_asd(i)%s ! 부 저장
pedi_po(j,2) = pedi_asd(i)%a ! 자손 저장
j = j + 1
END IF IF (pedi_asd(i)%d /= "0") THEN
pedi_po(j,1) = pedi_asd(i)%d ! 모 저장
pedi_po(j,2) = pedi_asd(i)%a ! 자손 저장
j = j + 1
END IF END DO !=============================================================================
! 정렬
! Sorting pedigree ...
!=============================================================================
CALL gnomesort_asd() !=============================================================================
! Sorting parent - offspring array ...
!============================================================================= CALL gnomesort_po()

!=============================================================================
! 부모가 반복되는 구간 알아내기
! Making parent range ...
!=============================================================================

! 유일 sire와 dam 개수 알아내기
no_of_usd = 1
DO i = 2, no_of_sd
IF (pedi_po(i-1,1) /= pedi_po(i,1)) THEN
no_of_usd = no_of_usd + 1
END IF
END DO

! 배열 할당
ALLOCATE(po_range(no_of_usd))
po_range(:)%f = 0 ! 어떤 부모가 pedi_po에서 어디에서 어디까지 반복되는지 시작과 끝을 기록
j = 1
po_range(j)%p = pedi_po(1,1)
po_range(j)%s = 1
DO i = 2, no_of_sd
IF (pedi_po(i-1,1) /= pedi_po(i,1)) THEN
po_range(j)%e = i - 1
j = j + 1
po_range(j)%p = pedi_po(i,1)
po_range(j)%s = i
END IF
END DO
po_range(j)%e = no_of_sd !=============================================================================
! 개체에 안 나타나는 부모를 개체에 올리기
!=============================================================================
! 부모가 개체에 나타나는지 조사
no_of_ubp = 0
DO i = 1, no_of_usd
! 부모가 개체에 안 나타나면 표시
IF (binary_search(pedi_asd(:)%a, po_range(i)%p) == 0) THEN
po_range(i)%f = 1
no_of_ubp = no_of_ubp + 1
END IF
END DO ! 개체에 안 나타나는 부모가 있다면
IF (no_of_ubp > 0) THEN

! pedi_asd와 grp_pedi의 크기 조정
DEALLOCATE(pedi_asd)
ALLOCATE(pedi_asd(no_of_pedi + no_of_ubp)) DEALLOCATE(grp_pedi)
ALLOCATE(grp_pedi(no_of_pedi + no_of_ubp)) ! 혈통 파일 다시 읽기
REWIND(20)
DO i = 1, no_of_pedi
READ(UNIT = 20, FMT = *, IOSTAT = status) pedi_asd(i)%a, pedi_asd(i)%s, pedi_asd(i)%d
IF (status /= 0) EXIT ! reach end of file
pedi_asd(i)%gn = 0
pedi_asd(i)%f = 0
END DO ! 개체에 안 나타나는 부모를 pedi_asd에 넣기
j = no_of_pedi + 1
DO i = 1, no_of_usd
IF (po_range(i)%f == 1) THEN
pedi_asd(j)%a = po_range(i)%p
pedi_asd(j)%s = '0'
pedi_asd(j)%d = '0'
pedi_asd(j)%gn = 0
pedi_asd(j)%f = 0
j = j + 1
END IF
END DO
END IF CLOSE(20) no_of_pedi = no_of_pedi + no_of_ubp ! 정렬
CALL gnomesort_asd()
!=============================================================================
! Grouping pedigree ...
!============================================================================= grp_number = 1
DO i = 1, no_of_pedi

! 개체가 grouping 되지 않았다면
IF (pedi_asd(i)%gn == 0) THEN

! grouping을 위한 임시 pedi에 저장
grp_pedi = ""
grp_pedi(1) = pedi_asd(i)%a
pedi_asd(i)%f = 1 ! grp_pedi에 추가했다고 표시

cur_pos = 1
ist_pos = 2

! grp_pedi를 처음부터 읽으며 부모와 자손을 grp_pedi에 추가
! 더 이상 추가가 안 될 때까지 추가
DO WHILE (cur_pos <= (ist_pos - 1))
! add parents
! 개체가 pedi_asd의 어디에 있는지 위치 저장
pos_a = binary_search(pedi_asd(:)%a, grp_pedi(cur_pos)) ! pedi_asd 배열에서 찾고, 그룹 할당되지 않았다면
IF (pedi_asd(pos_a)%gn == 0) THEN
! 그룹 할당
pedi_asd(pos_a)%gn = grp_number

! sire가 있을 때
IF (pedi_asd(pos_a)%s /= '0') THEN

pos_s = binary_search(pedi_asd(:)%a, pedi_asd(pos_a)%s) IF (pedi_asd(pos_s)%f == 0) THEN
! grp_pedi에 추가
grp_pedi(ist_pos) = pedi_asd(pos_a)%s
pedi_asd(pos_s)%f = 1

! 개체 추가 위치 변경
ist_pos = ist_pos + 1
END IF
END IF

! dam이 있을 때
IF (pedi_asd(pos_a)%d /= '0') THEN

pos_d = binary_search(pedi_asd(:)%a, pedi_asd(pos_a)%d) IF (pedi_asd(pos_d)%f == 0) THEN
! grp_pedi에 추가
grp_pedi(ist_pos) = pedi_asd(pos_a)%d
pedi_asd(pos_d)%f = 1

! 개체 추가 위치 변경
ist_pos = ist_pos + 1
END IF
END IF
END IF

! add_offspring
! 자손이 있는 개체인지 조사
pos_p = binary_search(po_range(:)%p, grp_pedi(cur_pos)) ! 자손이 있는 개체라면
IF (pos_p /= 0) THEN
rs = po_range(pos_p)%s
re = po_range(pos_p)%e

DO j = rs, re
pos_o = binary_search(pedi_asd(:)%a, pedi_po(j,2)) IF (pos_o /= 0 .AND. pedi_asd(pos_o)%f == 0) THEN
grp_pedi(ist_pos) = pedi_po(j,2)
pedi_asd(pos_o)%f = 1

ist_pos = ist_pos + 1
END IF
END DO
END IF

! grp_pedi의 현재 위치 증가
cur_pos = cur_pos + 1
END DO ! increase group number
grp_number = grp_number + 1 END IF END DO !=============================================================================
! Writing output pedigree file
!============================================================================= ! get output filename
CALL get_command_argument(2, ofn) ! get pedigree data filename from the command line ! open output file
OPEN(UNIT = 30, FILE = ofn, STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'output file open failed -- error message : ', error_msg
STOP
END IF ! Writing the result to file
DO i = 1, no_of_pedi
WRITE(30,*) pedi_asd(i)%a,pedi_asd(i)%s,pedi_asd(i)%d,pedi_asd(i)%gn
END DO

CLOSE(30) CONTAINS ! binary search function
FUNCTION binary_search(a, value)
INTEGER :: binary_search
CHARACTER(16), INTENT(in) :: a(:), value
INTEGER :: low, high, mid

low = 1
high = size(a) binary_search = 0

DO WHILE (low <= high)
mid = (low + high) / 2
IF (a(mid) > value) THEN
high = mid - 1
ELSE IF (a(mid) < value) THEN
low = mid + 1
ELSE
binary_search = mid
EXIT
END IF
END DO END FUNCTION binary_search

! see wikipedia for gnome sort
SUBROUTINE gnomesort_asd()
INTEGER :: i, j, max_pedi
TYPE(pedi_t) :: temp
i = 2
j = 3
max_pedi = no_of_pedi - 1
DO
IF (i > no_of_pedi) EXIT
IF (pedi_asd(i-1)%a <= pedi_asd(i)%a) THEN
i = j
j = j + 1
ELSE
temp = pedi_asd(i-1)
pedi_asd(i-1) = pedi_asd(i)
pedi_asd(i) = temp
i = i - 1
if (i == 1) i = 2
END IF
END DO

END SUBROUTINE

SUBROUTINE gnomesort_po()
INTEGER :: i, j, max_pedi
CHARACTER(16) :: temp(2)
i = 2
j = 3
max_pedi = no_of_sd - 1
DO
IF (i > no_of_sd) EXIT
IF (pedi_po(i-1,1) <= pedi_po(i,1)) THEN
i = j
j = j + 1
ELSE
temp(1) = pedi_po(i-1,1)
temp(2) = pedi_po(i-1,2)
pedi_po(i-1,1) = pedi_po(i,1)
pedi_po(i-1,2) = pedi_po(i,2)
pedi_po(i,1) = temp(1)
pedi_po(i,2) = temp(2)
i = i - 1
if (i == 1) i = 2
END IF
END DO

END SUBROUTINE

END PROGRAM

 

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

 

 

 

프로그램과 자료

1331192203_grouping.zip

다운로드

1331192203_grouping.zip