#include <stdio.h> #include <stdlib.h> #include <malloc.h> #include <time.h> unsigned long*Table,*ColSum,*RowSum,LowBound,UpBound,*TempTable; unsigned long CumColSum[2],CumRowSum[2]; unsigned long TTable[2][2],TColSum[2],TRowSum[2]; unsigned long Sum,temp,TSum,Size,Total; int TCol,TRow,Ind_C,Ind_R,Col,Row; double CHstats,Hstats,temp1; double Mvol,Hvol,Dvol,CHstats,Hstats; double Mvol_size,Hvol_size,Dvol_size; double Weight; unsigned long generate() { unsigned long int i; unsigned long int intu; intu=(65539*rand())%Total; for(i=0;i<1;i++) { intu=(65539*intu)%Total; } return (unsigned long)((double)intu/Total*(UpBound-LowBound+1)+LowBound); } 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; } //Evaluate the chi-square statistics double chisq(unsigned long*table) { int i,j; double chis2=0,tt; for(i=0;i<Row;i++) for(j=0;j<Col;j++) { tt=table[i*Col+j]-(double)RowSum[i]*ColSum[j]/(double)Sum; chis2+=tt*tt*(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) { int i,j; double ht=0,ht1=0,ht2=0; for(i=0;i<Col;i++) ht2+=ColSum[i]*ColSum[i]; for(i=0;i<Row;i++) { ht1+=RowSum[i]*RowSum[i]; for(j=0;j<Col;j++) { ht+=table[i*Col+j]*table[i*Col+j]; } } return ht-ht1*ht2/Sum/Sum; } //Build up the tree and Evaluate the statistics void FillIn() { if(TCol<Col-1&&TRow<Row-1) { if(TRow==0) { CumColSum[0]=0; CumColSum[1]=0; } else { for(Ind_R=0;Ind_R<TRow;Ind_R++) { if(Ind_R==0){CumColSum[0]=TempTable[TCol];} else{CumColSum[0]+=TempTable[Ind_R*Col+TCol];} temp=0; for(Ind_C=0;Ind_C<=TCol-1;Ind_C++) temp+=TempTable[Ind_R*Col+Ind_C]; if(Ind_R==0)CumColSum[1]=temp; else CumColSum[1]+=temp; if(Ind_R==0)TSum=RowSum[Ind_R]; else TSum+=RowSum[Ind_R]; } CumColSum[1]=TSum-CumColSum[1]; } for(Ind_C=TCol;Ind_C<Col;Ind_C++) { if(Ind_C==TCol) TColSum[1]=ColSum[Ind_C]; else TColSum[1]+=ColSum[Ind_C]; } TSum=TColSum[1]-CumColSum[1]; if(TCol==0) { CumRowSum[0]=0; } else { for(Ind_C=0;Ind_C<TCol;Ind_C++) { if(Ind_C==0){CumRowSum[0]=TempTable[TRow*Col];} else{CumRowSum[0]+=TempTable[TRow*Col+Ind_C];} } } TColSum[0]=ColSum[TCol]-CumColSum[0]; TRowSum[0]=RowSum[TRow]-CumRowSum[0]; TColSum[1]=TSum-TColSum[0]; TRowSum[1]=TSum-TRowSum[0]; LowBound=Max(0,TColSum[0]-TRowSum[1]); UpBound=Min(TColSum[0],TRowSum[0]); Weight=Weight*(UpBound-LowBound+1); TempTable[TRow*Col+TCol]=generate(); } else if(TRow<Row-1) { temp=0; for(Ind_C=0;Ind_C<TCol;Ind_C++) { temp+=TempTable[TRow*Col+Ind_C]; } LowBound=UpBound=RowSum[TRow]-temp; TempTable[TRow*Col+TCol]=LowBound; } else { temp=0; for(Ind_R=0;Ind_R<TRow;Ind_R++) { temp+=TempTable[Ind_R*Col+TCol]; } LowBound=UpBound=ColSum[TCol]-temp; TempTable[TRow*Col+TCol]=LowBound; if((TRow==Row-1)&&(TCol==Col-1)) { temp1=chisq(&TempTable[0]); /*M-Volumn Test*/ Mvol_size+=Weight; if(temp1<CHstats) Mvol+=Weight; /*D-Volumn Test*/ if((Table[0]-ColSum[0]*RowSum[0]/(double)Sum)*(TempTable[0]-ColSum[0]*RowSum[0]/(double)Sum)>=0) { Dvol_size+=Weight; if(temp1<CHstats) Dvol+=Weight; } /*H-Volumn Test*/ temp1=H_T(&TempTable[0]); if(temp1*Hstats>=0) { Hvol_size+=Weight; if(sign(temp1)*temp1<sign(Hstats)*Hstats) Hvol+=Weight; } } } } void evaluate() { int i,j,k,Limit; unsigned long l; /*Evaluate the Chi-square statistic of observed table*/ CHstats=chisq(&Table[0]); /*Evaluate the H statistics of observed table*/ Hstats=H_T(&Table[0]); srand((unsigned)time(NULL)); Limit=Min(Row-1,Col-1); for(l=0;l<Size;l++) { for(k=0;k<Limit;k++) { for(i=k;i<Col;i++) { if(i==0&&k==0) { TRowSum[0]=RowSum[0]; TColSum[0]=ColSum[0]; TRowSum[1]=Sum-RowSum[0]; TColSum[1]=Sum-ColSum[0]; LowBound=Max(0,TColSum[0]-TRowSum[1]); UpBound=Min(TColSum[0],TRowSum[0]); TRow=0; TCol=0; Weight=1; FillIn(); } else { TRow=k; TCol=i; FillIn(); } } for(j=k+1;j<Row;j++) { TRow=j; TCol=k; FillIn(); } } if(Row<Col) { TRow=Row-1; for(i=Row-1;i<Col;i++) { TCol=i; FillIn(); } } else { TCol=Col-1; for(i=Col-1;i<Row;i++) { TRow=i; FillIn(); } } } Mvol=Mvol/Mvol_size; Dvol=Dvol/Dvol_size; Hvol=sign(Hstats)*Hvol/Hvol_size; } main() { FILE*in,*out; char Input[256],Output[256]; int i,j,cell,length, flag; unsigned long control; /*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); } printf("\nHow many Monte Carlo samples?\n"); scanf("%ld",&Size); fscanf(in,"%d %d %d",&Row,&Col, &cell); for(i=0;i<Row;i++) { for(j=0;j<Col;j++) { fscanf(in,"%ld",&control); } } length=0; do { fscanf(in,"%ld",&control); if(control==0) length++; }while(control==0); length+=Row*Col; rewind(in); fprintf(out,"%s","Mvol Hvol\n"); Total=1; for(i=0;i<31;i++) Total=Total*2; while(!feof(in)) { Mvol=Hvol=Dvol=0; Mvol_size=Hvol_size=Dvol_size=0; Sum=0; /*Enter the dimensions information of the matrix*/ fscanf(in,"%d %d %d",&Row,&Col,&cell); Table=(unsigned long*)malloc(Row*Col*sizeof(unsigned long)); TempTable=(unsigned long*)malloc(Row*Col*sizeof(unsigned long)); ColSum=(unsigned long*)malloc(Col*sizeof(unsigned long)); RowSum=(unsigned long*)malloc(Row*sizeof(unsigned long)); /*Input the Table*/ for(i=0;i<2;i++) { CumRowSum[i]=0; CumColSum[i]=0; for(j=0;j<2;j++) TTable[i][j]=0; TColSum[i]=0; TRowSum[i]=0; } for(i=0;i<Row;i++) { for(j=0;j<Col;j++) { TempTable[i*Col+j]=0; fscanf(in,"%ld",&Table[i*Col+j]); Sum+=Table[i*Col+j]; if(j==0)RowSum[i]=Table[i*Col+j]; else RowSum[i]+=Table[i*Col+j]; if(i==0)ColSum[j]=Table[i*Col+j]; else ColSum[j]+=Table[i*Col+j]; } } for(i=0;i<length-cell;i++) fscanf(in,"%ld",&control); flag=0; for(i=0;i<Row;i++) { if (RowSum[i] == 0) { flag=1; break; } } if (flag == 0) { for(j=0;j<Col;j++) { if (ColSum[j] == 0) { flag=1; break; } } } if (flag == 1) { printf("Error message: Not all row and column sums are positive. Please remove rows and columns with zero sums.\n"); break; } evaluate(); fprintf(out,"%lf %lf\n",Mvol,Hvol); free(Table); free(TempTable); free(ColSum); free(RowSum); fscanf(in,"\n"); } fclose(out); fclose(in); }