SAS에서 외부 파일을 읽고, 외부 파일을 만들어야 할 때 의도하지 않은 폴더에서 자료를 읽으려고 하고, 쓸려고 해서 곤란한 적이 있다.

예를 들어 다음과 같은 자료를 다음과 같은 SAS 프로그램으로 읽으려 한다고 해보자. 현재 SAS 프로그램과 자료가 있는 폴더는 C:\Temp\SAS이다.

 

 

 

SAS 프로그램

data a1 ;

    infile 'hv1_192_finalReport.txt' delimiter='09'x MISSOVER DSD lrecl=32767 firstobs=11;
 
    informat snpname      $37. ;
    informat sampleid     $15. ;
    informat top1_hv1     $1. ;
    informat top2_hv1     $1. ;
    informat gcscore      $6. ;
    informat samplename   $1. ;
    informat forward1_hv1 $1. ;
    informat forward2_hv1 $1. ;
    informat ab1_hv1      $1. ;
    informat ab2_hv1      $1. ;
    informat chr          $2. ;
    informat position     best32. ;
    informat snp          $5. ;
    informat x            $5. ;
    informat y            $5. ;
 
    input
        snpname      $
        sampleid     $
        top1_hv1     $
        top2_hv1     $
        gcscore      $
        samplename   $
        forward1_hv1 $
        forward2_hv1 $
        ab1_hv1      $
        ab2_hv1      $
        chr          $
        position
        snp          $
        x            $
        y            $
    ;
run;

 

위 SAS 프로그램을 실행시키는 방법은 두 가지가 있다. 먼저 SAS 프로그램을 실행하고, 파일 -> 프로그램 열기 -> SAS 프로그램 선택하는 방법이 있다.

 

 

프로그램을 선택한 후 화면

 

 

그리고 프로그램을 실행한 결과

 

 

로그 윈도우의 메시지를 보면 그런 파일이 없다고 나온다. 메시지를 잘 보면 엉뚱한 폴더에서 데이터 파일을 찾고 있다. 이때 데이터 파일의 절대 경로를 적어서 자료를 읽을 수도 있다.

 

    infile 'D:\temp\SAS\hv1_192_finalReport.txt' delimiter='09'x MISSOVER DSD lrecl=32767 firstobs=11;

 

그리고 실행한 결과

 

 

실행이 잘 되었다. 이렇게 해결해도 되지만 작업 폴더를 지정할 수도 있다.

 

SAS 프로그램의 맨 아래쪽에 폴더가 있는 것을 볼 수 있다.

 

 

클릭하여 데이터 파일이 있는 폴더를 지정한다.

 

 

데이터 파일이 있는 폴더를 지정한 화면

 

 

아래쪽 폴더가 바뀐 것을 볼 수 있다. 위의 프로그램을 보면 데이터 파일 이름만 있고 경로는 없다. 실행한 결과.

 

 

실행이 잘 되었다. 이렇게 작업 폴더를 지정하면 파일 이름을 쓸 때, 절대 경로를 지정하지 않아도 된다. 그런데 이렇게 작업 폴더를 지정하지 않아도 된다.

SAS 프로그램을 실행시킬 때 윈도우 탐색기에서 실행할 SAS 프로그램을 마우스 오른쪽 버튼으로 클릭한다.

 

 

SAS 9.4(으)로 열기를 선택한다.

 

 

자동으로 현재 폴더가 작업 폴더로 지정된 것을 볼 수 있다. 이 상태에서 프로그램을 실행하면 데이터 파일을 잘 읽는다.

그런데 본 포스트를 작성하는 목적은 이렇게 클릭해서 작업 폴더를 지정하는 것이 아니라 프로그래밍적으로 작업 폴더를 지정하는 방법을 살펴보기 위해서이다. 윈도우 SAS에서는 클릭해서 작업 폴더를 지정하든, 프로그래밍적으로 작업 폴더를 지정하든 상관이 없으나, SAS에서 제공하는 인터넷 브라우저에서 사용하는 SAS Studio에서는 사정이 달라진다. SAS Studio는 좀 뒤에 설명한다.

다음은 SAS 프로그램을 실행한 화면

 

 

그리고 아래 프로그램을 입력하고 실행한다.

 

data _null_;
    wd = dlgcdir("D:\temp\SAS");
run;

 

다음은 실행한 화면

 

 

화면의 아래쪽 폴더를 보면 작업 폴더가 바뀐 것을 볼 수 있다.

이제 SAS Studio에 대해서 알아보자. SAS OnDemand for Academics 라는 것이 있다. 검색하면 많이 나온다. 간단히 설명하자면 비싼 SAS를 쓰기에 어려운(?) 학생들을 위하여 인터넷 브라우저에서 무료로 SAS를 사용할 수 있게 해 준거고, 그때의 SAS 프로그램 인터페이스를 SAS Studio라고 부르는 거 같다.

SAS OnDemand for Academics 접속 화면.

 

 

로그인하고 Access Now 클릭

 

 

다시 Sign In

 

 

Sign In 클릭한 화면

 

 

Launch 클릭

 

 

여기서는 파일이라는 폴더에 Microarray라는 폴더를 만들고 위에서 사용한 SAS 프로그램과 데이터 파일을 업로드한 상태이다. 그리고 SAS 프로그램을 연 상태이다. 여기서 사람이 뛰는 아이콘을 클릭하여 프로그램을 실행한다. 실행한 후 로그 화면이다.

 

 

역시나 작업 폴더가 지정이 되어 있지 않아 엉뚱한 폴더에서 데이터 파일을 찾고 있어 나타나는 오류이다. 현재 데이터 파일과 SAS 프로그램이 있는 폴더는 다음과 같다. Microarray 폴더를 오른쪽 버튼으로 클릭하고 속성을 클릭한다.

 

 

위치의 경로를 선택하여 복사한다.

다음을 프로그램의 앞에 복사한다.

 

/* 다음 디렉토리를 워킹 디렉토리(자료를 import하거나 export하는 디렉토리)로 설정 */
data _null_;
    wd = dlgcdir("/home/u63604154/Microarray");
run;

 

 

 

실행 화면

 

 

이런 방식으로 SAS Studio에서 손쉽게 자료를 읽거나 쓸 수 있다. 일부에서는 업로드한 데이터를 SAS dataset으로 변환 후, Libname을 적용하고 데이터 분석을 하라는 설명도 있으나 위의 설명처럼 작업 폴더를 지정하고 SAS 프로그램을 그대로 사용하는 것이 편리한 것으로 보인다. 아니면 코멘트 부탁.

 

