Why are there missing calls in a VCF file from exome sequencing?

Question

My data is a VCF file generated from an exome sequencing variant call pipeline. I'm not very familiar with the sequencing and variant calling process. I noticed that there are some missing genotypes, which are recorded as "./." at the GT field. From googling I learned that they're not homozygous reference genotype ("0/0"), but missing calls due to some sequencing failure.

In order to make my data "cleaner", I thought it would be better to filter out the loci with missing calls if there're not too many of them. I also checked the corresponding "DP" of the loci with and without missing calls. For example:

chr1:123 GT:DP 0/0:2 1/1:4 ./.
chr1:234 GT:DP 0/0:10 1/1:11 1/1:20
chr1:345 GT:DP 0/1:40 1/1:37 0/0:78
chr1:456 GT:DP 0/1:7 0/0:23 ./.
chr1:567 GT:DP 0/1:34 1/1:39 0/0:58

In the above toy example, there are 3 people and 5 loci. I checked the mean DP of all 5 loci and I found that the mean DP of loci with missing calls (1st and 4th loci) to be significantly lower than the mean DP of loci without missing calls (2nd, 3rd and 5th loci). This happens to my real data with a lot more loci and samples than this toy example. Is that a coincidence? Or are there any specific reason that the loci with missing calls have lower coverage than the "normal" loci? Thanks!

Answer 1

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses.

If you did your variant calling on samples separately and then merged them together, you'll have missing calls at loci called variant in one sample and homozygous reference in the other.

You can get over this by asking the variant caller to emit all sites within the target region, which would lead to a large VCF or do multi-sample calling, which is more computationally intensive.

One other way is to do variant calling on all samples separately, generating a BED file of all variant regions in the cohort and then asking the variant caller to genotype those loci across all samples before merging them together.

Answer 2

I've just been generating data like this, so can tell you about why/how missing calls are created in my dataset. There are two main reasons:

1. Sequencing failure

When reads don't map across the variant region, then it's impossible to accurately determine a genotype for that region. This will commonly happen just outside the borders of the selected regions for exome sequencing, but can also happen through natural random sampling of the genome, or through repeated sequences coupled with systematic error in either the assembled genome or the sequencer. In this case, it would be expected that called variant regions including missing data would be lower coverage. It may be helpful to plot the coverage within a particular region to work out if this is likely to be happening.

2. Merged non-variant data

Variant calling can be processed in two different modes, one which takes into account all samples when doing calling, and another which does the calling one sample at a time, followed by subsequent merging of the results. Multi-sample variant calling is more accurate, but computational limits (or experimental design) might mean that it is more appropriate for samples to be called one-by-one. It's fairly common (for the purpose of preserving space and computational power) for calling algorithms to only report variant information so if a single sample is identical to the reference, and that sample is the only one being called, then the variant and coverage information for that sample will be lost. When samples are merged after being separately called, the declared call for the "same as reference" samples will be set to missing. In this case, it would be generally expected that the other called samples would have high-coverage within the variant region.