Mercurial > repos > mmaiensc > ember

comparison GALAXY_FILES/tools/EMBER/Ember_Galaxy.c @ 3:037c3edda16e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mmaiensc
date Thu, 22 Mar 2012 13:49:52 -0400
parents 003f802d4c7d
children
comparison
equal deleted inserted replaced
2:1a84b8178b45 3:037c3edda16e
1 /*
2 * ember.c : Expectation Maximization of Binding and Expression pRofiles
3 * ----------------------------------------------------------------
4 * algorithm adapted from MEME algorithm: Bailey and Elkan, Proc Int Conf Intell Syst Mol Biol, 2:28-36 (1994).
5 * See Eqs. 3-13, plus "Implementation of MM" section for detailed description
6 * changes are as follows:
7 * background changes with column (j=1,...,numclasses), so it is a 2D histogram of the same size as the motif model; same for psuedocounts
8 * background model is determined by full list of genes (second input); z-values of genes in first input are subtracted
9 * update of background can be turned off by commenting out 8th line of z_update()
10 * in this case, also comment out index_genes at end of initialize(), as it will save some time
11 * sequences cannot overlap, since they are fixed per gene; this simplifies the enumeration of all possibilities
12 * a ZOOPS options is included (change in code, parameter zoops); this does not alter the algorithm,
13 * only adds a restriction on which genes are chosen as targets in choose_targets()
14 * general notes:
15 * the common dimensions are:
16 * numpeaks - number of peaks, size of peaklist, indexed by m
17 * peaklist[m].numgenes - number of genes matched to a peak, size of peaklist[m].genelist, indexed by n
18 * numclasses - number of comparisons made (i.e. motif width), indexed by j
19 * 5 - number of classifications in each comparison, assumed to be -2,-1,0,1,2, indexed by k
20 * 2 - number of models (motif=0 and background=1), indexed by l
21 * nmots - total number of motifs to find, indexed by i in main (always stored at l=0, as found motifs are erased)
22 * col_tog can be used to force the algorithm to ignore certain columns (behavior dimensions)
23 *
24 * things to add/consider:
25 * - best measure of motif quality score?
26 * - want to ensure that the same sequence can't get assigned to more than one motif?
27 */
28
29 #include "all.h"
30
31 typedef struct {
32 char id[50];
33 char name[500];
34 char chrom[50];
35 long start;
36 long end;
37 int sense;
38 double z;
39 double erase;
40 double score;
41 int index; // will refer to index of this id in bklist
42 int used; // set to 0 if a non-classified gene
43 int assigned; // set to 1 if assigned
44 int *values;
45 } gene;
46
47 typedef struct {
48 char info[5000];
49 char chrom[50];
50 long pkposn;
51 int numgenes;
52 int numused;
53 int assgngene; // set to # targets picked for this gene
54 gene *genelist;
55 } peak;
56
57 // FREE PARAMETERS TO CHOOSE IN COMMAND LINE
58 double mu; // scale background histogram by mu for psuedocounts
59 double tscale; // scales the threshold (to make choices more or less stringent)
60 int zoops; // set to 1 to enforce ZOOPS model (Zero Or One Per Strand)
61 int printall; // set to 1 to print all genes, not just targets
62 int mots; // number of motifs to search for
63 int repeat; // number of times to wean columns and repeat search
64 double rsca; // percent of average to set relative entropy cutoff
65 // END PARAMETERS TO CHOOSE
66
67 char *input, *background, *pref, *togs;
68 char *target_out, *list_out, *model_out, *annot_out;
69 long iseedval;
70 int numpeaks, numclasses, usedclasses, numbk, bkused, nmots;
71 peak *peaklist;
72 gene *bklist;
73 double ***hist, ***oldhist, lambda, t, **beta;
74 double *Beta;
75 int *col_tog, col_togl, set_togs, *col_tog_cpy;
76 int verbose;
77
78 void read_cmd_line( int argc, char **argv );
79 void print_usage();
80 void read_data( char *filename );
81 void read_bkgrnd( char *filename );
82 void initialize( long *seed );
83 void reinitialize( long *seed );
84 void beta_update();
85 void index_genes();
86
87 void hist_update();
88 void bkhist_update();
89 void lambda_update();
90 void toggle_update();
91 void toggle_reset();
92 void z_update();
93 void erase_update();
94 double lprob( int m, int n, int l );
95 double prob( int m, int n, int l );
96 void copy_hist();
97 double hist_diff();
98 double score( int m, int n );
99 void t_update();
100 double ElogL();
101 double ratio_score();
102 double rel_entropy();
103 double col_rel_entropy( int j );
104 double chisq();
105 int motif_count();
106
107 int assign_gene( int m, long *seed );
108 void score_genes();
109 void choose_targets();
110
111 void print_peak( FILE *fp, int i );
112 void print_gene( FILE *fp, int i, int j );
113 void print_state( FILE *fp, char *filename );
114 //void print_targets( FILE *fp, char *filename );
115 void print_model( FILE *fp, char *filename );
116 //void print_annotpeak( FILE *fp, char *filename, int motif );
117 char *nice_d( double val );
118
119 int main(int argc, char *argv[]){
120 int i,iter,r;
121 int pfreq;
122 double diff, conv = 1e-6;
123 long seed;
124 FILE *fout;
125
126 if(argc == 1 ){
127 print_usage();
128 }
129
130 /*
131 * read command line params
132 */
133 read_cmd_line( argc, argv );
134
135 /*
136 * allocate output names
137 */
138 // target_out = (char *) malloc((strlen(pref)+20)*sizeof(char));
139 // list_out = (char *) malloc((strlen(pref)+20)*sizeof(char));
140 // model_out = (char *) malloc((strlen(pref)+20)*sizeof(char));
141 // annot_out = (char *) malloc((strlen(pref)+20)*sizeof(char));
142
143 /*
144 * read in data, copy values from global to local variables
145 */
146 if( verbose != 0 ) printf("\nReading data and allocating\n");
147 read_data( input );
148 read_bkgrnd( background );
149 nmots = mots;
150 seed = iseedval;
151
152 /*
153 * allocate memory for other arrays
154 */
155 hist = (double ***) alloc3d(sizeof(double),2,numclasses,5);
156 oldhist = (double ***) alloc3d(sizeof(double),2,numclasses,5);
157 beta = (double **) alloc2d(sizeof(double),numclasses,5);
158 Beta = (double *) malloc(numclasses*sizeof(double));
159
160 /*
161 * find nmots motif models
162 */
163 for(i=1; i<= nmots; i++){
164 /*
165 * initialize, fill initial parameter values from z-guesses
166 */
167 if( verbose != 0 ) printf("Expression Pattern %i,\n Initializing\n", i);
168 if( i==1 ) initialize( &seed );
169 else reinitialize( &seed );
170 hist_update();
171 lambda_update();
172
173 /*
174 * now sample for awhile, repeating repeat times
175 * each repeat starts at the previously refined search (?)
176 */
177 for(r=0; r<= repeat; r++){
178 if( verbose != 0 ) printf(" Sampling\n");
179 if( verbose != 0 && repeat> 0 ) printf(" Repeat %i\n", r );
180 diff = 1;
181 pfreq = 10;
182 iter=1;
183 while( diff > conv ){
184 z_update(); // E step
185 copy_hist();
186 hist_update(); // M step
187 bkhist_update(); // M step
188 lambda_update(); // M step
189 diff = hist_diff();
190 if( iter%pfreq == 0 && verbose != 0 ) printf(" iteration %i, (normalized) ElogL %.2f, change %.2e\n", iter, ElogL(), diff);
191 iter++;
192 }
193 if( repeat > 0 && r< repeat ) toggle_update();
194 }
195
196 /*
197 * score and choose final targets
198 */
199 t_update();
200 if( verbose != 0 ) printf(" Done, scoring and choosing targets on threshold %.3f\n\n", t);
201 score_genes();
202 choose_targets();
203
204 /*
205 * print final targets
206 */
207 // sprintf( target_out,"%s-%i.targets", pref, i);
208 // sprintf( list_out,"%s-%i.list", pref, i);
209 // sprintf( model_out,"%s-%i.model", pref, i);
210 // sprintf( annot_out,"%s-%i.xls", pref, i);
211 print_state(fout, target_out);
212 print_model(fout, model_out);
213 // print_targets(fout, list_out);
214 // print_annotpeak( fout, annot_out, i);
215
216 /*
217 * erase found motifs (only if we're going around again)
218 */
219 if( i< nmots ) erase_update();
220 if( repeat > 0 ) toggle_reset();
221 }
222
223 }
224
225 /*==================== FUNCTIONS ======================*/
226
227 /*
228 * read command line inputs and options
229 */
230 void read_cmd_line( int argc, char **argv ){
231 int i, j;
232 int found_reqs=0;
233
234 // default values
235 mots = 1;
236 zoops = 0;
237 printall = 0;
238 mu = 0.1;
239 tscale = 1.0;
240 iseedval = make_seed();
241 set_togs = 0;
242 verbose = 1;
243 repeat = 0;
244 rsca = 1.0;
245
246 // parse options
247 for(i=1; i< argc; i++){
248 if( argv[i][0] != '-' ) print_usage();
249 i++;
250 if( i >= argc || argv[i][0] == '-' ) print_usage();
251 switch( argv[i-1][1] ){
252 case 'i':
253 input = (char *) malloc((strlen(argv[i])+5)*sizeof(char));
254 strcpy(input, argv[i]);
255 found_reqs++;
256 break;
257 case 'b':
258 background = (char *) malloc((strlen(argv[i])+5)*sizeof(char));
259 strcpy(background, argv[i]);
260 found_reqs++;
261 break;
262 case 'o':
263 target_out = (char *) malloc((strlen(argv[i])+5)*sizeof(char));
264 strcpy(target_out, argv[i]);
265 found_reqs++;
266 break;
267 case 'p':
268 model_out = (char *) malloc((strlen(argv[i])+5)*sizeof(char));
269 strcpy(model_out, argv[i]);
270 found_reqs++;
271 break;
272 case 'n':
273 mots = atoi(argv[i]);
274 if( mots <= 0 ){
275 printf("\nError: choose number of expression patterns to be at least 1\n");
276 print_usage();
277 }
278 break;
279 case 'z':
280 if( argv[i][0] == 'y' ) zoops = 1;
281 else if( argv[i][0] == 'n' ) zoops = 0;
282 else{
283 printf("\nError: if using -z option, choose y or n\n");
284 print_usage();
285 }
286 break;
287 case 'v':
288 if( argv[i][0] == 'y' ) verbose = 1;
289 else if( argv[i][0] == 'n' ) verbose = 0;
290 else{
291 printf("\nError: if using -v option, choose y or n\n");
292 print_usage();
293 }
294 break;
295 case 'm':
296 mu = atof(argv[i]);
297 if( mu<= 0 ){
298 printf("\nError: choose mu to be positive\n");
299 print_usage();
300 }
301 break;
302 case 't':
303 tscale = atof(argv[i]);
304 if( tscale<= 0 ){
305 printf("\nError: choose threshold to be positive\n");
306 print_usage();
307 }
308 break;
309 case 's':
310 iseedval = -1*atol(argv[i]);
311 if( iseedval>= 0 ){
312 printf("\nError: choose seed to be a large positive integer\n");
313 print_usage();
314 }
315 break;
316 case 'c':
317 set_togs = 1;
318 togs = (char *) malloc((strlen(argv[i])+5)*sizeof(char));
319 strcpy( togs, argv[i] );
320 break;
321 case 'a':
322 if( argv[i][0] == 'y' ) printall = 1;
323 else if( argv[i][0] == 'n' ) printall = 0;
324 else{
325 printf("\nError: if using -a option, choose y or n\n");
326 print_usage();
327 }
328 break;
329 case 'r':
330 repeat = atoi(argv[i]);
331 if( repeat< 0 ){
332 printf("\nError: choose repeat to be a non-negative integer\n");
333 print_usage();
334 }
335 break;
336 case 'R':
337 rsca = atof(argv[i]);
338 if( rsca <= 0 ){
339 printf("\nError: choose repeat scalar to be positive\n");
340 print_usage();
341 }
342 break;
343 default:
344 printf("\nError: unkown option: -%c\n", argv[i-1][1]);
345 print_usage();
346 }
347 }
348 // double check that the correct number of required inputs has been read in
349 if( found_reqs != 4 ){
350 printf("\nError: improper number of required inputs found\n");
351 print_usage();
352 }
353 }
354
355 /*
356 * usage info
357 */
358 void print_usage(){
359 printf("\nUsage: -i matched_data -b background_data -o target_output_name -p model_output_name [option value]\n");
360 printf("Options:\n");
361 printf(" -n number of expression patterns to find (positive integer, default 1)\n");
362 printf(" -z enforce ZOOPS model (y or n, default n)\n");
363 printf(" -m scalar for psuedocounts (positive real number, default 0.1)\n");
364 printf(" -t scalar for threshold (positive real number, default 1.0)\n");
365 printf(" -s seed (positive long integer, default chosen from date and time)\n");
366 printf(" -c column toggle file, used to ignore certain classifications (default: all on)\n");
367 printf(" -v verboseness, turn on/off printing to stdout (y or n, default: y)\n");
368 printf(" -a print all genes (and scores), rather than only those above the threshold (y or n, default: n)\n");
369 printf(" -r repeat search: after convergence of each expression pattern, turn off columns with lower relative entropy\n");
370 printf(" than average*scalar and repeat expression pattern search (non-negative integer, default 0)\n");
371 printf(" -R relative entropy scalar: if using the -r option, multiplies the average relative entropy\n");
372 printf(" by this value to set the threshold for turning off columns (positive real number, default 1.0)\n");
373 printf("\n");
374 exit(0);
375 }
376
377 /*
378 * function: read data from filename into the peaklist data structure
379 * notes on assumed format of data:
380 * peak lines start with "PEAK:" field
381 * gene lines start with "GENE:" field
382 * gene lines with only 1 classification are ignored
383 * assume all fields are <5000 characters
384 */
385 void read_data( char *filename ){
386 int j,m,n;
387 char val[5000], other[5000];
388 FILE *fp;
389
390 /*
391 * first count number of peaks, and allocate
392 */
393 fp = fopen(filename,"rt");
394 if( fp == NULL ){
395 printf("Error: can't open file %s\n", filename);
396 exit(0);
397 }
398 numpeaks = 0;
399 while( fscanf(fp, "%s", &val ) != EOF ){
400 if( contains( val, "PEAK" ) == 1 ){
401 numpeaks++;
402 }
403 }
404 fclose(fp);
405 peaklist = (peak *) malloc(numpeaks*sizeof(peak));
406
407 /*
408 * now count number of genes w/>1 classification, and allocate
409 * also count number of classifications, and allocate that at the end
410 * n tracks number of genes, j field; subtract genes that don't have full classifications
411 * m tracks peak index; record peak info
412 */
413 fp = fopen(filename,"rt");
414 n=-1;
415 m=-1;
416 j=0;
417 sprintf(other,"NA");
418 numclasses=0;
419 while( fscanf(fp, "%s", &val ) != EOF ){
420 if( contains( val, "PEAK" ) == 1 ){
421 if( j==7 ) ;
422 else{
423 if( j-6 > numclasses ) numclasses = j-6;
424 }
425 if( m>=0 ){
426 peaklist[m].numgenes = n;
427 peaklist[m].genelist = (gene *) malloc( peaklist[m].numgenes*sizeof(gene));
428 }
429 n=0;
430 j=0;
431 if( strcmp(other, "NA") != 0 ) strcpy(peaklist[m].info, other);
432 m++;
433 sprintf(other,"NA");
434 }
435 else if( contains( val, "GENE" ) == 1 ){
436 n++;
437 j=0;
438 }
439 else if( n==0 ){
440 if( j==0 ) strcpy(peaklist[m].chrom, val);
441 if( j==1 ) peaklist[m].pkposn = atol(val);
442 if( j> 1 ){
443 if( strcmp(other, "NA") == 0 ) strcpy(other, val);
444 else sprintf(other,"%s %s", other, val);
445 }
446 j++;
447 }
448 else{
449 j++;
450 }
451 }
452 peaklist[m].numgenes = n;
453 peaklist[m].genelist = (gene *) malloc( peaklist[m].numgenes*sizeof(gene));
454 strcpy( peaklist[m].info, other);
455 fclose(fp);
456
457 for(m=0; m< numpeaks; m++){
458 for(n=0; n< peaklist[m].numgenes; n++){
459 peaklist[m].genelist[n].values = (int *) malloc(numclasses*sizeof(int));
460 }
461 }
462
463 /*
464 * now read in gene info
465 */
466 fp = fopen(filename,"rt");
467 m=-1;
468 j=0;
469 while( fscanf(fp, "%s", &val ) != EOF ){
470 if( contains( val, "PEAK" ) == 1 ){
471 if( j==7 ) peaklist[m].genelist[n].used = 0;
472 if( j==6+numclasses) peaklist[m].genelist[n].used = 1;
473 n=-1;
474 j=0;
475 m++;
476 }
477 else if( contains( val, "GENE" ) == 1 ){
478 if( j==7 ) peaklist[m].genelist[n].used = 0;
479 if( j==6+numclasses ) peaklist[m].genelist[n].used = 1;
480 n++;
481 j=0;
482 }
483 else if( n==-1 ){
484 }
485 else {
486 if( j==0 ) strcpy( peaklist[m].genelist[n].id, val );
487 if( j==1 ) strcpy( peaklist[m].genelist[n].name, val );
488 if( j==2 ) strcpy( peaklist[m].genelist[n].chrom, val );
489 if( j==3 ) peaklist[m].genelist[n].start = atol(val);
490 if( j==4 ) peaklist[m].genelist[n].end = atol(val);
491 if( j==5 ) peaklist[m].genelist[n].sense = atoi(val);
492 if( j>=6 ) peaklist[m].genelist[n].values[j-6] = atoi(val);
493 j++;
494 }
495 }
496 if( j==7 ) peaklist[m].genelist[n].used = 0;
497 if( j==6+numclasses) peaklist[m].genelist[n].used = 1;
498 fclose(fp);
499
500 /*
501 * now count how many genes are to be used for each
502 */
503 for(m=0; m< numpeaks; m++){
504 peaklist[m].numused = 0;
505 for(n=0; n< peaklist[m].numgenes; n++){
506 peaklist[m].numused+= peaklist[m].genelist[n].used;
507 }
508 }
509 }
510
511 void read_bkgrnd( char *filename ){
512 int j,n;
513 char val[500];
514 FILE *fp;
515
516 /*
517 * first count number of genes, and allocate
518 */
519 numbk = count_lines( fp, filename );;
520 bklist = (gene *) malloc(numbk*sizeof(gene));
521
522 for(n=0; n< numbk; n++){
523 bklist[n].values = (int *) malloc(numclasses*sizeof(int));
524 }
525
526 /*
527 * now read in gene info
528 */
529 fp = fopen(filename,"rt");
530 n=-1;
531 j=0;
532 while( fscanf(fp, "%s", &val ) != EOF ){
533 if( contains( val, "_a" ) == 1 ){
534 if( j==7 ) bklist[n].used = 0;
535 if( j==6+numclasses ) bklist[n].used = 1;
536 j=1;
537 n++;
538 strcpy( bklist[n].id, val );
539 }
540 else{
541 if( j==1 ) strcpy( bklist[n].name, val );
542 if( j==2 ) strcpy( bklist[n].chrom, val );
543 if( j==3 ) bklist[n].start = atol(val);
544 if( j==4 ) bklist[n].end = atol(val);
545 if( j==5 ) bklist[n].sense = atoi(val);
546 if( j>=6 ) bklist[n].values[j-6] = atoi(val);
547 j++;
548 }
549 }
550 if( j==7 ) bklist[n].used = 0;
551 if( j==6+numclasses) bklist[n].used = 1;
552 fclose(fp);
553
554 /*
555 * now count how many genes are to be used for each
556 */
557 bkused = 0;
558 for(n=0; n< numbk; n++){
559 bkused+= bklist[n].used;
560 }
561 }
562
563 /*
564 * make initial guesses of z-values, set erase to 1
565 * read in column toggle file if provided
566 * also histogram background and scale into beta
567 */
568 void initialize( long *seed ){
569 int j,m,n;
570 int val;
571 FILE *fp;
572
573 for(m=0; m< numpeaks; m++){
574 for(n=0; n< peaklist[m].numgenes; n++){
575 if( peaklist[m].genelist[n].used == 1 ){
576 peaklist[m].genelist[n].erase = 1.0;
577 peaklist[m].genelist[n].z = 0.0;
578 }
579 }
580 if( peaklist[m].numused > 0 ){
581 n = assign_gene( m, seed );
582 peaklist[m].genelist[n].z = 1.0; // pick one random gene from each set to read into motif model
583 }
584 }
585
586 col_tog = (int *) malloc(numclasses*sizeof(int));
587 col_tog_cpy = (int *) malloc(numclasses*sizeof(int));
588 if( set_togs == 0 ){
589 for(j=0; j< numclasses; j++){
590 col_tog[j] = 1;
591 col_tog_cpy[j] = 1;
592 }
593 usedclasses = numclasses;
594 }
595 else{
596 // first count number of fields
597 col_togl = 0;
598 fp = fopen(togs,"rt");
599 if( fp == NULL ){
600 printf("Error: can't open file %s\n", togs);
601 exit(0);
602 }
603 j=0;
604 while( fscanf(fp,"%i", &val) != EOF ){
605 col_togl++;
606 }
607 fclose(fp);
608 if( col_togl != numclasses ){
609 printf("Error: number of fields in toggle file (%i) is different from number of classifications in data (%i)\n", col_togl, numclasses);
610 exit(0);
611 }
612 // now read in values
613 fp = fopen(togs,"rt");
614 j=0;
615 usedclasses = 0;
616 while( fscanf(fp,"%i", &val) != EOF ){
617 col_tog[j] = val;
618 col_tog_cpy[j] = val;
619 j++;
620 if( val == 1 ) usedclasses++;
621 if( val != 1 && val != 0 ){
622 printf("Error: choose the toggle values to be 0 or 1\n");
623 exit(0);
624 }
625 }
626 fclose(fp);
627 }
628
629 for(n=0; n< numbk; n++){
630 bklist[n].z = 1.0;
631 }
632 beta_update();
633 if( verbose != 0 ) printf(" Indexing genes in first pass\n");
634 index_genes();
635 }
636
637 /*
638 * make initial guesses of z-values for a second, etc motif
639 * recount number of used genes per peak to those with erase > 0.5
640 */
641 void reinitialize( long *seed ){
642 int m,n;
643 for(m=0; m< numpeaks; m++){
644 peaklist[m].numused = 0;
645 for(n=0; n< peaklist[m].numgenes; n++){
646 peaklist[m].genelist[n].z = 0.0;
647 if( peaklist[m].genelist[n].used == 1 && peaklist[m].genelist[n].erase> 0.5 ){
648 peaklist[m].numused++;
649 }
650 }
651 if( peaklist[m].numused > 0 ){
652 n = assign_gene( m, seed );
653 peaklist[m].genelist[n].z = 1.0; // pick one random gene from each set to read into motif model
654 }
655 }
656 for(n=0; n< numbk; n++){
657 bklist[n].z = 1.0;
658 }
659 }
660
661 /*
662 * histogram background into beta
663 */
664 void beta_update(){
665 int j,k,n;
666 double tot=0.0;
667 bkhist_update();
668 for(j=0; j< numclasses; j++){
669 // if( col_tog[j] == 1 ){
670 Beta[j] = 0.0;
671 for(k=0; k< 5; k++){
672 beta[j][k] = hist[1][j][k]*mu;
673 Beta[j]+= beta[j][k];
674 }
675 // }
676 }
677 }
678
679 /*
680 * for all the genes in peaklist[*].genelist[*], record the index of that gene in bklist
681 */
682 void index_genes(){
683 int m,n,nu;
684 int found;
685 for(m=0; m< numpeaks; m++){
686 for(n=0; n< peaklist[m].numgenes; n++){
687 if( peaklist[m].genelist[n].used == 1 ){
688 found = 0;
689 nu = 0;
690 while(nu< numbk && found == 0){
691 if( strcmp(peaklist[m].genelist[n].id, bklist[nu].id) == 0 ){
692 peaklist[m].genelist[n].index = nu;
693 found = 1;
694 }
695 nu++;
696 }
697 }
698 }
699 }
700 }
701
702 /*
703 * update motif model histogram
704 */
705 void hist_update(){
706 int j,k,m,n;
707 double tot=0;
708 for(j=0; j< numclasses; j++){
709 // if( col_tog[j] == 1 ){
710 for(k=0; k< 5; k++){
711 hist[0][j][k] = 0.0;
712 }
713 // }
714 }
715 for(m=0; m< numpeaks; m++){
716 for(n=0; n< peaklist[m].numgenes; n++){
717 if( peaklist[m].genelist[n].used == 1 ){
718 tot+=peaklist[m].genelist[n].erase*peaklist[m].genelist[n].z;
719 for(j=0; j< numclasses; j++){
720 // if( col_tog[j] == 1 ){
721 k = peaklist[m].genelist[n].values[j]+2;
722 hist[0][j][k]+= peaklist[m].genelist[n].erase*peaklist[m].genelist[n].z;
723 // }
724 }
725 }
726 }
727 }
728 for(j=0; j< numclasses; j++){
729 // if( col_tog[j] == 1 ){
730 for(k=0; k< 5; k++){
731 hist[0][j][k]+= beta[j][k];
732 hist[0][j][k]/= (tot+Beta[j]);;
733 }
734 // }
735 }
736 }
737
738 /*
739 * update background hist using z values from bklist
740 */
741 void bkhist_update(){
742 int j,k,n;
743 double tot=0.0;
744 for(j=0; j< numclasses; j++){
745 // if( col_tog[j] == 1 ){
746 for(k=0; k< 5; k++){
747 hist[1][j][k] = 0;
748 }
749 // }
750 }
751 for(n=0; n< numbk; n++){
752 if( bklist[n].used == 1 ){
753 tot+=bklist[n].z;
754 for(j=0; j< numclasses; j++){
755 // if( col_tog[j] == 1 ){
756 k = bklist[n].values[j]+2;
757 hist[1][j][k]+= bklist[n].z;
758 // }
759 }
760 }
761 }
762 for(j=0; j< numclasses; j++){
763 // if( col_tog[j] == 1 ){
764 for(k=0; k< 5; k++){
765 hist[1][j][k]/= tot;
766 }
767 // }
768 }
769 }
770
771 /*
772 * update lambda value
773 */
774 void lambda_update(){
775 int m,n;
776 double tot=0.0;
777 lambda=0;
778 for(m=0; m< numpeaks; m++){
779 for(n=0; n< peaklist[m].numgenes; n++){
780 if( peaklist[m].genelist[n].used == 1 ){
781 lambda+= peaklist[m].genelist[n].z;
782 tot++;
783 }
784 }
785 }
786 lambda/=tot;
787 }
788
789 /*
790 * update the toggle columns based on which columns have better than average relative entropy
791 */
792 void toggle_update(){
793 int j;
794 double avg=0.0;
795 avg = rsca*rel_entropy();
796 usedclasses=0;
797 for(j=0; j< numclasses; j++){
798 if( col_tog[j] == 1 ){
799 if( col_rel_entropy(j) < avg ) col_tog[j] = 0;
800 if( col_tog[j] == 1 ) usedclasses++;
801 }
802 }
803 }
804
805 /*
806 * reset col_tog values with their original ones
807 */
808 void toggle_reset(){
809 int j;
810 usedclasses=0;
811 for(j=0; j< numclasses; j++){
812 col_tog[j] = col_tog_cpy[j];
813 if( col_tog[j] == 1 ) usedclasses++;
814 }
815 }
816
817 /*
818 * update z
819 */
820 void z_update(){
821 int j,m,n;
822 double a, c, denom;
823 for(m=0; m< numpeaks; m++){
824 for(n=0; n< peaklist[m].numgenes; n++){
825 if( peaklist[m].genelist[n].used == 1 ){
826 // new z value for motif model
827 a = lprob(m,n,0)+log(lambda);
828 c = lprob(m,n,1)+log(1-lambda);
829 if( a> c ){
830 denom = a + log(1+exp(c-a));
831 }
832 else{
833 denom = c + log(1+exp(a-c));
834 }
835 peaklist[m].genelist[n].z = exp( lprob(m,n,0)+log(lambda) - denom );
836 // peaklist[m].genelist[n].z = prob(m,n,0)*lambda/( prob(m,n,0)*lambda + prob(m,n,1)*(1-lambda) );
837 // new z value for corresponding gene in background model
838 bklist[peaklist[m].genelist[n].index].z = 1.0-peaklist[m].genelist[n].z;
839 }
840 }
841 }
842 }
843
844 /*
845 * update erase parameters
846 */
847 void erase_update(){
848 int j,m,n;
849 for(m=0; m< numpeaks; m++){
850 for(n=0; n< peaklist[m].numgenes; n++){
851 if( peaklist[m].genelist[n].used == 1 ){
852 peaklist[m].genelist[n].erase*= (1-peaklist[m].genelist[n].z);
853 }
854 }
855 }
856 }
857
858 /*
859 * calculate (log) (multinomial) probability of sequence m,n in model l
860 */
861 double lprob( int m, int n, int l ){
862 int j,k;
863 double val=0.0;
864 for(j=0; j< numclasses; j++){
865 if( col_tog[j] == 1 ){
866 k = peaklist[m].genelist[n].values[j]+2;
867 val+= log(hist[l][j][k]);
868 }
869 }
870 return val;
871 }
872 double prob( int m, int n, int l ){
873 int j,k;
874 double val=1.0;
875 for(j=0; j< numclasses; j++){
876 if( col_tog[j] == 1 ){
877 k = peaklist[m].genelist[n].values[j]+2;
878 val*= hist[l][j][k];
879 }
880 }
881 return val;
882 }
883
884 /*
885 * copy hist into oldhist before updating
886 */
887 void copy_hist(){
888 int j,k,l;
889 for(l=0; l< 2; l++){
890 for(j=0; j< numclasses; j++){
891 // if( col_tog[j] == 1 ){
892 for(k=0; k< 5; k++){
893 oldhist[l][j][k] = hist[l][j][k];
894 }
895 // }
896 }
897 }
898 }
899
900 /*
901 * compute change in hist
902 */
903 double hist_diff(){
904 int j,k,l;
905 double val=0.0;
906 for(l=0; l< 2; l++){
907 for(j=0; j< numclasses; j++){
908 if( col_tog[j] == 1 ){
909 for(k=0; k< 5; k++){
910 val+= fabs( oldhist[l][j][k]-hist[l][j][k] );
911 }
912 }
913 }
914 }
915 return val;
916 }
917
918 /*
919 * calculate the score of a sequence
920 */
921 double score( int m, int n ){
922 int j,k;
923 double val=0.0;
924 for(j=0; j< numclasses; j++){
925 if( col_tog[j] == 1 ){
926 k = peaklist[m].genelist[n].values[j]+2;
927 val+= log( hist[0][j][k]/hist[1][j][k] );
928 }
929 }
930 return val;
931 }
932
933 /*
934 * calculate value of t (threshold)
935 */
936 void t_update(){
937 t = tscale*log( (1-lambda)/lambda );
938 }
939
940 /*
941 * quality scores of a motif, averaging over only those sequences picked as targets
942 * ElogL: log-likelihood of the combination of models, averaged by # genes
943 * ratio_score: log-likelihood ratio of motif vs background models
944 * rel_entropy: relative entropy of motif to background
945 * col_rel_entropy: same as rel_entropy, but for just one column (rel_entropy is average over all cols)
946 * chisq: chi-square statistic
947 */
948 double ElogL(){
949 int m,n;
950 double val=0.0, tot=0.0;
951 for(m=0; m< numpeaks; m++){
952 for(n=0; n< peaklist[m].numgenes; n++){
953 if( peaklist[m].genelist[n].used == 1 ){
954 val+= peaklist[m].genelist[n].z*lprob(m, n, 0 ) \
955 + (1-peaklist[m].genelist[n].z)*lprob(m, n, 1 )\
956 + peaklist[m].genelist[n].z*log(lambda) \
957 + (1-peaklist[m].genelist[n].z)*log(1-lambda);
958 tot++;
959 }
960 }
961 }
962 return val/tot;
963 }
964
965 double ratio_score(){
966 int m,n;
967 double val=0.0;
968 double tot=0.0;
969 for(m=0; m< numpeaks; m++){
970 for(n=0; n< peaklist[m].numgenes; n++){
971 if( peaklist[m].genelist[n].assigned == 1 ){
972 val+= lprob(m,n,0)-lprob(m,n,1);
973 tot++;
974 }
975 }
976 }
977 return val/tot;
978 }
979
980 double rel_entropy(){
981 int j,k;
982 double val=0.0;
983 for(j=0; j< numclasses; j++){
984 if( col_tog[j] == 1 ){
985 val+= col_rel_entropy( j );
986 }
987 }
988 return val/usedclasses;
989 }
990
991 double col_rel_entropy( int j ){
992 int k;
993 double val=0.0;
994 for(k=0; k< 5; k++){
995 if( hist[1][j][j] > 0 && hist[0][j][k] > 0 ){
996 // do in log base 2 for bits
997 val+= hist[0][j][k]*log(hist[0][j][k]/hist[1][j][k])/log(2);
998 }
999 }
1000 return val;
1001 }
1002
1003 double chisq(){
1004 int j, k;
1005 double val=0.0;
1006 for(j=0; j< numclasses; j++){
1007 if( col_tog[j] == 1 ){
1008 for(k=0; k< 5; k++){
1009 val+= pow(hist[0][j][k] - hist[1][j][k],2)/hist[1][j][k];
1010 }
1011 }
1012 }
1013 return val;
1014 }
1015
1016 /*
1017 * number of members of a motif model
1018 */
1019 int motif_count(){
1020 int m,n;
1021 int tot=0;
1022 for(m=0; m< numpeaks; m++){
1023 for(n=0; n< peaklist[m].numgenes; n++){
1024 if( peaklist[m].genelist[n].assigned == 1 ){
1025 tot++;
1026 }
1027 }
1028 }
1029 return tot;
1030 }
1031
1032 /*
1033 * assign a random used gene: one that is used and hasn't been erased
1034 */
1035 int assign_gene( int m, long *seed ){
1036 int n;
1037 n = floor( peaklist[m].numgenes*ran1(seed) );
1038 while( peaklist[m].genelist[n].used == 0 || peaklist[m].genelist[n].erase < 0.5 ){
1039 n = floor( peaklist[m].numgenes*ran1(seed) );
1040 }
1041 return n;
1042 }
1043
1044 /*
1045 * score all genes using the models
1046 */
1047 void score_genes(){
1048 int m,n;
1049 for(m=0; m< numpeaks; m++){
1050 for(n=0; n< peaklist[m].numgenes; n++){
1051 if( peaklist[m].genelist[n].used == 1 ){
1052 peaklist[m].genelist[n].score= score( m, n );
1053 }
1054 }
1055 }
1056 }
1057
1058 /*
1059 * choose targets using the scores, and ZOOPS setting
1060 * zoops = 1 -> 0 or 1 target/strand
1061 * else -> all above threshold
1062 */
1063 void choose_targets(){
1064 int m,n;
1065 int best;
1066 for(m=0; m< numpeaks; m++){
1067 peaklist[m].assgngene = 0;
1068 best=-1;
1069 for(n=0; n< peaklist[m].numgenes; n++){
1070 if( peaklist[m].genelist[n].used == 1 ){
1071 if( peaklist[m].genelist[n].score > t ){
1072 peaklist[m].genelist[n].assigned = 1;
1073 peaklist[m].assgngene++;
1074 }
1075 else peaklist[m].genelist[n].assigned = 0;
1076 if( best == -1 ) best = n;
1077 else if( peaklist[m].genelist[n].score > peaklist[m].genelist[best].score ) best = n;
1078 }
1079 }
1080 if( zoops == 1 ){
1081 for(n=0; n< peaklist[m].numgenes; n++){
1082 if( peaklist[m].genelist[n].used == 1 && n != best ) peaklist[m].genelist[n].assigned = 0;
1083 }
1084 if( peaklist[m].genelist[best].score > t ) peaklist[m].assgngene = 1;
1085 else peaklist[m].assgngene = 0;
1086 }
1087 }
1088 }
1089
1090 /*
1091 * functions to print to output (fp must already be opened for top two)
1092 */
1093 void print_peak( FILE *fp, int m ){
1094 int n;
1095 fprintf(fp,"PEAK: %s %ld %s\n", peaklist[m].chrom, peaklist[m].pkposn, peaklist[m].info);
1096 for(n=0; n< peaklist[m].numgenes; n++){
1097 if( printall == 0 ){
1098 if( peaklist[m].genelist[n].used == 1 && peaklist[m].genelist[n].assigned == 1 ){
1099 print_gene( fp, m, n );
1100 }
1101 }
1102 if( printall == 1 ){
1103 if( peaklist[m].genelist[n].used == 1 ){
1104 print_gene( fp, m, n );
1105 }
1106 }
1107 }
1108 }
1109
1110 void print_gene( FILE *fp, int m, int n ){
1111 int j;
1112 if( peaklist[m].genelist[n].used == 1 ){
1113 fprintf(fp, "TGENE: %0.3f %s %s %s %ld %ld %i", peaklist[m].genelist[n].score, peaklist[m].genelist[n].id, peaklist[m].genelist[n].name, peaklist[m].genelist[n].chrom, peaklist[m].genelist[n].start, peaklist[m].genelist[n].end, peaklist[m].genelist[n].sense);
1114 for(j=0; j< numclasses; j++){
1115 if( col_tog[j] == 0 ) fprintf(fp," %i ", peaklist[m].genelist[n].values[j]);
1116 else fprintf(fp," %i", peaklist[m].genelist[n].values[j]);
1117 }
1118 fprintf(fp,"\n");
1119 }
1120 }
1121
1122 void print_state(FILE *fp, char *filename){
1123 int m;
1124 fp = fopen(filename,"w");
1125 for(m=0; m< numpeaks; m++){
1126 print_peak(fp, m);
1127 }
1128 fclose(fp);
1129 }
1130
1131 /*
1132 * this prints a list of <region id> <name>
1133 */
1134 /****************
1135 void print_targets( FILE *fp, char *filename ){
1136 int m,n,nu;
1137 fp = fopen(filename,"w");
1138 for(m=0; m< numpeaks; m++){
1139 if( peaklist[m].assgngene > 0 ){
1140 fprintf(fp,"%s ", peaklist[m].regid);
1141 n=0;
1142 nu=0;
1143 while( n< peaklist[m].numgenes ){
1144 if( peaklist[m].genelist[n].used == 1 && peaklist[m].genelist[n].assigned == 1 ){
1145 if( nu< peaklist[m].assgngene-1 ) fprintf(fp,"%s,", peaklist[m].genelist[n].name);
1146 else fprintf(fp,"%s\n", peaklist[m].genelist[n].name);
1147 nu++;
1148 }
1149 n++;
1150 }
1151 }
1152 }
1153 fclose(fp);
1154 }
1155 *************/
1156
1157 /*
1158 * print the motif model and its score
1159 */
1160 void print_model( FILE *fp, char *filename ){
1161 int j,k;
1162 fp = fopen(filename,"w");
1163 fprintf(fp,"comparison\\class; avg relative entropy %0.3f, %i genes, threshold %0.3f, quality %0.1f\n", rel_entropy(), motif_count(), t, rel_entropy()*motif_count() );
1164 fprintf(fp," ++ + 0 - --\n");
1165 for(j=0; j< numclasses; j++){
1166 if( j< 10 ) fprintf(fp," ");
1167 if( col_tog[j] == 0 ) fprintf(fp,"#");
1168 else fprintf(fp," ");
1169 fprintf(fp,"%i,%.3f",j+1, col_rel_entropy(j) );
1170 for(k=4; k>=0; k--){
1171 if( hist[1][j][k] > 0 ) fprintf(fp,"%s", nice_d(log(hist[0][j][k]/hist[1][j][k])) );
1172 else fprintf(fp," NA");
1173 }
1174 fprintf(fp,"\n");
1175 }
1176 fclose(fp);
1177 }
1178
1179 /*
1180 * prints a double with 3 digits, aligned by spaces
1181 * longest would be -999.000, so 8 characters
1182 * shortest would be 1.000, 5 characters
1183 * up to 3 spaces
1184 */
1185 char *nice_d( double val ){
1186 char *returnval;
1187 char *dval;
1188 char *spaces;
1189 int i, numspaces=0;
1190
1191 returnval = (char *) malloc(10*sizeof(char));
1192 dval = (char *) malloc(10*sizeof(char));
1193 spaces = (char *) malloc(10*sizeof(char));
1194
1195 sprintf(dval,"%.3f", val);
1196
1197 if( val >= 0 ) numspaces++;
1198 if( fabs(val) < 100 ){
1199 if( fabs(val) < 10 ) numspaces+=2;
1200 else numspaces++;
1201 }
1202 if( numspaces == 0 ) sprintf(spaces,"");
1203 if( numspaces == 1 ) sprintf(spaces," ");
1204 if( numspaces == 2 ) sprintf(spaces," ");
1205 if( numspaces == 3 ) sprintf(spaces," ");
1206 sprintf(returnval,"%s%s", spaces, dval);
1207 return returnval;
1208 }
1209
1210 /*
1211 * print the peaks only in annotation format, for use in comppeaks, subpeaks, etc
1212 */
1213 /**************
1214 void print_annotpeak( FILE *fp, char *filename, int motif ){
1215 int m;
1216 fp = fopen(filename,"w");
1217 for(m=0; m< numpeaks; m++){
1218 if( peaklist[m].assgngene > 0 ){
1219 fprintf(fp,"%s-%i %s %s %ld %ld %ld ", peaklist[m].filename, motif, peaklist[m].regid, peaklist[m].chrom, peaklist[m].pkposn, peaklist[m].start, peaklist[m].end );
1220 fprintf(fp,"%f 1.0 1.0 1 1 1.0 ", peaklist[m].height );
1221 fprintf(fp,"1 1.0 1.0 1.0 1 -1.0 poop 1\n");
1222 }
1223 }
1224 close(fp);
1225 }
1226 *************/
1227