SAS OnDemand for Academics에서는 저장 공간 5Gb를 제공하는 것으로 보인다. 필요한 사람은 잘 써보자.

 

 

 

 

유전체 자료를 이용하여 실제 근교 계수 또는 혈연 계수를 구할 수 있다. 그러면 지금까지 혈통으로 구한 근교 계수 또는 혈연 계수는 실제가 아닌가? 그렇다. 혈통 정보에 기반하기 때문에 혈통이 없다면 또는 잘못되었다면 근교 계수 또는 혈연 계수가 0이 나오고, 혈통이 올바르다고 하여도 확률적인 근교 계수이다. 수 세대에 걸친 올바른 혈통이 있다면 확률적인 근교 계수와 실제 근교 계수의 차이가 크지 않을 것이다. 이 차이가 크다면 혈통이 올바르지 않거나 충분하지 않다고 생각할 수 있다. 여기서는 blupf90 family of programs의 하나인 pregsf90을 이용하여 유전체(실제) 근교 계수를 구해 본다. qcf90으로 구하면 좋겠지만 quality control과 근교 계수는 관계가 없는지 지금은 qcf90에서 구할 수 없다. pregsf90을 이용할 경우 renumf90 과정을 거쳐야 하고, 형식적이지만 자료 파일을 준비해야 한다는 귀찮음이 있다. 그래서 형식적인 자료로 단형질 개체 모형 예제의 자료를 이용한다.
 
혈통(animal, sire, dam)

K1 0 0
K2 0 0
K3 0 0
K4 K1 0
K5 K3 K2
K6 K1 K2
K7 K4 K5
K8 K3 K6

 
위 혈통 자료를 우리가 구하고자 하는 혈통에 붙인다.

 

 
자료(animal, sex, gain)

K4 1 4.5
K5 2 2.9
K6 2 3.9
K7 1 3.5
K8 1 5

 
파라미터 파일

# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
data.txt
TRAITS
3
FIELDS_PASSED TO OUTPUT

WEIGHT(S)

RESIDUAL_VARIANCE
40
EFFECT
2 cross numer
EFFECT
1 cross alpha
RANDOM
animal
FILE
pedigree.txt
FILE_POS
1 2 3
SNP_FILE
kbh_hv1.txt
PED_DEPTH
0
(CO)VARIANCES
20

 
snp 파일도 포함되어 있고, PED_DEPTH를 0으로 하여 혈통 파일의 모든 개체를 분석에 포함한다.
실행

renumf90 renumf90_stam.par | tee renumf90_stam_01.log

 

 
실행을 하면 여러 파일이 생긴다.
renf90.par : 새로운 파라미터 파일
renf90.inb : 근교계수
renf90.dat : 리넘버된 data 파일
renf90.fields : renf90.dat에서 각 필드에 대한 설명
renadd02.ped : 리넘버된 혈통
kbh_hv1.txt_XrefID : 유전체 자료 ID의 리넘버. cross reference id 등과 같은 파일이 생긴다.
 
renf90.par 파일을 수정하여 renf90_pregsf90.par 로 저장한다.

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)
 
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross 
 3     12629 cross 
RANDOM_RESIDUAL VALUES
   40.000    
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped                                                                    
(CO)VARIANCES
   20.000    
OPTION SNP_file kbh_hv1.txt
OPTION saveDiagHOrig
OPTION methodDiagH 1
OPTION saveGOrig

 
아래 세 줄을 추가한다. 
“OPTION methodDiagH 1” 이 부분은 “OPTION methodDiagH 2”로도 분석할 수 있다.

pregsf90 실행
SET OMP_STACKSIZE=64M 
pregsf90 renf90_pregsf90.par | tee pregsf90_01.log

 

 
diagH.txt

 
H matrix의 대각원소만 저장한 파일이다. 전 개체에 대해서 H matrix를 만들고 그 대각원소만 추출한 것이므로 이것들 중 유전체 자료를 가진 것만 추출하고 1을 빼면 유전체 근교 계수를 구할 수 있다.
OPTION methodDiagH 1을 OPTION methodDiagH 2로 변경하고 다시 실행해 보자.
실행

SET OMP_STACKSIZE=64M 
pregsf90 renf90_pregsf90.par | tee pregsf90_01.log
 

 
diagHdirect.txt

 

 

개체 ID, H matrix의 대각 원소, A matrix(혈통 기반 nemerator relationship matrix)의 대각 원소, 두 대각 원소의 차이.
A matrix의 대각 원소 값을 보면 혈통에 기반한 것이므로 1 또는 1.125등의 값을 보이지만, H matrix의 대각 원소는 유전체(실현) 근교 계수여서 매우 다양한 값을 볼 수 있다. 이 차이가 크다면 혈통 정보의 오류 또는 부족으로 볼 수 있다.
“OPTION methodDiagH 1”와 “OPTION methodDiagH 2”의 유전체 근교 계수를 비교해 보는 것도 좋을 것이다. 두 방법의 차이는 다음을 참조한다.

G_Orig.txt

 
G_Orig.txt 파일은 G matrix의 상삼각(또는 하삼각)을 보여준다. 이 값은 혈연 계수가 아니다. 혈연 계수는 근교 계수를 포함하여 계산한다. 하지만 근교 계수가 0이라고 가정한다면 혈연 계수라 볼 수 있다. 암튼 두 개체 사이의 관계 정도라 볼 수 있다. 혈통 정보가 부족할 경우 개체 사이의 혈연 관계 정도를 볼 수 있다.
유전체 자료는 있으나 혈통 정보가 부족한 소규모 집단의 계획 교배를 할 때 유전체 자료를 이용하여 개체의 근교 계수 또는 개체들 사이의 혈연 계수를 구하여 근친 교배를 피하고, 근교 퇴화를 예방할 수 있다.
 

육종을 하면서 집단에 있는 개체들의 근교계수를 계속 계산하는 것은 매우 중요한 일이다. 근친교배(inbreeding)가 일어나면 근교퇴화(inbreeding depression)가 발생하기 때문이다. 근교퇴화란 치사 유전자가 발현하여 죽어 버리거나 죽지 않더라도 경제 능력이 감소하는 현상을 말한다. 근교계수를 계산해 주는 수많은 프로그램들이 있지만 오늘은 inbupgf90을 이용하여 근교계수를 계산해 보자.
inbupgf90을 다운로드 하자.
http://nce.ads.uga.edu/html/projects/programs/Windows/64bit/

Index of /html/projects/programs/Windows/64bit

