Answers.Body: <p>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.</p>
  <p>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.</p>
  <p>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 <code>LowGQX</code> which means:</p>
  <blockquote>
  <p>The genotyping quality locus (GQX) is less than 30 or not present.</p>
  </blockquote>
  <p>And GQX is:</p>
  <blockquote>
  <p>The minimum of genotype quality assuming variant position and genotype quality assuming nonvariant position.</p>
  </blockquote>
  <p>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.</p>
  <p>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.</p>
  <p>Also, I urge you to have a look at tools like <a href="https://bedtools.readthedocs.io/en/latest/index.html" rel="nofollow noreferrer"><code>bedtools</code></a> and <a href="https://samtools.github.io/bcftools/bcftools.html" rel="nofollow noreferrer"><code>bcftools</code></a> that will help you filter your data.</p>
  <p>Finally, while your command will work, it is a bit inefficient, and you could simplify</p>
  <pre><code>gunzip -c individual.snp.vcf.gz | grep -v '^##' | 
    awk 'BEGIN {OFS=&quot;\t&quot;} {print $1, $2, $3, $4, $5, $6, $7}' | 
      awk 'BEGIN{OFS=&quot;\t&quot;} {print $0 , &quot;SNP&quot;}' | sed '1s/SNP/Type/'
  </code></pre>
  <p>to</p>
  <pre><code>zgrep -v '^##'  individual.snp.vcf.gz | 
    awk 'BEGIN {OFS=&quot;\t&quot;} {print $1, $2, $3, $4, $5, $6, $7, NR==1 ? &quot;Type&quot; : &quot;SNP&quot;}'
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-02T13:12:30.610Z
Answers.Text: Thank you for the answer! 
  I do have the fastq1 and fastq2 files. How should I proceed to filter by coverage in the various positions?
Answers.UserDisplayName: williantafsilva
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-02T13:30:52.530Z
Answers.Text: @williantafsilva you would have to align them to a genome to get a bam file and then you can filter the bam. If you need help with that, read [gatk best practices](https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels) and if you still need help ask a new question and we will try and help out.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-03T11:14:29.483Z
Answers.Text: Thank you! Very helpful!
Answers.UserDisplayName: williantafsilva
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-01T12:56:18.833Z
Answers.LastActivityDate: 2024-03-02T00:21:20.183Z
Answers.LastEditDate: 2024-03-02T00:21:20.183Z
Answers.UserDisplayName: terdon
Body: <p>I have a <code>individual.snp.vcf.gz</code> file of an individual genome and the <code>referencegenome.snp.vcf.gz</code> file of the reference genome.</p>
  <p>When I run the following code on the individual genome</p>
  <pre><code>gunzip -c individual.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS=&quot;\t&quot;} {print $1, $2, $3, $4, $5, $6, $7}' | awk 'BEGIN{OFS=&quot;\t&quot;} {print $0 , &quot;SNP&quot;}' | sed '1s/SNP/Type/'
  </code></pre>
  <p>I get this:</p>
  <pre><code>#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
  </code></pre>
  <p>When I run the equivalent code on the reference genome</p>
  <pre><code>gunzip -c referencegenome.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS=&quot;\t&quot;} {print $1, $2, $3, $4, $5, $6, $7}' | head -n100 | awk 'BEGIN{OFS=&quot;\t&quot;} {print $0 , &quot;SNP&quot;}' | sed '1s/SNP/Type/'
  </code></pre>
  <p>I get this:</p>
  <pre><code>#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
  </code></pre>
  <p>Now I want to get the genotype of the <em>individual genome</em> at all the sites listed in the <em>reference genome</em> table, so my question is: <strong>When I don't find a particular variant ID in the <em>individual genome</em> table (for example, <strong>rs1570391716</strong>), can I assume that the genotype of the individual is equal to the homozygous <em>REF</em> (<strong>AA</strong>, for <strong>rs1570391716</strong>) in the <em>reference genome</em> table?</strong> How certain can I be (assuming a low sequence reading error)?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-02-01T11:22:02.037Z
Comments.Text: What exactly is the `referencegenome.snp.vcf.gz` file? Is that a list of all known, reported, SNPs in this particular species? And how did you obtain the individual file? Note how all of the variants you show have a QUAL of 0 and are marked with [`LowGQX`](https://support.illumina.com/help/BaseSpace_App_WGS_v5_OLH_15050955_02/Content/Source/Informatics/Apps/VCFFileAnnotations_appENR_appISCWGS.htm), it is extremely unlikely that those are real calls.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-02-01T11:22:48.310Z
Comments.Text: Just FYI, there is no need for `gunzip`, you can just use `zgrep -v '^##' vcf.gz | awk ...`.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-02-01T12:07:23.723Z
Comments.Text: Yes, the `referencegenome.snp.vcf.gz` file is the list of all known, reported, SNPs in this particular species. And the individual file was obtained from the sequencing company (I don't know how exactly they obtained it) but I believe it was obtained through some standard variant calling methods.
  
  Regarding the QUAL of 0, let's ignore that for now.
Comments.UserDisplayName: williantafsilva
Tags: <vcf><genomics><genotyping>
Title: What does the absence of a variant in a VCF file mean?
ThreadText: # What does the absence of a variant in a VCF file mean?
  
  (Asked by: username_0 on 2024-02-01)
  
  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)?
  
  **Tags:** vcf, genomics, genotyping
  
  > **Comments:**
  > **username_1**: What exactly is the `referencegenome.snp.vcf.gz` file? Is that a list of all known, reported, SNPs in this particular species? And how did you obtain the individual file? Note how all of the variants you show have a QUAL of 0 and are marked with [`LowGQX`](https://support.illumina.com/help/BaseSpace_App_WGS_v5_OLH_15050955_02/Content/Source/Informatics/Apps/VCFFileAnnotations_appENR_appISCWGS.htm), it is extremely unlikely that those are real calls.  
  > **username_1**: Just FYI, there is no need for `gunzip`, you can just use `zgrep -v '^##' vcf.gz | awk ...`.  
  > **username_0**: Yes, the `referencegenome.snp.vcf.gz` file is the list of all known, reported, SNPs in this particular species. And the individual file was obtained from the sequencing company (I don't know how exactly they obtained it) but I believe it was obtained through some standard variant calling methods.  Regarding the QUAL of 0, let's ignore that for now.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-02-01. Score: 3 | Accepted answer
  
  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`](https://bedtools.readthedocs.io/en/latest/index.html) and [`bcftools`](https://samtools.github.io/bcftools/bcftools.html) 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"}'
  
  ```
  
  > **Comments:**
  > **username_0**: Thank you for the answer!  I do have the fastq1 and fastq2 files. How should I proceed to filter by coverage in the various positions?  
  > **username_1**: @username_0 you would have to align them to a genome to get a bam file and then you can filter the bam. If you need help with that, read [gatk best practices](https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels) and if you still need help ask a new question and we will try and help out.  
  > **username_0**: Thank you! Very helpful!  
  
  ---
  
  