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
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 프로그램에서 출력한 엑셀 자료를 이용한다.
'Animal Breeding > Axiom Microarray Genotyping' 카테고리의 다른 글
| Final Report 다운로드하기 (2) | 2024.11.06 |
|---|---|
| Final Report 다운로드 하기(개체 이름 템플릿 작성) (1) | 2024.11.06 |
| Genotyping Data를 다운로드하기(PLINK) (2) | 2024.11.05 |
| Genotyping Data를 다운로드하기(VCF) (0) | 2024.11.05 |
| Genotyping Data를 다운로드하기(TXT - Call Codes) (2) | 2024.11.05 |