nce.ads.uga.edu

 
c:\ 밑에 blupf90이라는 폴더를 만들고 다운로드 받은 inbupgf90을 복사한다.

그러나 이 파일만으로는 실행이 안된다. 라이브러리 파일이 필요하다. libiomp5md.dll이라는 파일이 필요하다.
검색을 해서 다운로드 받거나 다음 링크를 연결하여 다운로드 한다.
https://ko.dll-files.com/libiomp5md.dll.html

libiomp5md.dll 무료 다운로드 | DLL‑files.com

libiomp5md.dll, 파일 설명: Intel(R) OMP Runtime Library libiomp5md.dll와(과) 관련된 오류는 몇 가지 다른 이유로 발생할 수 있습니다. 예를 들어, 응용 프로그램에 결함이 있거나, libiomp5md.dll이(가) PC에 존재하

ko.dll-files.com

 
다운로드하고 압축을 푼 결과

 
명령 프롬프트 실행 : 윈도우 시작 -> 찾기 창에서 cmd 입력

 
cd .. 명령어로 상위 폴더로 이동
cd .. 명령어로 한번 더 상위 폴더로 이동
cd blupf90으로 blupf90 폴더로 이동
dir 명령어로 폴더 내의 파일 확인

 
inbupgf90 명령어로 프로그램 실행

설치가 잘 된 것을 확인할 수 있다.마지막으로 모든 폴더에서 실행이 가능하도록 path를 설정하자.윈도우 시작 -> 설정 -> 정보 -> 고급 시스템 설정

 
환경 변수 클릭

 
시스템 변수 -> Path -> 편집 클릭

 
새로 만들기 클릭

 
C:\blupf90 입력하고 확인명령 프롬프트를 다시 실행하면 어떤 폴더에서든지 inbupgf90을 실행할 수 있다.

지금까지 프로그램을 실행하기 위한 환경 설정이었다. 한가지 더. D: 드라이브로 이동하기 위해서는 D: 명령어를 입력한다. 

다시 기억해 보자. 위로 올라가는 명령어는 "cd .." 그리고 폴더로 이동하는 명령어는 "cd 폴더이름"이다.
 
이제 자료를 준비하고 프로그램을 실행하여 근교계수를 계산해 보자. 그리고 미래에 태어날 개체의 근교계수를 계산해 보자.
근교 계수를 구할 개체들의 혈통 자료 준비

 
여기서는 12422두를 준비하였다. 부모를 모를 경우 0으로 기입한다.
여기서 개체, 아비 및 어미 사이를 탭을 사용하여 구분하지 말고, 스페이스를 사용하여 구분한다. 엑셀에서 copy & paste를 할 경우 탭으로 구분된 경우가 있다. 이럴 경우 탭을 복사해서 찾아서 바꾸기를 할 수도 있고, 에디터에서 탭을 스페이스로 변환할 수 있다. 각자 편한 방법으로.
다음과 같은 간단한 명령어로 근교계수 계산이 끝난다.

inbupgf90 --pedfile kbh_pedi_ori_202305.txt --method 3 | tee inbupgf90_01.log

 
--pedfile 뒤에 혈통 파일의 이름을 적어 준다.
근교계수 계산 방법을 지정할 수 있다. 1번과 2번은 Aguilar & Misztal, 2008의 방법을 이용하고 3번은 Meuwissen & Luo 1992의 방법을 이용한다. --method 옵션을 생략하면 2번 방법을 이용하는데 메모리에서 계산하여 빠르지만 많은 메모리를 요구한다. 그러나 그렇게 많은 혈통이 있나 모르겠다. 그냥 생략해도 된다.
| tee inbupgf90_01.log 는 inbupgf90_01.log 파일에 화면에 나오는 정보를 저장하고 나중에 천천히 보기 위함이다.
그런데 이 tee 명령이 실행될려면 윈도우에 프로그램을 설치해야 하는데 그건 다음 포스팅을 참조한다. 복잡하다고 느낀다면 " | tee inbupgf90_01.log" 없이 실행한다.

2009.03.17 - [ETC] - GnuWin32의 CoreUtils 패키지 설치하기

GnuWin32의 CoreUtils 패키지 설치하기

윈도우즈에서 유닉스 명령어 사용하기GnuWin32가 있음http://gnuwin32.sourceforge.net/ GnuWin32에는 여러 패키지가 있는데 CoreUtils에 내가 필요한 sort 프로그램이 있음http://www.gnu.org/software/coreutils/ CoreUtils Dow

bhpark.tistory.com

 
실행 화면

25% 초과가 4마리 있다.
마지막에 근교계수가 저장된 파일의 이름이 나온다.
계산된 근교계수 확인

 
근교계수가 0이라고 해서 정말 근친교배가 아닌 것은 아니다. 혈통에 기반한 것이기 때문에 혈통이 부족하면 0이 나온다. 그래서 혈통 완전도 지수(pedigree complete index)를 보고 근교계수를 판단해야 한다.
2015.11.24 - [Animal Breeding/Animal Breeding Etc] - Pedigree Completeness Index

Pedigree Completeness Index

Pedig의 ngen에서 만들어진 결과 파일을 이용하여 혈통완성도지수(Pedigree Completeness Index)를 계산한다. Pedigree completeness index란 ○ 지정한 세대까지 혈통이 얼마나 완전한지를 나타내는 지수 ○ 혈통

bhpark.tistory.com

 
근교계수를 계산했는데 태어난 개체의 근교계수를 계산하면 뭐하나. 태어나기 전에 계산을 해서 근교계수가 높을 것 같으면 그런 교배를 회피해야 한다. 그래서 inbupgf90에서는 미래에 태어날 개체의 근교계수를 계산해 준다. 교배에 사용할 아비와 어미 목록을 준비하면 각각의 아비와 각각의 어미가 교배하여 태어날 개체의 근교계수를 알려준다.
먼저 교배에 사용할 아비 목록을 준비한다.

 
교배에 사용할 어미 목록을 준비한다.

 
다음과 같이 명령어를 실행한다.

inbupgf90 --pedfile kbh_pedi_ori_202305.txt --sire_file list_sire.txt --dam_file list_dam.txt --matings | tee inbupgf90_01.log

 
--sire_file 준비한 아비 목록이 담긴 파일 이름
--dam_file 준비한 어미 목록이 담긴 파일 이름
--matings 태어날 개체의 근교계수 계산
실행 화면

