VSQR resource format: which elements of vcf are needed?

grafalgrafal MPI TuebingenMember

Dear GATK team,

I am working with maize aDNA and would like to find SNPs called in aDNA samples that are at least as good as those in HapMap.

Am right to assume that variant recalibration is the correct tool for the job? My problem is that I can't find details about format of "resources" used for training. (I would like to see your hapmap.vcf but i can't log in to your ftp; it tells me that there can be only 25 users there).

From your documentation (https://www.broadinstitute.org/gatk/guide/article?id=2805 ; https://www.broadinstitute.org/gatk/guide/article?id=1259 and https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_variantrecalibration_VariantRecalibrator.php#--resource) I gather that the training set needs to be in vcf file.

Here is my question:
Which components are important in resource vcf file; positions? variants? annotations? are SNP calls for each individual important?
In other words, is the purpose of training "resource" to
1) find overlapping sites in my vcf and use annotations that are in my vcf or
2) is it the training resource annotations that are being used?
If the former is true, would it be possible to use bed files instead? If the latter is true, do I need individual calls? HapMap for my organism has > 1000 individuals and is very heavy; if possible I would like to avoid downloading it all.

Many thanks for your help,

Best wishes,
Rafal

Issue · Github
by Sheila

Issue Number
906
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
chandrans

Best Answers

Answers

  • grafalgrafal MPI TuebingenMember

    Hello again,

    Thought I would report that stripping hapmap.vcf of genotype calls causes problems with VariantRecalibrator:

    ERROR MESSAGE: The provided VCF file is malformed at approximately line number : there are 0 genotypes while the header requires that genotypes be present for all records at
    ERROR ------------------------------------------------------------------------------------------

    My workaround is to use hapmap.bed file to intersect de_novo_calls.vcf using bedtools intersect. Following that I am sorting intersected hapmap_de_novo_calls.vcf with picard SortVcf and feed it as resource (true positives) for training.

    If GATK ever comes to upgrading that, please allow bed files to select sites as positives for training. It is just way less redundant.

    Best,
    Rafal

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Rafal, can you tell me how you stripped the genotypes?

    We use a sites-only hapmap in our own work so I can guarantee it should work.

  • grafalgrafal MPI TuebingenMember

    Hi Geraldine,

    I wrote a perl script that checks if ref and alt alleles are congruent in hapmap.bed and de_novo_calls.vcf and outputs positions without columns coding for genotypes. Given your guarantee, I assume that I should have changed header as well in some way?

    Cheers,
    Rafal

Sign In or Register to comment.