metastats2.c

   1
   2 #include "metastats.h"
   3
   4 int metastat_main (char* outputFileName, int numRows, int numCols, double threshold, int numPermutations, double** data, int secondGroupingStart){
   5
   6   int size,c=0,i=0,j=0,k,counter=0, bflag=0;
   7   int B=numPermutations;
   8   int row = numRows;
   9   int col = numCols;
  10   int g = secondGroupingStart;
  11   double thresh=threshold;
  12   double placeholder=0,min=0;
  13
  14   char output[1024];
  15   strcpy(output, outputFileName);
  16   FILE *out;
  17
  18   if (g>=col || g<=0){
  19         printf("Check your g value\n");
  20   }
  21
  22   // Initialize the matrices
  23   size = row*col;
  24   double matrix[row][col];
  25   double pmatrix[size],pmatrix2[size],permuted[size];
  26   double storage[row][9];
  27
  28   for (i=0;i<row;i++){
  29         for (j =0;j<9;j++){
  30       storage[i][j]=0;
  31         }
  32   }
  33
  34   for(i=0; i<row; i++){
  35     for(j=0; j<col;j++){
  36       matrix[i][j]=data[i][j];
  37       pmatrix[c]=0; // initializing to zero
  38       permuted[c]=0;
  39       c++;
  40     }
  41   }
  42
  43   // Produces the sum of each column
  44   double total[col],total1=0,total2=0;
  45   double ratio[col];
  46
  47   for(i=0;i<col;i++){
  48     total[i]=0;
  49     ratio[i]=0; }
  50
  51   for(i=0; i<col; i++){
  52     for(j=0;j<row;j++){
  53       total[i]=total[i]+matrix[j][i];
  54     }
  55   }
  56
  57   for(i=0;i<g-1;i++){
  58         total1=total1+total[i];}
  59
  60   for(i=g-1;i<col;i++){
  61     total2=total2+total[i];}
  62
  63
  64   // Creates the ratios by first finding the minimum of totals
  65   min = total[0];
  66   if (col==2){
  67   if (total[0]<total[1]){
  68         min = total[1];}
  69   }
  70   if (col >2){
  71   for(i=1;i<col;i++){
  72     if (min > total[i]){
  73       min = total[i];}
  74   }
  75   }
  76   if (min<=0){
  77     printf("Error, the sum of one of the columns <= 0.");
  78     return 0;
  79   }
  80
  81
  82   // Ratio time...
  83   for(i=0;i<col;i++){
  84     ratio[i]=total[i]/min;
  85   }
  86
  87   // Change matrix into an array as received by R for compatibility.
  88
  89   c=0;
  90   for(i=0;i<col;i++){
  91     for(j=0;j<row;j++){
  92       pmatrix[c]=matrix[j][i];
  93       c++;
  94     }
  95   }
  96
  97   if(row == 1){
  98         for (i =0; i<col;i++){
  99                 pmatrix[i]=pmatrix[i]/ratio[i];
 100         }
 101   }
 102   else {
 103     counter = 0;
 104     j=-1;
 105     for (i=0; i<size; i++) {
 106       if (counter % row == 0) {
 107         j++;
 108       }
 109       pmatrix[i]=pmatrix[i]/ratio[j];
 110       counter++;
 111     }
 112   }
 113   // pass everything to the rest of the code using pointers. then
 114   // write to output file. below pointers for most of the values are
 115   // created to send everything by reference.
 116
 117   int ptt_size, *permutes,*nc,*nr,*gvalue;
 118
 119   nc = &col;
 120   nr = &row;
 121   gvalue = &g;
 122
 123   permutes = &B;
 124   ptt_size = B*row;
 125
 126   //changing ptt_size to row
 127   double permuted_ttests[row], pvalues[row], tinitial[row];
 128
 129   for(i=0;i<row;i++){
 130     permuted_ttests[i]=0;}
 131
 132   for(i=0;i<row;i++){
 133     pvalues[i]=0;
 134     tinitial[i]=0; }
 135
 136   // Find the initial values for the matrix.
 137   start(pmatrix,gvalue,nr,nc,tinitial,storage);
 138
 139   // Start the calculations.
 140
 141   if ( (col==2) || ((g-1)<8) || ((col-g+1) < 8) ){
 142
 143   double fish[row], fish2[row];
 144   for(i=0;i<row;i++){
 145         fish[i]=0;
 146         fish2[i]=0;}
 147
 148   for(i=0;i<row;i++){
 149
 150         for(j=0;j<g-1;j++){
 151                 fish[i]=fish[i]+matrix[i][j];
 152         }
 153
 154         for(j=g-1;j<col;j++){
 155                 fish2[i]=fish2[i]+matrix[i][j];
 156         }
 157
 158         double  f11,f12,f21,f22;
 159
 160         f11=fish[i];
 161         f12=fish2[i];
 162
 163         f21=total1-f11;
 164         f22=total2-f12;
 165
 166         double data[] = {f11, f12, f21, f22};
 167
 168         // CONTINGENGCY TABLE:
 169         //   f11   f12
 170         //   f21   f22
 171
 172         int *nr, *nc, *ldtabl, *work;
 173         int nrow=2, ncol=2, ldtable=2, workspace=100000;
 174         double *expect, *prc, *emin,*prt,*pre;
 175         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 176
 177         nr = &nrow;
 178         nc = &ncol;
 179         ldtabl=&ldtable;
 180         work = &workspace;
 181
 182         expect = &e;
 183         prc = &prc1;
 184         emin=&emin1;
 185         prt=&prt1;
 186         pre=&pre1;
 187
 188         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 189
 190         if (*pre>.999999999){
 191                 *pre=1;
 192         }
 193         storage[i][8] = *pre;
 194         pvalues[i]=*pre;
 195     }
 196   }
 197   else{
 198
 199   testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
 200
 201                // Checks to make sure the matrix isn't sparse.
 202   double sparse[row], sparse2[row];
 203   for(i=0;i<row;i++){
 204         sparse[i]=0;
 205         sparse2[i]=0;}
 206
 207   c=0;
 208   for(i=0;i<row;i++){
 209
 210         for(j=0;j<g-1;j++){
 211                 sparse[i]=sparse[i]+matrix[i][j];
 212         }
 213
 214         if(sparse[i] < (double)(g-1)){
 215                 c++;
 216         }
 217         for(j=g-1;j<col;j++){ // ?<= for col
 218                 sparse2[i]=sparse2[i]+matrix[i][j];
 219         }
 220
 221         if( (sparse2[i] <(double)(col-g+1))) {
 222                 c++;
 223         }
 224
 225         if (c==2){
 226                 c=0;
 227
 228         double  f11,f12,f21,f22;
 229
 230         f11=sparse[i];
 231         sparse[i]=0;
 232
 233         f12=sparse2[i];
 234         sparse2[i]=0;
 235
 236         f21=total1-f11;
 237         f22=total2-f12;
 238
 239         double data[] = {f11, f12, f21, f22};
 240
 241         int *nr, *nc, *ldtabl, *work;
 242         int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
 243         double *expect, *prc, *emin,*prt,*pre;
 244         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 245
 246         nr = &nrow;
 247         nc = &ncol;
 248         ldtabl=&ldtable;
 249         work = &workspace;
 250
 251         expect = &e;
 252         prc = &prc1;
 253         emin=&emin1;
 254         prt=&prt1;
 255         pre=&pre1;
 256
 257         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 258
 259         if (*pre>.999999999){
 260                 *pre=1;
 261         }
 262         storage[i][8] = *pre;
 263         pvalues[i]=*pre;
 264     }
 265   }
 266   // End of else statement
 267   bflag = 1;
 268   }
 269
 270   // Calculates the mean of counts (not normalized)
 271   double temp[row][2];
 272
 273   for(j=0;j<row;j++){
 274         for(i=0;i<2;i++){
 275                 temp[j][i]=0;
 276         }
 277   }
 278
 279   for (j=0;j<row;j++){
 280         for (i=1; i<=(g-1); i++){
 281                 temp[j][0]=temp[j][0]+matrix[j][i-1];
 282         }
 283         temp[j][0]= (double) temp[j][0]/(g-1);
 284         for(i=g;i<=col;i++){
 285                 temp[j][1]=temp[j][1]+matrix[j][i-1];
 286         }
 287         temp[j][1]= (double) temp[j][1]/(col-g+1);
 288   }
 289
 290   for(i=0;i<row;i++){
 291         storage[i][3]=temp[i][0];
 292         storage[i][7]=temp[i][1];
 293         storage[i][8]=pvalues[i];
 294   }
 295
 296 // BACKUP checks
 297
 298   for (i=0;i<row;i++){
 299     if(pvalues[i]<thresh){
 300         printf("Feature %d is significant, p = %.10lf \n",i+1,pvalues[i]);
 301         }
 302   }
 303
 304   // And now we write the files to a text file.
 305   struct tm *local;
 306   time_t t;
 307   t = time(NULL);
 308   local = localtime(&t);
 309
 310   out = fopen(output,"w");
 311
 312   fprintf(out,"Local time and date of test: %s\n", asctime(local));
 313   fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
 314   if (bflag == 1){
 315     fprintf(out,"%d permutations\n\n",B);
 316   }
 317
 318   //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file
 319   //storage 0 = meanGroup1 - line 529, 1 = varGroup1 - line 532, 2 = err rate1 - line 534, 3 = mean of counts group1?? - line 291, 4 = meanGroup2 - line 536, 5 = varGroup2 - line 539, 6 = err rate2 - line 541, 7 = mean of counts group2?? - line 292, 8 = pvalues - line 293
 320   fprintf(out, "OTU\tmean(group1)\tvariance(group1)\tstderr(group1)\tmean_of_counts(group1)\tmean(group2)\tvariance(group2)\tstderr(group2)\tmean_of_counts(group1)\tp-value\n");
 321
 322   for(i=0; i<row; i++){
 323         fprintf(out,"%d",(i+1));
 324
 325     for(j=0; j<9;j++){
 326       fprintf(out,"\t%.12lf",storage[i][j]);
 327     }
 328     fprintf(out,"\n");
 329   }
 330
 331   fprintf(out,"\n \n");
 332
 333  // fclose(jobj);
 334   fclose(out);
 335
 336   return 0;
 337 }
 338
 339 void testp(double *permuted_ttests,int *B,double *permuted,
 340            double *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double
 341            *ps) {
 342
 343   double Tvalues[*nr];
 344   int a, b, n, i, j,k=0;
 345
 346   a = *B;
 347   b = *nr;
 348   n = a*b;
 349
 350   double counter[b];
 351
 352   for(j=0;j<b;j++){
 353     counter[j]=0;
 354   }
 355
 356   for (j=1; j<=*B; j++){
 357     permute_matrix(Imatrix,nc,nr,permuted,g,Tvalues,Tinitial,counter);
 358    // for(i=0;i<*nr;i++){
 359    //   permuted_ttests[k]=fabs(Tvalues[i]);
 360         //    k++;
 361     }
 362
 363
 364   for(j=0;j<*nr;j++){
 365     ps[j]=((counter[j]+1)/(double)(a+1));
 366   }
 367 }
 368
 369 void permute_matrix(double *Imatrix,int *nc,int *nr,double *permuted,
 370                     int *g,double *trial_ts,double *Tinitial,double *counter1){
 371
 372   int i=0,j=0,n=0,a=0,b=0,f=0,c=0,k=0;
 373
 374   a = *nr; // number of rows
 375   b = *nc;
 376   n = a*b;
 377
 378   int y[b];
 379
 380   for (i=1; i<=*nc; i++){
 381     y[i-1] = i;
 382   }
 383
 384   permute_array(y, b);
 385
 386   for (i=0; i<*nc; i++){
 387     f = y[i]; //column number
 388     c=1;
 389     c*=(f-1);
 390     c*=a;
 391     if (f == 1){
 392       c = 0;
 393     } // starting value position in the Imatrix
 394     for(j=1; j<=*nr; j++){
 395       permuted[k] = Imatrix[c];
 396       c++;
 397       k++;
 398     }
 399   }
 400
 401   calc_twosample_ts(permuted,g,nc,nr,trial_ts,Tinitial,counter1);
 402 }
 403
 404 void permute_array(int *array, int n) {
 405   static int seeded = 0;
 406   int i;
 407
 408   if (! seeded) {
 409     seeded = 1;
 410     srandom(time(NULL));
 411   }
 412
 413   for (i = 0; i < n; i++) {
 414     int selection = random() % (n - i);
 415     int tmp = array[i + selection];
 416     array[i + selection] = array[i];
 417     array[i] = tmp;
 418   }
 419 }
 420
 421 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,
 422                        double *Ts,double *Tinitial,double *counter) {
 423   int i,a;
 424   a = *nr;
 425   a*=4;
 426
 427   double C1[*nr][3], C2[*nr][3], storage[a],tool[a];
 428   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 429
 430   nrows = (double) *nr;
 431   ncols = (double) *nc;
 432   gvalue= (double) *g;
 433
 434     meanvar(Pmatrix,g,nr,nc,storage);
 435     for(i=0;i<=a-1;i++){
 436       tool[i]=storage[i];
 437     }
 438     for (i=0; i<*nr;i++){
 439       C1[i][0]=tool[i];
 440       C1[i][1]=tool[i+*nr+*nr];
 441       C1[i][2]=C1[i][1]/(gvalue-1);
 442
 443       C2[i][0]=tool[i+*nr];
 444       C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
 445       C2[i][2]=C2[i][1]/(ncols-gvalue+1);
 446     }
 447
 448     for (i=0; i<*nr; i++){
 449       xbardiff = C1[i][0]-C2[i][0];
 450       denom = sqrt(C1[i][2]+C2[i][2]);
 451       Ts[i]=fabs(xbardiff/denom);
 452       if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
 453             counter[i]++;
 454       }
 455     }
 456 }
 457
 458 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
 459   double temp[*nr], temp2[*nr],var[*nr],var2[*nr],a,b;
 460
 461   int i,m,k,l,n;
 462
 463   a = (double) *g-1;
 464   b = (double) (*nc-a);
 465
 466   for (i = 0; i<*nr; i++){
 467     temp[i]=0;
 468     temp2[i]=0;
 469     var[i]=0;
 470     var2[i]=0;
 471   }
 472
 473   k = *nr; // number of rows
 474   l = *nc;
 475   n = k*l;
 476
 477     m=0;
 478     m=*g-1;
 479     k=*nr;
 480     m*=k; // m = g * nr now
 481     for (i=0;i<m;i++){
 482       temp[i%k]=temp[i%k]+pmatrix[i];
 483     }
 484     for (i=0;i<n;i++){
 485       temp2[i%k]=temp2[i%k]+pmatrix[i];
 486     }
 487     for (i=0;i<*nr;i++){
 488       temp2[i]=temp2[i]-temp[i];
 489     }
 490     for (i=0;i<=*nr-1;i++){
 491       store[i]=temp[i]/a;
 492       store[i+*nr]=temp2[i]/b;
 493     }
 494
 495     // That completes the mean calculations.
 496
 497     for (i=0;i<m;i++){
 498       var[i%k]=var[i%k]+pow((pmatrix[i]-store[i%k]),2);
 499     }
 500     for (i=m;i<n;i++){
 501       var2[i%k]=var2[i%k]+pow((pmatrix[i]-store[(i%k)+*nr]),2);
 502     }
 503
 504     for (i=0;i<=*nr-1;i++){
 505       store[i+2*k]=var[i]/(a-1);
 506       store[i+3*k]=var2[i]/(b-1);
 507     }
 508     // That completes var calculations.
 509 }
 510
 511 void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
 512                 double storage[][9]){
 513   int i, a = *nr;
 514   a*=4;
 515
 516   double store[a], tool[a], C1[*nr][3], C2[*nr][3];
 517   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 518
 519   nrows = (double) *nr;
 520   ncols = (double) *nc;
 521   gvalue= (double) *g;
 522
 523   meanvar(Imatrix,g,nr,nc,store);
 524
 525   for(i=0;i<=a-1;i++){
 526     tool[i]=store[i];
 527   }
 528   for (i=0; i<*nr;i++){
 529     C1[i][0]=tool[i]; //mean group 1
 530     storage[i][0]=C1[i][0];
 531     C1[i][1]=tool[i+*nr+*nr]; // var group 1
 532     storage[i][1]=C1[i][1];
 533     C1[i][2]=C1[i][1]/(gvalue-1);
 534     storage[i][2]=sqrt(C1[i][2]);
 535
 536     C2[i][0]=tool[i+*nr]; // mean group 2
 537     storage[i][4]=C2[i][0];
 538     C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
 539     storage[i][5]=C2[i][1];
 540     C2[i][2]=C2[i][1]/(ncols-gvalue+1);
 541     storage[i][6]=sqrt(C2[i][2]);
 542   }
 543   for (i=0; i<*nr; i++){
 544     xbardiff = C1[i][0]-C2[i][0];
 545     denom = sqrt(C1[i][2]+C2[i][2]);
 546     initial[i]=fabs(xbardiff/denom);
 547   }
 548 }
 549
 550
 551
 552