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=(row*sizeof(double *)+size*sizeof(double *));
 190         double *expect, *prc, *emin,*prt,*pre;
 191         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 192
 193         nr = &nrow;
 194         nc = &ncol;
 195         ldtabl=&ldtable;
 196         work = &workspace;
 197
 198         expect = &e;
 199         prc = &prc1;
 200         emin=&emin1;
 201         prt=&prt1;
 202         pre=&pre1;
 203
 204         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 205
 206         if (*pre>.999999999){
 207                 *pre=1;
 208         }
 209         storage[i][8] = *pre;
 210         pvalues[i]=*pre;
 211     }
 212   }
 213   else{
 214 printf("here before testp\n");
 215   testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
 216
 217                // Checks to make sure the matrix isn't sparse.
 218   double * sparse, * sparse2;
 219   sparse = (double *) malloc(size*sizeof(double *));
 220   sparse2 = (double *) malloc(size*sizeof(double *));
 221
 222   for(i=0;i<row;i++){
 223         sparse[i]=0;
 224         sparse2[i]=0;}
 225
 226   c=0;
 227   for(i=0;i<row;i++){
 228
 229         for(j=0;j<g-1;j++){
 230                 sparse[i]=sparse[i]+matrix[i][j];
 231         }
 232
 233         if(sparse[i] < (double)(g-1)){
 234                 c++;
 235         }
 236         for(j=g-1;j<col;j++){ // ?<= for col
 237                 sparse2[i]=sparse2[i]+matrix[i][j];
 238         }
 239
 240         if( (sparse2[i] <(double)(col-g+1))) {
 241                 c++;
 242         }
 243
 244         if (c==2){
 245                 c=0;
 246
 247         double  f11,f12,f21,f22;
 248
 249         f11=sparse[i];
 250         sparse[i]=0;
 251
 252         f12=sparse2[i];
 253         sparse2[i]=0;
 254
 255         f21=total1-f11;
 256         f22=total2-f12;
 257
 258         double data[] = {f11, f12, f21, f22};
 259
 260         int *nr, *nc, *ldtabl, *work;
 261         int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
 262         double *expect, *prc, *emin,*prt,*pre;
 263         double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 264
 265         nr = &nrow;
 266         nc = &ncol;
 267         ldtabl=&ldtable;
 268         work = &workspace;
 269
 270         expect = &e;
 271         prc = &prc1;
 272         emin=&emin1;
 273         prt=&prt1;
 274         pre=&pre1;
 275 printf("here before fexact2\n");
 276         fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 277
 278         if (*pre>.999999999){
 279                 *pre=1;
 280         }
 281         storage[i][8] = *pre;
 282         pvalues[i]=*pre;
 283     }
 284   }
 285   // End of else statement
 286   bflag = 1;
 287   }
 288
 289   // Calculates the mean of counts (not normalized)
 290   double ** temp;
 291   temp = malloc(row*sizeof(double *));
 292   for (i = 0; i<row; i++){
 293          temp[i] = malloc(col*sizeof(double));
 294   }
 295 printf("here\n");
 296   for(j=0;j<row;j++){
 297         for(i=0;i<2;i++){
 298                 temp[j][i]=0;
 299         }
 300   }
 301
 302   for (j=0;j<row;j++){
 303         for (i=1; i<=(g-1); i++){
 304                 temp[j][0]=temp[j][0]+matrix[j][i-1];
 305         }
 306         temp[j][0]= (double) temp[j][0]/(g-1);
 307         for(i=g;i<=col;i++){
 308                 temp[j][1]=temp[j][1]+matrix[j][i-1];
 309         }
 310         temp[j][1]= (double) temp[j][1]/(col-g+1);
 311   }
 312
 313   for(i=0;i<row;i++){
 314         storage[i][3]=temp[i][0];
 315         storage[i][7]=temp[i][1];
 316         storage[i][8]=pvalues[i];
 317   }
 318
 319 // BACKUP checks
 320
 321   for (i=0;i<row;i++){
 322     if(pvalues[i]<thresh){
 323         printf("Feature %d is significant, p = %.10lf \n",i+1,pvalues[i]);
 324         }
 325   }
 326
 327   // And now we write the files to a text file.
 328   struct tm *local;
 329   time_t t;
 330   t = time(NULL);
 331   local = localtime(&t);
 332
 333   out = fopen(output,"w");
 334
 335   fprintf(out,"Local time and date of test: %s\n", asctime(local));
 336   fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
 337   if (bflag == 1){
 338     fprintf(out,"%d permutations\n\n",B);
 339   }
 340
 341   //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file
 342   //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
 343   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");
 344
 345   for(i=0; i<row; i++){
 346         fprintf(out,"%d",(i+1));
 347
 348     for(j=0; j<9;j++){
 349       fprintf(out,"\t%.12lf",storage[i][j]);
 350     }
 351     fprintf(out,"\n");
 352   }
 353
 354   fprintf(out,"\n \n");
 355
 356  // fclose(jobj);
 357   fclose(out);
 358
 359   return 0;
 360 }
 361
 362 void testp(double *permuted_ttests,int *B,double *permuted,
 363            double *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double
 364            *ps) {
 365
 366   double Tvalues[*nr];
 367   int a, b, n, i, j,k=0;
 368
 369   a = *B;
 370   b = *nr;
 371   n = a*b;
 372
 373   double counter[b];
 374
 375   for(j=0;j<b;j++){
 376     counter[j]=0;
 377   }
 378
 379   for (j=1; j<=*B; j++){
 380     permute_matrix(Imatrix,nc,nr,permuted,g,Tvalues,Tinitial,counter);
 381    // for(i=0;i<*nr;i++){
 382    //   permuted_ttests[k]=fabs(Tvalues[i]);
 383         //    k++;
 384     }
 385
 386
 387   for(j=0;j<*nr;j++){
 388     ps[j]=((counter[j]+1)/(double)(a+1));
 389   }
 390 }
 391
 392 void permute_matrix(double *Imatrix,int *nc,int *nr,double *permuted,
 393                     int *g,double *trial_ts,double *Tinitial,double *counter1){
 394
 395   int i=0,j=0,n=0,a=0,b=0,f=0,c=0,k=0;
 396
 397   a = *nr; // number of rows
 398   b = *nc;
 399   n = a*b;
 400
 401   int y[b];
 402
 403   for (i=1; i<=*nc; i++){
 404     y[i-1] = i;
 405   }
 406
 407   permute_array(y, b);
 408
 409   for (i=0; i<*nc; i++){
 410     f = y[i]; //column number
 411     c=1;
 412     c*=(f-1);
 413     c*=a;
 414     if (f == 1){
 415       c = 0;
 416     } // starting value position in the Imatrix
 417     for(j=1; j<=*nr; j++){
 418       permuted[k] = Imatrix[c];
 419       c++;
 420       k++;
 421     }
 422   }
 423
 424   calc_twosample_ts(permuted,g,nc,nr,trial_ts,Tinitial,counter1);
 425 }
 426
 427 void permute_array(int *array, int n) {
 428   static int seeded = 0;
 429   int i;
 430
 431   if (! seeded) {
 432     seeded = 1;
 433     srand(time(NULL));
 434   }
 435
 436   for (i = 0; i < n; i++) {
 437     int selection = rand() % (n - i);
 438     int tmp = array[i + selection];
 439     array[i + selection] = array[i];
 440     array[i] = tmp;
 441   }
 442 }
 443
 444 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,
 445                        double *Ts,double *Tinitial,double *counter) {
 446   int i,a;
 447   a = *nr;
 448   a*=4;
 449
 450   double C1[*nr][3], C2[*nr][3], storage[a],tool[a];
 451   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 452
 453   nrows = (double) *nr;
 454   ncols = (double) *nc;
 455   gvalue= (double) *g;
 456
 457     meanvar(Pmatrix,g,nr,nc,storage);
 458     for(i=0;i<=a-1;i++){
 459       tool[i]=storage[i];
 460     }
 461     for (i=0; i<*nr;i++){
 462       C1[i][0]=tool[i];
 463       C1[i][1]=tool[i+*nr+*nr];
 464       C1[i][2]=C1[i][1]/(gvalue-1);
 465
 466       C2[i][0]=tool[i+*nr];
 467       C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
 468       C2[i][2]=C2[i][1]/(ncols-gvalue+1);
 469     }
 470
 471     for (i=0; i<*nr; i++){
 472       xbardiff = C1[i][0]-C2[i][0];
 473       denom = sqrt(C1[i][2]+C2[i][2]);
 474       Ts[i]=fabs(xbardiff/denom);
 475       if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
 476             counter[i]++;
 477       }
 478     }
 479 }
 480
 481 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
 482   double temp[*nr], temp2[*nr],var[*nr],var2[*nr],a,b;
 483
 484   int i,m,k,l,n;
 485
 486   a = (double) *g-1;
 487   b = (double) (*nc-a);
 488
 489   for (i = 0; i<*nr; i++){
 490     temp[i]=0;
 491     temp2[i]=0;
 492     var[i]=0;
 493     var2[i]=0;
 494   }
 495
 496   k = *nr; // number of rows
 497   l = *nc;
 498   n = k*l;
 499
 500     m=0;
 501     m=*g-1;
 502     k=*nr;
 503     m*=k; // m = g * nr now
 504     for (i=0;i<m;i++){
 505       temp[i%k]=temp[i%k]+pmatrix[i];
 506     }
 507     for (i=0;i<n;i++){
 508       temp2[i%k]=temp2[i%k]+pmatrix[i];
 509     }
 510     for (i=0;i<*nr;i++){
 511       temp2[i]=temp2[i]-temp[i];
 512     }
 513     for (i=0;i<=*nr-1;i++){
 514       store[i]=temp[i]/a;
 515       store[i+*nr]=temp2[i]/b;
 516     }
 517
 518     // That completes the mean calculations.
 519
 520     for (i=0;i<m;i++){
 521       var[i%k]=var[i%k]+pow((pmatrix[i]-store[i%k]),2);
 522     }
 523     for (i=m;i<n;i++){
 524       var2[i%k]=var2[i%k]+pow((pmatrix[i]-store[(i%k)+*nr]),2);
 525     }
 526
 527     for (i=0;i<=*nr-1;i++){
 528       store[i+2*k]=var[i]/(a-1);
 529       store[i+3*k]=var2[i]/(b-1);
 530     }
 531     // That completes var calculations.
 532 }
 533
 534 void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
 535                 double** storage){
 536   int i, a = *nr;
 537   a*=4;
 538
 539   double store[a], tool[a], C1[*nr][3], C2[*nr][3];
 540   double nrows,ncols,gvalue, xbardiff=0, denom=0;
 541
 542   nrows = (double) *nr;
 543   ncols = (double) *nc;
 544   gvalue= (double) *g;
 545
 546   meanvar(Imatrix,g,nr,nc,store);
 547
 548   for(i=0;i<=a-1;i++){
 549     tool[i]=store[i];
 550   }
 551
 552   for (i=0; i<*nr;i++){
 553     C1[i][0]=tool[i]; //mean group 1
 554     storage[i][0]=C1[i][0];
 555     C1[i][1]=tool[i+*nr+*nr]; // var group 1
 556     storage[i][1]=C1[i][1];
 557     C1[i][2]=C1[i][1]/(gvalue-1);
 558     storage[i][2]=sqrt(C1[i][2]);
 559     printf("here2\n");
 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