여러 파일이 생기는데 .matings 파일을 보자.

               DamID        SireId_MinInb      KOR002154583225     KOR002154584201     KOR002154584156     KOR002154584130     KOR002154583321     KOR002166240033     KOR002160500515     KOR002166240308     KOR002154584148     KOR002166240381
  
     KOR002149839834      KOR002154583225                0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000
     KOR002148174888      KOR002154584201                3.906               0.781               4.688               4.688               1.562               3.906               8.594               5.078               5.469              16.797
     KOR002155610552      KOR002154583225                0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000               0.000
     KOR002155610489      KOR002154584201                1.172               0.391               2.344               1.562               0.781               1.172               1.172               0.977               1.172               1.367
     KOR002148174958      KOR002154583225                0.000               0.391               0.781               0.195               0.000               0.391               0.391               0.195               0.000               0.195
     KOR002155610938      KOR002154584201                1.855               0.879               7.617               2.539               1.758               1.465               1.660               1.416               1.660               1.709
     KOR002149839699      KOR002154584201                1.172               0.391               2.344               1.562               0.781               1.172               1.172               0.977               1.172               1.367
     KOR002149839842      KOR002154584201                3.320               1.172               2.734               2.930               2.344               1.953               1.758               2.148               2.539               2.148
     KOR002158246735      KOR002154584201                1.562               0.781               7.031               1.953               1.562               1.172               1.367               1.172               1.172               1.367
     KOR002149839666      KOR002154584201                1.172               0.391               2.344               1.562               0.781               1.172               1.172               0.977               1.172               1.367

어미별로 각 아비와 교배했을 때 태어날 개체의 근교계수를 계산해 주고, 그 중 가장 근교계수를 보인 아비를 표시해 준다. 그러나 그 아비의 근교계수만 최소인 것은 아니다. 다른 아비와 교배해도 0인 경우가 많다. 보고 고르면 된다.
_relationships 파일은 아비 사이의 혈연 계수를 보여준다.

   KOR002154583225      KOR002154583225           1.047
   KOR002154583225      KOR002154584201           0.105
   KOR002154583225      KOR002154584156           0.086
   KOR002154583225      KOR002154584130           0.133
   KOR002154583225      KOR002154583321           0.180
   KOR002154583225      KOR002166240033           0.285
   KOR002154583225      KOR002160500515           0.090
   KOR002154583225      KOR002166240308           0.111
   KOR002154583225      KOR002154584148           0.129
   KOR002154583225      KOR002166240381           0.078
   KOR002154584201      KOR002154584201           1.000
   KOR002154584201      KOR002154584156           0.047
   KOR002154584201      KOR002154584130           0.080
   KOR002154584201      KOR002154583321           0.195
   KOR002154584201      KOR002166240033           0.125
   KOR002154584201      KOR002160500515           0.172
   KOR002154584201      KOR002166240308           0.273
   KOR002154584201      KOR002154584148           0.152
   KOR002154584201      KOR002166240381           0.086
   KOR002154584156      KOR002154584156           1.062
   KOR002154584156      KOR002154584130           0.125
   KOR002154584156      KOR002154583321           0.078
   KOR002154584156      KOR002166240033           0.078
   KOR002154584156      KOR002160500515           0.195
   KOR002154584156      KOR002166240308           0.070
   KOR002154584156      KOR002154584148           0.070
   KOR002154584156      KOR002166240381           0.141
   KOR002154584130      KOR002154584130           1.031
   KOR002154584130      KOR002154583321           0.156
   KOR002154584130      KOR002166240033           0.072
   KOR002154584130      KOR002160500515           0.090
   KOR002154584130      KOR002166240308           0.087
   KOR002154584130      KOR002154584148           0.305
   KOR002154584130      KOR002166240381           0.074
   KOR002154583321      KOR002154583321           1.000
   KOR002154583321      KOR002166240033           0.086
   KOR002154583321      KOR002160500515           0.117
   KOR002154583321      KOR002166240308           0.113
   KOR002154583321      KOR002154584148           0.086
   KOR002154583321      KOR002166240381           0.043
   KOR002166240033      KOR002166240033           1.016
   KOR002166240033      KOR002160500515           0.070
   KOR002166240033      KOR002166240308           0.117
   KOR002166240033      KOR002154584148           0.121
   KOR002166240033      KOR002166240381           0.082
   KOR002160500515      KOR002160500515           1.016
   KOR002160500515      KOR002166240308           0.133
   KOR002160500515      KOR002154584148           0.117
   KOR002160500515      KOR002166240381           0.293
   KOR002166240308      KOR002166240308           1.023
   KOR002166240308      KOR002154584148           0.154
   KOR002166240308      KOR002166240381           0.223
   KOR002154584148      KOR002154584148           1.016
   KOR002154584148      KOR002166240381           0.123
   KOR002166240381      KOR002166240381           1.055

 
_statistics 파일은 아비를 기준으로 평균, 최소, 최대 근교계수를 보여 준다.

   1 KOR002154583225           1.416     1.172     3.906 3
   2 KOR002154584201           0.518     0.391     1.172 7
   3 KOR002154584156           2.988     0.781     7.617 0
   4 KOR002154584130           1.699     0.195     4.688 0
   5 KOR002154583321           0.957     0.781     2.344 0
   6 KOR002166240033           1.240     0.391     3.906 0
   7 KOR002160500515           1.729     0.391     8.594 0
   8 KOR002166240308           1.294     0.195     5.078 0
   9 KOR002154584148           1.436     1.172     5.469 0
   10 KOR002166240381           2.632     0.195    16.797 0

마지막 열의 수는 _matings 파일에 추천 아비로 몇 번 되었는지를 나타내는 수이다.
교배계획을 어떻게 할 것인가? 근교계수가 최소로 상승하면서 최대의 능력을 발휘하는 교배가 좋은 교배이긴 하겠지만 그게 쉽지 않다. 어떤 아비가 집단의 어미와 근교도 안 되고 능력도 최고치로 올린다고 해서 모든 어미에 교배하면 그 다음 세대에 곤란한 일이 발생할 것이다. 이와 관련된 문제를 멋지게 풀어 보려고 Optimal Contribution Selection이란 개념도 생긴 거 같다. 
사용 가능한 아비를 골라 보고, 어미를 기준으로 어느 정도 이하의 근교계수를 보이는 아비들을 골라보고 그 아비 중 최대의 능력을 보이는 아비를 골라보자. 그리고 전체적으로 골고루 골라졌나 보자. 한 두 개의 아비로 몰렸다면 적당히 조정하자.

