2005-09-07

DNA配列のNucleotide Polymorphism (theta)とNucleotide Diversity (pi)

多型

長さk bp, 多型箇所数 s, 多型性配列数 nに対し、この領域がどのくらい多型に富んでいるか、塩基構成がどのくらい多彩かをそれぞれtheta, piで現す。thetaは長さkあたりに存在する多型箇所数 s(多型の密度s/n)と、多型性配列数nの関数となる。密度が大きいほど、thetaは大きく、その多型箇所の確認するのに要した多型性配列数が少ないほど大きくなる。他方、塩基構成の多彩さの度合いであるpiは、多型箇所の密度とそれをもたらす多型性配列数のみならず、個々の多型のアレル頻度の影響も加味した指標である。アレル頻度が50-50に近い多型は、アレル頻度が低い多型に較べて、piをより大きくする。



# thetaの計算

係数　a1 = sum(from 1 to n-1) 1/i;

係数　a2 = sum(from 1 to n-1) 1/i^2;

theta = (s/n)/a1;

Var(theta) = theta/(k*a1) + a2*theta^2/(a1^2);
# piの計算

係数　b1 = (n+1)/(3*(n-1));

係数　b2 = 2(n^2+n+3)/(9(n(n-1));

pi = sum(n配列の作る2配列ペア)(異なる多型箇所数)/(n配列の作る2配列ペア数*k);

Var(pi) = b1*pi/k+b2*pi^2;

これらを計算するJavaソースは以下の通り



//配列長 length, 多型性配列数 numhap, 多型箇所行列 haplotype[][]

	public static double[] NucDiversity(int length,int numhap,char[][] haplotype){

		double ret[];

		ret = new double[4];

		double degPoly=0;//theta

		double varDegPoly=0;//variance of theta

		double nucDiversity=0;//pi

		double varNucDiversity=0;//variance of pi

		

		int div[] = new int[haplotype[0].length];

		int numseg=0;

		

		for(int i=0;i<haplotype[0].length;i++){

			div[i]=0;

			for(int j=0;j<numhap;j++){

				for(int k=j+1;k<numhap;k++){

					if(haplotype[j][i]!=haplotype[k][i]){

						div[i]++;

					}

				}

			}

		}

		for(int i=0;i<haplotype[0].length;i++){

			if(div[i]>0){

				numseg++;

				nucDiversity+=div[i];

			}

		}

		

		double a1=0;

		double a2=0;

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

			a1+=(double)1/i;

			a2+=(double)1/(i*i);

		}

		

		degPoly = numseg/(length*a1);

		varDegPoly = degPoly/(double)(length*a1)+a2*Math.pow(degPoly,2.0)/Math.pow(a1,2.0);		

		nucDiversity/=(numhap*(numhap-1)/2*length);

		double b1=0;

		double b2=0;

		b1=(double)(numhap+1)/(3*(numhap-1));

		b2=(double)2*(numhap*numhap+numhap+3)/(9*numhap*(numhap-1));

		varNucDiversity = b1*nucDiversity/length+b2*nucDiversity*nucDiversity;

		ret[0]=degPoly;

		ret[1]=varDegPoly;

		ret[2]=nucDiversity;

		ret[3]=varNucDiversity;

		return ret;

	}

実行結果サンプル：長さ500塩基対、多型箇所 16、多型性配列５本のサンプルでの上記実行結果を示す。

theta 0.015360000000000002

var-theta 9.213050880000004E-5

pi 0.0158

var-pi 1.0733466666666668E-4

実行に用いた配列については以下の通り



	public static void main(String[] args) {

		int snp[]={132,142,162,192,198,201,207,240,246,351,354,372,375,405,417,483};

		int length = 500;

		int numhap=5;

		char seq[][];

		seq = new char[numhap][snp.length];

		

		char seq0[]={'T','C','T','A','C','C','T','C','C','T','C','G','G','T','T','A'};

		char seq1[]={'T','C','C','T','A','C','C','T','C','C','T','G','G','T','T','T'};

		char seq2[]={'C','T','C','C','C','C','C','T','C','T','T','T','G','C','T','A'};

		char seq3[]={'C','T','C','C','C','C','C','T','T','C','T','G','A','C','T','T'};

		char seq4[]={'C','T','C','C','C','T','C','T','T','T','T','G','G','C','C','A'};

		

		seq[0]=seq0;

		seq[1]=seq1;

		seq[2]=seq2;

		seq[3]=seq3;

		seq[4]=seq4;

		

		double diversity[] = PopGenetTools.NucDiversity(length,numhap,seq);

		

		System.out.println("theta var-theta pi var-pi " + diversity[0] + " "

				+ diversity[1] + " " + diversity[2] + " " + diversity[3]);

		

	}