To celebrate the release of GATK 4.0, we are giving away free credits for running the GATK4 Best Practices pipelines in FireCloud, our secure online analysis portal. It’s first come first serve, so sign up now to claim your free credits worth $250. Sponsored by Google Cloud. Learn more at https://software.broadinstitute.org/firecloud/documentation/freecredits

Variant recalibration (WQSR) after pooled calling

ClementKent2ClementKent2 Toronto, CanadaMember

Hi - we have whole genome samples consisting of non-barcoded DNA from 50 diploid individuals with 80-120x coverage. We use version 3.5-0-g36282e). We align, then use BQSR. For variant calling, we cannot use GATK - Haplotype caller doesn't support this and UG never finishes with such high ploidy levels (it's not just us; Huang et al 2015, "Evaluation of variant detection software for pooled next-generation sequence data", noted that UG failed to complete for large pools).

Thus, after BQSR we are calling variants with LoFreq. However, when we attempt to process the vcf file produced by LoFreq with VQSR, we get error messages. We suspect the absence of genotype columns is involved, although strictly speaking you can't call a genotype with a high ploidy pooled sample.

An example of the LoFreq vcf file is shown at the end of this discussion.

I'm quite willing to put together some awk scripts to add a non-haplotyped genotype column to our vcf, if this will be enough to make VQSR happy. Do you have any experience with this? A comment by delangel in 2013 stated the team had not at that time tried to do VQSR with pooled samples. Has this changed?

Example from LoFreq vcf:

INFO=<ID=DP4,Number=4,Type=Integer,Description="Counts for ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">

CHROM POS ID REF ALT QUAL FILTER INFO

1.1 1348 . C T 152 PASS DP=81;AF=0.345679;SB=8;DP4=27,26,9,19
1.1 1418 . G A 101 PASS DP=124;AF=0.233871;SB=2;DP4=45,50,16,13
1.1 1493 . T G 507 PASS DP=109;AF=0.522936;SB=9;DP4=30,21,24,33
1.1 1523 . A T 1640 PASS DP=86;AF=0.976744;SB=3;DP4=0,1,43,41

I hypothesize that if I produce something like the lines below this might help:

INFO=<ID=DP4,Number=4,Type=Integer,Description="Counts for ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samp1

1.1 1348 . C T 152 PASS DP=81;AF=0.345679;SB=8;DP4=27,26,9,19 GT:DP 1:81
1.1 1418 . G A 101 PASS DP=124;AF=0.233871;SB=2;DP4=45,50,16,13 GT:DP 1:124
1.1 1493 . T G 507 PASS DP=109;AF=0.522936;SB=9;DP4=30,21,24,33 GT:DP 1:109
1.1 1523 . A T 1640 PASS DP=86;AF=0.976744;SB=3;DP4=0,1,43,41 GT:DP 1:86

Is there some way forward here or is VQSR and high-ploidy pooled samples a no-go?

Thanks!