유전 평가하고 상관 없으나 혈통이 있을 때 부모 또는 자식으로 연결하였을 때 몇 개의 그룹으로 나누어질지 궁금할 때가 있다. 물론 그림을 그리면 되긴 하지만, 몇십 개 이상은 그림을 그리기가 힘들다. 물론 포트란 소스로 올린 적도 있으나 컴파일하기도 힘들어서 누가 쓸까 싶다. 그래서 간단히 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개마다 할 지 수정해야 한다.

 

 

 

seekparentf90을 이용한 친자감정과 부모 찾기는 이미 포스팅한 적이 있다.

2020.02.12 - [Animal Breeding/Genomic Selection] - SNP 유전체 자료를 이용한 아비 어미 찾기

 

SNP 유전체 자료를 이용한 아비 어미 찾기

imputation 또는 qc 과정에서, 사용하는 프로그램에 따라 친자감정 결과를 제시하거나 또는 진짜 부모라고 여겨지는 개체를 찾아 주기도 한다. 여기서는 아비, 어미를 본격적으로 찾아주는 프로그

bhpark.tistory.com

 

qcf90 또는 pregsf90 프로그램을 통해서 대부분의 SNP genotype quality control을 수행할 수 있다. 그러나 이런 프로그램들은 여러 버전의 마이크로 어레이를 동시에 이용하여 quality control을 할 수 없고, 하려 한다면 imputation을 해야 한다. imputation이 필요하다면 imputation을 해서 quality control을 하면 되나, raw data를 그대로 이용하여 친자감정을 하거나 부모를 찾으려 하는 경우도 있다. seekparentf90은 약간(?)의 자료 가공을 거쳐 imputation 하지 않고, 친자 감정하거나 부모를 찾을 수 있는 기능을 제공하고 있다.

 

먼저 BovineSNP50V2, BovineSNP50V3, Hanwoo50V1의 GenomeStudio에서 출력한 finalreport를 illumina2pregs 프로그램을 이용하여 coded하였다고 하자. 그럴 경우 각 버전의 genotype과 map file은 다음과 같다.

 

BovineSNP50V2

- genotype

 

- map file

 

BovineSNP50V3

- genotype

 

- map file

 

Hanwoo50V1

- genotype

 

- map file

 

먼저 map file을 합쳐야 한다.

 

세 개의 map file을 읽어서 seekparentf90을 위한 map file 만드는 SAS 프로그램

/*
file name : make_snpmapfile_for_seekparentf90.sas

seekparentf90 프로그램을 위한 map file 만들기

BovineSNP50V2, BovineSNP50V3, HanwooSNP50V1 map file을 하나로 합치기
*/

/* BovineSNP50V2 snp_map 파일 읽기 */
data a1 ;
    infile 'snp_map_bv2' firstobs = 2;

    informat name $37. ;
    informat bv2chrmo best32. ;
    informat bv2pos best32. ;
    informat bv2oindex best32. ;

    format name $37. ;
    format bv2chrmo best32. ;
    format bv2pos best32. ;
    format bv2oindex best32. ;

    input 
            name
            bv2chrmo
            bv2pos
            bv2oindex
    ;

    bv2index = _n_;
run;

/* BovineSNP50V3 snp_map 파일 읽기 */
data b1 ;
    infile 'snp_map_bv3' firstobs = 2;

    informat name $37. ;
    informat bv3chrmo best32. ;
    informat bv3pos best32. ;
    informat bv3oindex best32. ;

    format name $37. ;
    format bv3chrmo best32. ;
    format bv3pos best32. ;
    format bv3oindex best32. ;

    input 
            name
            bv3chrmo
            bv3pos
            bv3oindex
    ;

    bv3index = _n_;
run; 

/* HanwooSNP50V1 snp_map 파일 읽기 */
data c1 ;
    infile 'snp_map_hv1' firstobs = 2;

    informat name $37. ;
    informat hv1chrmo best32. ;
    informat hv1pos best32. ;
    informat hv1oindex best32. ;

    format name $37. ;
    format hv1chrmo best32. ;
    format hv1pos best32. ;
    format hv1oindex best32. ;

    input 
            name
            hv1chrmo
            hv1pos
            hv1oindex
    ;

    hv1index = _n_;
run;
 
/* bv2, bv3, hv1 머지하기 */
    /* sort */
proc sort data = a1;
    by name;
run;

proc sort data = b1;
    by name;
run;

proc sort data = c1;
    by name;
run;

    /* 머지하기 */
data t1;
    merge a1 b1 c1;
    by name;
run;

/* 이름은 같은데 염색체, 위치가 다른 SNP */
/* 있으면 안된다. 있다면 아마도 genome build가 칩마다 다른 것일 수도 */
data t1_e1;
    set t1;

	/* bv2와 bv3 비교 */
	if bv2chrmo ~= . and bv3chrmo ~= . and bv2chrmo ~= bv3chrmo then output;

	/* bv3와 hv1 비교 */
	if bv3chrmo ~= . and hv1chrmo ~= . and bv3chrmo ~= hv1chrmo then output;
run;

/* index, 염색체 번호, 위치 번호 정리 */
data t2;
    set t1;

    /* 없었던 SNP index는 0 */
    if bv2index = . then bv2index = 0;
    if bv3index = . then bv3index = 0;
    if hv1index = . then hv1index = 0;

    /* 각 SNP의 염색체 번호 */
    chr = bv2chrmo;
    if bv2chrmo   = . then chr = bv3chrmo;
    if bv2chrmo   = . and bv3chrmo = . then chr = hv1chrmo;

    /* 각 SNP의 위치 */
    pos = bv2pos;
    if bv2pos   = . then pos = bv3pos;
    if bv2pos   = . and bv3pos = . then pos = hv1pos;

run;

/* dataset을 열어서 확인 */
proc sort data = t2;
    by chr pos;
run;

data _null_;
    set t2;
    file 'mapfile_for_seekparentf90.txt';
    if _n_ = 1 then put "SNP_ID                                  Chr       Pos       chip1     chip2     chip3";
    put name @41 chr @51 pos @61 bv2index @71 bv3index @81 hv1index;
run;

 

결과적으로 만들어진 map file은 다음과 같다. 

 

각 genotype 파일에 chip number를 넣어야 한다.

BovineSNP50V2에 chip number 1을 넣는 awk 프로그램은 다음과 같다. insert_chip_number_bv2.awk로 저장한다.

BEGIN{
}
{
    print $1 " 1 " $2 > "kbh_bv2_02.txt"
}
END {
}

 

위 awk 프로그램을 실행하는 방법은 다음과 같다.

 

awk -f insert_chip_number_bv2.awk kbh_bv2.txt
pause

 

