R을 이용한 혈통 파일 처리
혈통 파일을 처리한다는 것은 혈통 파일의 오류를 찾아내고, 혈연 계수 행렬(Numerator Relationship Matrix, NRM)을 구하거나, NRM의 역행렬을 구하거나, 근교계수 등을 구하는 것을 말한다.
1) 혈통 파일의 처리는 전체 혈통에서 필요한 혈통을 추출하는데서 시작한다. 물론 전체 혈통 파일을 처리할 수도 있으나 관심이 가는 개체의 조상을 추적하면 다루어야 할 개체수를 줄일 수 있다.(Trace) 2) 다음으로 이렇게 추출한 개체들을 조상에서 자손까지 정렬하여야 한다. 생일이 있으면 쉽게 할 수 있으나 생일을 신뢰하기 어렵거나 생일이 아예 없는 경우도 있다. 생일이 없는 경우 선조의 세대를 개체의 세대보다 높게 매기고 세대를 정렬한다.(Stack) 3) 다음으로 계산의 편리를 위하여 개체의 이름을 일련번호로 바꾼다.(Renumber) 4) 이렇게 일련번호로 된 혈통을 가지고 NRM, NRM의 역행렬, 근교계수 등을 구한다.
자세한 이론에 대해서는 농촌진흥청에서 나온 “동물육종을 위한 혈통분석 알고리듬과 프로그래밍”을 참고하기 바란다.
http://lib.rda.go.kr/search/mediaView.do?mets_no=000000010979
■ 전체혈통
animal,sire,dam
a1,0,0
a2,0,0
a3,a1,a2
a6,a5,a2
a9,0,0
a10,0,0
a11,0,0
a14,0,0
a12,a11,a10
a15,a13,a12
a16,a15,a14
a4,a1,0
a5,a4,a3
a8,a7,a4
a7,a5,a6
a13,a9,a6
■ Key animals
key
a3
a4
a5
a6
■ Trace Pedigree
- R 프로그램
# trace pedigree
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# read all pedigree
pedi_01 <- read.table("pedi_all.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print pedigree
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
no_of_animals
# add flag column
pedi_02 = cbind(pedi_01, flag = rep(0, no_of_animals))
# print pedigree
dim(pedi_02)
head(pedi_02)
tail(pedi_02)
# read key animals
key_animals <- read.table("key_animals.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print key animal
dim(key_animals)
head(key_animals)
tail(key_animals)
# no_of_keys
no_of_keys = nrow(key_animals)
no_of_keys
# process each animal of key animals
for (i in 1 : no_of_keys) {
# key 파일의 개체를 temp에 넣기
temp = data.frame(animal = c(key[i, "key"]), stringsAsFactors = FALSE)
# temp를 읽어 내려가면서 처리하기
# j가 temp 끝까지 처리
j = 1
while (j <= nrow(temp)) {
# 개체의 위치 저장
index = which(pedi_02$animal == temp[j, "animal"])
# key에 있는 개체를 pedigree에서 찾았다면
if (length(index) != 0) {
# 개체를 아직 처리하지 않았다면
if (pedi_02[index, "flag"] == 0) {
# peid_02에 개체를 처리했다고 표시(flag)
pedi_02[index, "flag"] = 1
# 개체의 아비 찾기
sire = pedi_02[index, "sire"]
# 아비가 있는 건가 아닌가
if (sire != 0) {
# temp에 sire 추가
temp = rbind(temp, c(sire))
}
# 개체의 어미 찾기
dam = pedi_02[index, "dam"]
# 어미가 있는 건가 아닌가
if (dam != 0) {
# temp에 dam 추가
temp = rbind(temp, c(dam))
}
}
}
# key에 있는 개체를 pedigree에서 찾지 못했다면
else{
print(paste("개체 ", temp[j, "animal"], "을 혈통에서 찾지 못함"))
}
j = j + 1
}
#print(temp)
}
# pedi_02의 flag가 1인 것만 추출
pedi_03 = pedi_02[which(pedi_02$flag == 1),]
# print pedigree
dim(pedi_03)
head(pedi_03)
tail(pedi_03)
# animal, sire, dam만 추출
pedi_04 = pedi_03[,c("animal", "sire", "dam")]
# print pedigree
dim(pedi_04)
head(pedi_04)
tail(pedi_04)
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".trc"))
write.table(pedi_04, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 실행 결과
> # trace pedigree
>
> # set working directory
> setwd("D:/temp")
>
> # print working directory
> getwd()
[1] "D:/temp"
>
> # read all pedigree
> pedi_01 <- read.table("pedi_all.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
>
> # print pedigree
> dim(pedi_01)
[1] 16 3
> head(pedi_01)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
5 a9 0 0
6 a10 0 0
> tail(pedi_01)
animal sire dam
11 a16 a15 a14
12 a4 a1 0
13 a5 a4 a3
14 a8 a7 a4
15 a7 a5 a6
16 a13 a9 a6
>
> # number of animal
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 16
>
> # add flag column
> pedi_02 = cbind(pedi_01, flag = rep(0, no_of_animals))
>
> # print pedigree
> dim(pedi_02)
[1] 16 4
> head(pedi_02)
animal sire dam flag
1 a1 0 0 0
2 a2 0 0 0
3 a3 a1 a2 0
4 a6 a5 a2 0
5 a9 0 0 0
6 a10 0 0 0
> tail(pedi_02)
animal sire dam flag
11 a16 a15 a14 0
12 a4 a1 0 0
13 a5 a4 a3 0
14 a8 a7 a4 0
15 a7 a5 a6 0
16 a13 a9 a6 0
>
> # read key animals
> key_animals <- read.table("key_animals.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
>
> # print key animal
> dim(key_animals)
[1] 4 1
> head(key_animals)
key
1 a3
2 a4
3 a5
4 a6
> tail(key_animals)
key
1 a3
2 a4
3 a5
4 a6
>
> # no_of_keys
> no_of_keys = nrow(key_animals)
> no_of_keys
[1] 4
>
> # process each animal of key animals
> for (i in 1 : no_of_keys) {
+ # key 파일의 개체를 temp에 넣기
+ temp = data.frame(animal = c(key_animals[i, "key"]), stringsAsFactors = FALSE)
+
+ # temp를 읽어 내려가면서 처리하기
+ # j가 temp 끝까지 처리
+ j = 1
+ while (j <= nrow(temp)) {
+ # 개체의 위치 저장
+ index = which(pedi_02$animal == temp[j, "animal"])
+
+ # key에 있는 개체를 pedigree에서 찾았다면
+ if (length(index) != 0) {
+ # 개체를 아직 처리하지 않았다면
+ if (pedi_02[index, "flag"] == 0) {
+ # peid_02에 개체를 처리했다고 표시(flag)
+ pedi_02[index, "flag"] = 1
+
+ # 개체의 아비 찾기
+ sire = pedi_02[index, "sire"]
+
+ # 아비가 있는 건가 아닌가
+ if (sire != 0) {
+ # temp에 sire 추가
+ temp = rbind(temp, c(sire))
+ }
+
+ # 개체의 어미 찾기
+ dam = pedi_02[index, "dam"]
+
+ # 어미가 있는 건가 아닌가
+ if (dam != 0) {
+ # temp에 dam 추가
+ temp = rbind(temp, c(dam))
+ }
+ }
+ }
+ # key에 있는 개체를 pedigree에서 찾지 못했다면
+ else{
+ print(paste("개체 ", temp[j, "animal"], "을 혈통에서 찾지 못함"))
+ }
+
+ j = j + 1
+ }
+ #print(temp)
+ }
>
> # pedi_02의 flag가 1인 것만 추출
> pedi_03 = pedi_02[which(pedi_02$flag == 1),]
>
> # print pedigree
> dim(pedi_03)
[1] 6 4
> head(pedi_03)
animal sire dam flag
1 a1 0 0 1
2 a2 0 0 1
3 a3 a1 a2 1
4 a6 a5 a2 1
12 a4 a1 0 1
13 a5 a4 a3 1
> tail(pedi_03)
animal sire dam flag
1 a1 0 0 1
2 a2 0 0 1
3 a3 a1 a2 1
4 a6 a5 a2 1
12 a4 a1 0 1
13 a5 a4 a3 1
>
> # animal, sire, dam만 추출
> pedi_04 = pedi_03[,c("animal", "sire", "dam")]
>
> # print pedigree
> dim(pedi_04)
[1] 6 3
> head(pedi_04)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
12 a4 a1 0
13 a5 a4 a3
> tail(pedi_04)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
12 a4 a1 0
13 a5 a4 a3
>
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".trc"))
> write.table(pedi_04, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 결과 파일(pedi_20200410.trc)
animal,sire,dam
a1,0,0
a2,0,0
a3,a1,a2
a6,a5,a2
a4,a1,0
a5,a4,a3
■ Stack Pedigree
- R 프로그램
# stack pedigree and check errors
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.trc", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print pedigree
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# no_of_animals
no_of_animals = nrow(pedi_01)
no_of_animals
# copy pedigree
pedi_02 = pedi_01
# 아비 어미가 모두 개체에 나오게 하자
for (i in 1 : no_of_animals) {
# 아비가 0이 아니고, 개체에 등장하지 않으면
if (pedi_01[i, 2] != 0 && length(which(pedi_02$animal == pedi_01[i, 2])) == 0) {
# add sire to animal
pedi_02 = rbind(pedi_02, c(pedi_01[i, 2], 0, 0))
}
# 어미가 0이 아니고, 개체에 등장하지 않으면
if (pedi_01[i, 3] != 0 && length(which(pedi_02$animal == pedi_01[i, 3])) == 0) {
# add dam to animal
pedi_02 = rbind(pedi_02, c(pedi_01[i, 3], 0, 0))
}
}
# print pedigree
dim(pedi_02)
tail(pedi_02)
# 중복된 행 제거
pedi_03 = unique(pedi_02)
# number of animals in pedi_02 and pedi_03
nrow(pedi_02)
no_of_animals = nrow(pedi_03)
no_of_animals
# dup, sex, gen 컬럼 추가하기
pedi_04 = cbind(pedi_03, dup = rep(0, no_of_animals),
sex = rep(0, no_of_animals),
asd = rep(0, no_of_animals),
gen = rep(1, no_of_animals))
# print pedigree
dim(pedi_04)
head(pedi_04)
tail(pedi_04)
################################
# check duplication error
################################
# sorting
pedi_04 = pedi_04[order(pedi_04$animal), ]
head(pedi_04)
tail(pedi_04)
for (i in 2 : no_of_animals) {
# 이전 개체와 현재 개체가 같은 경우
if (pedi_04[i - 1, "animal"] == pedi_04[i, "animal"]){
# 아비가 다르거나 어미가 다르다면
if (pedi_04[i - 1, "sire"] != pedi_04[i, "sire"]
|| pedi_04[i - 1, "dam"] != pedi_04[i, "dam"]){
pedi_04[i - 1, "dup"] = 1
pedi_04[i , "dup"] = 1
print(paste("개체 ", pedi_04[i, "animal"] , "가 두 번 등장, 혈통이 다름"))
}
}
}
# print animals with duplication error
pedi_04[which(pedi_04$dup == 1), ]
##########################
# check sex error
##########################
for (i in 1 : no_of_animals) {
# sire와 dam을 저장
sire = pedi_04[i, "sire"]
dam = pedi_04[i, "dam"]
# sire가 있다면
if (sire != 0) {
# sire가 개체의 어디에 나오는가
index = which(pedi_04$animal == sire)
# sire의 sex가 결정되지 않았다면 수로 표시
if (pedi_04[index, "sex"] == 0) {
pedi_04[index, "sex"] = 2
}
# sire의 sex가 암으로 되어 있다면 error표시
else if (pedi_04[index, "sex"] == 1) {
pedi_04[index, "sex"] = 3
}
}
# dma이 있다면
if (dam != 0) {
# dam이 개체의 어디에 나오는가
index = which(pedi_04$animal == dam)
# dam의 sex가 결정되지 않았다면 암으로 표시
if (pedi_04[index, "sex"] == 0) {
pedi_04[index, "sex"] = 1
}
# dam의 sex가 수로 되어 있다면 error 표시
else if (pedi_04[index, "sex"] == 2) {
pedi_04[index, "sex"] = 3
}
}
}
# print animals with sex error
pedi_04[which(pedi_04$sex == 3), ]
########################################
# check animal = sire = dam error
########################################
# animal = sire인 개체의 위치 저장
index = which(pedi_04$animal == pedi_04$sire)
# animal = sire인 개체가 있다면 error 추가
if (length(index) != 0) {
pedi_04[index, "asd"] = pedi_04[index, "asd"] + 1
}
# animal = dam인 개체의 위치 저장
index = which(pedi_04$animal == pedi_04$dam)
# animal = dam인 개체가 있다면 error 추가
if (length(index) != 0) {
pedi_04[index, "asd"] = pedi_04[index, "asd"] + 2
}
# sire와 dam이 0이 아니고, sire = dam인 개체의 위치 저장
index = which(pedi_04$sire == pedi_04$dam & pedi_04$sire != 0 & pedi_04$dam != 0)
# sire = dam이면 error 추가
if (length(index) != 0) {
pedi_04[index, "asd"] = pedi_04[index, "asd"] + 4
}
# print animals with a = s = d error
pedi_04[which(pedi_04$asd != 0), ]
#########################
# check cycling error
#########################
for (i in 1 : 50) {
for (j in 1 : no_of_animals) {
sire = pedi_04[j, "sire"]
dam = pedi_04[j, "dam"]
# 아비가 0이 아니고, 아비의 세대수가 개체의 세대수보다 같거나 작다면
index = which(pedi_04$animal == sire)
if ((sire != 0) && (pedi_04[index, "gen"] <= pedi_04[j, "gen"])) {
pedi_04[index, "gen"] = pedi_04[j, "gen"] + 1
}
# 어미가 0이 아니고, 어미의 세대수가 개체의 세대수보다 같거나 작다면
index = which(pedi_04$animal == dam)
if ((dam != 0) && (pedi_04[index, "gen"] <= pedi_04[j, "gen"])) {
pedi_04[index, "gen"] = pedi_04[j, "gen"] + 1
}
}
}
# sort by generation
pedi_04 = pedi_04[order(-pedi_04$gen),]
# print pedigree : if generation is very large, there may be an error
head(pedi_04, 10)
tail(pedi_04, 10)
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".stk"))
write.table(pedi_04, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.trc", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # print pedigree
> dim(pedi_01)
[1] 6 3
> head(pedi_01)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
5 a4 a1 0
6 a5 a4 a3
> tail(pedi_01)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
5 a4 a1 0
6 a5 a4 a3
> # no_of_animals
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 6
> # copy pedigree
> pedi_02 = pedi_01
> # 아비 어미가 모두 개체에 나오게 하자
> for (i in 1 : no_of_animals) {
+
+ # 아비가 0이 아니고, 개체에 등장하지 않으면
+ if (pedi_01[i, 2] != 0 && length(which(pedi_02$animal == pedi_01[i, 2])) == 0) {
+ # add sire to animal
+ pedi_02 = rbind(pedi_02, c(pedi_01[i, 2], 0, 0))
+ }
+
+ # 어미가 0이 아니고, 개체에 등장하지 않으면
+ if (pedi_01[i, 3] != 0 && length(which(pedi_02$animal == pedi_01[i, 3])) == 0) {
+ # add dam to animal
+ pedi_02 = rbind(pedi_02, c(pedi_01[i, 3], 0, 0))
+ }
+ }
> # print pedigree
> dim(pedi_02)
[1] 6 3
> tail(pedi_02)
animal sire dam
1 a1 0 0
2 a2 0 0
3 a3 a1 a2
4 a6 a5 a2
5 a4 a1 0
6 a5 a4 a3
> # 중복된 행 제거
> pedi_03 = unique(pedi_02)
> # number of animals in pedi_02 and pedi_03
> nrow(pedi_02)
[1] 6
> no_of_animals = nrow(pedi_03)
> no_of_animals
[1] 6
> # dup, sex, gen 컬럼 추가하기
> pedi_04 = cbind(pedi_03, dup = rep(0, no_of_animals),
+ sex = rep(0, no_of_animals),
+ asd = rep(0, no_of_animals),
+ gen = rep(1, no_of_animals))
> # print pedigree
> dim(pedi_04)
[1] 6 7
> head(pedi_04)
animal sire dam dup sex asd gen
1 a1 0 0 0 0 0 1
2 a2 0 0 0 0 0 1
3 a3 a1 a2 0 0 0 1
4 a6 a5 a2 0 0 0 1
5 a4 a1 0 0 0 0 1
6 a5 a4 a3 0 0 0 1
> tail(pedi_04)
animal sire dam dup sex asd gen
1 a1 0 0 0 0 0 1
2 a2 0 0 0 0 0 1
3 a3 a1 a2 0 0 0 1
4 a6 a5 a2 0 0 0 1
5 a4 a1 0 0 0 0 1
6 a5 a4 a3 0 0 0 1
> # sorting
> pedi_04 = pedi_04[order(pedi_04$animal), ]
> head(pedi_04)
animal sire dam dup sex asd gen
1 a1 0 0 0 0 0 1
2 a2 0 0 0 0 0 1
3 a3 a1 a2 0 0 0 1
5 a4 a1 0 0 0 0 1
6 a5 a4 a3 0 0 0 1
4 a6 a5 a2 0 0 0 1
> tail(pedi_04)
animal sire dam dup sex asd gen
1 a1 0 0 0 0 0 1
2 a2 0 0 0 0 0 1
3 a3 a1 a2 0 0 0 1
5 a4 a1 0 0 0 0 1
6 a5 a4 a3 0 0 0 1
4 a6 a5 a2 0 0 0 1
> for (i in 2 : no_of_animals) {
+
+ # 이전 개체와 현재 개체가 같은 경우
+ if (pedi_04[i - 1, "animal"] == pedi_04[i, "animal"]){
+
+ # 아비가 다르거나 어미가 다르다면
+ if (pedi_04[i - 1, "sire"] != pedi_04[i, "sire"]
+ || pedi_04[i - 1, "dam"] != pedi_04[i, "dam"]){
+ pedi_04[i - 1, "dup"] = 1
+ pedi_04[i , "dup"] = 1
+
+ print(paste("개체 ", pedi_04[i, "animal"] , "가 두 번 등장, 혈통이 다름"))
+ }
+ }
+ }
> # print animals with duplication error
> pedi_04[which(pedi_04$dup == 1), ]
[1] animal sire dam dup sex asd gen
<0 행> <또는 row.names의 길이가 0입니다>
> ##########################
> # check sex error
> ##########################
> for (i in 1 : no_of_animals) {
+ # sire와 dam을 저장
+ sire = pedi_04[i, "sire"]
+ dam = pedi_04[i, "dam"]
+
+ # sire가 있다면
+ if (sire != 0) {
+ # sire가 개체의 어디에 나오는가
+ index = which(pedi_04$animal == sire)
+
+ # sire의 sex가 결정되지 않았다면 수로 표시
+ if (pedi_04[index, "sex"] == 0) {
+ pedi_04[index, "sex"] = 2
+ }
+ # sire의 sex가 암으로 되어 있다면 error표시
+ else if (pedi_04[index, "sex"] == 1) {
+ pedi_04[index, "sex"] = 3
+ }
+ }
+
+ # dma이 있다면
+ if (dam != 0) {
+ # dam이 개체의 어디에 나오는가
+ index = which(pedi_04$animal == dam)
+
+ # dam의 sex가 결정되지 않았다면 암으로 표시
+ if (pedi_04[index, "sex"] == 0) {
+ pedi_04[index, "sex"] = 1
+ }
+ # dam의 sex가 수로 되어 있다면 error 표시
+ else if (pedi_04[index, "sex"] == 2) {
+ pedi_04[index, "sex"] = 3
+ }
+ }
+ }
> # print animals with sex error
> pedi_04[which(pedi_04$sex == 3), ]
[1] animal sire dam dup sex asd gen
<0 행> <또는 row.names의 길이가 0입니다>
> ########################################
> # check animal = sire = dam error
> ########################################
> # animal = sire인 개체의 위치 저장
> index = which(pedi_04$animal == pedi_04$sire)
> # animal = sire인 개체가 있다면 error 추가
> if (length(index) != 0) {
+ pedi_04[index, "asd"] = pedi_04[index, "asd"] + 1
+ }
> # animal = dam인 개체의 위치 저장
> index = which(pedi_04$animal == pedi_04$dam)
> # animal = dam인 개체가 있다면 error 추가
> if (length(index) != 0) {
+ pedi_04[index, "asd"] = pedi_04[index, "asd"] + 2
+ }
> # sire와 dam이 0이 아니고, sire = dam인 개체의 위치 저장
> index = which(pedi_04$sire == pedi_04$dam & pedi_04$sire != 0 & pedi_04$dam != 0)
> # sire = dam이면 error 추가
> if (length(index) != 0) {
+ pedi_04[index, "asd"] = pedi_04[index, "asd"] + 4
+ }
> # print animals with a = s = d error
> pedi_04[which(pedi_04$asd != 0), ]
[1] animal sire dam dup sex asd gen
<0 행> <또는 row.names의 길이가 0입니다>
> #########################
> # check cycling error
> #########################
> for (i in 1 : 50) {
+ for (j in 1 : no_of_animals) {
+ sire = pedi_04[j, "sire"]
+ dam = pedi_04[j, "dam"]
+
+ # 아비가 0이 아니고, 아비의 세대수가 개체의 세대수보다 같거나 작다면
+ index = which(pedi_04$animal == sire)
+ if ((sire != 0) && (pedi_04[index, "gen"] <= pedi_04[j, "gen"])) {
+ pedi_04[index, "gen"] = pedi_04[j, "gen"] + 1
+ }
+
+ # 어미가 0이 아니고, 어미의 세대수가 개체의 세대수보다 같거나 작다면
+ index = which(pedi_04$animal == dam)
+ if ((dam != 0) && (pedi_04[index, "gen"] <= pedi_04[j, "gen"])) {
+ pedi_04[index, "gen"] = pedi_04[j, "gen"] + 1
+ }
+ }
+ }
> # sort by generation
> pedi_04 = pedi_04[order(-pedi_04$gen),]
> # print pedigree : if generation is very large, there may be an error
> head(pedi_04, 10)
animal sire dam dup sex asd gen
1 a1 0 0 0 2 0 4
2 a2 0 0 0 1 0 4
3 a3 a1 a2 0 1 0 3
5 a4 a1 0 0 2 0 3
6 a5 a4 a3 0 2 0 2
4 a6 a5 a2 0 0 0 1
> tail(pedi_04, 10)
animal sire dam dup sex asd gen
1 a1 0 0 0 2 0 4
2 a2 0 0 0 1 0 4
3 a3 a1 a2 0 1 0 3
5 a4 a1 0 0 2 0 3
6 a5 a4 a3 0 2 0 2
4 a6 a5 a2 0 0 0 1
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".stk"))
> write.table(pedi_04, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 결과 파일(pedi_20200410.stk)
animal,sire,dam,dup,sex,asd,gen
a1,0,0,0,2,0,4
a2,0,0,0,1,0,4
a3,a1,a2,0,1,0,3
a4,a1,0,0,2,0,3
a5,a4,a3,0,2,0,2
a6,a5,a2,0,0,0,1
■ Renumber Pedigree
- R 프로그램
# renumber pedigree
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.stk", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print pedigree
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# animal 기준 넘버링
animal_renum = as.data.frame(cbind(animal = pedi_01[, "animal"], arenum = seq(1, nrow(pedi_01))), stringsAsFactors = FALSE)
# 부모가 0인 경우 추가
animal_renum = rbind(animal_renum, c(0, 0))
# print animal renumber
dim(animal_renum)
head(animal_renum)
tail(animal_renum)
# 조인하기 위한 패키지 실행
#install.packages("plyr")
library(plyr)
# 개체에 번호 붙이기
pedi_02 = join(x = pedi_01, y = animal_renum, by = "animal", type = "left")
# print pedigree
head(pedi_02)
tail(pedi_02)
# animal에서 sire로 이름 변경
sire_renum = animal_renum
colnames(sire_renum) = c("sire", "srenum")
head(sire_renum)
tail(sire_renum)
# sire에 번호 붙이기
pedi_03 = join(x = pedi_02, y = sire_renum, by = "sire", type = "left")
# print pedigree
head(pedi_03)
tail(pedi_03)
# sire에서 dam으로 이름 변경
dam_renum = animal_renum
colnames(dam_renum) = c("dam", "drenum")
head(dam_renum)
# dam에 번호 붙이기
pedi_04 = join(x = pedi_03, y = dam_renum, by = "dam", type = "left")
# print pedigree
head(pedi_04)
tail(pedi_04)
# select arenum, srenum, drenum
pedi_05 = pedi_04[, c("arenum", "srenum", "drenum")]
head(pedi_05)
# renumbered pedigree 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".ren"))
write.table(pedi_05, output_filename, sep=",", row.names = FALSE, quote = FALSE)
# animal number 파일로 출력, 분리자는 ",", 따옴표 없고
arn_animal = pedi_04[,c("arenum", "animal", "sire", "dam")]
output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".arn"))
write.table(arn_animal, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.stk", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # print pedigree
> dim(pedi_01)
[1] 6 7
> head(pedi_01)
animal sire dam dup sex asd gen
1 a1 0 0 0 2 0 4
2 a2 0 0 0 1 0 4
3 a3 a1 a2 0 1 0 3
4 a4 a1 0 0 2 0 3
5 a5 a4 a3 0 2 0 2
6 a6 a5 a2 0 0 0 1
> tail(pedi_01)
animal sire dam dup sex asd gen
1 a1 0 0 0 2 0 4
2 a2 0 0 0 1 0 4
3 a3 a1 a2 0 1 0 3
4 a4 a1 0 0 2 0 3
5 a5 a4 a3 0 2 0 2
6 a6 a5 a2 0 0 0 1
> # animal 기준 넘버링
> animal_renum = as.data.frame(cbind(animal = pedi_01[, "animal"], arenum = seq(1, nrow(pedi_01))), stringsAsFactors = FALSE)
> # 부모가 0인 경우 추가
> animal_renum = rbind(animal_renum, c(0, 0))
> # print animal renumber
> dim(animal_renum)
[1] 7 2
> head(animal_renum)
animal arenum
1 a1 1
2 a2 2
3 a3 3
4 a4 4
5 a5 5
6 a6 6
> tail(animal_renum)
animal arenum
2 a2 2
3 a3 3
4 a4 4
5 a5 5
6 a6 6
7 0 0
> # 조인하기 위한 패키지 실행
> #install.packages("plyr")
> library(plyr)
> # 개체에 번호 붙이기
> pedi_02 = join(x = pedi_01, y = animal_renum, by = "animal", type = "left")
> # print pedigree
> head(pedi_02)
animal sire dam dup sex asd gen arenum
1 a1 0 0 0 2 0 4 1
2 a2 0 0 0 1 0 4 2
3 a3 a1 a2 0 1 0 3 3
4 a4 a1 0 0 2 0 3 4
5 a5 a4 a3 0 2 0 2 5
6 a6 a5 a2 0 0 0 1 6
> tail(pedi_02)
animal sire dam dup sex asd gen arenum
1 a1 0 0 0 2 0 4 1
2 a2 0 0 0 1 0 4 2
3 a3 a1 a2 0 1 0 3 3
4 a4 a1 0 0 2 0 3 4
5 a5 a4 a3 0 2 0 2 5
6 a6 a5 a2 0 0 0 1 6
> # animal에서 sire로 이름 변경
> sire_renum = animal_renum
> colnames(sire_renum) = c("sire", "srenum")
> head(sire_renum)
sire srenum
1 a1 1
2 a2 2
3 a3 3
4 a4 4
5 a5 5
6 a6 6
> tail(sire_renum)
sire srenum
2 a2 2
3 a3 3
4 a4 4
5 a5 5
6 a6 6
7 0 0
> # sire에 번호 붙이기
> pedi_03 = join(x = pedi_02, y = sire_renum, by = "sire", type = "left")
> # print pedigree
> head(pedi_03)
animal sire dam dup sex asd gen arenum srenum
1 a1 0 0 0 2 0 4 1 0
2 a2 0 0 0 1 0 4 2 0
3 a3 a1 a2 0 1 0 3 3 1
4 a4 a1 0 0 2 0 3 4 1
5 a5 a4 a3 0 2 0 2 5 4
6 a6 a5 a2 0 0 0 1 6 5
> tail(pedi_03)
animal sire dam dup sex asd gen arenum srenum
1 a1 0 0 0 2 0 4 1 0
2 a2 0 0 0 1 0 4 2 0
3 a3 a1 a2 0 1 0 3 3 1
4 a4 a1 0 0 2 0 3 4 1
5 a5 a4 a3 0 2 0 2 5 4
6 a6 a5 a2 0 0 0 1 6 5
> # sire에서 dam으로 이름 변경
> dam_renum = animal_renum
> colnames(dam_renum) = c("dam", "drenum")
> head(dam_renum)
dam drenum
1 a1 1
2 a2 2
3 a3 3
4 a4 4
5 a5 5
6 a6 6
> # dam에 번호 붙이기
> pedi_04 = join(x = pedi_03, y = dam_renum, by = "dam", type = "left")
> # print pedigree
> head(pedi_04)
animal sire dam dup sex asd gen arenum srenum drenum
1 a1 0 0 0 2 0 4 1 0 0
2 a2 0 0 0 1 0 4 2 0 0
3 a3 a1 a2 0 1 0 3 3 1 2
4 a4 a1 0 0 2 0 3 4 1 0
5 a5 a4 a3 0 2 0 2 5 4 3
6 a6 a5 a2 0 0 0 1 6 5 2
> tail(pedi_04)
animal sire dam dup sex asd gen arenum srenum drenum
1 a1 0 0 0 2 0 4 1 0 0
2 a2 0 0 0 1 0 4 2 0 0
3 a3 a1 a2 0 1 0 3 3 1 2
4 a4 a1 0 0 2 0 3 4 1 0
5 a5 a4 a3 0 2 0 2 5 4 3
6 a6 a5 a2 0 0 0 1 6 5 2
> # select arenum, srenum, drenum
> pedi_05 = pedi_04[, c("arenum", "srenum", "drenum")]
> head(pedi_05)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> # renumbered pedigree 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".ren"))
> write.table(pedi_05, output_filename, sep=",", row.names = FALSE, quote = FALSE)
> # animal number 파일로 출력, 분리자는 ",", 따옴표 없고
> arn_animal = pedi_04[,c("arenum", "animal", "sire", "dam")]
> output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".arn"))
> write.table(arn_animal, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 결과 파일
<pedi_20200410.ren>
arenum,srenum,drenum
1,0,0
2,0,0
3,1,2
4,1,0
5,4,3
6,5,2
<pedi_20200410.arn>
arenum,animal,sire,dam
1,a1,0,0
2,a2,0,0
3,a3,a1,a2
4,a4,a1,0
5,a5,a4,a3
6,a6,a5,a2
■ Numerator Relationship Matrix 만들기
- R 프로그램
# make numerator relationship matrix
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 혈통 파일 확인
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
no_of_animals
# empty NRM 생성
A = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
# empty NRM 확인
A
# 개체를 하나씩 읽으면서 NRM을 채워 나가기
for (i in 1:no_of_animals) {
animal = pedi_01[i, 1]
sire = pedi_01[i, 2]
dam = pedi_01[i, 3]
if (sire == 0 & dam == 0) {
# 부모를 모를 경우
A[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# 부만 알고 있을 경우
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, sire])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1
} else if (sire == 0 & dam != 0) {
# 모만 알고 있을 경우
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, dam])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1
} else {
# 부와 모를 알고 있을 경우
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, sire] + A[j, dam])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1 + 0.5 * A[sire, dam]
}
}
# NRM 출력
A
# install 'MASS' packages
#install.packages("MASS")
# activate package
library("MASS")
# inverse of A
ginv(A)
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 혈통 파일 확인
> dim(pedi_01)
[1] 6 3
> head(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> tail(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> # number of animal
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 6
> # empty NRM 생성
> A = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
> # empty NRM 확인
> A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 0 0 0 0 0 0
[5,] 0 0 0 0 0 0
[6,] 0 0 0 0 0 0
> # 개체를 하나씩 읽으면서 NRM을 채워 나가기
> for (i in 1:no_of_animals) {
+ animal = pedi_01[i, 1]
+ sire = pedi_01[i, 2]
+ dam = pedi_01[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 부모를 모를 경우
+ A[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # 부만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ A[j, animal] = 0.5 * (A[j, sire])
+ A[animal, j] = A[j, animal]
+ }
+ A[animal, animal] = 1
+ } else if (sire == 0 & dam != 0) {
+ # 모만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ A[j, animal] = 0.5 * (A[j, dam])
+ A[animal, j] = A[j, animal]
+ }
+ A[animal, animal] = 1
+ } else {
+ # 부와 모를 알고 있을 경우
+ for (j in 1:animal - 1) {
+ A[j, animal] = 0.5 * (A[j, sire] + A[j, dam])
+ A[animal, j] = A[j, animal]
+ }
+ A[animal, animal] = 1 + 0.5 * A[sire, dam]
+ }
+ }
> # NRM 출력
> A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00 0.000 0.5000 0.5000 0.5000 0.2500
[2,] 0.00 1.000 0.5000 0.0000 0.2500 0.6250
[3,] 0.50 0.500 1.0000 0.2500 0.6250 0.5625
[4,] 0.50 0.000 0.2500 1.0000 0.6250 0.3125
[5,] 0.50 0.250 0.6250 0.6250 1.1250 0.6875
[6,] 0.25 0.625 0.5625 0.3125 0.6875 1.1250
> # activate package
> library("MASS")
> # inverse of A
> ginv(A)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.833333e+00 5.000000e-01 -1.000000e+00 -6.666667e-01 -1.887379e-15 6.106227e-16
[2,] 5.000000e-01 2.033333e+00 -1.000000e+00 2.220446e-16 5.333333e-01 -1.066667e+00
[3,] -1.000000e+00 -1.000000e+00 2.500000e+00 5.000000e-01 -1.000000e+00 -8.881784e-16
[4,] -6.666667e-01 -4.440892e-16 5.000000e-01 1.833333e+00 -1.000000e+00 -7.771561e-16
[5,] -4.440892e-16 5.333333e-01 -1.000000e+00 -1.000000e+00 2.533333e+00 -1.066667e+00
[6,] -4.996004e-16 -1.066667e+00 -1.110223e-16 -1.110223e-16 -1.066667e+00 2.133333e+00
■ 근친교배를 무시한 NRM의 역행렬 만들기
- R 프로그램
# make inverse of numerator relationship matrix ignoring inbreeding
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 혈통 파일 확인
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
no_of_animals
# 혈연계수의 역행렬(A-1) 할당
Ainv = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
Ainv
# 혈연계수의 역행렬(A-1) 계산
for (i in 1 : no_of_animals) {
animal = pedi_01[i, 1]
sire = pedi_01[i, 2]
dam = pedi_01[i, 3]
if (sire == 0 & dam == 0) {
# 개체의 부모를 둘다 모를 경우
alpha = 1
Ainv[animal, animal] = alpha + Ainv[animal, animal]
} else if (sire != 0 & dam == 0) {
# 개체의 부만 알고 있을 경우
alpha = 4/3
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
} else if (sire == 0 & dam != 0) {
# 개체의 모만 알고 있을 경우
alpha = 4/3
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
} else {
# 개체의 부모를 모두 알고 있을 경우
alpha = 2
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
}
# 개체별 처리과정 출력 cat('Animal', i, '\n') print(Ainv)
}
# print A - Inverse
Ainv
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 혈통 파일 확인
> dim(pedi_01)
[1] 6 3
> head(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> tail(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> # number of animal
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 6
> # 혈연계수의 역행렬(A-1) 할당
> Ainv = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 0 0 0 0 0 0
[5,] 0 0 0 0 0 0
[6,] 0 0 0 0 0 0
> for (i in 1 : no_of_animals) {
+ animal = pedi_01[i, 1]
+ sire = pedi_01[i, 2]
+ dam = pedi_01[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 개체의 부모를 둘다 모를 경우
+ alpha = 1
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ } else if (sire != 0 & dam == 0) {
+ # 개체의 부만 알고 있을 경우
+ alpha = 4/3
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+ Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+ Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+ } else if (sire == 0 & dam != 0) {
+ # 개체의 모만 알고 있을 경우
+ alpha = 4/3
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+ Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+ Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+ } else {
+ # 개체의 부모를 모두 알고 있을 경우
+ alpha = 2
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+ Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+ Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+ Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+ Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+ Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+ Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+ Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+ }
+
+ # 개체별 처리과정 출력 cat('Animal', i, '\n') print(Ainv)
+ }
> # print A - Inverse
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.8333333 0.5 -1.0 -0.6666667 0.0 0
[2,] 0.5000000 2.0 -1.0 0.0000000 0.5 -1
[3,] -1.0000000 -1.0 2.5 0.5000000 -1.0 0
[4,] -0.6666667 0.0 0.5 1.8333333 -1.0 0
[5,] 0.0000000 0.5 -1.0 -1.0000000 2.5 -1
[6,] 0.0000000 -1.0 0.0 0.0000000 -1.0 2
■ Inbreeding Coefficient 계산하기
- R 프로그램
# calculate inbreeding coefficients using Meuwissen and Luo algorithm
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 혈통 파일 확인
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
no_of_animals
# 근교계수를 저장할 컬럼 추가
pedi_02 = cbind(pedi_01, inb = rep(0, no_of_animals))
dim(pedi_02)
head(pedi_02)
tail(pedi_02)
# Dii를 구하는 함수
getDii = function(animal) {
sire = pedi_02[animal, 2]
dam = pedi_02[animal, 3]
if (sire == 0 & dam == 0) {
return(1)
} else if (sire != 0 & dam == 0) {
return(0.75 - 0.25 * pedi_02[sire, 4])
} else if (sire == 0 & dam != 0) {
return(0.75 - 0.25 * pedi_02[dam, 4])
} else if (sire != 0 & dam != 0) {
return(0.5 - 0.25 * (pedi_02[sire, 4] + pedi_02[dam, 4]))
}
}
# 1열 즉 개체에 대하여 정렬
pedi_02 = pedi_02[order(pedi_02[, 1]), ]
head(pedi_02)
tail(pedi_02)
# 한 개체씩 근교계수 계산
for (i in 1 : no_of_animals) {
# id - t array
idt = data.frame(id = pedi_02[i, 1], t = 1)
# current row position of array
cur_pos = 1
# add position of sire or dam
add_pos = 2
# idt를 모두 처리할 때까지
while (cur_pos < add_pos) {
animal = idt[cur_pos, 1]
sire = pedi_02[animal, 2]
dam = pedi_02[animal, 3]
if (sire > 0) {
# sire가 있다면
idt = rbind(idt, c(sire, idt[cur_pos, 2] * 0.5))
add_pos = add_pos + 1
}
if (dam > 0) {
# dam이 있다면
idt = rbind(idt, c(dam, idt[cur_pos, 2] * 0.5))
add_pos = add_pos + 1
}
cur_pos = cur_pos + 1
}
# sorting
idt = idt[order(idt$id), ]
# 개체의 근교계수 계산
p_id = idt[1, 1]
t = idt[1, 2]
f = -1
if (nrow(idt) > 1) {
for (j in 2:nrow(idt)) {
if (p_id != idt[j, 1]) {
f = f + t * t * getDii(p_id)
p_id = idt[j, 1]
t = idt[j, 2]
} else {
t = t + idt[j, 2]
}
}
}
f = f + t * t * getDii(p_id)
pedi_02[i, 4] = f
}
# print inbreeding coefficient
pedi_02
# 혈통 자료 읽기
ori_pedi <- read.table("pedi_20200409.arn", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 파일 확인
dim(ori_pedi)
head(ori_pedi)
tail(ori_pedi)
# 조인하기 위한 패키지 실행
#install.packages("plyr")
library(plyr)
# add orignal pedigree
pedi_03 = join(x = pedi_02, y = ori_pedi, by = "arenum", type = "left")
# renumbered pedigree, inbreeding coefficient, original pedigree, sep = ",", no quote
output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".inb"))
write.table(pedi_03, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.ren", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 혈통 파일 확인
> dim(pedi_01)
[1] 6 3
> head(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> tail(pedi_01)
arenum srenum drenum
1 1 0 0
2 2 0 0
3 3 1 2
4 4 1 0
5 5 4 3
6 6 5 2
> # number of animal
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 6
> # 근교계수를 저장할 컬럼 추가
> pedi_02 = cbind(pedi_01, inb = rep(0, no_of_animals))
> dim(pedi_02)
[1] 6 4
> head(pedi_02)
arenum srenum drenum inb
1 1 0 0 0
2 2 0 0 0
3 3 1 2 0
4 4 1 0 0
5 5 4 3 0
6 6 5 2 0
> tail(pedi_02)
arenum srenum drenum inb
1 1 0 0 0
2 2 0 0 0
3 3 1 2 0
4 4 1 0 0
5 5 4 3 0
6 6 5 2 0
> # Dii를 구하는 함수
> getDii = function(animal) {
+ sire = pedi_02[animal, 2]
+ dam = pedi_02[animal, 3]
+
+ if (sire == 0 & dam == 0) {
+ return(1)
+ } else if (sire != 0 & dam == 0) {
+ return(0.75 - 0.25 * pedi_02[sire, 4])
+ } else if (sire == 0 & dam != 0) {
+ return(0.75 - 0.25 * pedi_02[dam, 4])
+ } else if (sire != 0 & dam != 0) {
+ return(0.5 - 0.25 * (pedi_02[sire, 4] + pedi_02[dam, 4]))
+ }
+ }
> # 1열 즉 개체에 대하여 정렬
> pedi_02 = pedi_02[order(pedi_02[, 1]), ]
> head(pedi_02)
arenum srenum drenum inb
1 1 0 0 0
2 2 0 0 0
3 3 1 2 0
4 4 1 0 0
5 5 4 3 0
6 6 5 2 0
> tail(pedi_02)
arenum srenum drenum inb
1 1 0 0 0
2 2 0 0 0
3 3 1 2 0
4 4 1 0 0
5 5 4 3 0
6 6 5 2 0
> # 한 개체씩 근교계수 계산
> for (i in 1 : no_of_animals) {
+ # id - t array
+ idt = data.frame(id = pedi_02[i, 1], t = 1)
+
+ # current row position of array
+ cur_pos = 1
+ # add position of sire or dam
+ add_pos = 2
+
+ # idt를 모두 처리할 때까지
+ while (cur_pos < add_pos) {
+ animal = idt[cur_pos, 1]
+ sire = pedi_02[animal, 2]
+ dam = pedi_02[animal, 3]
+
+ if (sire > 0) {
+ # sire가 있다면
+ idt = rbind(idt, c(sire, idt[cur_pos, 2] * 0.5))
+ add_pos = add_pos + 1
+ }
+
+ if (dam > 0) {
+ # dam이 있다면
+ idt = rbind(idt, c(dam, idt[cur_pos, 2] * 0.5))
+ add_pos = add_pos + 1
+ }
+ cur_pos = cur_pos + 1
+ }
+
+ # sorting
+ idt = idt[order(idt$id), ]
+
+ # 개체의 근교계수 계산
+ p_id = idt[1, 1]
+ t = idt[1, 2]
+ f = -1
+
+ if (nrow(idt) > 1) {
+ for (j in 2:nrow(idt)) {
+ if (p_id != idt[j, 1]) {
+ f = f + t * t * getDii(p_id)
+ p_id = idt[j, 1]
+ t = idt[j, 2]
+ } else {
+ t = t + idt[j, 2]
+ }
+ }
+ }
+ f = f + t * t * getDii(p_id)
+ pedi_02[i, 4] = f
+ }
> # print inbreeding coefficient
> pedi_02
arenum srenum drenum inb
1 1 0 0 0.000
2 2 0 0 0.000
3 3 1 2 0.000
4 4 1 0 0.000
5 5 4 3 0.125
6 6 5 2 0.125
> # 혈통 자료 읽기
> ori_pedi <- read.table("pedi_20200409.arn", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 파일 확인
> dim(ori_pedi)
[1] 6 4
> head(ori_pedi)
arenum animal sire dam
1 1 a1 0 0
2 2 a2 0 0
3 3 a3 a1 a2
4 4 a4 a1 0
5 5 a5 a4 a3
6 6 a6 a5 a2
> tail(ori_pedi)
arenum animal sire dam
1 1 a1 0 0
2 2 a2 0 0
3 3 a3 a1 a2
4 4 a4 a1 0
5 5 a5 a4 a3
6 6 a6 a5 a2
> # 조인하기 위한 패키지 실행
> #install.packages("plyr")
> library(plyr)
> # add orignal pedigree
> pedi_03 = join(x = pedi_02, y = ori_pedi, by = "arenum", type = "left")
> # renumbered pedigree, inbreeding coefficient, original pedigree, sep = ",", no quote
> output_filename = gsub("[ -]", "", paste("pedi_", Sys.Date(), ".inb"))
> write.table(pedi_03, output_filename, sep=",", row.names = FALSE, quote = FALSE)
- 결과 파일(pedi_20200410.inb)
arenum,srenum,drenum,inb,animal,sire,dam
1,0,0,0,a1,0,0
2,0,0,0,a2,0,0
3,1,2,0,a3,a1,a2
4,1,0,0,a4,a1,0
5,4,3,0.125,a5,a4,a3
6,5,2,0.125,a6,a5,a2
■ TDT’ 행렬 만들기
- R 프로그램
# make T, D matrix and then make A
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.inb", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 혈통 파일 확인
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
no_of_animals
# empty T 생성
T = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
# empty T 확인
T
# 개체를 하나씩 읽으면서 NRM을 채워 나가기
for (i in 1:no_of_animals) {
animal = pedi_01[i, 1]
sire = pedi_01[i, 2]
dam = pedi_01[i, 3]
if (sire == 0 & dam == 0) {
# 부모를 모를 경우
T[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# 부만 알고 있을 경우
for (j in 1:animal - 1) {
T[animal, j] = 0.5 * (T[sire, j])
}
T[animal, animal] = 1
} else if (sire == 0 & dam != 0) {
# 모만 알고 있을 경우
for (j in 1:animal - 1) {
T[animal, j] = 0.5 * (T[dam, j])
}
T[animal, animal] = 1
} else {
# 부와 모를 알고 있을 경우
for (j in 1:animal - 1) {
T[animal, j] = 0.5 * (T[sire, j] + T[dam, j])
}
T[animal, animal] = 1
}
}
# T 출력
T
# Assum that we know inbreeding coefficients
InbCo = pedi_01[, c("inb")]
InbCo
# empty D
D = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
D
for (i in 1:6) {
animal = pedi_01[i, 1]
sire = pedi_01[i, 2]
dam = pedi_01[i, 3]
if (sire == 0 & dam == 0) {
# 부모를 모를 경우
D[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# 부만 알고 있을 경우
D[animal, animal] = 0.75 - 0.25 * InbCo[sire]
} else if (sire == 0 & dam != 0) {
# 모만 알고 있을 경우
D[animal, animal] = 0.75 - 0.25 * InbCo[dam]
} else {
# 부와 모를 알고 있을 경우
D[animal, animal] = 0.5 - 0.25 * InbCo[sire] - 0.25 * InbCo[dam]
}
}
# print D
D
# print A
T %*% D %*% t(T)
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.inb", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 혈통 파일 확인
> dim(pedi_01)
[1] 6 7
> head(pedi_01)
arenum srenum drenum inb animal sire dam
1 1 0 0 0.000 a1 0 0
2 2 0 0 0.000 a2 0 0
3 3 1 2 0.000 a3 a1 a2
4 4 1 0 0.000 a4 a1 0
5 5 4 3 0.125 a5 a4 a3
6 6 5 2 0.125 a6 a5 a2
> tail(pedi_01)
arenum srenum drenum inb animal sire dam
1 1 0 0 0.000 a1 0 0
2 2 0 0 0.000 a2 0 0
3 3 1 2 0.000 a3 a1 a2
4 4 1 0 0.000 a4 a1 0
5 5 4 3 0.125 a5 a4 a3
6 6 5 2 0.125 a6 a5 a2
> # number of animal
> no_of_animals = nrow(pedi_01)
> no_of_animals
[1] 6
> # empty T 생성
> T = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
> # empty T 확인
> T
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 0 0 0 0 0 0
[5,] 0 0 0 0 0 0
[6,] 0 0 0 0 0 0
> # 개체를 하나씩 읽으면서 NRM을 채워 나가기
> for (i in 1:no_of_animals) {
+ animal = pedi_01[i, 1]
+ sire = pedi_01[i, 2]
+ dam = pedi_01[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 부모를 모를 경우
+ T[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # 부만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ T[animal, j] = 0.5 * (T[sire, j])
+ }
+ T[animal, animal] = 1
+ } else if (sire == 0 & dam != 0) {
+ # 모만 알고 있을 경우
+ for (j in 1:animal - 1) {
+ T[animal, j] = 0.5 * (T[dam, j])
+ }
+ T[animal, animal] = 1
+ } else {
+ # 부와 모를 알고 있을 경우
+ for (j in 1:animal - 1) {
+ T[animal, j] = 0.5 * (T[sire, j] + T[dam, j])
+ }
+ T[animal, animal] = 1
+ }
+ }
> # T 출력
> T
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00 0.000 0.00 0.00 0.0 0
[2,] 0.00 1.000 0.00 0.00 0.0 0
[3,] 0.50 0.500 1.00 0.00 0.0 0
[4,] 0.50 0.000 0.00 1.00 0.0 0
[5,] 0.50 0.250 0.50 0.50 1.0 0
[6,] 0.25 0.625 0.25 0.25 0.5 1
> # Assum that we know inbreeding coefficients
> InbCo = pedi_01[, c("inb")]
> InbCo
[1] 0.000 0.000 0.000 0.000 0.125 0.125
> # empty D
> D = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
> D
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 0 0 0 0 0 0
[5,] 0 0 0 0 0 0
[6,] 0 0 0 0 0 0
> for (i in 1:6) {
+ animal = pedi_01[i, 1]
+ sire = pedi_01[i, 2]
+ dam = pedi_01[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 부모를 모를 경우
+ D[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # 부만 알고 있을 경우
+ D[animal, animal] = 0.75 - 0.25 * InbCo[sire]
+ } else if (sire == 0 & dam != 0) {
+ # 모만 알고 있을 경우
+ D[animal, animal] = 0.75 - 0.25 * InbCo[dam]
+ } else {
+ # 부와 모를 알고 있을 경우
+ D[animal, animal] = 0.5 - 0.25 * InbCo[sire] - 0.25 * InbCo[dam]
+ }
+ }
> # print D
> D
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0.0 0.00 0.0 0.00000
[2,] 0 1 0.0 0.00 0.0 0.00000
[3,] 0 0 0.5 0.00 0.0 0.00000
[4,] 0 0 0.0 0.75 0.0 0.00000
[5,] 0 0 0.0 0.00 0.5 0.00000
[6,] 0 0 0.0 0.00 0.0 0.46875
> # print A
> T %*% D %*% t(T)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00 0.000 0.5000 0.5000 0.5000 0.2500
[2,] 0.00 1.000 0.5000 0.0000 0.2500 0.6250
[3,] 0.50 0.500 1.0000 0.2500 0.6250 0.5625
[4,] 0.50 0.000 0.2500 1.0000 0.6250 0.3125
[5,] 0.50 0.250 0.6250 0.6250 1.1250 0.6875
[6,] 0.25 0.625 0.5625 0.3125 0.6875 1.1250
>
■ 근친교배를 고려한 NRM의 역행렬 만들기
- R 프로그램
# make inverse of the numerator relationship matrix accounting for inbreedig
# set working directory
setwd("D:/temp")
# print working directory
getwd()
# 혈통 자료 읽기
pedi_01 <- read.table("pedi_20200410.inb", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# 입력한 혈통 파일 확인
dim(pedi_01)
head(pedi_01)
tail(pedi_01)
# number of animal
no_of_animals = nrow(pedi_01)
# 혈연계수의 역행렬(A-1) 할당
Ainv = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
# inbreeding coefficients
InbCo = pedi_01[, c("inb")]
InbCo
# 혈연계수의 역행렬(A-1) 계산
for (i in 1:no_of_animals) {
animal = pedi_01[i, 1]
sire = pedi_01[i, 2]
dam = pedi_01[i, 3]
if (sire == 0 & dam == 0) {
# 개체의 부모를 둘다 모를 경우
alpha = 1
Ainv[animal, animal] = alpha + Ainv[animal, animal]
} else if (sire != 0 & dam == 0) {
# 개체의 부만 알고 있을 경우
alpha = 1/(0.75 - 0.25 * InbCo[sire])
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
} else if (sire == 0 & dam != 0) {
# 개체의 모만 알고 있을 경우
alpha = 1/(0.75 - 0.25 * InbCo[dam])
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
} else {
# 개체의 부모를 모두 알고 있을 경우
alpha = 1/(0.5 - 0.25 * InbCo[sire] - 0.25 * InbCo[dam])
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
}
# 개체별 처리과정 출력 cat('Animal', i, '\n') print(Ainv)
}
# print A Inverse
Ainv
- 실행 결과
> # set working directory
> setwd("D:/temp")
> # print working directory
> getwd()
[1] "D:/temp"
> # 혈통 자료 읽기
> pedi_01 <- read.table("pedi_20200410.inb", header = TRUE, sep = ",", stringsAsFactors = FALSE)
> # 입력한 혈통 파일 확인
> dim(pedi_01)
[1] 6 7
> head(pedi_01)
arenum srenum drenum inb animal sire dam
1 1 0 0 0.000 a1 0 0
2 2 0 0 0.000 a2 0 0
3 3 1 2 0.000 a3 a1 a2
4 4 1 0 0.000 a4 a1 0
5 5 4 3 0.125 a5 a4 a3
6 6 5 2 0.125 a6 a5 a2
> tail(pedi_01)
arenum srenum drenum inb animal sire dam
1 1 0 0 0.000 a1 0 0
2 2 0 0 0.000 a2 0 0
3 3 1 2 0.000 a3 a1 a2
4 4 1 0 0.000 a4 a1 0
5 5 4 3 0.125 a5 a4 a3
6 6 5 2 0.125 a6 a5 a2
> # number of animal
> no_of_animals = nrow(pedi_01)
> # 혈연계수의 역행렬(A-1) 할당
> Ainv = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
> # inbreeding coefficients
> InbCo = pedi_01[, c("inb")]
> InbCo
[1] 0.000 0.000 0.000 0.000 0.125 0.125
> for (i in 1:no_of_animals) {
+ animal = pedi_01[i, 1]
+ sire = pedi_01[i, 2]
+ dam = pedi_01[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # 개체의 부모를 둘다 모를 경우
+ alpha = 1
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ } else if (sire != 0 & dam == 0) {
+ # 개체의 부만 알고 있을 경우
+ alpha = 1/(0.75 - 0.25 * InbCo[sire])
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+ Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+ Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+ } else if (sire == 0 & dam != 0) {
+ # 개체의 모만 알고 있을 경우
+ alpha = 1/(0.75 - 0.25 * InbCo[dam])
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+ Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+ Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+ } else {
+ # 개체의 부모를 모두 알고 있을 경우
+ alpha = 1/(0.5 - 0.25 * InbCo[sire] - 0.25 * InbCo[dam])
+ Ainv[animal, animal] = alpha + Ainv[animal, animal]
+ Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+ Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+ Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+ Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+ Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+ Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+ Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+ Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+ }
+
+ # 개체별 처리과정 출력 cat('Animal', i, '\n') print(Ainv)
+ }
> # print A Inverse
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.8333333 0.5000000 -1.0 -0.6666667 0.0000000 0.000000
[2,] 0.5000000 2.0333333 -1.0 0.0000000 0.5333333 -1.066667
[3,] -1.0000000 -1.0000000 2.5 0.5000000 -1.0000000 0.000000
[4,] -0.6666667 0.0000000 0.5 1.8333333 -1.0000000 0.000000
[5,] 0.0000000 0.5333333 -1.0 -1.0000000 2.5333333 -1.066667
[6,] 0.0000000 -1.0666667 0.0 0.0000000 -1.0666667 2.133333
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| R을 이용한 반복 모형(Repeatability Model)의 육종가 구하기 프로그래밍 (0) | 2020.06.14 |
|---|---|
| R을 이용하여 반복 모형(repeatability model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.06.13 |
| R을 이용하여 단형질 개체 모형(single trait animal model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.06.03 |
| R을 이용한 단형질 개체 모형(Single Trait Animal Model)의 육종가 구하기 프로그래밍 (0) | 2020.05.18 |
| R로 단형질 임의회귀모형의 해 구하기 (1) | 2011.06.09 |