8 void testp(double *permuted_ttests,int *B,double *permuted,double
9 *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double *ps);
10 void permute_matrix(double *Imatrix,int *nc,int *nr,double
11 *permuted,int *g,double *trial_ts,double *Tinitial,double
13 void permute_array(int *array, int n);
14 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,double
15 *Ts,double *Tinitial,double *counter1);
16 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *storage);
17 void start(double *Imatrix,int *g,int *nr,int *nc,double *testing,
20 int main (int argc, char *argv[]){
22 int col=-1, row=-1, size,g=0,c=0,i=0,j=0,k,counter=0,lines=0, B=10000;
23 double placeholder=0,min=0, thresh=0.05;
25 char arr[10000], str[51], a;
26 char location[41]="jobj.txt", output[41]="out.txt";
32 int u,rflag=0,cflag=0, bflag=0;
37 extern int optind, optopt, opterr;
39 while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
42 numbers = atoi(optarg);
43 printf("The number of features/rows is %d.\n", numbers);
48 numbers = atoi(optarg);
49 printf("The number of samples/columns is %d.\n", numbers);
54 numbers = atoi(optarg);
55 printf("Your g-value is %d.\n", numbers);
59 numbers = atoi(optarg);
60 printf("The number of permutations is %d\n", numbers);
65 printf("Threshold is is %lf\n", numb);
70 printf("filename input is %s\n", filename);
71 strcpy(location,filename);
75 printf("filename output %s\n", filename);
76 strcpy(output,filename);
79 printf("-%c without filename\n", optopt);
82 printf("unknown arg %c\n", optopt);
88 jobj=fopen(location,"r");
91 printf("Please don't forget to save your matrix in the active");
92 printf(" directory as \"%s\".\n",location);
96 // Gets the first line of samples names and checks for user error.
97 fgets(arr,10000,jobj);
105 printf("You didn't tell us how many subjects there are!\n");
106 printf("But we'll still do the analysis as if there are %d subjects.\n\n",col=counter-1);
109 if (col != counter-1){
110 printf("We would expect %d subjects, but you said %d.\n",counter-1,col);
114 while((a = fgetc(jobj)) != EOF){
120 printf("You didn't specify the number of features!\n");
121 printf("But we'll still do the analysis assuming %d features.\n\n", row=lines-1);
124 if ( lines-1 != row ){
125 printf("We would expect %d features, but you said %d.\n",lines-1,row);
130 printf("Check your g value\n");
133 // Initialize the matrices
135 double matrix[row][col];
136 double pmatrix[size],pmatrix2[size],permuted[size];
137 double storage[row][9];
144 // Reset file below and create a separate matrix.
146 fgets(arr,10000,jobj);
148 for(i=0; i<row; i++){
149 fscanf(jobj,"%s",str);
151 fscanf(jobj,"%lf",&placeholder);
152 matrix[i][j]=placeholder;
153 if(isalnum(placeholder)!=0){ // check for ""
154 printf("Your matrix isn't set up properly!\n");
157 pmatrix[c]=0; // initializing to zero
166 // Produces the sum of each column
167 double total[col],total1=0,total2=0;
174 for(i=0; i<col; i++){
176 total[i]=total[i]+matrix[j][i];
181 total1=total1+total[i];}
183 for(i=g-1;i<col;i++){
184 total2=total2+total[i];}
187 // Creates the ratios by first finding the minimum of totals
190 if (total[0]<total[1]){
200 printf("Error, the sum of one of the columns <= 0.");
207 ratio[i]=total[i]/min;
210 // Change matrix into an array as received by R for compatibility.
215 pmatrix[c]=matrix[j][i];
221 for (i =0; i<col;i++){
222 pmatrix[i]=pmatrix[i]/ratio[i];
228 for (i=0; i<size; i++) {
229 if (counter % row == 0) {
232 pmatrix[i]=pmatrix[i]/ratio[j];
236 // pass everything to the rest of the code using pointers. then
237 // write to output file. below pointers for most of the values are
238 // created to send everything by reference.
240 int ptt_size, *permutes,*nc,*nr,*gvalue;
249 //changing ptt_size to row
250 double permuted_ttests[row], pvalues[row], tinitial[row];
253 permuted_ttests[i]=0;}
259 // Find the initial values for the matrix.
260 start(pmatrix,gvalue,nr,nc,tinitial,storage);
262 // Start the calculations.
264 if ( (col==2) || ((g-1)<8) || ((col-g+1) < 8) ){
266 double fish[row], fish2[row];
274 fish[i]=fish[i]+matrix[i][j];
277 for(j=g-1;j<col;j++){
278 fish2[i]=fish2[i]+matrix[i][j];
281 double f11,f12,f21,f22;
289 double data[] = {f11, f12, f21, f22};
291 // CONTINGENGCY TABLE:
295 int *nr, *nc, *ldtabl, *work;
296 int nrow=2, ncol=2, ldtable=2, workspace=100000;
297 double *expect, *prc, *emin,*prt,*pre;
298 double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
311 fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
313 if (*pre>.999999999){
316 storage[i][8] = *pre;
322 testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
324 // Checks to make sure the matrix isn't sparse.
325 double sparse[row], sparse2[row];
334 sparse[i]=sparse[i]+matrix[i][j];
337 if(sparse[i] < (double)(g-1)){
340 for(j=g-1;j<col;j++){ // ?<= for col
341 sparse2[i]=sparse2[i]+matrix[i][j];
344 if( (sparse2[i] <(double)(col-g+1))) {
351 double f11,f12,f21,f22;
362 double data[] = {f11, f12, f21, f22};
364 int *nr, *nc, *ldtabl, *work;
365 int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
366 double *expect, *prc, *emin,*prt,*pre;
367 double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
380 fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
382 if (*pre>.999999999){
385 storage[i][8] = *pre;
389 // End of else statement
393 // Calculates the mean of counts (not normalized)
403 for (i=1; i<=(g-1); i++){
404 temp[j][0]=temp[j][0]+matrix[j][i-1];
406 temp[j][0]= (double) temp[j][0]/(g-1);
408 temp[j][1]=temp[j][1]+matrix[j][i-1];
410 temp[j][1]= (double) temp[j][1]/(col-g+1);
414 storage[i][3]=temp[i][0];
415 storage[i][7]=temp[i][1];
416 storage[i][8]=pvalues[i];
422 if(pvalues[i]<thresh){
423 printf("Feature %d is significant, p = %.10lf \n",i+1,pvalues[i]);
427 // And now we write the files to a text file.
431 local = localtime(&t);
433 jobj= fopen(location,"r");
434 fgets(arr,10000,jobj);
436 out = fopen(output,"a+");
438 fprintf(out,"Local time and date of test: %s\n", asctime(local));
439 fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
441 fprintf(out,"%d permutations\n\n",B);
443 for(i=0; i<row; i++){
444 fscanf(jobj,"%s",str);
445 fprintf(out,"%s",str);
447 fscanf(jobj,"%*lf",&placeholder);
450 fprintf(out,"\t%.12lf",storage[i][j]);
455 fprintf(out,"\n \n");
463 void testp(double *permuted_ttests,int *B,double *permuted,
464 double *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double
468 int a, b, n, i, j,k=0;
480 for (j=1; j<=*B; j++){
481 permute_matrix(Imatrix,nc,nr,permuted,g,Tvalues,Tinitial,counter);
482 // for(i=0;i<*nr;i++){
483 // permuted_ttests[k]=fabs(Tvalues[i]);
489 ps[j]=((counter[j]+1)/(double)(a+1));
493 void permute_matrix(double *Imatrix,int *nc,int *nr,double *permuted,
494 int *g,double *trial_ts,double *Tinitial,double *counter1){
496 int i=0,j=0,n=0,a=0,b=0,f=0,c=0,k=0;
498 a = *nr; // number of rows
504 for (i=1; i<=*nc; i++){
510 for (i=0; i<*nc; i++){
511 f = y[i]; //column number
517 } // starting value position in the Imatrix
518 for(j=1; j<=*nr; j++){
519 permuted[k] = Imatrix[c];
525 calc_twosample_ts(permuted,g,nc,nr,trial_ts,Tinitial,counter1);
528 void permute_array(int *array, int n) {
529 static int seeded = 0;
537 for (i = 0; i < n; i++) {
538 int selection = random() % (n - i);
539 int tmp = array[i + selection];
540 array[i + selection] = array[i];
545 void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,
546 double *Ts,double *Tinitial,double *counter) {
551 double C1[*nr][3], C2[*nr][3], storage[a],tool[a];
552 double nrows,ncols,gvalue, xbardiff=0, denom=0;
554 nrows = (double) *nr;
555 ncols = (double) *nc;
558 meanvar(Pmatrix,g,nr,nc,storage);
562 for (i=0; i<*nr;i++){
564 C1[i][1]=tool[i+*nr+*nr];
565 C1[i][2]=C1[i][1]/(gvalue-1);
567 C2[i][0]=tool[i+*nr];
568 C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
569 C2[i][2]=C2[i][1]/(ncols-gvalue+1);
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 Ts[i]=fabs(xbardiff/denom);
576 if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
582 void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
583 double temp[*nr], temp2[*nr],var[*nr],var2[*nr],a,b;
588 b = (double) (*nc-a);
590 for (i = 0; i<*nr; i++){
597 k = *nr; // number of rows
604 m*=k; // m = g * nr now
606 temp[i%k]=temp[i%k]+pmatrix[i];
609 temp2[i%k]=temp2[i%k]+pmatrix[i];
612 temp2[i]=temp2[i]-temp[i];
614 for (i=0;i<=*nr-1;i++){
616 store[i+*nr]=temp2[i]/b;
619 // That completes the mean calculations.
622 var[i%k]=var[i%k]+pow((pmatrix[i]-store[i%k]),2);
625 var2[i%k]=var2[i%k]+pow((pmatrix[i]-store[(i%k)+*nr]),2);
628 for (i=0;i<=*nr-1;i++){
629 store[i+2*k]=var[i]/(a-1);
630 store[i+3*k]=var2[i]/(b-1);
632 // That completes var calculations.
635 void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
636 double storage[][9]){
640 double store[a], tool[a], C1[*nr][3], C2[*nr][3];
641 double nrows,ncols,gvalue, xbardiff=0, denom=0;
643 nrows = (double) *nr;
644 ncols = (double) *nc;
647 meanvar(Imatrix,g,nr,nc,store);
652 for (i=0; i<*nr;i++){
653 C1[i][0]=tool[i]; //mean group 1
654 storage[i][0]=C1[i][0];
655 C1[i][1]=tool[i+*nr+*nr]; // var group 1
656 storage[i][1]=C1[i][1];
657 C1[i][2]=C1[i][1]/(gvalue-1);
658 storage[i][2]=sqrt(C1[i][2]);
660 C2[i][0]=tool[i+*nr]; // mean group 2
661 storage[i][4]=C2[i][0];
662 C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
663 storage[i][5]=C2[i][1];
664 C2[i][2]=C2[i][1]/(ncols-gvalue+1);
665 storage[i][6]=sqrt(C2[i][2]);
667 for (i=0; i<*nr; i++){
668 xbardiff = C1[i][0]-C2[i][0];
669 denom = sqrt(C1[i][2]+C2[i][2]);
670 initial[i]=fabs(xbardiff/denom);