Answers.Body: <p><strong>Possible solution</strong>:</p>
  <pre class="lang-py prettyprint-override"><code>for record in vcf_reader:
      if not record.is_snp or len(record.ALT) != 1:
          continue  # skip if not a biallelic SNP
  
      gt = record.genotype('Sample')['GT']
  
      if gt == '0/1' or gt == '1/0':
          obs_hets += 1
      elif gt == '0/0' or gt == '1/1':
          obs_hom += 1
  
  num_variants = obs_hom + obs_hets
  exp_hets = num_variants * 0.5  # under the assumption of Hardy-Weinberg Equilibrium
  
  f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient
  </code></pre>
  <p>Example output: &quot;-0.25&quot;</p>
  <p>Interpretation: We observed more heterozygosity than we expected.</p>
  <p>I don't know if this method is valid, however.</p>
  <p><strong>A different way to calculate it</strong>:</p>
  <pre class="lang-py prettyprint-override"><code>for record in vcf_reader:
      if not record.is_snp or len(record.ALT) != 1:
          continue  # skip if not a biallelic SNP
  
      gt = record.genotype('Sample')['GT']
  
      if gt == '0/1' or gt == '1/0':
          obs_hets += 1
          allele_counts['A'] += 1
          allele_counts['a'] += 1
      elif gt == '0/0':
          allele_counts['A'] += 2
      elif gt == '1/1':
          allele_counts['a'] += 2
  
  total_alleles = sum(allele_counts.values())
  
  # calculate allele frequencies
  p = allele_counts['A'] / total_alleles
  q = allele_counts['a'] / total_alleles
  
  # calculate expected heterozygosity
  exp_hets = 2 * p * q * total_alleles / 2  # the total number of variants
  
  f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient
  </code></pre>
  <p>Example output: &quot;-0.45&quot;</p>
  <p>Interpretation: We observed more heterozygosity than we expected. I don't have a good idea of the meaning of this number.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-06-23T05:45:57.203Z
Answers.LastActivityDate: 2024-03-19T15:37:25.970Z
Answers.LastEditDate: 2024-03-19T15:37:25.970Z
Answers.UserDisplayName: BigMistake
Body: <p>I would like to know the level of inbreeding for a single sample VCF file. I'm agnostic to the specific method of measuring inbreeding; any should work. However, I'd like to avoid having to specify a &quot;background&quot; or &quot;reference&quot; population. Ideally, the only input should be the single VCF file containing one human genome sample (WGS, hg38). (I will run this multiple times on different samples, but I'd like to run them individually.)</p>
  <p>The existing tools seem to all require multiple VCFs.</p>
  <p>Relatedly, is runs of homozygosity an acceptable metric for evaluating inbreeding in this context? And are there any special considerations to consider for various ethnic groups, like Ashkenazi Jews, that might influence the interpretation?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-06-23T08:50:55.470Z
Comments.Text: https://samtools.github.io/bcftools/howtos/roh-calling.html will output data, but I don't know how to get its output into a usable format. Maybe the autozygosity is just the number of lines in the ST section divided by the number of lines processed?
Comments.UserDisplayName: BigMistake
Tags: <vcf><genome>
Title: How do I calculate inbreeding on a single VCF file?
ThreadText: # How do I calculate inbreeding on a single VCF file?
  
  (Asked by: username_0 on 2023-06-23)
  
  I would like to know the level of inbreeding for a single sample VCF file. I'm agnostic to the specific method of measuring inbreeding; any should work. However, I'd like to avoid having to specify a "background" or "reference" population. Ideally, the only input should be the single VCF file containing one human genome sample (WGS, hg38). (I will run this multiple times on different samples, but I'd like to run them individually.)
  
  The existing tools seem to all require multiple VCFs.
  
  Relatedly, is runs of homozygosity an acceptable metric for evaluating inbreeding in this context? And are there any special considerations to consider for various ethnic groups, like Ashkenazi Jews, that might influence the interpretation?
  
  **Tags:** vcf, genome
  
  > **Comments:**
  > **username_0**: https://samtools.github.io/bcftools/howtos/roh-calling.html will output data, but I don't know how to get its output into a usable format. Maybe the autozygosity is just the number of lines in the ST section divided by the number of lines processed?  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2023-06-23. Score: 1
  
  **Possible solution**:
  
  ```py
  for record in vcf_reader:
      if not record.is_snp or len(record.ALT) != 1:
          continue  # skip if not a biallelic SNP
  
      gt = record.genotype('Sample')['GT']
  
      if gt == '0/1' or gt == '1/0':
          obs_hets += 1
      elif gt == '0/0' or gt == '1/1':
          obs_hom += 1
  
  num_variants = obs_hom + obs_hets
  exp_hets = num_variants * 0.5  # under the assumption of Hardy-Weinberg Equilibrium
  
  f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient
  
  ```
  
  Example output: "-0.25"
  
  Interpretation: We observed more heterozygosity than we expected.
  
  I don't know if this method is valid, however.
  
  **A different way to calculate it**:
  
  ```py
  for record in vcf_reader:
      if not record.is_snp or len(record.ALT) != 1:
          continue  # skip if not a biallelic SNP
  
      gt = record.genotype('Sample')['GT']
  
      if gt == '0/1' or gt == '1/0':
          obs_hets += 1
          allele_counts['A'] += 1
          allele_counts['a'] += 1
      elif gt == '0/0':
          allele_counts['A'] += 2
      elif gt == '1/1':
          allele_counts['a'] += 2
  
  total_alleles = sum(allele_counts.values())
  
  # calculate allele frequencies
  p = allele_counts['A'] / total_alleles
  q = allele_counts['a'] / total_alleles
  
  # calculate expected heterozygosity
  exp_hets = 2 * p * q * total_alleles / 2  # the total number of variants
  
  f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient
  
  ```
  
  Example output: "-0.45"
  
  Interpretation: We observed more heterozygosity than we expected. I don't have a good idea of the meaning of this number.
  
  
  ---
  
  