how to filter a gvcf

Hi GATK team,

I generated a gvcf file to get sequencing information for both variant and non-variant positions. I have couple of questions about how to filter out the unreliable callings.

1, in a gvcf file, I found records with all values of DP, GQ, MIN_DP and PL are 0 (see below), is it mean that this region is missing from the data (no reads were mapped to this region of reference genome)?
Scaffold10021 1 . A . . END=2870 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0

2, for case, with high depth but GQ=0. Can I trust the calling, or should set it to missing.
Scaffold10012 1189 . C . . END=1198 GT:DP:GQ:MIN_DP:PL 0/0:401:0:370:0,0,0
Scaffold10012 904 . A . . END=905 GT:DP:GQ:MIN_DP:PL 0/0:415:0:415:0,0,0

3,Which parameters were usually used to filter a gvcf file?

4, To accept a calling, what is minimum values for GQ and PL?

Thanks very much!

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @baosheng.wang1gmail.com
    Hello,

    Is there a reason you are trying to filter a GVCF? We don't recommend using the GVCF for final analysis; it is meant to be an intermediate file. You should use GenotypeGVCFs to generate a final VCF which you can use in the filtering step.

    As for filtering, you have two options: VQSR or hard filtering.

    -Sheila

  • Hi Shelia,

    Thanks very much for your reply. My reason to use gcvf instead of vcf file is:
    I calls SNPs (VCF file) in several data sets that each used in different project. To determine the ancestral status of SNPs in each data set, I aligned an outgroup species' reads to reference genome and generated a gvcf file. To avoid re-call SNPs in each data sets by combining in-group samples (500+ samples) and the out-group, I want to extract callings of both variant and non-variant position from gvcf of the out-group, and then annotate the ancestral status for SNPs called in different datasets. Is that make sense, or I should re-call the SNPs again?

    By the way, could you answer my previous questions 1-2 about the gvcf file? Here attached again. Thanks a lot for your time.
    1, in a gvcf file, I found records with all values of DP, GQ, MIN_DP and PL are 0 (see below), is it mean that this region is missing from the data (no reads were mapped to this region of reference genome)?
    Scaffold10021 1 . A . . END=2870 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0

    2, for case, with high depth but GQ=0. Can I trust the calling, or should set it to missing.
    Scaffold10012 1189 . C . . END=1198 GT:DP:GQ:MIN_DP:PL 0/0:401:0:370:0,0,0
    Scaffold10012 904 . A . . END=905 GT:DP:GQ:MIN_DP:PL 0/0:415:0:415:0,0,0

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @baosheng.wang1gmail.com
    Hi,

    Can you explain a little more about how you are planning to "annotate the ancestral status for SNPs called in different datasets"?

    As for your last two questions, can you please post the original bam file and bamout file for those regions?

    Thanks,
    Sheila

  • Hi,
    I have a reference genome, and sequencing data of 400, 300 and 500 samples of the same species (same species of the reference genome) from three different projects. Since I do not want to use all of these 1200 samples together, I called SNPs for each of the three cohorts by using the same reference genome, separately. Now, for a specific SNP with two alleles, I want know which one is the ancestral allele and which one is the derived. For this purpose, I sequenced a relative species, and called SNPs by using the same reference genome. For a SNP, the allele share by two species will be identified as the ancestral one. For SNP that neither reference nor alternative allele was found in the relative species, or both of them are found in relative species, we cannot determine the ancestral status, and these SNP will not be used.
    As I talked before, I have called SNPs (vcf file) for each of the three groups, and a gvcf file (using -ERC BP_RESOLUTION) for the single relative species. If I can filter the gvcf file and extract all varied and unvaried loci for the relative species, I would have a list of the allele of the relative species on each locus (those can be used explained above), and I could know the ancestral status of an allele by comparing this allele list and vcf file. Thus, I would not have to combine the gvcf of relative species with those of the three groups of samples and call SNPs by using GenotypeGVCFs again.
    Do you have any suggestions to identify ancestral status?

    For the last two questions,
    1, the first one "Scaffold10021 1 . A . . END=2870 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0", I found that there is no reads.
    2, send one I still not know why with high depth but GQ=0.
    Scaffold10012 1189 . C . . . GT:AD:DP:GQ:PL 0/0:408,0:408:0:0,0,0
    Scaffold10012 1190 . T . . . GT:AD:DP:GQ:PL 0/0:405,0:405:0:0,0,0

    3, I have a new question, heterozygous loci was called as a homozygous.
    Scaffold10012 15 . C . . . GT:AD:DP:GQ:PL 0/0:49,68:117:0:0,0,0
    Scaffold10012 16 . G . . . GT:AD:DP:GQ:PL 0/0:71,50:121:0:0,0,868

    bam and bamout files for region of Scaffold10012 were attached. Thanks very much for your help.
    I could not upload bam file, thus I added a txt extension to each file, you can simply remove that.

    Best regards,
    Baosheng

Sign In or Register to comment.