How do I calculate inbreeding on a single VCF file?

Question

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?

Answer 1

Possible solution:

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:

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.