metastats2.c

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