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