Let's Study Animal Breeding

개체이름으로, 일부의 SNP만, 일루미나 SNP 이름으로 Final Report를 다운로드하면 할 일이 끝난 줄 알았으나 끝이 아니었다. Final Report와 SNP Map 파일의 SNP 순서가 달라서 일관된 순서 즉 SNP 이름 순으로 정렬을 해야 한다. 그래야 다음 프로그램(illumina2pregs)이 동작을 한다. Long Format Exprt Tool도 끝까지 사람을 귀찮게 한다. 어쩌랴, 하는 수 밖에 없다. 

Final Report는 개체, SNP 이름순으로 정렬된 것처럼 보이나, SNP Map 파일은 SNP 이름순으로 정렬된 것으로 보이지 않는다. Final Report도 확실히 정렬되어 있는지 알 수가 없으므로 그냥 정렬한다.정렬 프로그램이 문제다. 윈도우즈의 정렬 프로그램은 생각처럼 동작하지 않는다. 유닉스(리눅스)의 정렬 명령어를 사용하기 바란다. 윈도우에서 리눅스 명령어를 사용하는 가장 좋은 방법은 마이크로소프트의 권고대로 WSL을 사용하는 것이다. WSL을 사용하기 싫다면 리눅스 컴퓨터로 자료를 옮긴 후 정렬하기 바란다. 파일의 앞 뒤 내용을 확인하고, 행의 개수를 세고, 분리하고, 정렬하고, 합치는데 불과 간단한 몇 개의 리눅스 명령어로 해결이 가능하다.

wc -l output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt
head output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt
tail output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt

head -n 1 output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt > abv3_snpmap_header.txt
head abv3_snpmap_header.txt

tail -n +2 output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt > abv3_snpmap_body.txt
wc -l abv3_snpmap_body.txt

sort abv3_snpmap_body.txt > abv3_snpmap_body_sorted.txt
head abv3_snpmap_body_sorted.txt
tail abv3_snpmap_body_sorted.txt

cat abv3_snpmap_header.txt abv3_snpmap_body_sorted.txt > abv3_snpmap.txt

wc -l output_01_Axiom_BovMDv3.r2_newsnpname.txt
head -n 15 output_01_Axiom_BovMDv3.r2_newsnpname.txt
tail -n 15 output_01_Axiom_BovMDv3.r2_newsnpname.txt

head -n 11 output_01_Axiom_BovMDv3.r2_newsnpname.txt > abv3_finalreport_header.txt
cat abv3_finalreport_header.txt

tail -n +12 output_01_Axiom_BovMDv3.r2_newsnpname.txt > abv3_finalreport_body.txt
head abv3_finalreport_body.txt
tail abv3_finalreport_body.txt

sort -k2,2 -k1,1 abv3_finalreport_body.txt > abv3_finalreport_body_sorted.txt
head abv3_finalreport_body_sorted.txt
tail abv3_finalreport_body_sorted.txt

cat abv3_finalreport_header.txt abv3_finalreport_body_sorted.txt > abv3_finalreport.txt

illumina2pregs --snpfile abv3_finalreport.txt --mapfile  abv3_snpmap.txt --sortmap --snpfile_out snp_abv3.txt --mapfile_out snpmap_abv3.txt | tee illumia2pregs_01.log

wc -l output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt

행의 개수

head output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt

파일 앞 부분 내용 확인

tail output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt

파일 뒷 부분 내용 확인

head -n 1 output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt > abv3_snpmap_header.txt

한 줄 분리하여 저장

head abv3_snpmap_header.txt

분리한 내용 확인

tail -n +2 output_01_Axiom_BovMDv3.r2_SNPMap_newsnpname.txt > abv3_snpmap_body.txt

2줄 이후 분리하여 저장

wc -l abv3_snpmap_body.txt

행의 개수 확인

sort abv3_snpmap_body.txt > abv3_snpmap_body_sorted.txt

정렬

head abv3_snpmap_body_sorted.txt

파일의 앞 부분 확인

tail abv3_snpmap_body_sorted.txt

파일의 뒷 부분 확인

cat abv3_snpmap_header.txt abv3_snpmap_body_sorted.txt > abv3_snpmap.txt

합쳐서 저장

wc -l output_01_Axiom_BovMDv3.r2_newsnpname.txt

행의 개수 확인

head -n 15 output_01_Axiom_BovMDv3.r2_newsnpname.txt

파일의 앞 부분 확인

tail -n 15 output_01_Axiom_BovMDv3.r2_newsnpname.txt

