Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

many `LowQual` flags in output .vcf when using `HaplotypeCaller` with `--alleles` arg

luyitianluyitian ShenZhenMember

Hi GATK team,
I am working on a pipeline for exome sequencing variant calling. And I am only interested in the genotype for some specific positions so I used GENOTYPE_GIVEN_ALLELES mode with given vcf file. I only have "chromosome position ID REF ALT" info so I keeped other columns as . . I found that this mode is super slow, compared with discovery mode. And when it finished I found most of my records are marked with LowQual, but they have good sequencing depth and quality score. Can you tell me why they are marked with LowQual flag. Here are some output:

chr10 101421279 rs35229854 G A 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=80;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;SOR=0.407 GT:AD:DP:GQ:PL 0/0:80,0:80:99:0,241,2781
chr10 101421288 rs147747082 C A 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=81;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;SOR=0.332 GT:AD:DP:GQ:PL 0/0:81,0:81:99:0,244,2812
chr10 101421324 rs150267092 A G 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=64;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;SOR=0.173 GT:AD:DP:GQ:PL 0/0:64,0:64:99:0,193,2224
chr10 101421366 rs370286436 C T 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=59;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;SOR=0.101 GT:AD:DP:GQ:PL 0/0:59,0:59:99:0,178,2052
chr10 101421367 rs61729539 G A 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=58;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;SOR=0.105 GT:AD:DP:GQ:PL 0/0:58,0:58:99:0,175,2014
chr10 101451259 rs11190245 T C 605.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-0.278;ClippingRankSum=2.143;DB;DP=32;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.119;QD=18.93;ReadPosRankSum=1.904;SOR=0.507 GT:AD:DP:GQ:PL 0/1:11,21:32:99:634,0,286

I used default setting for other step. references are from gatk-bundle-2.8-hg19

