metastats2.c

   1 #include <stdio.h>
   2 #include <string.h>
   3 #include <stdlib.h>
   4 #include <time.h>
   5 #include <math.h>
   6 #include "fisher2.h"
   7
   8 void testp(double *permuted_ttests,int *B,double *permuted,double
   9            *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double *ps);
  10 void permute_matrix(double *Imatrix,int *nc,int *nr,double
  11                     *permuted,int *g,double *trial_ts,double *Tinitial,double
  12                     *counter);
  13 void permute_array(int *array, int n);
  14 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,double
  15                        *Ts,double *Tinitial,double *counter1);
  16 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *storage);
  17 void start(double *Imatrix,int *g,int *nr,int *nc,double *testing,
  18                         double storage[][9]);
  19
  20 int main (int argc, char *argv[]){
  21
  22   int col=-1, row=-1, size,g=0,c=0,i=0,j=0,k,counter=0,lines=0, B=10000;
  23   double placeholder=0,min=0, thresh=0.05;
  24
  25   char arr[10000], str[51], a;
  26   char location[41]="jobj.txt", output[41]="out.txt";
  27
  28   for(i=0;i<10000;i++){
  29         arr[i]='q';
  30   }
  31
  32 int u,rflag=0,cflag=0, bflag=0;
  33 char *filename;
  34 int numbers;
  35 double numb;
  36 extern char *optarg;
  37 extern int optind, optopt, opterr;
  38
  39 while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
  40     switch(u) {
  41     case 'r':
  42         numbers = atoi(optarg);
  43         printf("The number of features/rows is %d.\n", numbers);
  44         row = numbers;
  45         rflag = 1;
  46         break;
  47     case 'c':
  48         numbers = atoi(optarg);
  49         printf("The number of samples/columns is %d.\n", numbers);
  50         col = numbers;
  51         cflag = 1;
  52         break;
  53     case 'g':
  54         numbers = atoi(optarg);
  55         printf("Your g-value is %d.\n", numbers);
  56         g = numbers;
  57         break;
  58     case 'b':
  59         numbers = atoi(optarg);
  60         printf("The number of permutations is %d\n", numbers);
  61         B = numbers;
  62         break;
  63     case 't':
  64         numb = atof(optarg);
  65         printf("Threshold is is %lf\n", numb);
  66         thresh = numb;
  67         break;
  68     case 'f':
  69         filename = optarg;
  70         printf("filename input is %s\n", filename);
  71         strcpy(location,filename);
  72         break;
  73     case 'o':
  74         filename = optarg;
  75         printf("filename output %s\n", filename);
  76         strcpy(output,filename);
  77         break;
  78     case ':':
  79         printf("-%c without filename\n", optopt);
  80         break;
  81     case '?':
  82         printf("unknown arg %c\n", optopt);
  83         break;
  84     }
  85 }
  86
  87   FILE *jobj, *out;
  88   jobj=fopen(location,"r");
  89
  90   if(jobj == NULL){
  91     printf("Please don't forget to save your matrix in the active");
  92     printf(" directory as \"%s\".\n",location);
  93     return 0;
  94   }
  95
  96   // Gets the first line of samples names and checks for user error.
  97   fgets(arr,10000,jobj);
  98
  99   for(i=0;i<10000;i++){
 100     if(isspace(arr[i])){
 101       counter++; }
 102   }
 103
 104   if (cflag == 0) {
 105       printf("You didn't tell us how many subjects there are!\n");
 106       printf("But we'll still do the analysis as if there are %d subjects.\n\n",col=counter-1);
 107   }
 108   if (cflag == 1) {
 109         if (col != counter-1){
 110           printf("We would expect %d subjects, but you said %d.\n",counter-1,col);
 111     }
 112   }
 113
 114   while((a = fgetc(jobj)) != EOF){
 115     if(a == '\n'){
 116       lines++; }
 117   }
 118
 119   if (rflag == 0) {
 120       printf("You didn't specify the number of features!\n");
 121       printf("But we'll still do the analysis assuming %d features.\n\n", row=lines-1);
 122   }
 123   if (rflag == 1) {
 124     if ( lines-1 != row ){
 125      printf("We would expect %d features, but you said %d.\n",lines-1,row);
 126     }
 127   }
 128
 129   if (g>=col || g<=0){
 130         printf("Check your g value\n");
 131   }
 132
 133   // Initialize the matrices
 134   size = row*col;
 135   double matrix[row][col];
 136   double pmatrix[size],pmatrix2[size],permuted[size];
 137   double storage[row][9];
 138
 139   for (i=0;i<row;i++){
 140         for (j =0;j<9;j++){
 141       storage[i][j]=0;
 142         }
 143   }
 144   // Reset file below and create a separate matrix.
 145   rewind(jobj);
 146   fgets(arr,10000,jobj);
 147
 148   for(i=0; i<row; i++){
 149     fscanf(jobj,"%s",str);
 150     for(j=0; j<col;j++){
 151       fscanf(jobj,"%lf",&placeholder);
 152       matrix[i][j]=placeholder;
 153       if(isalnum(placeholder)!=0){ // check for ""
 154         printf("Your matrix isn't set up properly!\n");
 155         return 0;
 156                 }
 157       pmatrix[c]=0; // initializing to zero
 158       permuted[c]=0;
 159       c++;
 160     }
 161   }
 162
 163   fclose(jobj);
 164
 165
 166   // Produces the sum of each column
 167   double total[col],total1=0,total2=0;
 168   double ratio[col];
 169
 170   for(i=0;i<col;i++){
 171     total[i]=0;
 172     ratio[i]=0; }
 173
 174   for(i=0; i<col; i++){
 175     for(j=0;j<row;j++){
 176       total[i]=total[i]+matrix[j][i];
 177     }
 178   }
 179
 180   for(i=0;i<g-1;i++){
 181         total1=total1+total[i];}
 182
 183   for(i=g-1;i<col;i++){
 184     total2=total2+total[i];}
 185
 186
 187   // Creates the ratios by first finding the minimum of totals
 188   min = total[0];
 189   if (col==2){
 190   if (total[0]<total[1]){
 191         min = total[1];}
 192   }
 193   if (col >2){
 194   for(i=1;i<col;i++){
 195     if (min > total[i]){
 196       min = total[i];}
 197   }
 198   }
 199   if (min<=0){
 200     printf("Error, the sum of one of the columns <= 0.");
 201     return 0;
 202   }
 203
 204
 205   // Ratio time...
 206   for(i=0;i<col;i++){
 207     ratio[i]=total[i]/min;
 208   }
 209
 210   // Change matrix into an array as received by R for compatibility.
 211
 212   c=0;
 213   for(i=0;i<col;i++){
 214     for(j=0;j<row;j++){
 215       pmatrix[c]=matrix[j][i];
 216       c++;
 217     }
 218   }
 219
 220   if(row == 1){
 221         for (i =0; i<col;i++){
 222                 pmatrix[i]=pmatrix[i]/ratio[i];
 223         }
 224   }
 225   else {
 226     counter = 0;
 227     j=-1;
 228     for (i=0; i<size; i++) {
 229       if (counter % row == 0) {
 230         j++;
 231       }
 232       pmatrix[i]=pmatrix[i]/ratio[j];
 233       counter++;
 234     }
 235   }
 236   // pass everything to the rest of the code using pointers. then
 237   // write to output file. below pointers for most of the values are
 238   // created to send everything by reference.
 239
 240   int ptt_size, *permutes,*nc,*nr,*gvalue;
 241
 242   nc = &col;
 243   nr = &row;
 244   gvalue = &g;
 245
 246   permutes = &B;
 247   ptt_size = B*row;
 248
 249   //changing ptt_size to row
 250   double permuted_ttests[row], pvalues[row], tinitial[row];
 251
 252   for(i=0;i<row;i++){
 253     permuted_ttests[i]=0;}
 254
 255   for(i=0;i<row;i++){
 256     pvalues[i]=0;
 257     tinitial[i]=0; }
 258
 259   // Find the initial values for the matrix.
 260   start(pmatrix,gvalue,nr,nc,tinitial,storage);
 261
 262   // Start the calculations.
 263
 264   if ( (col==2) || ((g-1)<8) || ((col-g+1) < 8) ){
 265
 266   double fish[row], fish2[row];
 267   for(i=0;i<row;i++){
 268         fish[i]=0;
 269         fish2[i]=0;}
 270
 271   for(i=0;i<row;i++){
 272
 273         for(j=0;j<g-1;j++){
 274                 fish[i]=fish[i]+matrix[i][j];
 275         }
 276
 277         for(j=g-1;j<col;j++){
 278                 fish2[i]=fish2[i]+matrix[i][j];
 279         }
 280
 281         double  f11,f12,f21,f22;
 282
 283         f11=fish[i];
 284         f12=fish2[i];
 285
 286         f21=total1-f11;
 287         f22=total2-f12;
 288
 289         double data[] = {f11, f12, f21, f22};
 290
 291         // CONTINGENGCY TABLE:
 292         //   f11   f12
 293         //   f21   f22
 294
 295         int *nr, *nc, *ldtabl, *work;
 296         int nrow=2, ncol=2, ldtable=2, workspace=100000;
 297         double *expect, *prc, *emin,*prt,*pre;
 298         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 299
 300         nr = &nrow;
 301         nc = &ncol;
 302         ldtabl=&ldtable;
 303         work = &workspace;
 304
 305         expect = &e;
 306         prc = &prc1;
 307         emin=&emin1;
 308         prt=&prt1;
 309         pre=&pre1;
 310
 311         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 312
 313         if (*pre>.999999999){
 314                 *pre=1;
 315         }
 316         storage[i][8] = *pre;
 317         pvalues[i]=*pre;
 318     }
 319   }
 320   else{
 321
 322   testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
 323
 324                // Checks to make sure the matrix isn't sparse.
 325   double sparse[row], sparse2[row];
 326   for(i=0;i<row;i++){
 327         sparse[i]=0;
 328         sparse2[i]=0;}
 329
 330   c=0;
 331   for(i=0;i<row;i++){
 332
 333         for(j=0;j<g-1;j++){
 334                 sparse[i]=sparse[i]+matrix[i][j];
 335         }
 336
 337         if(sparse[i] < (double)(g-1)){
 338                 c++;
 339         }
 340         for(j=g-1;j<col;j++){ // ?<= for col
 341                 sparse2[i]=sparse2[i]+matrix[i][j];
 342         }
 343
 344         if( (sparse2[i] <(double)(col-g+1))) {
 345                 c++;
 346         }
 347
 348         if (c==2){
 349                 c=0;
 350
 351         double  f11,f12,f21,f22;
 352
 353         f11=sparse[i];
 354         sparse[i]=0;
 355
 356         f12=sparse2[i];
 357         sparse2[i]=0;
 358
 359         f21=total1-f11;
 360         f22=total2-f12;
 361
 362         double data[] = {f11, f12, f21, f22};
 363
 364         int *nr, *nc, *ldtabl, *work;
 365         int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
 366         double *expect, *prc, *emin,*prt,*pre;
 367         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 368
 369         nr = &nrow;
 370         nc = &ncol;
 371         ldtabl=&ldtable;
 372         work = &workspace;
 373
 374         expect = &e;
 375         prc = &prc1;
 376         emin=&emin1;
 377         prt=&prt1;
 378         pre=&pre1;
 379
 380         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 381
 382         if (*pre>.999999999){
 383                 *pre=1;
 384         }
 385         storage[i][8] = *pre;
 386         pvalues[i]=*pre;
 387     }
 388   }
 389   // End of else statement
 390   bflag = 1;
 391   }
 392
 393   // Calculates the mean of counts (not normalized)
 394   double temp[row][2];
 395
 396   for(j=0;j<row;j++){
 397         for(i=0;i<2;i++){
 398                 temp[j][i]=0;
 399         }
 400   }
 401
 402   for (j=0;j<row;j++){
 403         for (i=1; i<=(g-1); i++){
 404                 temp[j][0]=temp[j][0]+matrix[j][i-1];
 405         }
 406         temp[j][0]= (double) temp[j][0]/(g-1);
 407         for(i=g;i<=col;i++){
 408                 temp[j][1]=temp[j][1]+matrix[j][i-1];
 409         }
 410         temp[j][1]= (double) temp[j][1]/(col-g+1);
 411   }
 412
 413   for(i=0;i<row;i++){
 414         storage[i][3]=temp[i][0];
 415         storage[i][7]=temp[i][1];
 416         storage[i][8]=pvalues[i];
 417   }
 418
 419 // BACKUP checks
 420
 421   for (i=0;i<row;i++){
 422     if(pvalues[i]<thresh){
 423         printf("Feature %d is significant, p = %.10lf \n",i+1,pvalues[i]);
 424         }
 425   }
 426
 427   // And now we write the files to a text file.
 428   struct tm *local;
 429   time_t t;
 430   t = time(NULL);
 431   local = localtime(&t);
 432
 433   jobj= fopen(location,"r");
 434   fgets(arr,10000,jobj);
 435
 436   out = fopen(output,"a+");
 437
 438   fprintf(out,"Local time and date of test: %s\n", asctime(local));
 439   fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
 440   if (bflag == 1){
 441     fprintf(out,"%d permutations\n\n",B);
 442   }
 443   for(i=0; i<row; i++){
 444     fscanf(jobj,"%s",str);
 445         fprintf(out,"%s",str);
 446         for(k=0;k<col;k++){
 447                 fscanf(jobj,"%*lf",&placeholder);
 448         }
 449     for(j=0; j<9;j++){
 450       fprintf(out,"\t%.12lf",storage[i][j]);
 451     }
 452     fprintf(out,"\n");
 453   }
 454
 455   fprintf(out,"\n \n");
 456
 457   fclose(jobj);
 458   fclose(out);
 459
 460   return 0;
 461 }
 462
 463 void testp(double *permuted_ttests,int *B,double *permuted,
 464            double *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double
 465            *ps) {
 466
 467   double Tvalues[*nr];
 468   int a, b, n, i, j,k=0;
 469
 470   a = *B;
 471   b = *nr;
 472   n = a*b;
 473
 474   double counter[b];
 475
 476   for(j=0;j<b;j++){
 477     counter[j]=0;
 478   }
 479
 480   for (j=1; j<=*B; j++){
 481     permute_matrix(Imatrix,nc,nr,permuted,g,Tvalues,Tinitial,counter);
 482    // for(i=0;i<*nr;i++){
 483    //   permuted_ttests[k]=fabs(Tvalues[i]);
 484         //    k++;
 485     }
 486
 487
 488   for(j=0;j<*nr;j++){
 489     ps[j]=((counter[j]+1)/(double)(a+1));
 490   }
 491 }
 492
 493 void permute_matrix(double *Imatrix,int *nc,int *nr,double *permuted,
 494                     int *g,double *trial_ts,double *Tinitial,double *counter1){
 495
 496   int i=0,j=0,n=0,a=0,b=0,f=0,c=0,k=0;
 497
 498   a = *nr; // number of rows
 499   b = *nc;
 500   n = a*b;
 501
 502   int y[b];
 503
 504   for (i=1; i<=*nc; i++){
 505     y[i-1] = i;
 506   }
 507
 508   permute_array(y, b);
 509
 510   for (i=0; i<*nc; i++){
 511     f = y[i]; //column number
 512     c=1;
 513     c*=(f-1);
 514     c*=a;
 515     if (f == 1){
 516       c = 0;
 517     } // starting value position in the Imatrix
 518     for(j=1; j<=*nr; j++){
 519       permuted[k] = Imatrix[c];
 520       c++;
 521       k++;
 522     }
 523   }
 524
 525   calc_twosample_ts(permuted,g,nc,nr,trial_ts,Tinitial,counter1);
 526 }
 527
 528 void permute_array(int *array, int n) {
 529   static int seeded = 0;
 530   int i;
 531
 532   if (! seeded) {
 533     seeded = 1;
 534     srandom(time(NULL));
 535   }
 536
 537   for (i = 0; i < n; i++) {
 538     int selection = random() % (n - i);
 539     int tmp = array[i + selection];
 540     array[i + selection] = array[i];
 541     array[i] = tmp;
 542   }
 543 }
 544
 545 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,
 546                        double *Ts,double *Tinitial,double *counter) {
 547   int i,a;
 548   a = *nr;
 549   a*=4;
 550
 551   double C1[*nr][3], C2[*nr][3], storage[a],tool[a];
 552   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 553
 554   nrows = (double) *nr;
 555   ncols = (double) *nc;
 556   gvalue= (double) *g;
 557
 558     meanvar(Pmatrix,g,nr,nc,storage);
 559     for(i=0;i<=a-1;i++){
 560       tool[i]=storage[i];
 561     }
 562     for (i=0; i<*nr;i++){
 563       C1[i][0]=tool[i];
 564       C1[i][1]=tool[i+*nr+*nr];
 565       C1[i][2]=C1[i][1]/(gvalue-1);
 566
 567       C2[i][0]=tool[i+*nr];
 568       C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
 569       C2[i][2]=C2[i][1]/(ncols-gvalue+1);
 570     }
 571
 572     for (i=0; i<*nr; i++){
 573       xbardiff = C1[i][0]-C2[i][0];
 574       denom = sqrt(C1[i][2]+C2[i][2]);
 575       Ts[i]=fabs(xbardiff/denom);
 576       if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
 577             counter[i]++;
 578       }
 579     }
 580 }
 581
 582 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
 583   double temp[*nr], temp2[*nr],var[*nr],var2[*nr],a,b;
 584
 585   int i,m,k,l,n;
 586
 587   a = (double) *g-1;
 588   b = (double) (*nc-a);
 589
 590   for (i = 0; i<*nr; i++){
 591     temp[i]=0;
 592     temp2[i]=0;
 593     var[i]=0;
 594     var2[i]=0;
 595   }
 596
 597   k = *nr; // number of rows
 598   l = *nc;
 599   n = k*l;
 600
 601     m=0;
 602     m=*g-1;
 603     k=*nr;
 604     m*=k; // m = g * nr now
 605     for (i=0;i<m;i++){
 606       temp[i%k]=temp[i%k]+pmatrix[i];
 607     }
 608     for (i=0;i<n;i++){
 609       temp2[i%k]=temp2[i%k]+pmatrix[i];
 610     }
 611     for (i=0;i<*nr;i++){
 612       temp2[i]=temp2[i]-temp[i];
 613     }
 614     for (i=0;i<=*nr-1;i++){
 615       store[i]=temp[i]/a;
 616       store[i+*nr]=temp2[i]/b;
 617     }
 618
 619     // That completes the mean calculations.
 620
 621     for (i=0;i<m;i++){
 622       var[i%k]=var[i%k]+pow((pmatrix[i]-store[i%k]),2);
 623     }
 624     for (i=m;i<n;i++){
 625       var2[i%k]=var2[i%k]+pow((pmatrix[i]-store[(i%k)+*nr]),2);
 626     }
 627
 628     for (i=0;i<=*nr-1;i++){
 629       store[i+2*k]=var[i]/(a-1);
 630       store[i+3*k]=var2[i]/(b-1);
 631     }
 632     // That completes var calculations.
 633 }
 634
 635 void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
 636                 double storage[][9]){
 637   int i, a = *nr;
 638   a*=4;
 639
 640   double store[a], tool[a], C1[*nr][3], C2[*nr][3];
 641   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 642
 643   nrows = (double) *nr;
 644   ncols = (double) *nc;
 645   gvalue= (double) *g;
 646
 647   meanvar(Imatrix,g,nr,nc,store);
 648
 649   for(i=0;i<=a-1;i++){
 650     tool[i]=store[i];
 651   }
 652   for (i=0; i<*nr;i++){
 653     C1[i][0]=tool[i]; //mean group 1
 654     storage[i][0]=C1[i][0];
 655     C1[i][1]=tool[i+*nr+*nr]; // var group 1
 656     storage[i][1]=C1[i][1];
 657     C1[i][2]=C1[i][1]/(gvalue-1);
 658     storage[i][2]=sqrt(C1[i][2]);
 659
 660     C2[i][0]=tool[i+*nr]; // mean group 2
 661     storage[i][4]=C2[i][0];
 662     C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
 663     storage[i][5]=C2[i][1];
 664     C2[i][2]=C2[i][1]/(ncols-gvalue+1);
 665     storage[i][6]=sqrt(C2[i][2]);
 666   }
 667   for (i=0; i<*nr; i++){
 668     xbardiff = C1[i][0]-C2[i][0];
 669     denom = sqrt(C1[i][2]+C2[i][2]);
 670     initial[i]=fabs(xbardiff/denom);
 671   }
 672 }