위 내용을 run_insert_chip_number_bv2.awk.bat에 저장하면 더블 클릭하여 실행할 수 있다.

 

결과적으로 생긴 파일은 다음과 같다.

 

 

BovineSNPV3 genotype에 chip number 2를 넣는 awk 프로그램

 

BEGIN{
}
{
    print $1 " 2 " $2 > "kbh_bv3_02.txt"
}
END {
}

 

위 소스를 insert_chip_number_bv3.awk으로 저장

 

위 awk 프로그램을 실행하는 방법

awk -f insert_chip_number_bv3.awk kbh_bv3.txt
pause

 

위 명령어를 run_insert_chip_number_bv3.awk.bat에 저장하면 더블 클릭하여 실행할 수 있다.

 

결과적으로 생긴 파일은 다음과 같다. 

 

HanwooSNPV1 genotype에 chip number 3을 넣는 awk 프로그램

 

BEGIN{
}
{
    print $1 " 3 " $2 > "kbh_hv1_02.txt"
}
END {
}

 

위 소스를 insert_chip_number_hv1.awk으로 저장

 

위 awk 프로그램을 실행하는 방법

awk -f insert_chip_number_hv1.awk kbh_hv1.txt
pause

 

위 명령어를 run_insert_chip_number_hv1.awk.bat에 저장하면 더블 클릭하여 실행할 수 있다.

 

결과적으로 생긴 파일은 다음과 같다.

 

 

위 세 genotype 파일을 하나로 합치는 명령어

cat kbh_bv2_02.txt kbh_bv3_02.txt kbh_hv1_02.txt > kbh_genotype.txt
pause

 

위 명령어를 make_one_genotype_file.bat으로 저장하면 더블 클릭하여 실행할 수 있다.

 

결과적으로 생긴 파일은 다음과 같다.

 

 

친자감정을 하는 것이므로 개체의 아비와 어미가 누구인지를 나타내는 혈통파일이 있어야 한다. 준비한 혈통 파일은 다음과 같다(kbh_pedigree_yob.txt)

 

이제 seekparentf90으로 친자감정을 하고, 친자감정 불일치일 경우 아비 또는 어미를 찾아 보자.

 

seekparentf90을 실행하는 명령어는 다음과 같다.

seekparentf90 --pedfile kbh_pedigree_yob.txt --snpfile kbh_genotype.txt --chips mapfile_for_seekparentf90.txt --yob --seeksire_in_ped --seekdam_in_ped --exclude_chr 30 31 32 33 --seektype 1 --duplicate --full_log_checks | tee seekparentf90_01.log
pause

 

위 명령어를 run_seekparentf90.bat으로 저장하면 더블클릭으로 실행할 수 있다.

 

명령어 설명

 

seekparentf90 명령어

--pedfile kbh_pedigree_yob.txt 혈통 파일. 혈통 파일은 header 없음.

--snpfile kbh_genotype.txt 위에서 준비한 genotype 파일

--chips mapfile_for_seekparentf90.txt 위에서 준비한 map file

--yob 혈통의 넷째 컬럼에 출생연도가 있음

--seeksire_in_ped 혈통에서 아비 찾기

--seekdam_in_ped 혈통에서 어미 찾기

--exclude_chr 30 31 32 33 30에서 33번 염색체는 제외

--seektype 1 불일치인 경우만 아비 및 어미 찾기

--duplicate 유전자형이 중복인 개체들이 있는지 검사

--full_log_checks 많은 로그 기록

| tee seekparentf90_01.log 화면에 프린트되는 것을 파일에도 기록

 

로그는 다음과 같다.

  ---------------------------------------------------- 
 |                  SeekParentf90                     |
 |              Parent-Offspring tests                |
 |            and detection of paternity              |
 |               based on SNP markers                 |
 |                                                    |
 |                2013 - Version 1.47                 |
 |           (Last update:  Dec 15, 2022)             |
 |                                                    |
 |                                                    |
 |             INIA Las Brujas, Uruguay               |
 |        University of Georgia, Athens, US           |
  ---------------------------------------------------- 

 Options
    --pedfile kbh_pedigree_yob.txt
    --snpfile kbh_genotype.txt
    --yob
    --seeksire_in_ped
    --seekdam_in_ped
    --seektype 1
    --chips mapfile_for_seekparentf90.txt
    --exclude_chr 30 31 32 33
    --full_log_checks

 Using Jenkins hash function

 Maximum length for alphanumeric fields: 20

 Fields to read:  animal, sire, dam, yob (4)

 Elapsed time for reading and hash: .250

 Pedigree file "kbh_pedigree_yob.txt": 12412 records

 Total number of animals in pedigre: 12412
 
 Statistics for Year of Birth
  Minimun:         1990
  Maximun:         2023

 Read multiple Chip SNP file: "kbh_genotype.txt"

 Number of SNP from snp_info file: 60961

 ***** WARNING *****  Number of duplicate SNP by chr_pos: 93
 ***** List of duplicate SNP: "duplicate_snp_chr_pos"

 Exclude Chr(s)
   Chr:  30  1190
   Chr:  31    11
   Chr:  32    13
   Chr:  33   775

 Chip identification
   Chip   1: CHIP1
   Chip   2: CHIP2
   Chip   3: CHIP3
 

 Number of SNP by Chip
   Chip   1:    54609   50908   47514
   Chip   2:        0   53218   49760
   Chip   3:        0       0   53866
 
 Exclude SNP: 1989
 Maximum Available SNP for calculations: 58972
 

 Number of effective SNP by Chip for parent-progeny conflicts
   Chip   1:    52886   49233   46109
   Chip   2:        0   51278   48090
   Chip   3:        0       0   52195
    Column position in file for the first marker: 19
    Format to read SNP file: (18x,800000i1)                                                                                      
    Number of SNPs: 60961

 Threshold for exclusions - Maximun number of SNP: 58972
 Number of SNP > 130, set the threshold probability to exclude= 1.000
 Number of SNP > 130, set the threshold probability to assign = .500
 Threshold based on percent of conflicts to exclude: 1.000
 Threshold based on percent of conflicts to assign: .500

 Number of records in SNP file: 1884
 Number of records present in pedigree: 1883 1883

 Number of records by Chip
    Chip 1: 21
    Chip 2: 252
    Chip 3: 1611

 Call Rate
   Min CR:   5346 .90
   Max CR:     42 1.00
   AVG CR:    668 .99
   Number of samples with low call rate (  0.90 ):      0
 
 Samples with low call rate for effective SNP           0

 Check all pedigree for parent-progeny conflicts

 Pedigree records with genotype: 1883

 Records not tested (low call rate <  0.90 ): 0

 Records tested: 1883
 
 Pair parent/progeny tested: 1403
 Pair with conflicts: 19   1.4 %
   NOT tested (parent low call rate): 0
 
     Sire-progeny tested: 715
     Sire-progeny with conflicts: 11   1.5 %
       NOT tested (parent low call rate): 0
 
     Dam-progeny tested: 688
     Dam-progeny with conflicts: 8   1.2 %
       NOT tested (parent low call rate): 0

 Seek parent for SIRES
   Get list of parents: detected from pedigree        
     Number of Sires: 171

   Seek method: Only for animals with conflicts    

   Total number of animals: 11
      Found: 7
      Not Found: 4
      Found but cannot assign: 0

 Seek parent for DAMS
   Get list of parents: detected from pedigree        
     Number of Dams : 666

   Seek method: Only for animals with conflicts    

   Total number of animals: 8
      Found: 3
      Not Found: 5
      Found but cannot assign: 0

 Output files
     Pedigree file with removed parents for animals with conflicts: "Check_kbh_pedigree_yob.txt"

 Full reports
     "Parent_Progeny_Conflicts.txt"
     "Parent_Progeny_Conflicts_Summary.txt"
     "Seek_Sire.txt"
     "Seek_Dam.txt"
 

 

