앞선 포스팅의 두가지 과정을 거치고 나서 생성된 VCF 파일을 이용하면 드디어 분석 가능한 변이들을 확인할 수 있습니다. 그러나 실제로 이 데이터를 열어보면, 지저분하고 활용하기 위해서는 어느 정도 가공이 필요합니다. 그래서 이번 포스팅은 VCF Filter를 적용하여 분석을 위한 변이들을 정제하고, 분석에 참조하기 위한 Annotation 작업을 위한 Code까지 정리해보도록 하겠습니다.
관련 포스팅 보기>
NGS 분석 파이프 라인의 이해: GATK Best Practice
NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam
NGS DNA-seq pipeline: GATK Best Practice Code – Part2. Bam to Vcf
[일반적인 VCF 파일의 구조]
IV. Variant Filtration and Annotation : Genotype Refinement Workflow
VCF call을 하고 나면, 많은 변이 데이터가 생산되는데, 실제로 데이터의 퀄리티와 오류를 보정하는 작업이 필요합니다. 이 부분은 아직까지 확립된 최적의 Work Flow는 없지만, GATK에서는 Genotype Refinement Workflow를 제공하고 있습니다. 여기서는 Call된 genotype이 인구집단이나 가계도 정보에 근거하여 적절한 Call인지를 판단하게 됩니다. 가계도의 경우에는 양쪽 부모 정보가 모두 (Trio) 필요합니다. 아래에서는 이 Workflow를 따라 Code를 작성해보도록 하겠습니다.
1. CalculateGenotypePosteriors
VCF call을 통해 생성된 genotype을 기반으로 이번에는 거꾸로 Genotype의 posterior를 계산해줍니다. 이 결과를 토대로 변이 call을 판정해서, Genotype의 quality를 계산합니다. 이 과정에는 인구 집단이나 가계도의 Prior 정보를 이용하는데, 가계도 정보를 이용하기 위해서는 PED 파일 생성이 필요합니다.
PED 파일 생성 정보는 아래 링크를 참고하시기 바랍니다.
http://zzz.bwh.harvard.edu/plink/data.shtml
gatk --java-options 'Xmx16g' 'Xms8g' CalculateGenotypePosteriors -V [cohort.recal.vcf] \ -ped family.ped --skip-population-priors -O [cohort.GP.vcf]2. Filter Low Quality Genotypes
위의 과정에서 계산된 GQ 값을 기준으로 Filter 조건을 걸어서 지저분한 변이 call들을 제외 시킬 수 있습니다. 아직까지 표준화된 Filter 조건은 없기 때문에, 사용자가 적절하게 필터 조건을 설정해주어야 합니다. 여기서 선언해주는 조건은 흔히 JEXL expression이라고 부르는 방법으로 작성을 해야한다고 하는데, 여기서는 널리 쓰이는 Filter option인 GQ < 20 을 걸도록 하겠습니다. 필요에 따라 수치를 변경시키거나, 조건을 바꾸면 됩니다. 자세한 내용은 아래에서 찾아보시기 바랍니다.
article about using JEXL expressions
마지막으로 filter조건에 부합하는 변이만 남기기 위해, GATK의 SelectVariants 이용하도록 하겠습니다.
gatk --java-options 'Xmx16g' 'Xms8g' VariantFiltration -R [hg19_reference.fa] -V [cohort.GP.vcf] \ --genotype-filter-expression "GQ<20" --genotype-filter-name "lowGQ" -O [cohort.filter.vcf] gatk --java-options 'Xmx16g' 'Xms8g' SelectVariants -V [cohort.filter.vcf] --exclude-filtered true --max-nocoll-fraction 0.1 -O [cohort.filtered.vcf]3. VCF annotation using Annovar
이제 어느 정도 정제가 된 변이 데이터셋이 생성되었습니다. 이제 불러온 변이들에 맞는 annotation 작업을 진행하면 분석을 진행할 수가 있습니다. 최근에 GATK에서 Funcotator를 개발하여, annotation이 가능해졌는데 아직까지는 Annovar를 이용하는게 대세이므로, Annotation 작업은 Annovar를 이용하도록 하겠습니다.
관련 포스팅 보기>
Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용
Annovar 홈페이지를 통해 다운로드 신청을 하면, 등록된 메일로 Kai Wang에게서 자동으로 다운로드 url을 받을 수 있습니다. Annotation을 원하는 Database의 크기가 크기 때문에, 처음에는 DB 다운로드에 시간이 많이 소요가 될 수 있습니다.
Annovar를 설치하고, 아래 부위에 원하는 DB의 이름을 찾아서 넣으면 다운로드가 진행됩니다.
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar [refGene] humandbAnnotation은 아래 code를 통해서 진행합니다. 여기서는 refGene, cytoBand, avsnp150, dbnsfp35a, exac03, gnomad_exome, clinvar_20190305 DB를 이용하여 annotation 작업을 진행하도록 하겠습니다. 각각의 DB는 Gene-based, Region-based, Filter-based annotation에 해당하는 카테고리가 있고, 이를 operation 옵션에 g, r, f를 이용하여 나타내줍니다.
table_annovar.pl [cohort.filtered.vcf] humandb/ -buildver hg19 -out [myanno] -protocol refGene,cytoBand,avsnp150,dbnsfp35a,exac03,gnomad_exome,clinvar_20190305 -operation g,r,f,f,f,f,f -nastring . --vcfinput --remove --thread 8[References]
VCFtools를 이용한 VCF filtering: 위에서는 GATK를 이용하였지만, 과거에 많이 사용하였던 VCFtools를 이용하는 방법도 있습니다. 자세한 내용은 아래를 참조하시기 바랍니다.
http://www.ddocent.com/filtering/
GriffithLab에서 소개하고 있는 VCF filter 과정: 해당 내용은 GATK4로 넘어오면서 코드가 약간 달라졌습니다. 다만 전반적인 과정은 유사하므로, 참고하시기 바랍니다.
https://pmbio.org/module-04-germline/0004/02/02/Germline_SnvIndel_FilteringAnnotationReview/
두마디 정밀의료 