-#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;
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;}
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;}
// 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;
emin=&emin1;
prt=&prt1;
pre=&pre1;
-
+
fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
if (*pre>.999999999){
}
}
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;}
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;
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){
}
// 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;
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]);
}
fprintf(out,"\n \n");
- fclose(jobj);
+ // fclose(jobj);
fclose(out);
return 0;
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;
}
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);
- }
+ }
}
+
+
+
+
projects / mothur.git / blobdiff