X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=metastats2.c;h=a6424f5fa95086083d4c09658f7c9e2338d2381c;hb=8173238f9f94af9baab8471de58bed7c8830948d;hp=c4d1f46feaff17f90a6e6668126a0f0a1ecd96ac;hpb=2bf3df7736ef2a17286d99394e211f51751d6829;p=mothur.git

diff --git a/metastats2.c b/metastats2.c
index c4d1f46..a6424f5 100644
--- a/metastats2.c
+++ b/metastats2.c
@@ -1,167 +1,52 @@
-#include <stdio.h>
-#include <string.h>
-#include <stdlib.h>
-#include <time.h>
-#include <math.h>
-#include "fisher2.h"
 
-void testp(double *permuted_ttests,int *B,double *permuted,double 
-	   *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double *ps);
-void permute_matrix(double *Imatrix,int *nc,int *nr,double 
-		    *permuted,int *g,double *trial_ts,double *Tinitial,double 
-		    *counter);
-void permute_array(int *array, int n);
-void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,double 
-		       *Ts,double *Tinitial,double *counter1);
-void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *storage);
-void start(double *Imatrix,int *g,int *nr,int *nc,double *testing,
-			double storage[][9]);
+#include "metastats.h"
 
-int main (int argc, char *argv[]){
-  
-  int col=-1, row=-1, size,g=0,c=0,i=0,j=0,k,counter=0,lines=0, B=10000;
-  double placeholder=0,min=0, thresh=0.05;
-
-  char arr[10000], str[51], a;
-  char location[41]="jobj.txt", output[41]="out.txt";
-  
-  for(i=0;i<10000;i++){
-  	arr[i]='q';
-  }
-  
-int u,rflag=0,cflag=0, bflag=0;
-char *filename;
-int numbers;
-double numb;
-extern char *optarg;
-extern int optind, optopt, opterr;
+//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).
 
-while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
-    switch(u) {
-    case 'r':
-        numbers = atoi(optarg);
-        printf("The number of features/rows is %d.\n", numbers);
-        row = numbers;
-        rflag = 1;
-        break;
-    case 'c':
-        numbers = atoi(optarg);
-        printf("The number of samples/columns is %d.\n", numbers);
-        col = numbers;
-        cflag = 1;
-        break;
-    case 'g':
-        numbers = atoi(optarg);
-        printf("Your g-value is %d.\n", numbers);
-        g = numbers;
-        break;
-    case 'b':
-        numbers = atoi(optarg);
-        printf("The number of permutations is %d\n", numbers);
-        B = numbers;
-        break;
-    case 't':
-        numb = atof(optarg);
-        printf("Threshold is is %lf\n", numb);
-        thresh = numb;
-        break;        
-    case 'f':
-        filename = optarg;
-        printf("filename input is %s\n", filename);
-        strcpy(location,filename);
-        break;
-    case 'o':
-        filename = optarg;
-        printf("filename output %s\n", filename);
-        strcpy(output,filename);
-        break;
-    case ':':
-        printf("-%c without filename\n", optopt);
-        break;
-    case '?':
-        printf("unknown arg %c\n", optopt);
-        break;
-    }
-}
-  
-  FILE *jobj, *out;
-  jobj=fopen(location,"r");
+int metastat_main (char* outputFileName, int numRows, int numCols, double threshold, int numPermutations, double** data, int secondGroupingStart){
   
-  if(jobj == NULL){
-    printf("Please don't forget to save your matrix in the active");
-    printf(" directory as \"%s\".\n",location);
-    return 0;
-  }
- 
-  // Gets the first line of samples names and checks for user error.
-  fgets(arr,10000,jobj);
-  
-  for(i=0;i<10000;i++){
-    if(isspace(arr[i])){
-      counter++; }
-  }
-
-  if (cflag == 0) {
-      printf("You didn't tell us how many subjects there are!\n");
-      printf("But we'll still do the analysis as if there are %d subjects.\n\n",col=counter-1);
-  }
-  if (cflag == 1) {
-  	if (col != counter-1){
-  	  printf("We would expect %d subjects, but you said %d.\n",counter-1,col);
-    } 
-  }
+  int size,c=0,i=0,j=0,k,counter=0, bflag=0; 
+  int B=numPermutations;
+  int row = numRows;
+  int col = numCols;
+  int g = secondGroupingStart;
+  double thresh=threshold;
+  double placeholder=0,min=0; 
   
-  while((a = fgetc(jobj)) != EOF){
-    if(a == '\n'){
-      lines++; }
-  }
+  char output[1024];
+  strcpy(output, outputFileName);
+  FILE *out;
   
-  if (rflag == 0) {
-      printf("You didn't specify the number of features!\n");
-      printf("But we'll still do the analysis assuming %d features.\n\n", row=lines-1);   
-  }
-  if (rflag == 1) {
-    if ( lines-1 != row ){
-     printf("We would expect %d features, but you said %d.\n",lines-1,row);
-    }
-  }
- 
   if (g>=col || g<=0){
   	printf("Check your g value\n");	
   }
  
   // Initialize the matrices
-  size = row*col;
-  double matrix[row][col];
-  double pmatrix[size],pmatrix2[size],permuted[size];  
-  double storage[row][9];
-  
-  for (i=0;i<row;i++){
-  	for (j =0;j<9;j++){
-      storage[i][j]=0; 		
-	}	
-  }
-  // Reset file below and create a separate matrix.
-  rewind(jobj);
-  fgets(arr,10000,jobj);
-  
+	size = row*col;
+	
+	double ** matrix, * pmatrix, * permuted, ** storage;
+	matrix = malloc(row*sizeof(double *));
+	storage =  malloc(row*sizeof(double *));
+	for (i = 0; i<row; i++){
+		matrix[i] = malloc(col*sizeof(double));
+	}
+	for (i = 0; i<row;i++){
+		storage[i] = malloc(9*sizeof(double));
+	}
+	
+	pmatrix = (double *) malloc(size*sizeof(double *));
+	permuted = (double *) malloc(size*sizeof(double *));
+	
+		
+	
   for(i=0; i<row; i++){
-    fscanf(jobj,"%s",str); 	
     for(j=0; j<col;j++){
-      fscanf(jobj,"%lf",&placeholder);
-      matrix[i][j]=placeholder;
-      if(isalnum(placeholder)!=0){ // check for ""
-      	printf("Your matrix isn't set up properly!\n");
-      	return 0;
-		}
+      matrix[i][j]=data[i][j];
       pmatrix[c]=0; // initializing to zero
       permuted[c]=0;
       c++;
     }
   }
-  
-  fclose(jobj);
-
 
   // Produces the sum of each column
   double total[col],total1=0,total2=0;
@@ -245,9 +130,12 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
   
   permutes = &B;
   ptt_size = B*row;
-  
+
   //changing ptt_size to row
-  double permuted_ttests[row], pvalues[row], tinitial[row];
+  double * permuted_ttests, * pvalues, * tinitial;
+  permuted_ttests = (double *) malloc(size*sizeof(double *));
+  pvalues = (double *) malloc(size*sizeof(double *));
+  tinitial = (double *) malloc(size*sizeof(double *));
   
   for(i=0;i<row;i++){
     permuted_ttests[i]=0;}
@@ -257,13 +145,17 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
     tinitial[i]=0; }
   
   // Find the initial values for the matrix.
