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