Comments

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Hi there,

    You're correct that GATK's germline callers won't do well on high ploidy pools, and indeed given a ploidy like yours they would never complete at all unless you limit the lowest frequency you're willing to discover to a large enough fraction. We have done quite a bit of work since that comment by Guillermo DelAngel, so the ploidy that is callable is higher than it used to be, but for your use case I agree you are probably better off using a different tool.

    That being said, the VariantRecalibrator should work just fine on a vcf that does not have genotypes. It doesn't use genotype information at all, and in fact one of the ways we improve speed of processing in production for large cohorts is to train the model on a sites-only vcf. Is the failure you encounter happening in VariantRecalibrator or ApplyRecalibration? If you post the text of the error we'll be in a better position to diagnose and solve your problem.
  • ClementKent2ClementKent2 Toronto, CanadaMember

    Hi Geraldine. The problem occurs in VariantRecalibrator. The full error message is attached. The relevant line:

    ERROR MESSAGE: The provided VCF file is malformed at approximately line number 96106: there are 117 genotypes while the header requires that 125 genotypes be present for all records at 3.5:252897

    I won't attach the full 142MB vcf file. However, the text below is the header, the area around line 96106, and the area around 3.5:252897.
    Note that tabs are in fact present in the file but aren't in some cases showing in the web page here.

    Header:

    fileformat=VCFv4.0

    fileDate=20170116

    source=lofreq call --call-indels -f /data3/01_ngm_BQSR_haplotypecaller_VQSR/am45new.fasta --no-default-filter -r 1.1:1-1382403 -o /tmp/lofreq2_call_parallelbJZhTh/0.vcf.gz /data2/ON/fastq/pilot_run1_rerun/filtered_data/4080-T.BQSR.bam

    reference=/data3/01_ngm_BQSR_haplotypecaller_VQSR/am45new.fasta

    INFO=<ID=DP,Number=1,Type=Integer,Description="Raw Depth">

    INFO=<ID=AF,Number=1,Type=Float,Description="Allele Frequency">

    INFO=<ID=SB,Number=1,Type=Integer,Description="Phred-scaled strand bias at this position">

    INFO=<ID=DP4,Number=4,Type=Integer,Description="Counts for ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">

    INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">

    INFO=<ID=CONSVAR,Number=0,Type=Flag,Description="Indicates that the variant is a consensus variant (as opposed to a low frequency variant).">

    INFO=<ID=HRUN,Number=1,Type=Integer,Description="Homopolymer length to the right of report indel position">

    FILTER=<ID=min_snvqual_79,Description="Minimum SNV Quality (Phred) 79">

    FILTER=<ID=min_indelqual_66,Description="Minimum Indel Quality (Phred) 66">

    FILTER=<ID=min_dp_10,Description="Minimum Coverage 10">

    FILTER=<ID=sb_fdr,Description="Strand-Bias Multiple Testing Correction: fdr corr. pvalue > 0.001000">

    FILTER=<ID=min_snvqual_101,Description="Minimum SNV Quality (Phred) 101">

    FILTER=<ID=min_indelqual_88,Description="Minimum Indel Quality (Phred) 88">

    CHROM POS ID REF ALT QUAL FILTER INFO

    1.1 1089 . T A 806 PASS DP=246;AF=0.353659;SB=4;DP4=71,87,45,42
    1.1 1152 . G A 225 PASS DP=197;AF=0.263959;SB=0;DP4=49,48,26,26

    near 96106:
    1.16 93429 . G A 1257 PASS DP=170;AF=0.488235;SB=0;DP4=46,41,44,39
    1.16 93518 . G A 581 PASS DP=175;AF=0.274286;SB=0;DP4=59,68,22,26
    1.16 93614 . A G 141 PASS DP=163;AF=0.110429;SB=2;DP4=74,71,8,10
    1.16 93630 . T C 1284 PASS DP=163;AF=0.521472;SB=1;DP4=37,41,43,42
    1.16 93632 . T A 109 PASS DP=163;AF=0.104294;SB=0;DP4=71,74,8,9 <--- line 96106
    1.16 93634 . T C 102 PASS DP=158;AF=0.107595;SB=0;DP4=70,71,8,9
    1.16 93686 . G A 1072 PASS DP=151;AF=0.509934;SB=0;DP4=38,35,39,38
    1.16 93703 . G A 1101 PASS DP=141;AF=0.531915;SB=1;DP4=27,26,41,34
    1.16 93789 . T C 1693 PASS DP=151;AF=0.695364;SB=0;DP4=22,23,50,56

    around 3.5:252897:
    3.5 252851 . G A 133 PASS DP=149;AF=0.436242;SB=3;DP4=44,40,29,36
    3.5 252853 . C T 544 PASS DP=150;AF=0.433333;SB=3;DP4=44,40,29,36
    3.5 252873 . A G 451 PASS DP=144;AF=0.937500;SB=0;DP4=1,2,66,69
    3.5 252897 . G A 869 PASS DP=146;AF=0.513699;SB=4;DP4=29,42,37,38 <---3.5 252897
    3.5 252919 . G A 694 PASS DP=144;AF=0.437500;SB=0;DP4=36,45,29,34
    3.5 253056 . T C 2533 PASS DP=154;AF=0.876623;SB=2;DP4=12,6,77,58
    3.5 253368 . T C 348 PASS DP=156;AF=0.205128;SB=7;DP4=68,56,13,19
    3.5 253423 . G A 1877 PASS DP=147;AF=0.727891;SB=8;DP4=17,23,61,46

    If you need the whole vcf file let me know where to send it.

    Issue · Github
    by Sheila

    Issue Number
    1676
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @ClementKent2
    Hi,

    Can you post the entire VCF record for position 3.5:252897? Please also post the entire VCF record for the very first site in your VCF. I need to see the FORMAT section as well, because it looks like that is where the issue is occurring.

    Thanks,
    Sheila

  • ClementKent2ClementKent2 Toronto, CanadaMember

    Hi Sheila. The lines I posted are in fact the entire VCF records at those places. As noted in my original post:

    "after BQSR we are calling variants with LoFreq. However, when we attempt to process the vcf file produced by LoFreq with VQSR, we get error messages. We suspect the absence of genotype columns is involved, although strictly speaking you can't call a genotype with a high ploidy pooled sample."

    I repeat: the vcf file we get from LoFreq does not have a FORMAT and sample column, and hence does not have a GT tag.

    To which Geraldine replied:

    " the VariantRecalibrator should work just fine on a vcf that does not have genotypes. It doesn't use genotype information at all, and in fact one of the ways we improve speed of processing in production for large cohorts is to train the model on a sites-only vcf."

    So, the entire line at position 3.5:252897 is:
    3.5 252897 . G A 869 PASS DP=146;AF=0.513699;SB=4;DP4=29,42,37,38

    Thanks,
    Clement

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Hi Clement, the fact that error message states that there are genotypes at all when clearly there are none, and the odd conflation of a line number and record position that have nothing to do with each other, both suggest a rather severe parsing issue. This could either be due to a corrupted file copy or a transient system error. Is the error reproducible? If yes, have you validated the copy of the file you are running on?
  • ClementKent2ClementKent2 Toronto, CanadaMember

    Hi Geraldine, yes the error is reproducible. We ran VariantRecalibrator with 4 completely different LoFreq vcf files, and got the same error each time. Although the line number and snp location reported differed for each file, curiously the message part "there are 117 genotypes while the header requires that 125 genotypes be present" was the same in all 4 runs. I concur that it sounds like some kind of parsing error might be involved.

    You asked "have you validated the copy of the file you are running on? ". Is there a specific validation process you are referring to?

  • ClementKent2ClementKent2 Toronto, CanadaMember

    Here is a log file of two runs, each with a different input vcf file from the one first reported.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    There's a GATK tool called ValidateVariants that you can use for this purpose.
Sign In or Register to comment.