+ 
   start(pmatrix,gvalue,nr,nc,tinitial,storage);
   
   // Start the calculations.
  
   if ( (col==2) || ((g-1)<8) || ((col-g+1) < 8) ){  
 
-  double fish[row], fish2[row];
+  double * fish, *fish2;
+  fish = (double *) malloc(size*sizeof(double *));
+  fish2 = (double *) malloc(size*sizeof(double *));
+	  
   for(i=0;i<row;i++){
   	fish[i]=0;
   	fish2[i]=0;}
@@ -293,9 +185,13 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
 	//   f21   f22
 	
 	int *nr, *nc, *ldtabl, *work;
-	int nrow=2, ncol=2, ldtable=2, workspace=100000;
+	int nrow=2, ncol=2, ldtable=2;
+	int workspace = 6*(row*col*sizeof(double *)); 
 	double *expect, *prc, *emin,*prt,*pre;
 	double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
+	  
+	prt = (double *) malloc(size*sizeof(double *));
+	prc = (double *) malloc(size*sizeof(double *));
 
 	nr = &nrow;
 	nc = &ncol;
@@ -307,7 +203,7 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
 	emin=&emin1;
 	prt=&prt1;
 	pre=&pre1;
-	
+
 	fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 	
 	if (*pre>.999999999){
@@ -318,11 +214,14 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
     }
   }
   else{
-  	 
+printf("here before testp\n");	  	 
   testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
   	 
   	       // Checks to make sure the matrix isn't sparse.
-  double sparse[row], sparse2[row];
+  double * sparse, * sparse2;
+  sparse = (double *) malloc(size*sizeof(double *));
+  sparse2 = (double *) malloc(size*sizeof(double *));
+	  
   for(i=0;i<row;i++){
   	sparse[i]=0;
   	sparse2[i]=0;}
@@ -362,7 +261,7 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
 	double data[] = {f11, f12, f21, f22};
 
 	int *nr, *nc, *ldtabl, *work;
-	int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
+	int nrow=2, ncol=2, ldtable=2, workspace=INT_MAX; // I added two zeros for larger data sets
 	double *expect, *prc, *emin,*prt,*pre;
 	double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
 
@@ -376,7 +275,7 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
 	emin=&emin1;
 	prt=&prt1;
 	pre=&pre1;
-	
+printf("here before fexact2\n");		
 	fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
 	
 	if (*pre>.999999999){
@@ -391,8 +290,12 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
   }
   
   // Calculates the mean of counts (not normalized)
-  double temp[row][2];
-  
+  double ** temp;
+  temp = malloc(row*sizeof(double *));
+  for (i = 0; i<row; i++){
+	 temp[i] = malloc(col*sizeof(double));
+  }
+printf("here\n");	  
   for(j=0;j<row;j++){
   	for(i=0;i<2;i++){
   		temp[j][i]=0;
@@ -430,22 +333,21 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
   t = time(NULL);
   local = localtime(&t);
   
-  jobj= fopen(location,"r");
-  fgets(arr,10000,jobj);
-  
-  out = fopen(output,"a+");
+  out = fopen(output,"w");
   
   fprintf(out,"Local time and date of test: %s\n", asctime(local));
   fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
   if (bflag == 1){
     fprintf(out,"%d permutations\n\n",B);	
   }
+  
+  //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file
+  //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
+  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");
+    
   for(i=0; i<row; i++){
-    fscanf(jobj,"%s",str); 	
-	fprintf(out,"%s",str);
-	for(k=0;k<col;k++){
-		fscanf(jobj,"%*lf",&placeholder);
-	}
+	fprintf(out,"%d",(i+1));
+	
     for(j=0; j<9;j++){
       fprintf(out,"\t%.12lf",storage[i][j]);
     }
@@ -454,7 +356,7 @@ while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
   
   fprintf(out,"\n \n");
   
-  fclose(jobj);
+ // fclose(jobj);
   fclose(out);
   
   return 0;
@@ -531,11 +433,11 @@ void permute_array(int *array, int n) {
   
   if (! seeded) {
     seeded = 1;
-    srandom(time(NULL));
+    srand(time(NULL));
   }
   
   for (i = 0; i < n; i++) {
-    int selection = random() % (n - i);
+    int selection = rand() % (n - i);
     int tmp = array[i + selection];
     array[i + selection] = array[i];
     array[i] = tmp;
@@ -633,40 +535,47 @@ void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
 }
 
 void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
-		double storage[][9]){
+		double** storage){
   int i, a = *nr;
   a*=4;
-  
+	
   double store[a], tool[a], C1[*nr][3], C2[*nr][3];
   double nrows,ncols,gvalue, xbardiff=0, denom=0;
-  
+ 
   nrows = (double) *nr;
   ncols = (double) *nc;
   gvalue= (double) *g;
-  
+   
   meanvar(Imatrix,g,nr,nc,store);
-  
+   
   for(i=0;i<=a-1;i++){
     tool[i]=store[i];
   }
-  for (i=0; i<*nr;i++){
+
+  for (i=0; i<*nr;i++){ 
     C1[i][0]=tool[i]; //mean group 1
     storage[i][0]=C1[i][0];
     C1[i][1]=tool[i+*nr+*nr]; // var group 1
     storage[i][1]=C1[i][1];
     C1[i][2]=C1[i][1]/(gvalue-1);
     storage[i][2]=sqrt(C1[i][2]);
-    
+    printf("here2\n"); 
     C2[i][0]=tool[i+*nr]; // mean group 2
-    storage[i][4]=C2[i][0];    
+    storage[i][4]=C2[i][0]; 
     C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2 
-    storage[i][5]=C2[i][1];        
-    C2[i][2]=C2[i][1]/(ncols-gvalue+1);
+    storage[i][5]=C2[i][1];  
+	C2[i][2]=C2[i][1]/(ncols-gvalue+1);
     storage[i][6]=sqrt(C2[i][2]);   
+	 
   }
+ 
   for (i=0; i<*nr; i++){
     xbardiff = C1[i][0]-C2[i][0];
     denom = sqrt(C1[i][2]+C2[i][2]);
     initial[i]=fabs(xbardiff/denom);
-  }											
+  }	 
 }
+
+
+
+