SNP Effect를 이용한 DGV(or GEBV) 추정

single step GBLUP의 장점은 SNP genotype이 있는 개체와 없는 개체를 동시에 분석할 수 있다는 것이다. 그러나 유전평가를 한 직후에 몇 마리의 genotype이 확보되고 그 개체의 DGV를 계산할 필요성이 생겼다면 다시 전체 유전평가를 하기란 쉬운 일이 아니다. 그래서 ssGBLUP 후에 PostGSF90을 이용하여 SNP effect를 추정하였고, 여기서는 이들 SNP effect를 이용하여 새로이 genotype을 확보한 개체의 DGV를 추정할 것이다. 여기서는 PREDF90 프로그램을 이용한다.

 

준비할 것은 PostGSF90이 출력한 snp_pred 파일과 새로운 개체의 genotype 파일이다. 그런데 새로운 개체의 genotype 파일을 준비하는데 몇 가지 과정을 거쳐야 한다. 새로운 개체가 몇 마리이든 기존 개체와 함께 imputation을 해야 한다. imputation 과정은 이미 https://bhpark.tistory.com/198?category=866943 으로 설명한 바 있다. 여기서는 함께 imputation한 개체 중 DGV를 추정할 개체만 뽑아내야 한다.

 

imputation된 결과 파일을 50k_for_gs.txt 이라 하자.

 

 

추출할 개체의 개체 번호 파일을 작성하고 cow_list.txt 이라하자.

 

 

전체 파일에서 일부 개체를 추출하여 저장하는 프로그램은 다음과 같다. sort, join은 유닉스 명령어이다. 자세한 사항은 https://bhpark.tistory.com/58 을 참고한다.

 

sort -k 1,1 50k_for_gs.txt > 50k_for_gs_02.txt sort -k 1,1 cow_list.txt > cow_list_02.txt join 50k_for_gs_02.txt cow_list_02.txt > selected_cows_for_predf90.txt pause

 

 

추출한 개체의 유전자형이다.

 

 

여기서 주의해야 할 점은 유전평가할 때와 동일한 imputation을 해야 한다는 점이다. 예를 들어 유전평가할 때는 BovineSNP50v3로 했는데 DGV를 추정할 때는 HanwooSNP50v1으로 imputation하면 안 된다는 점이다. 그리고 imputation에 동원된 개체들의 genotype에 따라서 최종 imputation된 SNP가 달라질 수 있다. 예를 들어 유전평가할 때는 특정 SNP가 모두 missing(no call)이어서 imputation 결과에 특정 SNP가 빠졌는데, 후에 DGV를 추정하기 위하여 한 imputation에서는 해당 SNP의 genotype이 있어 imputation 에 포함될 수 있다. 그래서 유전평가할 때 imputation 된 SNP 수와 DGV 추정할 때 imputation 된 SNP 수를 비교하여야 하고, 다르다면 원인을 파악해야 한다.

 

다음으로 고려해야 할 것이 imputation 된 SNP의 수와 QC 후에 남아 있는 (그리고 effect를 추정한) SNP의 개수가 다르다는 점이다. 그래서 imputation된 SNP 중에서 effect를 추정한 SNP만 추출해야 한다. SNP를 추출하는 프로그램을 만들어야 하는데, 매번 effect를 가지고 있는 SNP가 달라지므로, SNP를 추출하는 프로그램을 만드는 프로그램을 작성해야 한다. 그건 SNP Map 파일에서 시작한다.

 

QC를 하고 난 SNP Map 파일의 이름은 snp_map_for_gs.txt_clean이다.

 

 

- 1열은 일련번호, 2열은 염색체 번호, 3열을 위치, 4열을 이름, 5열은 원래 일련번호

- 5열의 원래 일련번호를 이용해서 SNP 추출

 

SNP를 추출하는 프로그램을 작성하는 프로그램을 작성하는데 파일 이름을make_select_snp_awk_program.awk 로 한다.

BEGIN { printf "{print $1, " } { printf "substr($2,%d,1) ", $5 } END { printf "}\n" }

 

 

 

## 2023.3.20.

QC한 후 출력하는 map file 형식이 바뀌었다. 

 

header가 들어갔고, 1열이 사라졌다. 그래서 다음과 같이 AWK 프로그램이 바뀌어야 한다.

