Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.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.

GVCF from HaplotypeCaller using BAM with or without IndelRealignment+BQSR

mmokrejsmmokrejs Czech RepublicMember
edited May 2017 in Ask the GATK team

Hi,
I am studying why PDF figures from VQSR recalibration (https://software.broadinstitute.org/gatk/documentation/topic?name=methods#methods39) of my human exomes seem messy (http://gatkforums.broadinstitute.org/gatk/discussion/comment/38386/#Comment_38386). I went to use other BAM files of same samples for the same HaplotypeCaller step. Specifically, so far I used mysample.bwa.postalt.removed_duplicates.realignedtogether.BQSR.namesorted.fixmate.calmd.sorted.bam files as input but now I went for less derived BAM file in mysample.bwa.postalt.removed_duplicates.realignedtogether.bam. I wonder whether the content is actually recognized by VariantRecalibrator:

Here I create the GVCF file from either of the two input BAMs:

"$java_path" $java_opts -Xmx${ask_for_memory}g -jar "$gatk_binpath"GenomeAnalysisTK.jar -T HaplotypeCaller --genotyping_mode DISCOVERY --useNewAFCalculator --emitRefConfidence GVCF --pcr_indel_model "$pcr_indel_model" --sample_name "$sample" -R "$reference_flatfile" $ranges_arg -I "$sample_calmd_bam" --dbsnp "$dbsnp_file_Broad_style" -o "$sample"."$aligner"gatk.HaplotypeCaller.g.vcf

Here I create merge the 35 exome sample GVCFs:

$java_path $java_opts -Xmx${ask_for_memory}g -jar ${gatk_binpath}GenomeAnalysisTK.jar -T GenotypeGVCFs -nt ${threads} -R "$reference_flatfile" --dbsnp "$dbsnp_file_Broad_style" --useNewAFCalculator -o "$prefix".raw.vcf --variant mysample1.g.vcf --variant ...

Finally here I build the training models.

$java_path $java_opts -Xmx${ask_for_memory}g -jar ${gatk_binpath}GenomeAnalysisTK.jar -T VariantRecalibrator -nt 1 -R "$reference_flatfile" --input CR-MGUS.raw.vcf -resource:hapmap,known=false,training=true,truth=true,prior=15.0 "$hapmap_file" -resource:omni,known=false,training=true,truth=true,prior=12.0 "$omni_file" -resource:1000G,known=false,training=true,truth=false,prior=10.0 "$KG_file" -resource:NA12877,known=false,training=true,truth=true,prior=15.0 "$NA12877" -resource:NA12878,known=false,training=true,truth=true,prior=15.0 "$NA12878" -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 "$dbsnp_file_Broad_style" -an QD -an FS -an SOR -an MQRankSum -an ReadPosRankSum -an InbreedingCoeff --target_titv 3.2 -mode SNP -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 -recalFile "$prefix".recalibrate_SNP.recal -tranchesFile "$prefix".recalibrate_SNP.tranches -rscriptFile "$prefix".recalibrate_SNP_plots.R

I had above also -an MQ -MQCap 70 as an argument for VariantRecalibrator but dropped that hoping the models will improve. From the logs it turned out my MQ is uniformative. I thought it is because IndelRealignment and BQSR increased the MAPQ values by 10 (https://software.broadinstitute.org/gatk/blog?id=7847) and because I do not know how to restore the original MAPQ0 values I went a few BAM files backwards in my pipeline. But now I see the GVCFs differ in the TAGs annotated as well and I wonder whether that is the primary reason for VQSR failing.

$ grep 942242 mysample-18_18-PB.bwa.postalt.removed_duplicates.realignedtogether.BQSR.namesorted.fixmate.calmd.sorted.g.vcf
chr1       942242  .       C       <NON_REF>       .       .       END=942242      GT:DP:GQ:MIN_DP:PL      0/0:9:0:9:0,0,200
chr1       27942224        .       T       <NON_REF>       .       .       END=27942242    GT:DP:GQ:MIN_DP:PL      0/0:5:12:4:0,12,108
chr1       155942241       .       G       <NON_REF>       .       .       END=155942242   GT:DP:GQ:MIN_DP:PL      0/0:35:96:35:0,96,1440
$ grep 942242 mysample-18_18-PB.bwa.postalt.removed_duplicates.realignedtogether.g.vcf
chr1    942242  .       C       A,<NON_REF>     17.60   .       BaseQRankSum=-1.242;ClippingRankSum=0.000;DP=8;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQ=28800.00;ReadPosRankSum=1.242 GT:AD:DP:GQ:PL:SB       0/1:6,2,0:8:25:25,0,171,43,176,220:2,4,1,1
chr1    27942224        .       T       <NON_REF>       .       .       END=27942242    GT:DP:GQ:MIN_DP:PL      0/0:5:12:4:0,12,111
chr1    155942241       .       G       <NON_REF>       .       .       END=155942242   GT:DP:GQ:MIN_DP:PL      0/0:35:96:35:0,96,1440

I somewhat suspect the former GVCF is not recognized properly by GenotypeGVCFs but I did not receive any error during parsing. Where is the MQ value for example? There is only MQRankSum in GVCF files produced after IndelRealigner step but is missing in files after BQSR postprocessed using bamsort (from biobambam2) which calculated the MD/NM tags and did the sorting+indexing in one go. Please document which tags are needed by GenotypeGVCFs and which safety measures are in the code to ensure the GVCF file provided by user meets the specs (check for minimum relative incidence of required annotation tags?).

I use GATK-3.7.

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mmokrejs
    Hi,

    We have ValidateVariants with --validateGVCF option to validate GVCFs. I have asked the team to look into your issue. We will get back to you soon.

    Sheila

  • mmokrejsmmokrejs Czech RepublicMember

    No errors reported for any of my g.vcf files by ValidateVariants tool.

    No AS_InbreedingCoefficient records are in any of them ( http://gatkforums.broadinstitute.org/gatk/discussion/qna/accept?commentid=38676 ). Could that be a problem?

    However, I would like to see in this particular thread what are the required fields and automated checks being improved if input GVCF files do not comply to the conditions. Also, highest MQ being reported to ease decision what the -MQCap 70 value should actually be for a given dataset.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mmokrejs We aren't able to support use cases where you're processing the data through some other software that's not part of our Best Practices. If you follow the full pipeline without modification, we can guarantee that all necessary annotations should be present. Otherwise, there's no guarantee and it's up to you to check that. I understand that it would be convenient if the annotation outputs and interdependencies were documented, and we'll look into making that happen, but we're just not able to devote resources to that right now.

Sign In or Register to comment.