2007-11-12

CNV検定

CNV多型についての論文メモと、CNV多型のケースコントロール関連検定についての現状をまとめた(記事はこちら)。
おおまかに言って、６、７年前のSNP２ｘ３分割表検定の頃と同じ事情のようだ。
とりあえず、ケースコントロールデータを発表するにあたり、ひとまず、おもいつく分割表検定を行う、というイメージ。
具体的には、個人のジェノタイプを、保有するコピー数に還元(SNPの場合で言うと、２ｘ２分割表検定を作り直す)して検定している。コピー数の多寡で２分しているのも、その一方法であるし、コピー数についてMann-Whitneyをしているのもそう。
SNPについてトレンド検定が適用されるのと同じ考え方(あくまでもディプロタイプ観測数をよりどころにする、ジェノタイプには、適切な重みをつける)というやり方でやってみると同時に、そのような方法と現状の方法との差異についての検討をすることを考慮して、以下のソースを実装。
Mann-Whitneyの検算は群馬大青木先生のサイトのデータで。

package CNV;

public class CNVcalculator {
	/*
	 * アレル単位での最大コピー数をK
	 * 最小コピー数をkとしたとき
	 * Na=K-k+1
	 * Ng=Na*(Na+1)/2
	 * とすること！
	 */

	public int Na; //No. allele
	public int Ng; //No. genotype
	public int[][] t;//contingency table
	//public double[] w;//weight vector
	public double N;//No.samples
	public double R;//No.cases
	public double S;//No.controls
	public double[] n;//number vector
	int df=1;//自由度
	public org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
		new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df);
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ
		int[][] t = {{1,2,3,4,5,6},{2,3,4,5,6,7}};
		//int[][] t = {{5,10,3},{8,9,2}};
		CNVcalculator cnv = new CNVcalculator(t);

		double[] w = {1,2,3,4,5,6};
		double[] cnvTrend = cnv.test(w);
		double[] cnvCNpool = cnv.testCNpool();
		boolean cor=false;
		double[] cnvMWW = cnv.testMWW(cor);
	}
	
	public CNVcalculator(int[][] t_){
		t = StatUtilsX.MiscUtil.DeepCopyInt2(t_);
		Ng=t[0].length;
		
		Na=(int)(Math.round((Math.sqrt(1+8*(double)Ng)-1)/2));
		System.out.println(""+Na);
		N=0;
		R=0;
		S=0;
		n = new double[t[0].length];
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				N+=t[i][j];
				if(i==0){
					R+=t[i][j];
				}else if(i==1){
					S+=t[i][j];
				}
				n[j]+=t[i][j];
			}
		}
		//w = new double[Ng];
	}
	public double[] test(double[] w){
		double[] ret = new double[2];
		if(N==0){
			return ret;
		}
		
		
		double[] p = new double[Ng];
		double P = R/N;
		for(int i=0;i<p.length;i++){
			if(n[i]!=0){
				p[i]=t[0][i]/n[i];
			}
			
		}
		double[] wn = new double[Ng];
		double W = 0;
		for(int i=0;i<wn.length;i++){
			wn[i]=n[i]*w[i];
			W+=wn[i];
		}
		double Wav=W/N;
		double[] wcor=new double[Ng];
		for(int i=0;i<wcor.length;i++){
			wcor[i]=w[i]-Wav;
		}
		double[] rsqn=new double[Ng];
		double Rsqn=0;
		for(int i=0;i<rsqn.length;i++){
			rsqn[i]=t[0][i]*t[0][i]/n[i];
			Rsqn+=rsqn[i];
		}
		double[] numer=new double[Ng];
		double Numer=0;
		for(int i=0;i<numer.length;i++){
			numer[i]=n[i]*(p[i]-P)*wcor[i];
			Numer+=numer[i];
		}
		double[] denom=new double[Ng];
		double Denom=0;
		for(int i=0;i<denom.length;i++){
			denom[i]=n[i]*wcor[i]*wcor[i];
			Denom+=denom[i];
		}
		ret[0]=Numer*Numer/(Denom*P*(1-P));
		/*
		ret[0]=(N*Math.pow((N*(t[0][1]+2*t[0][2])-R*(n[1]+2*n[2])),2))/
		(R*S*(N*(n[1]+4*n[2])-Math.pow((n[1]+2*n[2]),2)));
	*/
		try{
			ret[1]=1-chidf1.cumulativeProbability(ret[0]);
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		//double s=Math.pow(3,2);
		System.out.println(ret[0]+"\t"+ret[1]);
		return ret;
	}
	public double[] testMWW(boolean correct){
		/*
		 * Mann-Whiteney-Willcoxon
		 */
		double[] ret = new double[2];
		int maxCN=(Na-1)*2;
		
		int[][] t2= new int[2][maxCN+1];
		int a1=Na-1;
		int a2=Na-1;
		
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			a2--;
			if(a2<0){
				
				a1--;
				a2=a1;
			}
		}
		int count=0;
		int counted=1;
		double countedNum=S;
		if(R>S){
			count=1;
			counted=0;
			countedNum=R;
		}
		
		double U=0;
		System.out.println("countedNum="+countedNum);
		for(int i=0;i<t2[count].length;i++){
			countedNum-=t2[counted][i];
			System.out.println("i="+i+"\t"+countedNum+"\t"+t2[count][i]);
			U+=t2[count][i]*(countedNum+(double)t2[counted][i]/2);
		}
		double tmpU=R*S-U;
		U=Math.min(U,tmpU);
		double[] tie=new double[t2[0].length];
		double Tsum=0;
		for(int i=0;i<t2[0].length;i++){
			tie[i]=t2[0][i]+t2[1][i];
			Tsum+=tie[i]*tie[i]*tie[i]-tie[i];
		}
		double V=R*S*(N*N*N-N-Tsum)/12/(N*N-N);
		double E=R*S/2;
		double cor = 0;
		if(correct){
			cor=0.5;
		}else{
			
		}
		double Z=(Math.abs(U-E)-cor)/Math.sqrt(V);
		ret[0]=Z;
		System.out.println("U=\t"+U);
		System.out.println("V=\t"+V);
		System.out.println("E=\t"+E);
		org.apache.commons.math.distribution.NormalDistributionImpl norm = 
			new org.apache.commons.math.distribution.NormalDistributionImpl();
		try{
			ret[1]= (1-norm.cumulativeProbability(Z))*2;
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		System.out.println("Z\t"+Z+"\t"+"p\t"+ret[1]);
		return ret;
	}
	public double[] testCNpool(){
		
		int maxCN=(Na-1)*2;
		
		int[][] t2= new int[2][maxCN+1];
		int a1=Na-1;
		int a2=Na-1;
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			a2--;
			if(a2<0){
				
				a1--;
				a2=a1;
			}
		}
		for(int i=0;i<t2[0].length;i++){
			System.out.println(t2[0][i]+"\t"+t2[1][i]);
		}
		CNVcalculator cnv=new CNVcalculator(t2);
		double[] w= new double[t2[0].length];
		w[0]=2*(Na-1);
		for(int i=1;i<w.length;i++){
			w[i]=w[i-1]-1;
		}
		double[] ret = cnv.test(w);
		
		
		
		return ret;
	}

}