BEGIN {
    printf "{print $1, "
}
{
    if (NR >= 2){
        printf "substr($2,%d,1) ", $4
    }
}
END {
    printf "}\n"
}

NR >= 2라는 문장을 통해서 2행부터 뭔가 실행하는 것으로 한다. 그리고 $4 4열을 printf 문의 %d에 넣는다.

 

결과적으로 생긴 SNP를 추출하는 생성된 AWK 프로그램은 다음과 같다.

{print $1, substr($2,1,1) substr($2,3,1) substr($2,5,1) substr($2,8,1) 
substr($2,9,1) substr($2,11,1) substr($2,12,1) substr($2,14,1) 
substr($2,15,1) substr($2,17,1) substr($2,18,1) substr($2,20,1) 
substr($2,21,1) substr($2,22,1) substr($2,23,1) substr($2,24,1) 
substr($2,25,1) substr($2,26,1) substr($2,27,1) substr($2,29,1) 
...
}

 

2023.3.20. ##

실행은 다음과 같이 한다.

awk -f make_select_snp_awk_program.awk snp_map_for_gs.txt_clean > select_snp_for_predf90.awk

 

실행화면이다.

 

 

실행결과 만들어진 select_snp.awk 프로그램은 다음과 같다.

{print $1, substr($2,1,1) substr($2,2,1) substr($2,4,1) substr($2,6,1) substr($2,7,1) substr($2,11,1) substr($2,12,1) substr($2,13,1) substr($2,15,1) substr($2,16,1) substr($2,18,1) substr($2,20,1) substr($2,21,1) substr($2,22,1) substr($2,23,1) substr($2,24,1) substr($2,25,1) ... }

1열은 그대로 출력하고, 2열의 1번째 문자, 2열의 2번째 문자, 2열의 4번째 문자 등을 출력한다.

 

SNP를 추출하는 프로그램을 다음과 같이 실행한다.

awk -f select_snp_for_predf90.awk selected_cows_for_predf90.txt > selected_cows_snps_for_predf90.txt

 

 

SNP를 추출한 결과는 다음과 같다.

 

 

한 행의 길이가 42577 = 개체번호(15) + 공백(1) + SNP(46561) 임을 알 수 있다.

 

입력 파일 준비가 완료되었으므로 DGV를 계산하자.

 

준비한 파일은 다음과 같다.

- selected_cows_snps_for_predf90.txt

- snp_pred : postgsf90을 실행한 결과로 나온 파일

 

 

임의 효과(형질과 관련 효과의 수)에 대한 정보, 42561개 SNP의 gene frequency, solutions of SNP effects가 포함되어 있다. 여기서는 1개 형질만 유전평가 했으므로 1개 형질에 대한 SNP 효과가 있다.

 

프로그램 실행은 다음과 같다.

predf90 --snpfile selected_cows_snps_for_predf90.txt | tee predf90_01.log

 

 

실행 결과로 만들어진 파일은 snp_predictions이다.

 

 

1열은 개체번호, 2열은 call rate(여기서는 imputation 후이므로 모두 1.0), 형질의 DGV이다.

 

같은 개체를 ssGBLUP으로 GEBV를 구하고, PREDF90으로 DGV를 구했을 경우 값이 차이가 나는데, 일정한 값이 차이가 나므로 사용에 문제는 없을 것으로 보인다. 순위 상관은 98 ~ 99% 정도 되는 것으로 나온다. 일정한 차이 이외에 나타나는 약간의 차이는 GEBV를 구할 때 GRM과 NRM을 섞지만, PREDF90에서는 오로지 genotype만 이용하기 때문인 것으로 보인다. 이러한 차이가 있으나 한 농장의 어린 가축들의 DGV를 구하고 사용하는데는 문제가 없을 것으로 보인다.

 

EBV(estimated breeding value), GEBV(genomic estimated breeding value), DGV(direct genomic value)와 같은 용어들이 사용되고 있는데 EBV는 혈통(pedigree)과 검정자료(phenotype)를 이용한 추정치, GEBV는 혈통, 검정자료와 유전자자료(SNP genotype)를 이용한 추정치, DGV는 SNP effect를 이용하여 SNP genotype으로부터 구한 추정치를 의미한다. 물론 사용자들마다 이것을 혼용해 사용하기도 한다.