혈통 그루핑

유전 평가하고 상관 없으나 혈통이 있을 때 부모 또는 자식으로 연결하였을 때 몇 개의 그룹으로 나누어질지 궁금할 때가 있다. 물론 그림을 그리면 되긴 하지만, 몇십 개 이상은 그림을 그리기가 힘들다. 물론 포트란 소스로 올린 적도 있으나 컴파일하기도 힘들어서 누가 쓸까 싶다. 그래서 간단히 R로 만들어 보았다. 포트란보다 많이 느리겠지만 ...

먼저 혈통 자료

animal,sire,dam
a1,0,0
a2,0,0
a3,0,0
a4,0,0
a9,0,0
a10,0,0
a11,0,0
a14,0,0
a5,a1,a2
a6,a3,a4
a7,a5,a6
a12,a9,a10
a13,a12,a11
a8,a5,a7

그림을 그려 보면 세 개의 그룹으로 이루어진 것을 알 수 있다.

R 소스

#
# group.R
#

# 혈통 자료 읽기
pedi <-
    read.table(
        "pedi_group.csv",
        header = TRUE,
        sep = ",",
        stringsAsFactors = FALSE
    )
pedi

pedi2 = cbind(pedi, flag = rep(0, nrow(pedi)), group = rep(0, nrow(pedi)))
pedi2

group_num = 0
# pedi2 파일의 개체 하나씩 처리하기
for (i in 1:nrow(pedi2)) {

    # 처리한 개체가 아니라면
    if (pedi2[i, "flag"] == 0) {
        
        #group_num 증가
        group_num = group_num + 1
        
        # pedi2 파일의 개체를 temp에 넣기
        temp = data.frame(animal = pedi2[i, "animal"], stringsAsFactors = FALSE)
        
        # temp를 읽어 내려가면 처리하기
        j = 1
        while (j <= nrow(temp)) {
            
            # temp의 개체가 pedi에 존재하는가 그리고 처리를 안 했는가
            index = which(pedi2$animal == temp[j, "animal"])
            if (length(index) != 0 &&
                pedi2[index, "flag"] == 0) {
                # pedid2에 개체를 처리했다고 표시(flag)
                pedi2[index, "flag"] = 1
                
                # pedid2에 group_num 입력
                pedi2[index, "group"] = group_num
                
                # 개체의 아비 찾기
                sire = pedi2[index, "sire"]
                
                # 아비가 있다면
                if (sire != 0) {
                    # temp에 sire 추가
                    temp = rbind(temp, sire)
                }
                
                # 개체의 어미 찾기
                dam = pedi2[index, "dam"]
                
                # 어미가 있다면
                if (dam != 0) {
                    # temp에 dam 추가
                    temp = rbind(temp, dam)
                }
                
                # 아비라고 생각하고 자식 찾기
                index = which(pedi2$sire == temp[j, "animal"] & pedi2[which(pedi2$sire == temp[j, "animal"]), "flag"] == 0)
                progeny = as.data.frame(pedi2[index, "animal"])
                names(progeny) = names(temp)                                        
                temp = rbind(temp, progeny)

                # 어미라고 생각하고 자식 찾기
                index = which(pedi2$dam == temp[j, "animal"] & pedi2[which(pedi2$dam == temp[j, "animal"]), "flag"] == 0)
                progeny = as.data.frame(pedi2[index, "animal"])
                names(progeny) = names(temp)                                        
                temp = rbind(temp, progeny)
                
                # 1000개마다 처리 중임을 표시
                if( (j %% 1000) == 1 ){
                    print(paste("Group " , group_num , ": ", j, " processing"))
                }

            }
            j = j + 1
        }
    }
}

pedi2 = pedi2[order(pedi2$group), ]
pedi2

 

실행 결과

> # 혈통 자료 읽기
> pedi <-
+     read.table(
+         "pedi_group.csv",
+         header = TRUE,
+         sep = ",",
+         stringsAsFactors = FALSE
+     )
> pedi
   animal sire dam
1      a1    0   0
2      a2    0   0
3      a3    0   0
4      a4    0   0
5      a9    0   0
6     a10    0   0
7     a11    0   0
8     a14    0   0
9      a5   a1  a2
10     a6   a3  a4
11     a7   a5  a6
12    a12   a9 a10
13    a13  a12 a11
14     a8   a5  a7
> 
> pedi2 = cbind(pedi, flag = rep(0, nrow(pedi)), group = rep(0, nrow(pedi)))
> pedi2
   animal sire dam flag group