파일의 뒷 부분 확인

head -n 11 output_01_Axiom_BovMDv3.r2_newsnpname.txt > abv3_finalreport_header.txt

파일의 앞 11줄 분리하여 저장

cat abv3_finalreport_header.txt

분리한 파일 확인

tail -n +12 output_01_Axiom_BovMDv3.r2_newsnpname.txt > abv3_finalreport_body.txt

파일의 12줄 이후 분리하여 저장

head abv3_finalreport_body.txt

파일의 앞 부분 확인

tail abv3_finalreport_body.txt

파일의 뒷 부분 확인

sort -k2,2 -k1,1 abv3_finalreport_body.txt > abv3_finalreport_body_sorted.txt

2열, 1일 순으로 정렬하여 저장

head abv3_finalreport_body_sorted.txt

정렬한 파일의 앞 부분 확인

tail abv3_finalreport_body_sorted.txt

정렬한 파일의 뒷 부분 확인

cat abv3_finalreport_header.txt abv3_finalreport_body_sorted.txt > abv3_finalreport.txt

합쳐서 저장

illumina2pregs --snpfile abv3_finalreport.txt --mapfile  abv3_snpmap.txt --sortmap --snpfile_out snp_abv3.txt --mapfile_out snpmap_abv3.txt | tee illumia2pregs_01.log

illumin2pregs 실행

(드디어 Axiom 칩으로 실험을 하였을 경우 Final Report로 출력하는 방법을 모두 설명하였다. 이걸 알아 내는데 (나 빼고) 많은 사람이  고생하였다. 나는 그저 어렵게 알아낸 걸 배워서 블로그에 올릴 뿐이다. 써모피셔의 도움을 기대하였으나 그건 오산이었다. 시원한 대답을 몇 번 못 들어봤다. 선택은 여러분의 몫이다.)

필요한 일부의 SNP만을 출력하였는데 출력하고 보니 일루미나의 SNP 이름과 Axiom의 SNP이 이름이 서로 달랐다. 이름, Chromsome, Position이 같아야 imputation을 위한 snp_info 파일을 만들기 편하다. 그래서 여기서는 SNP를 다운로드할 때 Axiom의 SNP 이름이 아니라 일루미나의 SNP 이름으로 다운로드하는 방법에 대해서 알아본다. Axiom의 어느 SNP 이름을 일루미나의 어느 SNP 이름으로 바꿀지 미리 준비하여야 하는데 그것은 이전 글을 참고하기 바란다. 이전 글에서 준비한 엑셀파일에서 Axiom SNP 이름, 일루미나 SNP 이름 순으로 준비한다.

Long Format Export Tool에는 재밌는 기능이 있다. 내가 원하는 일부 SNP만 다운로드하는 기능이다. 왜 일부 SNP만 다운로드할 필요가 있을까? 예를 들어 일루미나의 BovineSNP50v3와 써모피셔의 Axiom_BovMDv3.r2를 같이 분석하는데 일루미나의 BovineSNP50v3로 imputation해서 자료 분석을 한다고 가정해 보자. imputation을 하면 기준칩과 공통되는 snp만 이용하여 나머지 snp를 imputation하게 된다.(물론 양쪽으로 imputation 하는 경우도 있는 것 같지만 여기서는 한 쪽으로 imputation하는 것을 가정한다.) 즉 Axiom_BovMDv3.r2의 모든 snp가 필요한 것이 아니라 BovineSNP50v3와 공통된 snp만 필요하다는 뜻이다. 이런 의미에서 좀 더 빠르게 다운로드하고 자료 처리량을 줄이려면 일부 snp만 다운로드하는 것도 필요하게 된다. 그러면 일부 snp라고 했는데 일부 snp를 고르는 방법은? 그건 다음 SAS 프로그램으로 처리한다.

Illumina BovineSNP50v3 SNP Map

Thermofisher Axiom_BovMDv3.r2 SNP Map

/*
Illumina Bovine v3 칩과 Axiom BovMDv3.r2 칩의 공통 SNP 찾기
 - 1에서 29번 염색체 사이에서
 - chromosome, position 중복 SNP 제외 등
*/

