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.

Variant quality score recalibration

natenate Member
edited July 2012 in Ask the GATK team

Dear team,

I have a problem when running the vqsr. When running the apply recalibration I get the error that the recal file that is generated by the variant recalibration step is malformet. If I take a look at this file it looks like this:

fileformat=VCFv4.1

CHROM POS ID REF ALT QUAL FILTER INFO

1 14574 . N . . END=14574;VQSLOD=-11.1678;culprit=MQ
1 14907 . N . . END=14907;VQSLOD=-9.1812;culprit=MQ
1 14930 . N . . END=14930;VQSLOD=-9.5518;culprit=MQ
1 14933 . N . . END=14933;VQSLOD=-13.0475;culprit=MQ
1 14948 . N . . END=14948;VQSLOD=-13.8117;culprit=MQ
1 15118 . N . . END=15118;VQSLOD=-7.5026;culprit=MQ
1 16571 . N . . END=16571;VQSLOD=-12.1727;culprit=MQ
1 16580 . N . . END=16580;VQSLOD=-15.5872;culprit=MQ
1 63671 . N . . END=63671;VQSLOD=-7.1377;culprit=MQ
1 63697 . N . . END=63697;VQSLOD=-8.1332;culprit=MQ
1 69270 . N . . END=69270;VQSLOD=-10.5032;culprit=MQ
1 136048 . N . . END=136048;VQSLOD=-15.2755;culprit=MQ
1 568214 . N . . END=568214;VQSLOD=-6.0049;culprit=MQ
1 663097 . N . . END=663097;VQSLOD=-8.7168;culprit=MQ
1 752894 . N . . END=752894;VQSLOD=1.4287;culprit=MQ
1 753269 . N . . END=753269;VQSLOD=-2.3944;culprit=MQ

The input file is a multi-sample vcf.

This is the commant I use for the variant recalibration

java -Xmx4g -jar $GATK \
-T VariantRecalibrator -R $REF \
-input $VCFFILE \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /panfs1/breastcancer/naspeh/tools/GATK/GenomeAnalysisTK-1.0.5506/hapmap_3.3.b37.sites.vcf \
-resource:omni,known=false,training=true,truth=false,prior=12.0 /panfs1/breastcancer/naspeh/tools/GATK/GenomeAnalysisTK-1.0.5506/1000G_omni2.5.b37.sites.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=8.0 /home/projects/pr_46298/naspeh/resources/dbsnp_135.b37.vcf \
-an QD \
-an HaplotypeScore \
-an MQRankSum \
-an ReadPosRankSum \
-an MQ \
-recalFile $RECALFILE \
-tranchesFile $TRANCESFILE \
-rscriptFile $RSCRIPTFILE \
--percentBadVariants 0.05 \
--maxGaussians 4

I cannot see what I am doing wrong here.

Best,
N

Best Answer

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MA admin
    Accepted Answer

    As for your problem, this looks like a user error. There must be something wrong with your inputs, all the reference bases in the file excerpt you show are N. That's not something the GATK would generate. You need to validate the input vcf.