1      a1    0   0    0     0
2      a2    0   0    0     0
3      a3    0   0    0     0
4      a4    0   0    0     0
5      a9    0   0    0     0
6     a10    0   0    0     0
7     a11    0   0    0     0
8     a14    0   0    0     0
9      a5   a1  a2    0     0
10     a6   a3  a4    0     0
11     a7   a5  a6    0     0
12    a12   a9 a10    0     0
13    a13  a12 a11    0     0
14     a8   a5  a7    0     0
> 
> group_num = 0
> # pedi2 파일의 개체 하나씩 처리하기
> for (i in 1:nrow(pedi2)) {
+ 
+     # 처리한 개체가 아니라면
+     if (pedi2[i, "flag"] == 0) {
+         
+         #group_num 증가
+         group_num = group_num + 1
+         
+         # pedi2 파일의 개체를 temp에 넣기
+         temp = data.frame(animal = pedi2[i, "animal"], stringsAsFactors = FALSE)
+         
+         # temp를 읽어 내려가면 처리하기
+         j = 1
+         while (j <= nrow(temp)) {
+             
+             # temp의 개체가 pedi에 존재하는가 그리고 처리를 안 했는가
+             index = which(pedi2$animal == temp[j, "animal"])
+             if (length(index) != 0 &&
+                 pedi2[index, "flag"] == 0) {
+                 # pedid2에 개체를 처리했다고 표시(flag)
+                 pedi2[index, "flag"] = 1
+                 
+                 # pedid2에 group_num 입력
+                 pedi2[index, "group"] = group_num
+                 
+                 # 개체의 아비 찾기
+                 sire = pedi2[index, "sire"]
+                 
+                 # 아비가 있다면
+                 if (sire != 0) {
+                     # temp에 sire 추가
+                     temp = rbind(temp, sire)
+                 }
+                 
+                 # 개체의 어미 찾기
+                 dam = pedi2[index, "dam"]
+                 
+                 # 어미가 있다면
+                 if (dam != 0) {
+                     # temp에 dam 추가
+                     temp = rbind(temp, dam)
+                 }
+                 
+                 # 아비라고 생각하고 자식 찾기
+                 index = which(pedi2$sire == temp[j, "animal"] & pedi2[which(pedi2$sire == temp[j, "animal"]), "flag"] == 0)
+                 progeny = as.data.frame(pedi2[index, "animal"])
+                 names(progeny) = names(temp)                                        
+                 temp = rbind(temp, progeny)
+ 
+                 # 어미라고 생각하고 자식 찾기
+                 index = which(pedi2$dam == temp[j, "animal"] & pedi2[which(pedi2$dam == temp[j, "animal"]), "flag"] == 0)
+                 progeny = as.data.frame(pedi2[index, "animal"])
+                 names(progeny) = names(temp)                                        
+                 temp = rbind(temp, progeny)
+                 
+                 # 1000개마다 처리 중임을 표시
+                 if( (j %% 1000) == 1 ){
+                     print(paste("Group " , group_num , ": ", j, " processing"))
+                 }
+ 
+             }
+             j = j + 1
+         }
+     }
+ }
[1] "Group  1 :  1  processing"
[1] "Group  2 :  1  processing"
[1] "Group  3 :  1  processing"
> 
> pedi2 = pedi2[order(pedi2$group), ]
> pedi2
   animal sire dam flag group
1      a1    0   0    1     1
2      a2    0   0    1     1
3      a3    0   0    1     1
4      a4    0   0    1     1
9      a5   a1  a2    1     1
10     a6   a3  a4    1     1
11     a7   a5  a6    1     1
14     a8   a5  a7    1     1
5      a9    0   0    1     2
6     a10    0   0    1     2
7     a11    0   0    1     2
12    a12   a9 a10    1     2
13    a13  a12 a11    1     2
8     a14    0   0    1     3
>

 

실제 자료일 경우 출력을 하면 안 될 것이고, 대신 head, tail, nrow 등으로 파악해야 할 것이다. 그리고 처리 중임을 표시하는 부분도 1,000개마다 할 지, 10,000개마다 할 지 수정해야 한다.