Thanks for your help : )

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @luyitian
    Hi,

    Can you post the exact command you ran and let us know which version of GATK you are using? Can you also try running without GENOTYPE_GIVEN_ALLELES mode, and instead input the VCF with -L?

    It seems like the LowQual sites are the homozygous reference sites only. Can you confirm that is the case?

    Thanks,
    Sheila

  • luyitianluyitian ShenZhenMember

    Hi Sheila,

    I checked my result. All homozygous(0/0) are assigned as LowQual. And the majority of others(0/1, 1/1) are normal, just few of them are LowQual.

    the version info:

    The Genome Analysis Toolkit (GATK) v3.4-0-g7e26428, Compiled 2015/05/15 03:25:41

    this is the structure of my command line:

    "java -jar {GATK_path} -T HaplotypeCaller -R {genome} -L {bed_path} -D {snp_path} --alleles {known_snp} -nct 10 -I {input_path} --genotyping_mode GENOTYPE_GIVEN_ALLELES -stand_emit_conf 10 -stand_call_conf 30 -o {output_path}"

    the data I analysed is from exome sequencing so I added -L for exome interval.

    I tried to remove GENOTYPE_GIVEN_ALLELES and use -L to put my .vcf file. The LowQual disappears but I can only get heterozygote(0/1 or 1/1) genotype. The program does not report 0/0. But I hope they will report all genotypes of my interested SNP, includes 0/0. It is also the reason that I choose GENOTYPE_GIVEN_ALLELES mode with given SNP.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @luyitian
    Hi,

    The default HaplotypeCaller output is variant sites only.

    The best thing to do is run HaplotypeCaller with just -L "your VCF" in GVCF mode, then use GenotypeGVCFs with -allSites to get all the variant and non-variant sites in your final VCF.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Actually I think you can run this with HaplotypeCaller in normal mode and --allSites output instead of going through the GVCF workflow. You may need to add some interval padding to get the calls emitted correctly, using e.g. -ip 100. This will cause more sites to be output to your VCF than you really want, but you can easily select out the sites you want using SelectVariants with -L your_sites.vcf.

    This replaces the GENOTYPE_GIVEN_ALLELES functionality which seems to be a bit temperamental in HC.

  • luyitluyit ShenZhenMember

    @Geraldine_VdAuwera said:
    Actually I think you can run this with HaplotypeCaller in normal mode and --allSites output instead of going through the GVCF workflow. You may need to add some interval padding to get the calls emitted correctly, using e.g. -ip 100. This will cause more sites to be output to your VCF than you really want, but you can easily select out the sites you want using SelectVariants with -L your_sites.vcf.

    This replaces the GENOTYPE_GIVEN_ALLELES functionality which seems to be a bit temperamental in HC.

    still does not work. --allSites just report a few SNPs that overlap with my interest positions.

    in my older result with GENOTYPE_GIVEN_ALLELES arg, there are some homozygous with DP>100 but they disappear with --allSites args. So I am confused now. If they have DP>100 why they cannot be discovered by --allSites argument.

    My requirement is quite simple. There are some chromosome positions that I interested:
    chr1 10001
    chr1 20002
    chr2 3434

    I want to get genotype information at these positions, report all results.

    How can I use gatk to achieve this?

    Thanks for all your help and I wish you have a great Christmas holiday.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Can you post the exact commands that you ran?

  • luyitluyit ShenZhenMember

    @Geraldine_VdAuwera said:
    Can you post the exact commands that you ran?

    OK. I tried two mode. --output_mode EMIT_ALL_SITES and --genotyping_mode GENOTYPE_GIVEN_ALLELES with --alleles. these are my commands:
    normal mode and output all sites. There is no --allSites arg in HaplotypeCaller so I chose --output_mode EMIT_ALL_SITES

    java -jar /store/tianly/requirement/tools/GenomeAnalysisTK.jar -T HaplotypeCaller -R /store/tianly/requirement/genome/gatk-bundle-2.8-hg19/ucsc.hg19.fasta -D /store/tianly/requirement/genome/dbsnp_138.hg19.vcf -nct 10 --output_mode EMIT_ALL_SITES -I /store/tianly/autism_child/tmp/autism_child_bqsr.bam -stand_emit_conf 10 -stand_call_conf 30 -o /store/tianly/autism_child/autism_child_array_allsite.vcf

    this is the GENOTYPE_GIVEN_ALLELES. All the homozygous are marked as LowQual as I mentioned before.

    java -jar /store/tianly/requirement/tools/GenomeAnalysisTK.jar -T HaplotypeCaller -R /store/tianly/requirement/genome/gatk-bundle-2.8-hg19/ucsc.hg19.fasta -D /store/tianly/requirement/genome/dbsnp_138.hg19.vcf --alleles /store/tianly/requirement/snp_db/known_snp.vcf -nct 10 -I /store/tianly/autism_child/tmp/autism_child_bqsr.bam --genotyping_mode GENOTYPE_GIVEN_ALLELES -stand_emit_conf 10 -stand_call_conf 30 -o /store/tianly/autism_child/autism_child_array.vcf

    ../snp_db/known_snp.vcf contains some GWAS SNPs that I interested. The version of my GATK is 3.4-0-g7e26428.

    Besides the large number of homozygous that is not reported by EMIT_ALL_SITES, I also found some discripency between these two mode:

    a few heterozygous are reported by EMIT_ALL_SITES but not normal mode with EMIT_ALL_SITES

    in output .vcf
    chr3 128199918 rs200597976 A G 98.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=2.997;ClippingRankSum=-1.468;DB;DP=46;FS=7.832;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=-1.009;QD=2.15;ReadPosRankSum=-1.743;SOR=2.694 GT:AD:DP:GQ:PL 0/1:39,7:46:99:127,0,1236

    mpileup result:
    chr3 128199918 A 56 GGG.G.,TG,.GG..GG.,..,.........,,.,...,.,......,..,...,^]. ;9BBBBeBBbBGBBBBBGgTXgT\Z\QZZZ\ccX[WX_b`cafeeecccccc^RRE

    I guess such sites get filtered out by normal mode because of strand bias.

    Some discrepancy is arisen because GENOTYPE_GIVEN_ALLELES only checks the genotype in the given .vcf files.

    position: chr4 187004217
    known_snp.vcf : chr4 187004217 rs3775290 C A . . .
    output by normal mode with EMIT_ALL_SITES: 'C', 'T', '1/1'
    output by GENOTYPE_GIVEN_ALLELES mode: 'C', 'A', '0/0'
    mpileup result:
    chr4 187004217 C 51 ttTtTTtTtT.TTtTtTTtTtTTTtTTTtTTTTtTTTTTTTTTTttT.TTT 8ceebgdih_eiiiiihiifhiiiiiiihiiiihhiiii[hiiieeY^RBB

    So it seems GENOTYPE_GIVEN_ALLELES won't perform any filter and sometimes report false result.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @luyit
    Hi,

    Can you try running HaplotypeCaller in GVCF mode and then running GenotypeGVCFs like I suggested above? I think HaplotypeCaller's EMIT_ALL_SITES mode does not work properly. You can add -ip 100 to your command as well.

    -Sheila

  • aklakl Member

    I Ran HaplotypeCaller in GVCF and then running GenotypeGVCF did not help for indel. For example, when I run

    java -jar "E:\GenomeAnalysisTK.jar" -T HaplotypeCaller -R E:\genome.fa -I ML4delEX1_EX7-31031_S2.bam -o ERC/ ML4delEX1_EX7-31031_S2.bam.alleles.g.vcf --downsampling_type NONE -ERC GVCF --output_mode EMIT_ALL_SITES -L E:\ROI.V1.Sorted.vcf.gz --alleles E:\ROI.V1.Sorted.vcf.gz
    Got the following records
    chr1 11863157 . C . . END=11863158 GT:DP:GQ:MIN_DP:PL 0/0:2380:99:2356:0,120,1800

    My ROI.V1.Sorted.vcf.gz has
    chr1 11863157 . CT C . PASS None

    Then running GenotypeGVCF, I did get that site 11863157 but I also got 11863158
    chr1 11863157 . CT . 145.23 . AN=2;DP=2356;MLEAC=.;MLEAF=. GT:DP:RGQ 0/0:2380:99
    chr1 11863158 . T 148.23 . AC=0;AF=0.00;AN=2;DP=2356;MLEAC=0;MLEAF=0.00 GT:AD:DP:RGQ 0/0:2356,0:2356:120

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @akl
    Hi,

    Setting stand_call_conf 0 will allow you to be more sensitive, but you will also decrease specificity.

    Are you running GenotypeGVCFs with just the one sample or with multiple samples? Also, can you try removing --downsampling_type NONE --output_mode EMIT_ALL_SITES --alleles E:\ROI.V1.Sorted.vcf.gz from your HaplotypeCaller command?

    Thanks,
    Sheila

Sign In or Register to comment.