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

```
# thetaの計算

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

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;
}

```

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]);

}```