/* Illumina Bovine v3 snp map 읽기 */
data a1    ;
    infile 'SNP_Map_ibv3.txt' delimiter='09'x MISSOVER DSD lrecl=32767 firstobs=2 ;
    
    informat Index_ibv3 best32. ;
    informat Name_ibv3 $40. ;
    informat chr $2. ;
    informat pos best32. ;
    informat GenTrain_Score best32. ;
    informat SNP $5. ;
    informat ILMN_Strand $3. ;
    informat Customer_Strand $3. ;
    informat NormID best32. ;
    
    input
             Index_ibv3
             Name_ibv3  $
             chr $
             pos
             GenTrain_Score
             SNP  $
             ILMN_Strand  $
             Customer_Strand  $
             NormID
    ;
run;
/*
no of data : 53218
*/

proc freq data = a1;
    table chr;
run;


/* 1에서 29 염색체 추출, 문자 -> 숫자 변환  */
data a2;
    set a1;

    if chr = '0' or chr = 'X' or chr = 'Y' or chr = 'MT' then delete;
    
    chr2 = input(chr, 2.);

    keep name_ibv3 chr2 pos;
run;
/*
no of data : 51278
*/

/* 염색체 1에서 29 추출, 문자 -> 숫자 변환 후 빈도 */
proc freq data = a2;
    table chr2;
run;

/* chromosome, position이 중복되는 snp 추출 */
proc sort data = a2;
    by chr2 pos;
run;

data a2_e1;
    set a2;
    if lag(chr2) = chr2 and lag(pos) = pos then output;
    keep chr2 pos;
run;
/*
no of data : 29
*/

/* chromosome, position이 중복되는 snp 쌍 제거 */
data a3;
    merge a2 (in = x1) a2_e1 (in = x2) ;
    by chr2 pos;
    if x2 then delete;
run;
/*
no of data : 51220
*/

/* 염색체 별 SNP 빈도 */
proc freq data = a3;
    table chr2;
run;

/* Axiom BovMDv3.r2 SNP Map 읽기 */
data b1    ;
    infile 'output_01_Axiom_BovMDv3.r2_SNPMap.txt' delimiter='09'x MISSOVER DSD lrecl=32767 firstobs=2 ;
 
    informat Index $1. ;
    informat Name_abv3 $12. ;
    informat chr $13. ;
    informat pos best32. ;
    informat GenTrain_Score $1. ;
    informat SNP $5. ;
    informat ILMN_Strand $1. ;
    informat Customer_Strand $1. ;
    informat NormID $1. ;
 
    input
             Index  $
             Name_abv3  $
             chr  $
             pos
             GenTrain_Score  $
             SNP  $
             ILMN_Strand  $
             Customer_Strand  $
             NormID  $
    ;
run;
/*
no of data : 63655
*/

proc freq data = b1;
    table chr;
run;
/* 다양한 염색체 이름 때문에 당황 */

/* 에러가 나지만 문자 -> 숫자 변환. 숫자로 변환된 것 중 염색체 1 ~ 29만 추출 */
data b2;
    set b1;
    chr2 = input(chr, 2.);
    if chr2 = .  then delete;
    keep name_abv3 chr2 pos;
run;
/*
no of data : 60598
*/

proc freq data = b2;
   table chr2;
run;

/* chromosome, position이 동일한 snp 추출 */
proc sort data = b2;
    by chr2 pos;
run;

data b2_e1;
    set b2;
    if lag(chr2) = chr2 and lag(pos) = pos then output;
    keep chr2 pos;
run;
/*
no of data : 29
*/

/* unique chromosome, position 추출 */
data b2_e2;
   set b2_e1;
   by chr2 pos;
   if first.pos then output;
run;
/*
no of data : 28
*/

/* chromosome, position이 동일한 snp 쌍 제거  */
data b3;
    merge b2 (in = x1) b2_e2 (in = x2);
    by chr2 pos;
    if x2 then delete;
run;
/*
no of data : 60541
*/

/* illumina bovine v3와 axiom bovine md v3.r2 snp map 머지 */
proc sort data = a3;
    by chr2 pos;
run;

proc sort data = b3;
    by chr2 pos;
run;

data c1;
    merge a3 b3;
    by chr2 pos;
    if name_ibv3 ~= '' and name_abv3 ~= '' then output;
run;
/*
no of data : 40422
*/

/* 엑셀로 출력 */
PROC EXPORT DATA = C1 
            OUTFILE = "common_snp_between_ibv3_and_abv3r2.xlsx" 
            DBMS =EXCEL REPLACE;
     SHEET="common_snp"; 
RUN;

위 SAS 프로그램에서 출력한 엑셀 자료를 이용한다.

