1
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ 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? $\endgroup$
    – BigMistake
    Jun 23, 2023 at 8:50

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.