What does the absence of a variant in a VCF file mean?

Question

I have a individual.snp.vcf.gz file of an individual genome and the referencegenome.snp.vcf.gz file of the reference genome.

When I run the following code on the individual genome

gunzip -c individual.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

I get this:

#CHROM  POS ID  REF ALT QUAL    FILTER  Type
chr1    10067   .   T   A   0   LowGQX  SNP
chr1    10070   .   C   A   0   LowGQX  SNP
chr1    10073   .   T   C   0   LowGQX  SNP
chr1    10105   .   A   C   0   LowGQX  SNP
chr1    10108   rs62651026  C   T   0   LowGQX  SNP
chr1    10157   .   T   C   0   LowGQX  SNP
chr1    10177   rs201752861 A   C   17  LowGQX  SNP
chr1    10180   rs201694901 T   C   0   LowGQX  SNP
chr1    10181   rs1246412344    A   T   0   LowGQX  SNP
chr1    10250   rs199706086 A   C   0   LowGQX  SNP
chr1    10257   rs111200574 A   C   0   LowGQX  SNP
chr1    10285   rs866375379 T   C   67  LowGQX  SNP
chr1    10321   rs1002315756    C   T   0   LowGQX  SNP
chr1    10327   rs112750067 T   C   36  LowGQX  SNP
chr1    10332   rs1175748383    C   A   0   LowGQX  SNP

When I run the equivalent code on the reference genome

gunzip -c referencegenome.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | head -n100 | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

I get this:

#CHROM  POS ID  REF ALT QUAL    FILTER  Type
1   10001   rs1570391677    T   A   .   .   SNP
1   10002   rs1570391692    A   C   .   .   SNP
1   10003   rs1570391694    A   C   .   .   SNP
1   10008   rs1570391698    A   G   .   .   SNP
1   10009   rs1570391702    A   G   .   .   SNP
1   10015   rs1570391706    A   G   .   .   SNP
1   10020   rs1570391708    A   C   .   .   SNP
1   10021   rs1570391710    A   G   .   .   SNP
1   10026   rs1570391712    A   C   .   .   SNP
1   10027   rs1570391716    A   C,G .   .   SNP
1   10032   rs1570391720    A   C   .   .   SNP
1   10033   rs1570391722    A   G   .   .   SNP
1   10039   rs978760828 A   C   .   .   SNP
1   10043   rs1008829651    T   A   .   .   SNP
1   10045   rs1570391729    A   C,G .   .   SNP
1   10051   rs1052373574    A   C,G .   .   SNP
1   10055   rs892501864 T   A   .   .   SNP
1   10056   rs1570391738    A   C   .   .   SNP

Now I want to get the genotype of the individual genome at all the sites listed in the reference genome table, so my question is: When I don't find a particular variant ID in the individual genome table (for example, rs1570391716), can I assume that the genotype of the individual is equal to the homozygous REF (AA, for rs1570391716) in the reference genome table? How certain can I be (assuming a low sequence reading error)?

Answer 1

The absence of a variant means the variant was not found. However, that does not mean you can assume the position in the sequenced individual was homozygous for what is in the reference. Perhaps the position was simply not sequenced. Or the coverage was too low. So your first step will have to be to focus only on the regions covered by whatever sequencing assay was used to sequence your data, and only those regions with sufficient coverage for a variant to be called.

There is no magic number here, the threshold will depend on how the data were generated. For instance, a minimum coverage of 8 is a reasonable lower bound, but these days, we often have coverage values of 60x or more, and in the case of amplicon assays, we can have coverage in the thousands. So the threshold is something you will need to decide on.

Next, you really need to take the quality and filters into account. You haven't shown us the full VCF line, but my bet is that all the variants you show in your individual file are most likely not genotyped at all. These are not actually present. They all have a score of 0 and are marked with LowGQX which means:

The genotyping quality locus (GQX) is less than 30 or not present.

And GQX is:

The minimum of genotype quality assuming variant position and genotype quality assuming nonvariant position.

Variants with a QUAL score of 0 and low genotype quality are almost certainly not variants at all and you really shouldn't consider them.

So, since you apparently don't have the raw data, neither fastq nor bam, that means you cannot easily filter by coverage in the various positions. However, at the very least, you need to find out what sequencing assay and method was used to sequence your sample. Then, you can get a bed file with the regions of interest, and then limit your analysis to variants in those regions.

Also, I urge you to have a look at tools like bedtools and bcftools that will help you filter your data.

Finally, while your command will work, it is a bit inefficient, and you could simplify

gunzip -c individual.snp.vcf.gz | grep -v '^##' | 
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | 
    awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

to

zgrep -v '^##'  individual.snp.vcf.gz | 
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7, NR==1 ? "Type" : "SNP"}'