드디어 Final Report를 다운로드하는 단계에 왔다. Final Report를 다운로드하는데 반드시 Export Complete를 확인해야 한다. SNP 6만개, 30두 정도는 몇 분 만에 끝나지만 SNP 70만개, 5000두는 며칠이 걸릴지 모른다. 며칠이 걸려도 끝나기만 하면 좋겠는데 너무 용량이 커서 끝나지 못하는 경우도 있다. PC의 메모리 영향을 받는 것 같다. 경험상 16G 램에서 SNP 40만개라면 1000두 정도, SNP 6만개라면 5000두 정도를 처리하는 것으로 보인다. 행의 수로는 3 ~ 4억 정도의 행을 출력하는 것으로 보인다. 자신의 컴퓨터에 알맞은 SNP 수, 샘플 두수를 정하기를 바란다. 그리고 해당 컴퓨터에서 다른 일을 하고 있다면 다운로드할 수 있는 행의 수는 아마도 더 줄어들 수도 있을 것이다.

snp map 파일과 final report가 만들어 졌으면 잘 만들어 졌는지 확인해야 하는데 snp map 파일이야 5만에서 70만 행이므로 열어볼 수 있겠지만 final report가 3억 개의 행을 가지고 있다면 열어 보기가 쉽지 않다. 그럴 경우 Unix 명령어 head, tail, wc 등을 이용하여 확인하여야 한다.

다음은 head, tail  명령어로 파일의 앞 뒤를 확인하는 화면이다.

다음은 wc 명령어로 행의 수를 확인하는 그림이다.

SNP Map 파일의 행의 수는 63656이었다. header를 제외하면 SNP 개수는 63655개임을 알 수 있다. Final Report의 행의 수는 1,909,661개인데 이것은 30두 x 63655 SNP 개수 + 11 행 header = 1,909,661이다. 잘 다운로드된 것을 확인할 수 있다. 경험에서 우러나오는 얘기를 하자면 다운로드된 파일을 믿으면 안된다. 일루미나의 지놈스튜이도, 써모피셔의 Axiom Analysis Suite도 오류를 일으킨다. 자기사 실수했다고 실토하면 좋겠지만 그렇지 않은 경우가 많다. 결국 다운로드가 잘 되었는지 안 되었는지 확인해야 하는 것은 사용자의 몫이다. 항상 하는 얘기지만 위와 같은 확인을 하는 것이 정신 건강상 좋고, 결국 시간을 절약하는 길이다.

(위 화면은 우분투가 설치된 윈도우즈 WSL2 환경에서 명령어를 입력한 화면이다. 윈도우즈에 WSL2를 설치하면 유닉스-리눅스 명령어를 사용할 수 있다. 무료로 설치할 수 있다. 예전에는 WSL2 설치가 복잡했는데 요즘에는 무척이나 간단해 졌다. 윈도우즈에서 리눅스 명령어를 사용하기 싫다면 리눅스가 설치된 컴퓨터로 파일을 보내서 확인해야 한다.)

지금까지 Genotype Data를 다운로드하는 형태는 기본적으로 사각형 방식이었다. 기본적으로 열(또는 행)에 개체가 있고, 행(또는 열)에 SNP가 있는 형태였다. 그런데 일루미나 사의 지놈스튜디오는 한 개체의 한 SNP가 한 행을 차지하는 출력형태를 지원했다. 지놈스튜디오에서는 Final Report라 불려졌다.아마도 Axiom Analysis Suite는 이러한 Final Report 형태의 출력을 처음부터 지원한 것은 아닌 거 같다. 그러니 Axiom Analysis Suite 프로그램에 Final Report 출력 기능이 포함된 것이 아니라 Axiom Long Export Tool 이라는 외부 프로그램에 이런 기능이 있는 것으로 보인다. 그런데 이렇게 되니 불편한 점이 생겼다. 이미 Axiom Analysis Suite에 Alternate Sample Name을 넣어 놨는데, Axiom Long Export Tool에서는 이것을 이용하지 못해 다시 템플릿을 다시 만들어 주어야 한다. 짜증이 나지만 어쩔 수 없으니 Axiom Long Export Tool에서 사용할 개체 이름 템플릿을 만들어 보자.

(cel 파일 이름에 개체 이름이 있을 경우 개체 이름을 추출하는 방법에 대해서는 맨 아래에 있는 링크 참조)

Axiom Analysis Suite에 개체 이름 입력하기, cel 파일 이름에서 개체 이름 추출하기

2024.10.31 - [Animal Breeding/Axiom Microarray Genotyping] - 개체 이름으로 Genotype 다운로드 하기