2006-09-29

HWEの正確検定

わけあって、J.E. Wigginton, G.R. AbecasisらによるHWE正確検定をJava化する必要が出た。このソースは、著者らにより、C/C++,R,Fortranにて公開されている。

N人、2N本染色体、アレル(A/B)、 $n_A+n_B=N$ , $n_A:No.allele A$ , $n_B:No.allele B$ 、 $n_{AA}+n_{AB}+n_{BB}=N$ , $n_{XY}:No.individulas with genotype XY$ 、と表すこととすると、

$n_{AA}=￥frac{n_A-n_{AB}}{2}$ , $n_{BB}=￥frac{n_{B}-n_{AB}}{2}$ と表せる。

また、2N本のアレルの分け方は全部で $￥frac{(2N)!}{n_A!n_B!}$ であり、N人の3ジェノタイプへの分け方は全部で $￥frac{N!}{n_{AA}!n_{AB}!n_{BB}!}$ ある。

今、ヘテロの人数 $n_{AB}$ について、その２アレルの取り方はそれぞれ２通りあるので、N人について、 $n_A$ 本のAアレルを観測したときに、 $n_{AB}$ 人のヘテロを観測する確率は

$P(N_{AB}=n_{AB}|N,n_A)=2^{n_{AB}}￥times ￥frac{￥frac{N!}{n_{AA}!n_{AB}!n_{BB}!}}{￥frac{(2N)!}{n_A!n_B!}}$

$n_{A}$ が奇数のときは、 $n_{B}$ も奇数で、 $n_{AB}$ は奇数。 $n_{AB}$ が偶数である確率は０。逆に、 $n_{A}$ が偶数のときは、 $n_{B}$ も偶数で、 $n_{AB}$ も偶数。 $n_{AB}$ が奇数となる確率は０となることに留意する。

別の説明はこちら
文献はこちら
ExactHWEPを算出するエクセルはこちらの記事で紹介。
ソース



public class Fisher {

	

	public static double hweFisherAbecasis(int[] da){

		int sum = da[0]+da[1] +da[2];

		int allele = da[0]*2 +da[1];

		double p = (double)allele/((double)sum*2);

		double[] genp = new double[3];

		

		genp[0] = p*p;

		genp[1] = 2*p*(1-p);

		genp[2] = (1-p)*(1-p);

		//System.out.println("sum " + sum);

		//System.out.println("allele " + allele);

		//System.out.println("p " + p);

		//System.out.println("genp " + genp[0] + " " + genp[1] + " " + genp[2]);

		double exp0 = genp[0] *sum;

		//System.out.println("exp0 " + exp0);

		boolean less = true;

		if(da[0]>exp0){

			less = false;

		}

		double ret =0 ;

		double total = 0;

		

		

		double combin0 = (logfact(sum)+logfact(allele)+logfact(sum*2-allele))

		-(logfact(da[0]) + logfact(da[1]) + logfact(da[2]) +logfact(sum*2));

		double P0=Math.exp(combin0 +  da[1] * Math.log(2) );

		//System.out.println("P0 " + P0);

		for(int i=0;i<=sum;i++){

			int hetero = allele-i*2;

			if(hetero<0){

				continue;

			}

			int neghomo = sum-i-hetero;

			if(neghomo<0){

				continue;

			}

			//System.out.println("i hetero neghomo " + i + " " + hetero + " " + neghomo);

			double combin = (logfact(sum)+logfact(allele)+logfact(sum*2-allele))

					-(logfact(i) + logfact(hetero) + logfact(neghomo) +logfact(sum*2));

			//System.out.println("combin "+ combin);

			double P1 = Math.exp(combin +  hetero * Math.log(2) );

			if(P0>=P1){

				ret += P1;

			}

			/*

			if(less){

				if(i<=da[0]){

					ret += Math.exp(combin +  hetero * Math.log(2) );

					total +=  Math.exp(combin +  hetero * Math.log(2) );

					

				}else if(i<=exp0){

					total +=  Math.exp(combin +  hetero * Math.log(2) );

					

				}else{

					//break;

				}

				

			}

			if(!less){

				if(i>=da[0]){

					ret += Math.exp(combin +  hetero * Math.log(2) );

					total +=  Math.exp(combin +  hetero * Math.log(2) );

					

				}else if(i>=exp0){

					total +=  Math.exp(combin +  hetero * Math.log(2) );

					

				}else{

					//break;

				}

					

			}

			*/

			//total +=  Math.exp(combin + i * Math.log(genp[0]) + hetero * Math.log(genp[1]) + neghomo * Math.log(genp[2]));

		}

		//System.out.println("ret total "+ ret + " " + total);

		//ret /= total;

		//ret = Math.exp(ret);

		return ret;

	}


		

	public static double logfact(int a){

		double ret =0;

		for(int i=1;i<=a;i++){

			ret += Math.log(i);

		}

		

		return ret;

	}

	public static double logP(int[] da){

		double ret = 0;

		double num = logPnumerator(da);

		double denom = logPdenominator(da);

		ret = num-denom;

		return ret;

	}

	public static double P(int[] da){

		double ret =0;

		double logP = logP(da);

		ret = Math.exp(logP);

		return ret;

	}

	public static double logPdenominator(int[] da){

		double ret =0;

		for(int i=0;i<da.length;i++){

			ret += logfact(da[i]);

		}

		

		

		return ret;

	}

	public static double logPnumerator(int[] da){

		double ret = 0;

		int e = da[0] + da[1];

		int f = da[2] + da[3];

		int g = da[0] + da[2];

		int h = da[1] + da[3];

		int n = e + f;

		ret = logfact(e)+logfact(f)+logfact(g)+logfact(h)-logfact(n);

		

		return ret;

		

	}

}