Holiday Notice:
The Frontline Support team will be slow to respond December 17-18 due to an institute-wide retreat and offline December 22- January 1, while the institute is closed. Thank you for your patience during these next few weeks. Happy Holidays!

Combine GVCF files problem

mahyarheymahyarhey BostonMember

I used the following command to combine 3 VCF files which are outputs of HaplotypeCaller:

java -jar data/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar \
-R data/ucsc.hg19.fasta \
-T CombineGVCFs \
--variant data/47V_post.ERC.vcf \
--variant data/48V_post.ERC.vcf \
--variant data/49V_post.ERC.vcf \
--out data/Combined_3files.vcf

However, after combined all 3 files, in the output final VCF, I can only see ./. genotypes. What is the problem? how I can to fix this?
Thanks

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mahyarhey
    Hi,

    You don't need to worry about this. It is the expected output of CombineGVCFs. Once you run GenotypeGVCFs on those combined gvcfs, you will get the genotypes.

    -Sheila

  • mahyarheymahyarhey BostonMember

    Hi Sheila,
    Thanks for your prompt reply!
    Can I use "genotypeGVCFs" in the above script instead of "CombineGVCFs" for saving time or it is necessary first run combineGVCFs and then run genotypeGVCFs? thanks

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mahyarhey
    Hi,

    You can use GenotypeGVCFs instead of CombineGVCFs. CombineGVCFs is useful if you have more than 200 gvcfs. For 3 gvcfs, you can simply run GenotypeGVCFs.

    -Sheila

  • mahyarheymahyarhey BostonMember

    Can I extract SNPs "rsID" first from each VCF file and then run GenotypeGVCFs? I have 34 VCF files
    because I just need the "SNP" information, the rest of the VCF is "."

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mahyarhey
    What do you mean? When you run GenotypeGVCFs, the rsid information will be carried over.

    If you are only interested in SNPs, you can extract them using SelectVariants.

    -Sheila

  • mahyarheymahyarhey BostonMember

    OK, let me to explain more. Supposed I got 3 VCF files after running HaplotypeCaller. Then Using "grep" command in UNIX I extracted only "rsID" from each VCF. The number of SNPs for each file is different. for example for file1=186,766 file2=257,392, and for file3=231,039. So, after merging these 3 files using "GenotypeGVCFs" what's happen to the uncommon SNPs. The program put ./. for those uncommon variants? Because the number is different in each file what would be expect the number of SNPs in final VCF file?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mahyarhey
    Hi,

    All variant sites in all of the sample GVCFs will be emitted in the final vcf from GenotypeGVCFs. If only one sample has a variant at a site, but the other two samples do not, the site will be emitted in the final vcf, but only the sample that has a variant will have a variant genotype. The other two samples will be homozygous reference.

    The number of SNPs in the final vcf should be the greatest number of SNPs in any of the samples. In your case, sample 2 has the most SNPs, so the final vcf will have 257,392.

    I hope this helps.

    -Sheila

  • cdrurycdrury United StatesMember

    Hi, after running combineGVCFs I'm also getting many sites with ./. as the genotype call. Earlier you mentioned that these should go away after running GenotypeGVCFs, but they remain, even for samples with read data. Even if the number of reads for a sample is low, after running genoypteGVCFs it reamins as ./. for the GT field.

    I also have examples where the same (low) number of reads does received at GT field call.

    Ex:

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @cdrury
    Hi,

    It looks like your example got cut off.

    After CombineGVCFs, all genotypes will be ./.

    If the ./. genotype is present in the ouput of GenotypeGVCFs, that means we don't have enough information to make a call there. I suspect the ./. samples have a GQ of 0. Again, this means the tool is not confident in the genotype.

    -Sheila

  • cdrurycdrury United StatesMember

    @Sheila

    Thanks again, I see what you mean where all sites are ./. after combiningGVCFs.

    However, after genotyping I have examples of same number of reads that are called and recieve a ./. and others that recieve a 1/1 or 0/0. I suppose this is based on the population level data, quality scores, etc. that the tool is using to establish confidence and not just on the read depth?

    -Ford

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @cdrury
    Hi Ford,

    Yes, HaplotypeCaller takes into account many other factors, not just read depth. Have a look at this document for more information. If two or more genotypes are equally likely, then the genotype is a no-call (./.).

    -Sheila

Sign In or Register to comment.