Hardy Weinberg Equilibrium exact test

Documentation

According to the comments:

This code implements an exact SNP test of Hardy-Weinberg Equilibrium as described in Wigginton, JE, Cutler, DJ, and Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg Equilibrium. American Journal of Human Genetics. 76: 000 - 000

Written by Jan Wigginton

Converted to python from the C/C++ code from: http://www.sph.umich.edu/csg/abecasis/Exact/c_instruct.html

This test is used by plink when applying the --hardy option (http://pngu.mgh.harvard.edu/~purcell/plink/summary.shtml#hardy) and from Beagle utilities (http://faculty.washington.edu/browning/beagle_utilities/utilities.html) to compute Hardy-Weinberg Equilibrium.

Parameters

Observed heterozygosity:
Observed AA homozygosity:
Observed aa homozygosity:

<inputs>
<param name="obs_hets" type="eval" value="" label="Observed heterozygosity: "/>
<param name="obs_hom1" type="eval" value="" label="Observed AA homozygosity: "/>
<param name="obs_hom2" type="eval" value="" label="Observed aa homozygosity: "/>
</inputs>

Return

A float value with the Hardy Weinberg Equilibrium exact test

See also

• Similar articles: Hardy_Weinberg_Equilibrium_asymptotic_test

Code

def Hardy_Weinberg_Equilibrium_exact_test(obs_hets, obs_hom1, obs_hom2):
        if obs_hom1 < 0 or obs_hom2 < 0 or obs_hets < 0:
                raise Exception("FATAL ERROR - SNP-HWE: Current genotype configuration (%s  %s %s) includes negative count" % (obs_hets, obs_hom1, obs_hom2))

        obs_homc = obs_hom2 if obs_hom1 < obs_hom2 else obs_hom1
        obs_homr = obs_hom1 if obs_hom1 < obs_hom2 else obs_hom2

        rare_copies = 2 * obs_homr + obs_hets
        genotypes   = obs_hets + obs_homc + obs_homr

        het_probs = [0.0] * (rare_copies + 1)

        #start at midpoint
        mid = rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes)

        #check to ensure that midpoint and rare alleles have same parity
        if (rare_copies & 1) ^ (mid & 1):
                mid += 1

        curr_hets = mid
        curr_homr = (rare_copies - mid) / 2
        curr_homc = genotypes - curr_hets - curr_homr

        het_probs[mid] = 1.0
        sum = float(het_probs[mid])

        for curr_hets in xrange(mid, 1, -2):

                het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0) / (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0))

                sum += het_probs[curr_hets - 2];

                # 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote
                curr_homr += 1
                curr_homc += 1

        curr_hets = mid
        curr_homr = (rare_copies - mid) / 2
        curr_homc = genotypes - curr_hets - curr_homr

        for curr_hets in xrange(mid, rare_copies - 1, 2):

                het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc / ((curr_hets + 2.0) * (curr_hets + 1.0))

                sum += het_probs[curr_hets + 2]

                #add 2 heterozygotes for next iteration -> subtract one rare, one common homozygote
                curr_homr -= 1
                curr_homc -= 1

        for i in xrange(0, rare_copies + 1):
                het_probs[i] /= sum

        #alternate p-value calculation for p_hi/p_lo
        p_hi = float(het_probs[obs_hets])
        for i in xrange(obs_hets, rare_copies+1):
                p_hi += het_probs[i]

        p_lo = float(het_probs[obs_hets])
        for i in xrange(obs_hets-1, -1, -1):
                p_lo += het_probs[i]

        p_hi_lo = 2.0 * p_hi if p_hi < p_lo else 2.0 * p_lo

        p_hwe = 0.0
        #  p-value calculation for p_hwe
        for i in xrange(0, rare_copies + 1):
                if het_probs[i] > het_probs[obs_hets]:
                        continue;
                p_hwe += het_probs[i]

        p_hwe = 1.0 if p_hwe > 1.0 else p_hwe

        return p_hwe

Unit Tests

def uni1():
	if int(10000 * Hardy_Weinberg_Equilibrium_exact_test(762, 269, 618)) != 1973:
		return "Failed to compute HWE of 762, 269, 618"

	return True

Development Code

def Hardy_Weinberg_Equilibrium_exact_test():
	pass

Permissions

Documentation Permissions

Kantale

Code Permissions

Kantale

Unit Tests Permissions

Kantale

Permissions Permissions

Kantale