중간에

 

***** WARNING ***** Number of duplicate SNP by chr_pos: 93

***** List of duplicate SNP: "duplicate_snp_chr_pos“

 

로그가 있는데 duplicate_snp_chr_pos 파일을 열어 보면 다음과 같다,

 
Duplicate SNP chr_pos: 1_59409838 2
Duplicate SNP chr_pos: 13_25606469 2
Duplicate SNP chr_pos: 3_58040470 2
Duplicate SNP chr_pos: 14_27271835 2
Duplicate SNP chr_pos: 29_28647816 2
Duplicate SNP chr_pos: 15_15014275 2
Duplicate SNP chr_pos: 2_50689222 2
Duplicate SNP chr_pos: 11_33265117 2
Duplicate SNP chr_pos: 7_1071389 2
Duplicate SNP chr_pos: 15_79531511 2
Duplicate SNP chr_pos: 9_45729853 2
Duplicate SNP chr_pos: 25_10359385 2
Duplicate SNP chr_pos: 22_56526462 2
Duplicate SNP chr_pos: 20_676757 2
Duplicate SNP chr_pos: 22_8308367 2
Duplicate SNP chr_pos: 20_17837675 2
Duplicate SNP chr_pos: 14_56761589 2
Duplicate SNP chr_pos: 3_76797698 2
Duplicate SNP chr_pos: 25_26191380 2
Duplicate SNP chr_pos: 19_43264699 2
Duplicate SNP chr_pos: 21_65198296 2
Duplicate SNP chr_pos: 8_106174871 2
Duplicate SNP chr_pos: 15_15367990 2
Duplicate SNP chr_pos: 2_111155237 2
Duplicate SNP chr_pos: 14_76524093 2
Duplicate SNP chr_pos: 13_71212300 2
Duplicate SNP chr_pos: 14_25401722 2
Duplicate SNP chr_pos: 12_80629629 2
Duplicate SNP chr_pos: 10_49791638 2
Duplicate SNP chr_pos: 6_54524018 2
Duplicate SNP chr_pos: 3_40654331 2
Duplicate SNP chr_pos: 15_15595454 2
Duplicate SNP chr_pos: 7_21298998 2
Duplicate SNP chr_pos: 15_15738972 2
Duplicate SNP chr_pos: 7_18454636 2
Duplicate SNP chr_pos: 7_80042191 2
Duplicate SNP chr_pos: 26_38233337 2
Duplicate SNP chr_pos: 28_35331560 2
Duplicate SNP chr_pos: 15_77675440 2
Duplicate SNP chr_pos: 15_15908105 2
Duplicate SNP chr_pos: 26_8221270 2
Duplicate SNP chr_pos: 18_1839733 2
Duplicate SNP chr_pos: 15_15697628 2
Duplicate SNP chr_pos: 11_1703612 2
Duplicate SNP chr_pos: 15_21207529 2
Duplicate SNP chr_pos: 15_79187295 2
Duplicate SNP chr_pos: 22_31841994 2
Duplicate SNP chr_pos: 20_51449833 2
Duplicate SNP chr_pos: 15_15544958 2
Duplicate SNP chr_pos: 15_15520367 2
Duplicate SNP chr_pos: 15_15656347 2
Duplicate SNP chr_pos: 15_15305071 2
Duplicate SNP chr_pos: 6_61079420 2
Duplicate SNP chr_pos: 15_15213319 2
Duplicate SNP chr_pos: 16_38116988 2
Duplicate SNP chr_pos: 28_44261945 2
Duplicate SNP chr_pos: 9_10783460 2
Duplicate SNP chr_pos: 14_48380429 2
Duplicate SNP chr_pos: 29_13695183 2
Duplicate SNP chr_pos: 1_136603695 2
Duplicate SNP chr_pos: 15_15162470 2
Duplicate SNP chr_pos: 15_15245842 2
Duplicate SNP chr_pos: 3_116448759 2
Duplicate SNP chr_pos: 15_38078775 2
Duplicate SNP chr_pos: 15_15190242 2
Duplicate SNP chr_pos: 8_95410507 2
Duplicate SNP chr_pos: 3_91910014 2
Duplicate SNP chr_pos: 17_60069612 2
Duplicate SNP chr_pos: 9_98483346 2
Duplicate SNP chr_pos: 12_79289169 2
Duplicate SNP chr_pos: 1_151349514 2
Duplicate SNP chr_pos: 4_17200594 2
Duplicate SNP chr_pos: 19_47734925 2
Duplicate SNP chr_pos: 26_47781927 2
Duplicate SNP chr_pos: 8_88974063 2
Duplicate SNP chr_pos: 14_29987025 2
Duplicate SNP chr_pos: 14_31322421 2
Duplicate SNP chr_pos: 14_21343649 2
Duplicate SNP chr_pos: 1_1277227 2
Duplicate SNP chr_pos: 24_4118163 2
Duplicate SNP chr_pos: 13_49053256 2
Duplicate SNP chr_pos: 10_55611885 2
Duplicate SNP chr_pos: 28_8508619 2
Duplicate SNP chr_pos: 6_55756520 2
Duplicate SNP chr_pos: 15_15949175 2
Duplicate SNP chr_pos: 4_94176209 2
Duplicate SNP chr_pos: 14_23893220 2
Duplicate SNP chr_pos: 25_26198573 2
Duplicate SNP chr_pos: 16_33671664 2
Duplicate SNP chr_pos: 3_71860062 2
Duplicate SNP chr_pos: 5_63150400 2
Duplicate SNP chr_pos: 7_37537822 2
Duplicate SNP chr_pos: 9_32549297 2

 

