KoGES imputation (1)

 

한국인 유전체 역학조사사업 (KoGES) genotype data를 imputation하는 과정입니다.

 

KoGES는 KARE (안산안성), CAVAS, HEXA cohort로 구성되어 있습니다.

 

가장 처음 genotype 된 KARE는 8840명이고, 이것은 1000 G로 imputation된 데이터를 같이 줍니다.

 

(정확히는 같이 신청할 수 있습니다)

 

그리고 KCHIP ver 1.0으로 만들어진 7만여명의 데이터는 5000여명의 KARE, 그리고 HEXA 5만여명, CAVAS 만여명으로

 

구성되어 있습니다.

 

처음에 받으면 map, ped 파일로 되어있습니다.

 

그런데 1000 genome으로 imputation 된 데이터 imputation이 굉장히 오래되어서 좀 이상한 부분이 많고

 

KCHIP은 imputation이 되어있지 않다는 것이 문제입니다.

 

 

이것을 imputation 하는 방법입니다.

 

 

우선 AXID를 rsID로 변경해야 합니다.

 

https://nephro.tistory.com/9

Kchip의 AX-ID -> RS-ID로 변경

 

 

참조하셔서 변경하세요.

 

그리고 map, ped를 binary 형식은 bed, bim, ped로 변경합니다. (plink는 사용하실 줄 안다고 가정하였습니다.)

 

Kchip_RSID.bed, bim, ped 라고 이름 정하겠습니다.

 

원래 genotype 데이터는 QC를 해야하는데, 찾아보니 이미 QC한 자료를 배포한 것이라고 해서

 

QC도 넘어가도록 하겠습니다

 

 

imputation은 예전에는 impute2가 gold standard였는데 지금은 너무 커져서 도저히 impute2로 할 수가 없습니다.

 

그래서 UKBB도 Jonathan Marchini 선생님(impute 만드신 분)이 impute2를 recoding해서 imupe4를 만들어서 했다고

 

되어있습니다.

 

 

impute5로 오면서, 너무 느린 문제를 해결하기 위하여, 알고리즘을 바꾸었다고 해서

 

저는 그냥 beagle을 사용하기로 하였습니다.

 

왜냐하면 imputation 전에 prephasing을 해야하는데, impute는 shapeit이라고 하는 프로그램을 써야 합니다만

 

beagle은 prephasing, imputation 모두 하나 프로그램으로 가능하고

 

또 속도도 빠른 것 같습니다.

 

http://faculty.washington.edu/browning/beagle/beagle.html

Beagle 5.4

 

 

 

우선 beagle은 human genetic map, reference file까지 모두 만들어서 올려놓았기 때문에

 

사용하기 간편한 것 같습니다.

 

리눅스가 필요하고, 어느 정도 사용법에 익숙하셔야 합니다.

 

(1) 

#!/bin/bash

## name of file

FILE_NAME=Kchip_RSID

 

우선 bash임을 선언하고,

 

FILE_NAME을 변수로 만듭니다.

 

(2) 다음으로 chromosome별로 파일을 sepration 해야합니다.

## (1) seperation by chromosome

for chr in {1..22}; do

  plink --bfile ${FILE_NAME} --chr $chr --make-bed --recode vcf --out prephased/${FILE_NAME}_chr${chr} ;

done

 

for 문을 이용해서

 

chromosome 1부터 22까지 쪼개줍니다.

 

vcf 파일로 만들어야 해서, --recode vcf 명령어를 사용하였습니다.

 

 

(3) 다음으로 prephasing 하기 전에 conformation이라는 과정을 거쳐야 합니다.

 

chromosome은 2개이니, 이 유전 block이 어디에 붙어야 하는지 결정하는 과정 같습니다.

 

그냥 필요하다 하고 너머가시면 될 것 같습니다.

## (2) conform

BEAGLE=/home/user_name/bin/beagle5.4/

REFERENCE=/home/user_name/driveD/ref_gene/1K_phase3_b37_vcf/EAS/

THREAD=18

for chr in {1..22}; do \

  java -jar ${BEAGLE}conform-gt.24May16.cee.jar \

  gt=prephased/${FILE_NAME}_chr${chr}.vcf \

  ref=${REFERENCE}chr${chr}.1kg.phase3.v5a.EAS.vcf.gz \

  chrom=${chr} \

  out=prephased/${FILE_NAME}_chr${chr}_conformed;

done

 

BEAGLE 프로그램 위치를 저장해주고

REFERENCE 위치도 저장해줍니다.

 

위에 beagle 홈페이지에 가시면 vcf 파일로 만들어져 있습니다.

 

1K genome reference이고, east asian subpopulation 503명만 이용한 것입니다.

 

THREAD는 사용할 process 수입니다. 저는 20개짜리라, 2개는 다른 일 할 때 쓰려고 (R같은...)

 

남겨두었습니다.

 

 

for 문을 이용해서 실행되는데

 

gt는 위에서 만들어진 sepration된 vcf 파일 넣어주면 됩니다.

 

 

(4) 다음으로 phasing 입니다.

 

## (3) phasing

MAP=/home/user_name/driveD/ref_gene/1K_phase3_b37_bref3/

for chr in {1..22}; do \

  java -Xmx500g -jar ${BEAGLE}beagle.05May22.33a.jar \

  gt=prephased/${FILE_NAME}_chr${chr}_conformed.vcf.gz \

  ref=${REFERENCE}chr${chr}.1kg.phase3.v5a.EAS.vcf.gz \

  map=${MAP}plink.chr${chr}.GRCh37.map \

  chrom=${chr} \

  out=prephased/${FILE_NAME}_chr${chr}_1K_beagle_b37_EAS \

  impute=false \

  nthreads=${THREAD};

done

 

여기서는 MAP이라는 파일들이 사용되는데

 

역시나 beagle 홈페이지에 다 만들어져 있습니다.

 

 

(5) 다음으로 imputation 직접 실행입니다.

## (4) imputation

for chr in {1..22}; do \

  java -Xmx500g -jar ${BEAGLE}beagle.05May22.33a.jar \

  gt=prephased/${FILE_NAME}_chr${chr}_1K_beagle_b37_EAS.vcf.gz \

  ref=${REFERENCE}chr${chr}.1kg.phase3.v5a.EAS.vcf.gz \

  map=${MAP}plink.chr${chr}.GRCh37.map \

  chrom=${chr} \

  out=imputed/${FILE_NAME}_chr${chr}_1K_imputed_beagle_b37_EAS \

  impute=true \

  window=40.0 \

  nthreads=${THREAD};

done

 

phasing과 차이라면, impute=true로 되어있는 것, 그리고 window 입니다.

 

얼마 단위로 imputation 할 거냐라는 이야기인데, 저는 이걸 오래전에 setting해놔서

 

왜 40이었는지는 가물가물 합니다. 

 

아마 저자들의 recommendation이었던 것 같습니다.

 

 

이걸 한 bash 파일에 집어넣어서 실행하시면 되겠습니다.

 

이제 실행시키고 나면, 1~2주 쯤 결과가 나와있을 겁니다.

 

 

그리고 이제 나온 결과를 이용해서, filtering을 해야합니다.

 

filtering은 다음 기회에...