잃어버린 부모 찾기 - VCF 파일을 PLINK 포맷으로 변경하기

즐거운 분석 놀이 - 역시나 제일 재미있는 것은 뭐랄까? 데이터를 분석하고 거기에서 의미를 찾아내는 것이 아닐까?   ^_____________^   마치 보물찾기와 같은...그럼 오늘은 저번 galaxy를 이용한 분석의 연장선상으로 발굴한 유전변이를 가지고 다양한 분석을  해보도록 하겠습니다.

그런데 안타깝게도 오늘 사용할 주재료인 VCF 파일은 galaxy에서 아직 지원을 안하고 있습니다. VCF (Variant Call Format)는 유전변이를 나타내는 표준 파일 형식인데요. NGS를 통해 발굴된 유전변이를 VCF 포맷으로 만들어야 하건만, 이때 Picard나 GATK와 같은 툴들이 사용되는데, 아직 galaxy에서는 이 툴들을 지원하지 않고 있습니다. 뭐 조만간 지원할 것이라고 보는데요. 그건 테스트 galaxy 서버에는 해당 툴들이 있기 때문입니다. 그러니까 곧 정식 galaxy에도 도입된다는 이야기가 되겠죠. 암튼 NGS를 통해 나온 유전변이를 VCF 포맷의 파일로 만들었다고 가정하고 진행하도록 하겠습니다.

한국인 데이터가 없다구요?!!
현재 국내에서 가장 많은 한국인 NGS 데이터를 공개하고 있는 곳은 KPGP(Korean Personal Genome Project)입니다. 이 프로젝트에서는 현재 42명의 한국인 데이터를 공개하고 있습니다.  

	남	여	합
한국인	20	16	36
한국인 혼혈 (East Asian/Caucasian)	2	0	2
한국인 일란성 쌍둥이	0	2	2
한국인 이란성 쌍둥이	2	0	2
합계	24	18	42

일란성 쌍둥이 데이터 다운로드
KPGP 데이터는 opengenome.net 을 통해  다운로드가 가능하며, 공개되는 데이터는 FASTQ 포맷의 raw data, BAM 포맷의 mapping 데이터, 각종  유전변이 데이터가 존재합니다. 여기서는 VCF 포맷의 SNV 데이터를 다운로드 하도록 하겠습니다. gzip으로 압축된 형태로 약 80 MB의 용량 (압축 해제시 약 400 MB)을 차지합니다. KPGP 웹 페이지를 보시면 일란성 쌍둥이 (Monozygotic Twin) 샘플이 있습니다.  일란성 쌍둥이의 경우 SNV 데이터는 정확히 일치한다고 알려져 있습니다. 물론 그래야만 하구요 ^^;; KPGP_00088, KPGP_00089가 일란성 쌍둥이 샘플명입니다. 각각을 클릭하여 SNV 데이터를 다운로드 합니다.

VCF 파일을 PLINK 포맷으로 변경하기
다양한 genome 분석을 위해서는 VCF 포맷의 데이터를 PLINK 포맷으로 변경해야 합니다. 이때 사용되는 툴이 VCFTools 입니다. VCFTools의 기능 중에는 PLINK에서 사용 가능한 PED, MAP 파일로 변경해주는 옵션이 있습니다.

$ vcftools --vcf KPGP89_G_110915_HiSeq_EastAsian_Kor_F.vcf.SNV --plink

이렇게 명령을 주면, out.map 파일과 out.ped 파일이 생성됩니다. 그럼 이 파일을 PLINK를 이용하여 Binary 형태로 변경합니다. 그런데 주의하셔야 할것은 기본적으로 KPGP에서 생성한 VCF 파일에는 염색체를 나타내는 첫번째 컬럼이 "chr1"과 같이 chr이라는 접두사가 붙습니다. 이렇게 되면 PLINK에서 인식하지 못하기 때문에 그냥 "1"이라고 변경해주어야 합니다. vi 등에서 문자열 치환을 통해 "chr"을 모두 삭제합니다. 또 하나 주의하여야 할것이 있는데요. KPGP의 VCF 포맷에는 ID 즉 해당 변이의 식별자에 해당하는 세번째 컬럼 또한 모두 "."으로 되어 있습니다. 이 부분 또한 임의의 ID를 부여해야합니다. 저는 "염색체_포지션" 즉 "chr1_112000"의 형태로 임의의 ID를 부여했습니다. 이상의 두 가지에 대해서 변경을 하신 후에 VCFtools를 이용하여 PLINK포맷으로 변경합니다. 

$ plink --file out --make-bed --out KPGP_89

위의 명령을 실행하면, out.map과 out.ped 파일은 KPGP_89.bed, KPGP_89.bim, KPGP.fam의 세개의 파일로 변경됩니다. 이 과정을 두 쌍둥이에 대해서 각각 수행합니다. 이제는 각각의 PLINK 파일을 병합합니다.

$ plink --bfile KPGP_89 --bmerge  KPGP_88.bed KPGP_88.bim KPGP_88.fam --make-bed --out TWINS

그럼 TWINS.bed, TWINS.bim, TWINS.fam 파일에 두명의 쌍둥이 데이터가 들어가게 됩니다. 이러한 방법으로 VCF파일을 PLINK로 변환하여 association 분석 등을 수행할 수 있습니다.

혈연관계 분석하기
그럼 과연 이 쌍둥이 데이터가 정말 혈연관계가 있는건지 확인해 보도록 하겠습니다. PLINK에서는 IBS를 이용하여 혈연관계를 분석할 수 있습니다만, 여기서는 KING(Kinship-based INference for Gwas) 이라는 프로그램으로 혈연관계를 분석하도록 하겠습니다. KING은 pairwise로 사람들간의 혈연지수 (kinship coefficient)를 계산해줍니다. kinship coefficient는 다음과 같습니다.

Relationship	Kinship cofficient	Coeffcient of relatedness
자기자신과 비교, SELF	0.5	1.0
일란성 쌍둥이, Monozygotic twins	0.5	1.0
부모-자식, Parent-child	0.25	0.5
친형제, Full siblings	0.25	0.5
배다른 형제, Half siblings	0.125	0.25
친사촌, First cousins	0.0625	0.125
관계없음, Unrelated	0	0

$ king -b TWINS.bed --kinship

위의 결과로 king.kin0 파일이 생성되며 여기에 비교하고자 하는 샘플들간의 Kinship cofficient가 있습니다. 아래와 같이 두 쌍둥이간에 0.4988이 나온것을 확인하실 수 있는데요. 바로 일란성 쌍둥이이거나 혹은 같은 사람이라는 것을 의미합니다. 

 
가족 찾기 문제
자 이제 문제 갑니다. 아까 말한 OPEN KPGP의 공개된 42명의 데이터에는 서로 가족인 사람이 있습니다. 부모와 자식2이 있습니다. 한번 분석해 보시면 간단하게 어떠한 샘플들이 서로 혈연 관계인지를 확인하 실 수 있을 겁니다. 힌트 나갑니다. 부(KPGP_00009) 모(KPGP_00010) 입니다. 자 그럼 여러분들이 잃어버린 자식을 한번 찾아 주시기 바랍니다. ^________________^ 

덧) KPGP 데이터를 다운로드 하실경우 olleh ucloud라고 되어 있는 곳의 RAW DATA를 다운로드 하실 경우 클라우드를 통해 빠르게 다운로드 하실 수 있습니다. 현재 KPGP데이터는 일반 FTP와 클라우드에 나뉘어 제공되니 참고 하시기 바랍니다.