[유전체] Imputation

마이크로 어레이를 이용한 genotyped 데이터는 기술적 한계나 depth 등으로 SNP의 결측치가 발생한다.

이러한 데이터를 보완하기 위해 결측치를 보정하는 것을 imputation이라고 하고, 

impute2, shapeit, minimac, MaCH, beagle 등 많은 보정 프로그램이 있다. 

 

결측치를 확인하고 보완하는 과정은

계산의 반복을 줄이고 효율적으로 진행하기 위해 보통 pre-phase와 phase 단계로 나누어 진행하는 것을 추천하고 있다. 

 

Imputation을 하기 위해, 

우선 해당 array data의 QC를 하고,

AffyMetrix chip인 경우에는 RS ID를 바꾸는 과정이 필요했다.

 

이후 염색체나 위치 정보를 바꿔주는데, affymetrix 사에서 제공하는 annotation 파일을 사용하거나

필요한 경우 liftover(https://genome.ucsc.edu/cgi-bin/hgLiftOver)를 이용해 build 버전을 GRCh/hg를 맞춰준다.

그리고 1000genome이나 HapMap panel을 이용할 경우, 보두 forward(+) strand로 맞춰줘야 하기 때문에,

이 또한 annotation file을 이용해서 strand 목록을 만들어 strand 정보를 교정한다.

이는 Plink에서

--flip <strand list.txt>

옵션을 사용해 할 수 있다.

strand list는 'RS ID' 'RS _ID'를 tab을 구분자로 하여 작성한다. 

(annotation에서 strand가 '-'로 되어 있는 SNP를 추출하여 만들었음)   

 

 

혼자 시작하는 분석이라 시행착오도 많고 여기저기 검색과 논문에 의존하여 진행하고 있다. 

참고하는 사이트는 어래와 같다. 

plink  zzz.bwh.harvard.edu/plink/

shapeit  mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html

biostars  www.biostars.org/ 

 

이렇게 데이터가 준비되면

phasing을 하기 위해 shapeit2를 서버(우분투 16)에 다운 받는다.

shapeit은 염색체별로 쪼개지 않으면 position error가 쉽게 나는 것 같다. 

1KG 데이터를 통째로 돌렸더니 previous position과 current position의 duplicate site 에러가 났고, 

검색해 보니 염색체 별로 쪼개서 진행해야 한다고 한다. 

 

일단 리눅스 커맨드 창에서 plink loop을 실행한다.

리눅스에 익숙하지 않고, 바이오 데이터를 다루면서 거의 처음 사용하는데, 

sh 스크립트로 짜야하는 줄 알았는데 커맨드에서 바로 실행 가능하더라.

for chr in $(seq 1 23) ; \
do ./plink --noweb --file data/1KG/1KG --chr $chr --make-bed \
--out 1KG_${chr}; \done

이 과정을 실행하면 모든 염색체가 들어있는 1개의 파일에서 

염색체 1~22번까지의 bim, bed, fam 파일이 각각 생성된다. 

 

이렇게 생성된 파일을 input으로 하여,  

이번에는 각 염색체 별로 shapeit을 이용해 haps 파일과 sample 파일을 생성한다.

for chr in $(seq 1 22); do ./shapeit \
--input-bed 1KG_${chr}.bed 1KG_${chr}.bim 1KG_1KG_${chr}.fam \
--output-max result/1KG_${chr}.phased.haps result/1KG_${chr}.phased.sample \
--thread 8; done

이 과정에서 1시간 당 샘플 500개 정도 처리 시간이 걸린다. 

 

hapgen 옵션을 따로 선택하지 않으면 0.0045가 걸리는 듯.