Fork me on GitHub 단맛만좋아요 :: Samtools를 이용한 genotype likelihoods 구하기

Samtools를 이용한 genotype likelihoods 구하기
2014/10/31 13:32 | 바이오인포매틱스
samtools가 버전 0.1.19를 마지막으로 major 번호가 올라갔습니다. 바로 1.0 버전대가 탄생한것이죠. 홈페이지도 이제는 www.htslib.org를 사용합니다. 흔히 우리가 말하는 Samtools는 Samtools, BCFtools, HTSlib 3개로 구성되어 있습니다.  

Samtools는 SAM/BAM/CRAM 포맷의 파일을 읽고/쓰고/편집하고/인덱싱하고/볼 수 있는 툴입니다. BCFtools는 BCF2/VCF/gVCF 포맷의 파일을 읽고/쓸 수 있으며 SNP나 short indel의 sequence variants를 calling/filtering/summarising 할 수 있는 툴입니다. 그리고 그 기반은 high-throughput sequencing 데이터를 다루는 바로  C로 작성된 HTSlib이라는 라이브러리를 기반으로 하고 있습니다. Samtools를 이용한 Variant Calling에서도 다음과 같은 명령으로 변경되었습니다. 
samtools mpileup -ugf ref.fa sample1.bam sample2.bam sample3.bam | bcftools call -vmO z -o study.vcf.gz
잠깐 살펴보면 예전에는 bcftools의 view가 call이라는 명령어로 변경된 것을 확인 할 수 있을 겁니다. 그럼 간단하게 Samtools를 이용한 variants를 identifying하는 방법에 대해서 알아보도록 하겠습니다. mpipleup 명령은 aligned read 즉 BAM 포맷의 파일로 부터 genotype likelihood를 계산하는 명령입니다. 해당 샘플의 모든 site에 대해서 genotype likelihood를 계산하게 됩니다.

The Variant Call Format (VCF) Version 4.1 Specification을 보면 INFO 필드의 몇몇 키워드에 대한 설명이 나오는데요. GL 필드와 이와 관련한 PL, GP 필드에 대해서 우선 알아보겠습니다. GL 필드는 앞선 genotype likelihood에 대한 정보를 담고 있는 필드 (Floats)로 genotype likelihoods comprised of comma separated floating point log10-scaled likelihood for all possible genotypes given the set of alleles defined in the REF and ALT fields.

예를 들어 아래와 같이 reference allele가 T라고 할때 alternative allele가 C와 G 두 개가 있다고 합시다 (T=0, C=1, G=2) . 그리고 PL 필드를 보면 콤마(,)로 구분된 총 6개의 값이 보입니다. diploid인 경우 AA, AB, BB, AC, BC, CC 6개의 genotype에 대해서 likelihood 값을 계산하게 됩니다. REF는 0 , ALT는 1,2로 표현되는 경우 (0,0), (0,1), (1,1), (0,2), (1,2), (2,2)로 표시할 수 있겠습니다. 그럼 여기서 TG genotype은 몇 번째에 있는 값인가요?라고 물어보면 F(j,k) = (k*(k+1)/2)+j 즉, F(0, 2) = 3 이므로 'TG'는 genotype likelihood는 4번째에 존재하는 99가 genotype likelihood 값이 됩니다.
chr1 12887054 . T C,G 194 . DP=253;VDB=0.423501;SGB=-0.693147;RPB=2.7979e-06;MQB=2.50009e-15;MQSB=3.79278e-05;BQB=0.860152;MQ0F=0.41502;AC=1,1;AN=2;DP4=60,0,147,20;MQ=15 GT:PL 1/2: 255,100,99,213,0,231
 Genotype Phread Scaled Score   Likelihood p-value
AA, (0,0), TT 255
AB, (0,1), TC 100
BB, (1,1), CC 99
AC, (0,2), TG 213  
BC, (1,2), CG 0  
CC, (2,2), GG  231  




glf 포맷,,, glfsingle, samtools-hybrid,,, 등등등 추가해야함,,, 우선 skip 추후에 지속 업데이트하기로 한다. 

아이고 의미없다.  
저작자 표시 비영리 동일 조건 변경 허락
hongiiv (Changbum Hong) is the software enginner of GenomeCloud. He covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014/10/31 13:32
Currently 댓글이 없습니다. comments want to say something now?

티스토리 툴바