Hardy-Weinberg平衡检验计算方法.doc

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148/* Copyright (c) 2009 The Broad Institute* Permission is hereby granted, free of charge, to any person* obtaining a copy of this software and associated documentation* files (the Software), to deal in the Software without* restriction, including without limitation the rights to use,* copy, modify, merge, publish, distribute, sublicense, and/or sell* copies of the Software, and to permit persons to whom the* Software is furnished to do so, subject to the following* conditions:* The above copyright notice and this permission notice shall be* included in all copies or substantial portions of the Software.* THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND,* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR* OTHER DEALINGS IN THE SOFTWARE.*/package org.broad.tribble.util.popgen;/* The Broad Institute* SOFTWARE COPYRIGHT NOTICE AGREEMENT* This software and its documentation are copyright 2004 by the* Broad Institute/Massachusetts Institute of Technology. All rights are reserved.* This software is supplied without any warranty or guaranteed support whatsoever. Neither* the Broad Institute nor MIT can be responsible for its use, misuse, or functionality.*/* This class calculates a HardyWeinberg p-value given three values representing* the observed frequences of homozygous and heterozygous genotypes in the* test population.* author Bob Handsaker*/public final class HardyWeinbergCalculation /* * This class is not instantiable. */ private HardyWeinbergCalculation() /* * Calculates exact two-sided hardy-weinberg p-value. Parameters * are number of genotypes, number of rare alleles observed and * number of heterozygotes observed. * * (c) 2003 Jan Wigginton, Goncalo Abecasis () */ public static double hwCalculate(int obsAA, int obsAB, int obsBB) int diplotypes = obsAA + obsAB + obsBB; int rare = (obsAA * 2) + obsAB; int hets = obsAB; /make sure rare allele is really the rare allele if (rare diplotypes) rare = (2 * diplotypes) - rare; /make sure numbers arent screwy if (hets rare) return -1; double tailProbs = new doublerare + 1; for (int z = 0; z 1; het -= 2) tailProbshet - 2 = (tailProbshet * het * (het - 1.0) / (4.0 * (hom_r + 1.0) * (hom_c + 1.0); sum += tailProbshet - 2; /2 fewer hets for next iteration - add one rare and one common homozygote hom_r+; hom_c+; het = mid; hom_r = (rare - mid) / 2; hom_c = diplotypes - het - hom_r; for (het = mid; het subtract one rare and one common homozygote hom_r-; hom_c-; for (int z = 0; z tailProbs.length; z+) tailProbsz /= sum; double top = tailPr