Answers

  • SophiaSophia Member

    there is no # at the beginning of the fileformat and CHROM lines?

  • natenate Member

    Yes, there is,

    fileformat=VCFv4.1

    CHROM POS ID REF ALT QUAL FILTER INFO

    1 14574 . N . . END=14574;VQSLOD=-11.1678;culprit=MQ
    1 14907 . N . . END=14907;VQSLOD=-9.1812;culprit=MQ
    1 14930 . N . . END=14930;VQSLOD=-9.5518;culprit=MQ
    1 14933 . N . . END=14933;VQSLOD=-13.0475;culprit=MQ

  • natenate Member

    But the format change when I post in this forum

  • Mark_DePristoMark_DePristo Broad InstituteMember admin

    Try posting surrounded with "pre" html tags

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Either use <pre> tags or you can use Markdown (see link below the comment box) to tell the forum not to reformat. The best way is to have 4 spaces or a tab at the start of each line (you can do this in one go for the whole block of text with any good text editor). The result will look like this:

    ##fileformat=VCFv4.1
    ##fileDate=20090805
    ##source=myImputationProgramV3.1
    ##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
    ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
    ##phasing=partial
    ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
    ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
    ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
    ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
    ##FILTER=<ID=q10,Description="Quality below 10">
    ##FILTER=<ID=s50,Description="Less than 50% of samples have data">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
    ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
    #CHROM POS     ID        REF    ALT     QUAL FILTER INFO                              FORMAT      NA00001        NA00002        NA00003
    20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
    20     17330   .         T      A       3    q10    NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3   0/0:41:3
    20     1110696 rs6040355 A      G,T     67   PASS   NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2   2/2:35:4
    20     1230237 .         T      .       47   PASS   NS=3;DP=13;AA=T                   GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
    20     1234567 microsat1 GTC    G,GTCT  50   PASS   NS=3;DP=9;AA=G                    GT:GQ:DP    0/1:35:4       0/2:17:2       1/1:40:3
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    Accepted Answer

    As for your problem, this looks like a user error. There must be something wrong with your inputs, all the reference bases in the file excerpt you show are N. That's not something the GATK would generate. You need to validate the input vcf.

  • natenate Member

    Hi again,

    I just ran the commands once again changing the version of gatk. With version GenomeAnalysisTK-1.5-30-g27e7e17 I do not have any problems and the .recal file looks fine like this:

    7,100551008,100551008,-3260.1158,MQ
    7,100551004,100551004,-2884.8729,MQ
    7,100551003,100551003,-2841.6938,ReadPosRankSum
    7,100551012,100551012,-2556.0079,ReadPosRankSum
    7,100550808,100550808,-2331.9863,QD
    7,100550995,100550995,-1688.9692,MQ
    7,100550806,100550806,-1537.2323,MQ
    7,100551015,100551015,-1442.7138,HaplotypeScore
    11,1016918,1016918,-1223.9386,HaplotypeScore
    11,1016919,1016919,-1079.2547,HaplotypeScore
    11,1017893,1017893,-951.2670,HaplotypeScore
    7,100550837,100550837,-946.2444,HaplotypeScore
    

    But with version GenomeAnalysisTK-1.6-2-gc2b74ec it looks like this

    ##fileformat=VCFv4.1
    #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    
    1       14574   .       N         .       .       END=14574;VQSLOD=-11.1678;culprit=MQ
    1       14907   .       N         .       .       END=14907;VQSLOD=-9.1812;culprit=MQ
    1       14930   .       N         .       .       END=14930;VQSLOD=-9.5518;culprit=MQ
    1       14933   .       N         .       .       END=14933;VQSLOD=-13.0475;culprit=MQ
    1       14948   .       N         .       .       END=14948;VQSLOD=-13.8117;culprit=MQ
    1       15118   .       N         .       .       END=15118;VQSLOD=-7.5026;culprit=MQ
    1       16571   .       N         .       .       END=16571;VQSLOD=-12.1727;culprit=MQ
    1       16580   .       N         .       .       END=16580;VQSLOD=-15.5872;culprit=MQ
    1       63671   .       N         .       .       END=63671;VQSLOD=-7.1377;culprit=MQ
    1       63697   .       N         .       .       END=63697;VQSLOD=-8.1332;culprit=MQ
    1       69270   .       N         .       .       END=69270;VQSLOD=-10.5032;culprit=MQ
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Nate,

    Both of those are retired versions of the software that we no longer support. Please update to GATK 2.0, run your analysis again and let us know if you still have problems.

  • natenate Member

    Ok, sorry for that. I have been looking at your old pages and did not realize there is a new version ready.

  • rpoplinrpoplin Member ✭✭✭

    Yeah, it looks like the problem might be from mixing different version of the GATK when the format of the recal file changed in the middle on you.

  • natenate Member

    Yes, that might the reason...Thanks a lot!

Sign In or Register to comment.