일부 SNP에 대해서 Final Report 다운로드하기

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

SNP_Map_ibv3.txt

3.04MB

Thermofisher Axiom_BovMDv3.r2 SNP Map

output_01_Axiom_BovMDv3.r2_SNPMap.txt

2.20MB

 

/*
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 프로그램에서 출력한 엑셀 자료를 이용한다.