93개의 염색체 번호와 위치에서 염색체 번호와 위치는 같은데 서로 다른 SNP 이름을 갖고 있다는 에러를 뿌려줌. 유전체 선발에서는 제외를 하겠지만 여기서는 제외를 해야 올바른 것인지 잘 모르겠다. 분석에서 93쌍 즉 186개 SNP를 이용했는지 안 했는지 모르겠다.

 

Pair parent/progeny tested: 1403

Pair with conflicts: 19 1.4 %

 

1403개 부모/자식을 테스트해서 19개 관계가 친자감정 불일치

 

Sire-progeny tested: 715

Sire-progeny with conflicts: 11 1.5 %

 

715개 아비/자식을 테스트해서 11개 관계가 불일치

 

Dam-progeny tested: 688

Dam-progeny with conflicts: 8 1.2 %

 

688개 어미/자식을 테스트해서 8개 관계가 불일치

 

parent_progeny_conflicts.txt

 

개체와 부모의 친자감정 결과를 보여준다. 그러나 너무 자세해서 보기 어렵다.

 

정렬을 해서 animal-dam 그리고 animal-sire 부분을 엑셀로 복사하고, D열을 내림차순 정렬하면 친자감정 불일치 개체 19두를 찾을 수 있다.

 

친자감정 불일치가 일어나는 경우는 세 가지이다. 1) 부(모)의 유전자형이 잘못 되었다. 2) 자손의 유전자형이 잘못되었다. 3) 혈통이 잘못되었다.

 

친자감정 불일치가 일어날 경우 1)과 2)를 생각안하고 무턱대고 3)을 생각할 경우 심각한 문제를 초래할 수 있다. 그러면 유전자형이 정말 제대로 되었는지는 어떻게 알 수 있는가?

 

parent_progeny_conflicts_summary.txt를 보자.

 

KOR002085189555 개체의 경우는 아비로서 9번 친자감정을 했고 그중 1번만 불일치다. 그러면 8번은 맞은 경우이므로 자기 자신의 유전자형이라고 확신할 수 있다.

 

그러나 KOR002072558437 는 개체로서 2번 친자감정 되었는데 2번다 친자감정 불일치가 났다. 이런 경우 자신의 유전자형이라고 보기 매우 힘들다. (아니라는 확신은 못한다.) 다시 genotyping을 할 수 있다면 하고, 다른 친자감정 결과 예를 들어 Microsatellite 친자감정 결과가 있다면 참고하여 결정하여야 한다.

 

조직을 샘플링할 때, 실험할 때, 결과를 처리할 때 등등 자신의 유전자형이 아닐 확률은 너무나 많다. 위와 같이 친자감정 일치가 하나라도 나왔다면 자신의 유전자형이라고 확신할 수 있지만 친자감정 일치가 없다면 자신의 유전자형인지 강한 의심을 하고, 그다음으로 혈통을 고민해야 한다.

 

A라는 개체가 B의 유전자형을 들고 우수한 개체로 판명이 나고, 그 결과를 이용하여 부모로 자손을 생산한다면 얼마나 큰 실수인가.

아래는 자신의 유전자형을 확신할 때 보는 부분이다. 즉 자신의 유전자형이 맞으므로 혈통을 수정해야 할 경우 보는 부분이다.

 

Seek method: Only for animals with conflicts

 

Total number of animals: 11

Found: 7

Not Found: 4

Found but cannot assign: 0

 

아비 불일치 11개 중 7개를 찾았다.

 

Seek method: Only for animals with conflicts

 

Total number of animals: 8

Found: 3

Not Found: 5

 

어미 불일치 8개 중 3개를 찾고, 5개는 찾지 못했다.

 

Output files

Pedigree file with removed parents for animals with conflicts: "Check_kbh_pedigree_yob.txt“

 

불일치인 개체의 부 또는 모 정보를 제거한 혈통 파일

 

아비 찾기 결과 seek_sire.txt

 

 

위에선 7개의 아비를 찾았다고 했으나 8개를 찾은 결과가 나옴. 이건 뭐지?

 

seek_dam.txt

 

 

3개의 어미를 찾았다.

SNP 정보를 기반으로 아비와 어미를 찾아주었다고 하여 무조건 믿는 것도 안된다. 같은 농장인지, 같은 농장이 아니라면 인공수정이나 수정란 이식인지, 부모의 나이가 많은지 등등을 종합적으로 살펴보고 수정하여야 한다.
 

이렇게 자기 자신의 올바른 유전자형과 올바른 혈통이 있다면 imputation 정확도가 매우 높아진다. 

 

요약

1) 파일 준비
 * genotype file
  - kbh_bv2.txt
  - kbh_bv3.txt
  - kbh_hv1.txt
 
 * chip 번호 넣는 awk 파일과 awk을 실행시키는 batch 파일
  - insert_chip_number_bv2.awk, run_insert_chip_number_bv2.awk.bat
  - insert_chip_number_bv3.awk, run_insert_chip_number_bv3.awk.bat
  - insert_chip_number_hv1.awk, run_insert_chip_number_hv1.awk.bat

 * chip 번호 넣은 genotype file을 합치는 파일
  - make_one_genotype_file.bat

 * map file
  - snp_map_bv2
  - snp_map_bv3
  - snp_map_hv1
  
 * 통합 map file을 만드는 SAS 프로그램
  - make_snpmapfile_for_seekparentf90.sas 

 * 혈통 파일
  - kbh_pedigree_yob.txt
  
 * seekparentf90 실행 파일
  - run_seekparentf90.bat

2) 실행
 * genotyp file 준비
  - run_insert_chip_number_bv2.awk.bat
  - run_insert_chip_number_bv3.awk.bat
  - run_insert_chip_number_hv1.awk.bat
  - make_one_genotype_file.bat
  
  * map file 준비
   - make_snpmapfile_for_seekparentf90.sas (새로운 칩이 추가 되지 않는다면 기존 결과(mapfile_for_seekparentf90.txt) 사용 가능)
   
  *  seekparentf90 실행
   - run_seekparentf90.bat
   

seekparentf90_for_multiple_chips.zip
0.00MB

+ Recent posts