#include <stdio.h> #include <math.h> #define dim 2 /*Return the minimum value*/ unsigned long max(long x, long y) { if(x<y)return (unsigned long)y; else return (unsigned long)x; } /*Return the maximum value*/ unsigned long min(unsigned long x,unsigned long y) { if(x<y)return x; else return y; } /*Chisquare statistics*/ double chisq(unsigned long*table,unsigned long*ColSum,unsigned long*RowSum,unsigned long sum) { int i,j; double chis2=0; for(i=0;i<dim;i++) for(j=0;j<dim;j++) { chis2+=pow(table[i*dim+j]-(double)RowSum[i]*ColSum[j]/(double)sum,2)*(double)sum/(RowSum[i]*ColSum[j]); } return chis2; } /*Sign function*/ int sign(double x) { if(x>0)return 1; else return -1; } /*Excess of homozygosity*/ double H_T(unsigned long*table,unsigned long*ColSum,unsigned long*RowSum,unsigned long sum) { int i,j; double ht=0,ht1=0,ht2=0; for(i=0;i<dim;i++) { ht1+=RowSum[i]*RowSum[i]; ht2+=ColSum[i]*ColSum[i]; for(j=0;j<dim;j++) { ht+=table[i*dim+j]*table[i*dim+j]; } } return ht-ht1*ht2/sum/sum; } main() { unsigned long Table[dim][dim],ColSum[dim],RowSum[dim],LowBound,UpBound,n; unsigned long TTable[dim][dim],sum; double Mvol=0,Dvol=0,Hvol=0,CHstats,Hstats,temp; double Mvol_size=0,Dvol_size=0,Hvol_size=0; int i,j, flag; FILE*in,*out; char Input[256],Output[256]; /*Input the Table*/ printf("\nPlease enter the data file name (including the folders,i.e. c:\\user\\table.txt):\n"); scanf("%s",Input); printf("Please enter the output file name (including the folders,i.e. c:\\user\\output.txt):\n"); scanf("%s",Output); while((in=fopen(Input,"r"))==NULL) { printf("\nPlease enter the data file name (including the folders,i.e. c:\\user\\table.txt):\n"); scanf("%s",Input); } while((out=fopen(Output,"w"))==NULL) { printf("Please enter the output file name (including the folders,i.e. c:\\user\\output.txt):\n"); scanf("%s",Output); } fprintf(out,"%s","Dvol Mvol Hvol\n"); while(!feof(in)) { Mvol=Dvol=Hvol=0; Mvol_size=Dvol_size=Hvol_size=0; sum=0; for(i=0;i<dim;i++) { for(j=0;j<dim;j++) { fscanf(in,"%ld",&Table[i][j]); sum+=Table[i][j]; if(j==0)RowSum[i]=Table[i][j]; else RowSum[i]+=Table[i][j]; if(i==0)ColSum[j]=Table[i][j]; else ColSum[j]+=Table[i][j]; } } flag=0; for(i=0;i<dim;i++) { if (RowSum[i] == 0) { flag=1; break; } } if (flag == 0) { for(j=0;j<dim;j++) { if (ColSum[j] == 0) { flag=1; break; } } } if (flag == 1) { printf("Error message: Now all row and column sums are positive. Please remove rows and columns with zero sums.\n"); break; } LowBound=max(0,ColSum[0]-RowSum[1]); UpBound=min(ColSum[0],RowSum[0]); /*Evaluate the Chi-square statistic of observed table*/ CHstats=chisq(&Table[0][0],&ColSum[0],&RowSum[0],sum); /*Evaluate the H statistics of observed table*/ Hstats=H_T(&Table[0][0],&ColSum[0],&RowSum[0],sum); for(n=LowBound;n<=UpBound;n++) { for(i=0;i<dim;i++) for(j=0;j<dim;j++) { if(i==0&&j==0)TTable[i][j]=n; else if(i==0)TTable[i][j]=RowSum[i]-n; else if(i==1)TTable[i][j]=ColSum[j]-TTable[i-1][j]; } temp=chisq(&TTable[0][0],&ColSum[0],&RowSum[0],sum); /*M-Volumn Test*/ if(temp<CHstats) Mvol+=1; /*D-Volumn Test*/ if((Table[0][0]-ColSum[0]*RowSum[0]/(double)sum)*(TTable[0][0]-ColSum[0]*RowSum[0]/(double)sum)>=0) { Dvol_size+=1; if(temp<CHstats) Dvol+=1; } /*H-Volumn Test*/ temp=H_T(&TTable[0][0],&ColSum[0],&RowSum[0],sum); if(temp*Hstats>=0) { Hvol_size+=1; if(sign(temp)*temp<sign(Hstats)*Hstats) Hvol+=1; } } Mvol_size=UpBound-LowBound+1; Mvol=Mvol/Mvol_size; Dvol=Dvol/Dvol_size; Hvol=sign(Hstats)*Hvol/Hvol_size; fprintf(out,"%lf %lf %lf\n",Dvol,Mvol,Hvol); fscanf(in,"\n"); } fclose(out); fclose(in); }