1 /* Code by Ethan Wolf 1995. Graphics code by David Wilson. */
2 /*** See dimer_trimer_like.c, function calculate_likelihood to see how */
3 /*** to compute likelihood line number before call DimerTrimerScore() */
7 #define XCOORD_SCORED_BY 430
8 #define XCOORD_SEQSCORES 0
9 #define YCOORD_SEQSCORES 55 /* Each char takes "15" so put us below printing */
10 #define WIDTH_ADD_ON 30 /** Empty space after seq in main window. **/
11 #define xmin 120 /* Want this far enough over so have room for seqscoreb */
12 /** seqscoreb is 108 wide. xmin used to be 60 **/
16 #include "sc2seq.h" /* This also includes scconst.h, sc.h, interface.h */
18 #include <X11/keysym.h>
24 Window window, zoom, printb, quit, next, prev, autob, zoomb, tableb, paramb, helpb, help, seqscoreb;
25 GC gc, gcgray, gcdash;
29 #include <sys/types.h>
37 #include "debug_printout_flag.h"
40 /* FILE *debug_out; */
43 /*************************** MULTICOIL LIKELIHOOD STUFF. *********************/
44 double both_determinants[2][MAX_CLASS_NUMBER];
45 double *determinants= &both_determinants[0][0];
47 double both_inverse_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
48 double *inverse_covars= &both_inverse_covars[0][0][0][0];
50 double both_means_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
51 double *means_submatrix= &both_means_submatrix[0][0][0];
53 double init_class_prob[MAX_CLASS_NUMBER];
55 int combine_dist[MAX_TABLE_NUMBER];
56 int number_score_dim=0;
57 int num_dist[MAX_TABLE_NUMBER];
58 /* Used to detrmine the number of different library */
59 /* distance values to do gaussian param over. */
63 FILE *ftotal_like=NULL; /* Extern for now. ***/
64 double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER];
65 /* When ftotal_like is read from get_defaults(), the */
66 /* i,j entry will hold the total gaussian like when */
67 /* assuming dimers have init_prob = */
68 /* MIN_DIMER_PROB + i*DIMER_GRID_STEP */
69 /* and trimers have init_prob = */
70 /* MIN_TRIMER_PROB + j*TRIMER_GRID_STEP */
71 /* These constants are defined in scconst.h. */
72 /* The value is updated during "output_seq()". */
75 /************************** General Scoring parameters. ********************/
76 int functnum=0, pair_functnum=0;
78 double log_bound = 0; /*In PRN_MODE only scores > bound output to logfile.*/
80 int WINDOW = PAIRWINDOW;
81 int window_length[MAX_TABLE_NUMBER];
83 double SCALE0 = DEFAULT_SCALE0;
84 double scale0s[MAX_TABLE_NUMBER]; /** Both used as bonnie's scale0 and **/
85 double scale0p[MAX_TABLE_NUMBER]; /** in computing posterior prob. from **/
86 /** prior as sigma^2(prior)=scale0*sigma^2 */
87 /** where these are the standard deviation *
88 /** of the n*(prior_prob) and the distib */
89 /** of number of hits when sampling */
91 int log_offset_to_use;
92 int offset_to_use=7; /* Default offset in DimerTrimerScore is all. */
93 int avg_max=1; /* Default means use max coil scores in log file. */
94 /** Value of 0 means do average, 1 means do max. **/
95 /** 2 means do both. */
96 int method = MultiCoil;
97 int main_method = MultiCoil;
98 int main_preprocessor_method = MultiCoil;
99 int preprocessor_method = MultiCoil;
104 /*********************** Parameters involved in Coil Scores ***********/
105 int by_coil_or_seq =0; /*** To determe what type of coil_out file to do. **/
106 int start_coil_at_reg_shift =0; /* For averaging scores over coil length. */
107 int Coil_Score =0; /* 0 means for MultiCoil do residues score, 1 means */
108 /* average scores over coils, 2 means take max coil */
114 char *methodname[] = {"MultiCoil", "Coil","NEWCOILS", "HMMCoil","ActualCoil",
116 "PairCoilDiffAvg"}; /* "Coil" is PairCoil or */
125 double seq_scores[MAX_CLASS_NUMBER];
127 double *scores; /** The space these pointers points to is alocated by **/
128 /** using static variable arrays in ScoreSeq */
129 /** all_scores is scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM +1]*/
130 /** scores is scores[current_score_algor][MAXSEQLEN][POSNUM+1] **/
132 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
133 /** Since preprocess_likelihood is used everywhere, just */
134 /** make it an external variable. */
135 /** This is given a value in ScoreSeq and sometime in */
137 double *preprocess_like;
138 /* This will point to the current table version */
139 /* of all_preprocess_like */
140 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
144 /************************ XWindow stuff. ***********************************/
146 int ps_res_per_line=0; /* 0 is a flag that it wasn't set in get_defaults */
147 /* so should just print it using nominal_length */
148 /* as done on the screen. If it is set to -1 */
149 /* then should print it all on one line. */
154 " *********************************",
156 " ** MultiCoil (version 1.0) **",
159 " *********************************",
161 "MultiCoil is described in:",
162 " Ethan Wolf, Peter S. Kim, Bonnie Berger,",
163 " `MultiCoil: A Program for Predicting Two- and",
164 " Three-stranded Coiled Coils'.",
166 "Email inquiries to: multicoil-help@theory.lcs.mit.edu.",
171 char *disclaimer[] = {
172 " *********************************",
174 " * MultiCoil Release 1.0 Notes: *",
176 " *********************************",
178 " The scores obtained using the PairCoil sub-method",
179 " will not match the scores of the official PairCoil",
180 " release. Parameters were optimized based on the",
181 " performance of MultiCoil. However, any differences",
182 " in the PairCoil predictions should be minor.",
188 "help: Clicking on help displays this message.",
190 "next: Displays the the next sequence.",
191 " Same as pressing 'n' or down-arrow.",
193 "prev: Displays the previous sequence.",
194 " Same as pressing 'p' or up-arrow.",
196 "auto: Automatically steps through the sequences",
199 "zoom: Gives more information on the highest",
200 " scoring residue. Same as pressing 'z'.",
201 " Clicking on a specific region will give",
202 " more information on that region.",
204 "print: Prints the sequence.",
206 "'table:': Switches the structure scores displayed.",
207 " and For MultiCoil: 1= dimeric, 2= trimeric",
208 "key 't': 3= total coiled-coil prob.",
209 "key 'i': Same as 't' for the second prediction method.",
211 "key 'm': Switches the scoring method.",
212 "key 'w': Switches scoring method for second predictor.",
214 "key 'Up' or 'Down': moves the bound up or down by .01",
215 "key 'Return': Rescores the sequence and coil scores after",
216 " changing the bound.",
218 "key 'o': Toggle between oligomerization ratio and",
219 " coiled-coil probability",
221 "key 'c': Switch between MultiCoil residue score,",
222 " average score over coiled coils, and coiled",
223 " coil score corresponding to the residue with",
224 " the maximum total coiled-coil prob over the",
225 " length of the coiled coil.",
227 "key '7': MultiCoil max score is obtained by taking tuple",
228 " of each dimension's max score (combo register).",
229 #ifdef DEBUG_PRINTOUT
230 "key 'x': MultiCoil computes scores for each offset, taking",
231 " max of those scores to be the max score (max reg).",
234 "quit: Quits the program and completes log file.",
235 " 'q' quits without completing the log.",
236 " Pressing 'q' in any window also closes ",
243 ShowTitle (Window win)
247 for (i = 0; motd[i][0]; ++i)
248 XDrawString(display, win, gc, 5, 17*i+20, motd[i], strlen(motd[i]));
249 XDrawString(display,win,gcgray,100,17*i+20,
250 "Press Any Key To Continue",25);
253 ShowDisclaimer(Window win)
257 XClearWindow(display,win);
258 for (i = 0; disclaimer[i][0]; ++i)
259 XDrawString(display, win, gc, 5, 17*i+20, disclaimer[i], strlen(disclaimer[i]));
260 XDrawString(display,win,gcgray,100,17*i+20,
261 "Press Any Key To Continue",25);
264 ShowHelp (Window win)
268 /* Places hmotd text in win */
269 for (i=0; hmotd[i][0]; ++i)
270 XDrawString(display,win,gc,5,17*i+20, hmotd[i],strlen(hmotd[i]));
275 static int seqnum=0, sequp=0, seqlow=1;
278 main (int argc, char *argv[])
281 struct itimerval interval;
285 char buff[MAXLINE+MAXSEQLEN];
289 char *command_line = NULL;
291 FILE *fgin=NULL, *fout=NULL, *fout_coils=NULL, *flog=NULL;
292 FILE *fpin[MAX_TABLE_NUMBER];
294 char pir_name[MAXLINE];
295 char likelihoods[MAX_TABLE_NUMBER][MAXLINE];
296 /* Files of likelihood lines. */
297 char printfile[500]; /* The filename to output postscript print to. */
298 char *print= printfile;
300 char class_sc_filenames[2][MAX_CLASS_NUMBER][MAXLINE];
301 double both_class_covars[2][MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX]
302 [NUM_DIM_IN_ORIG_MATRIX];
303 double both_class_means[2][MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX];
305 char gauss_param[2][MAXLINE];
306 /** 0 entry is for combo offset, 1 entry is for **/
307 /** max of other registers. ***/
308 extern int functnum; /* What is used in scscore.c. Takes value from */
309 /* pair_functnum and dimension stuff for multicoil. */
311 char lib[MAXFUNCTNUM]; /* Use MAXFUNCTUM=7 here since are 7 distances */
312 /* in coil, so we know AT MOST 7 lib functions. */
313 extern int pair_functnum; /* The number of distances in lib[].*/
315 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM];
316 extern int num_dist[MAX_TABLE_NUMBER];
320 int mode = 0; /* If the i-th option is set, mode gets 1 in i-th bit. */
323 double m[MAX_TABLE_NUMBER], m_single[MAX_TABLE_NUMBER],
324 b[MAX_TABLE_NUMBER], b_single[MAX_TABLE_NUMBER];
325 /*The likelihood line is mx+b for paircoil. */
329 int number_tables =0; /** Currently not used (replaced by NUMBER_TABLES */
331 int table_to_remove_from = -1;
332 /* Used if want to remove ALL input seq */
333 /* from that table at outset. **/
334 char *command_line_config = NULL;
336 /*** INITIALIZATION ***/
340 /*FILE *test = fopen ("/tmp/multicoil/test.txt", "wb");
341 fprintf(test, "test!\n");
345 printf("writing environment\n");
352 printf("%s\n",equality);
355 initialize(&num_dist[0], window_length, scale0s,
357 likelihoods, &pir_name[0], print, gauss_param,
358 &mode, init_class_prob);
361 printf("Printing arguments\n");
362 for (i=0;i<argc;i++){
363 printf("%s\n",argv[i]);
366 read_command_line(argc, argv, &command_line, combine_dist,
367 &mode, num_dist, multi_lib, &command_line_config);
369 fprintf(stderr,"\n\nRead commandline, command_line=%s",command_line);
371 /* Get the default options in the .paircoil file!!! */
372 get_defaults(command_line_config, argv[1], &mode, &log_bound, &fgin, fpin,
374 &flog, &fout_coils, &by_coil_or_seq, likelihoods,
375 pir_name, lib, multi_lib,
376 &pair_functnum, combine_dist,
377 print, &main_method, &main_preprocessor_method,
378 &main_table, &offset_to_use,
379 &avg_max, &Coil_Score,
382 gauss_param, num_dist,
384 &table_to_remove_from,
385 command_line, window_length, scale0s, scale0p);
387 if ((mode & RAW_OUT) && flog) {
388 fprintf(stderr,"\n Log file will not be written when doing raw scores.");
392 if ((mode & RAW_OUT) && fout_coils) {
394 "\n Sequence scores will not be written when doing raw scores.");
401 fprintf(stderr,"\nTable %d, num_dist = %d, combine_dist =%d", i,
402 num_dist[i], combine_dist[i]);
404 fprintf(stderr,"\nOut of get_def");
410 if (ftotal_like) { /** Want to output total gauss like over classes*/
411 for (i=0; i<NUMBER_CLASSES; i++) /* for various init_prob settings. */
412 init_class_prob[i]=1; /** All 1 values is flag to just output */
413 init_totals(total_gauss_like); /** gauss values, not like values in */
414 } /** convert_raw_scores_to_gauss_like. */
416 /** 7*(1-combine_dist) + 1*combine_dist gives the number of **/
417 /** dimensions started with, and **/
418 /** num_dist*(1-combine_dist) + 1*combine_dist gives the new number. **/
419 number_score_dim = num_table_dimensions(num_dist[0],combine_dist[0]) +
420 num_table_dimensions(num_dist[1],combine_dist[1]);
423 if (gauss_param[i][0] != ',') {
424 get_gauss_param_all_classes2(class_sc_filenames[i],
425 num_table_dimensions(7,combine_dist[0]) +
426 num_table_dimensions(7,combine_dist[1]),
427 both_class_means[i], both_class_covars[i],
430 /******Note we store the submatrices in the original matrices....****/
431 compute_submatrices2(both_class_means[i],both_means_submatrix[i],
432 both_class_covars[i],
433 both_inverse_covars[i],
434 both_determinants[i], multi_lib,
435 num_table_dimensions(7,combine_dist[0]),
436 num_table_dimensions(7,combine_dist[1]),
437 num_table_dimensions(num_dist[0],combine_dist[0]),
438 num_table_dimensions(num_dist[1],combine_dist[1]));
439 #ifdef DEBUG_PRINTOUT
440 for (class=0; class<NUMBER_CLASSES; class++) {
441 fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
442 print_matrix(both_means_submatrix[i][class],1,
443 number_score_dim,stderr);
444 fprintf(stderr,"\n\nInver Covar Submatrix %d:\n", class);
445 print_matrix(both_inverse_covars[i][class],
446 number_score_dim, number_score_dim,stderr);
450 fprintf(stderr,"\nDet0 =%lf, det1= %lf, det2= %lf\n",
451 both_determinants[i][0],both_determinants[i][1],
452 both_determinants[i][2]);
454 fprintf(stderr,"\nnew_num_dim_tab0=%d, new_num_dim_tab1=%d, old_num_dim_tab0=%d, old_num_dim_tab1=%d, number_score_dim=%d",
455 num_table_dimensions(num_dist[0],combine_dist[0]),
456 num_table_dimensions(num_dist[1],combine_dist[1]),
457 num_table_dimensions(7,combine_dist[0]),
458 num_table_dimensions(7,combine_dist[1]),
465 switch_gauss_param(offset_to_use, main_method);
470 log_offset_to_use = offset_to_use;
471 method = main_method;
472 preprocessor_method = main_preprocessor_method;
476 if ( (mode & PRN_MODE) && flog ){
477 log_file_header(flog, mode,
478 argc, argv, avg_max, log_bound,
479 preprocessor_method, preproc_table, log_offset_to_use,
480 start_coil_at_reg_shift, main_method,
481 main_table, lib, multi_lib);
484 if (mode & POS_MODE) ActualCoil_Method = available;
486 /****** Now if didn't get locations from .paircoil, get from ENV variable */
487 if (!fgin) { /* Didn't get fgin in */
488 fprintf(stderr,"\n\n **********************");
489 fprintf(stderr,"\nDid not get genbnk location from config file.");
490 fprintf(stderr,"\n Will instead use uniform distrib. as underlying distrib.");
491 fprintf(stderr,"\n**********************\n\n");
493 /**** exit(-1); ****/ /* get_defaults. */
495 if (!fpin[0]) { /* Didn't get fpin in */
496 fprintf(stderr,"\nDid not get positive data set from config file.");
497 exit(-1); /* get_defaults. */
499 /******************************************/
501 pfreqp = malloc(sizeof(long) * MAX_TABLE_NUMBER*AANUM*AANUM*POSNUM*POSNUM);
503 PairCoilData(fgin,fpin, mode & PROLINE_FREE_MODE,
504 table_to_remove_from, argv[1], mode, pfreqp);
505 /* If fpin[1] != NULL this will call */
506 /* TrimerDimerDiff to set up table diff */
508 /* If it is NULL then table is average */
509 /* of probabilities from 1 and 2. */
513 hmm = CCHMM(fgin,fpin);
517 if (mode & USE_LIKE_LINE) {
518 /** Don't want to compute likelihood line if */
520 functnum = pair_functnum;
521 get_likelihood_line(likelihoods, m, b, m_single, b_single,
522 lib, pair_functnum,pir_name,mode,NUMBER_TABLES);
525 #ifdef DEBUG_PRINTOUT
526 for (i=0; i<NUMBER_TABLES; i++)
527 fprintf(stderr, "\n m[%d]= %lf, b[%d]= %lf\n",i,m[i],i,b[i]);
532 while (i< NUMBER_TABLES) {
536 if (fgin) sclose(fgin);
542 fin = sopen(argc>1?argv[1]:"-","r");
544 fprintf(stderr, "\nfin is opened to %s",argc>1?argv[1]:"-");
547 /*** Note: If can't open xdisplay then OpenScreen will write log and */
548 /*** out files if those options are set. */
550 OpenScreen(mode,log_bound,fin,flog, fout_coils,
551 fout, m,b,m_single, b_single,
553 main_method,main_preprocessor_method, &seqnum, pfreqp);
556 /* Let user read message while initializing */
562 interval.it_value.tv_sec = 2;
563 interval.it_value.tv_usec = 0;
564 interval.it_interval.tv_sec = 2;
565 interval.it_interval.tv_usec = 0;
568 XSelectInput(display, window, KeyPressMask | ButtonPressMask | ExposureMask);
570 if (autoadvance) do {
571 FD_SET(ConnectionNumber(display),&fds);
572 if (!select(FD_SETSIZE,&fds,&fdnone,&fdnone,&interval))
573 NextSeq(fin,lib,multi_lib,m,b,m_single,b_single,
574 mode,bound,flog,fout_coils,fout);
576 while ( (!XPending(display)) && (autoadvance));
577 /* Autoadvance until button is hit or end of seq file. */
578 XNextEvent(display,&report);
579 switch (report.type) {
581 if (report.xany.window == zoom)
583 else if (report.xany.window == help)
585 else if (sequence.seqlen)
586 ShowScores(scores,&sequence,NULL,bound, mode);
589 #ifdef SHOW_DISCLAIMER
590 XNextEvent(display, &report);
591 ShowDisclaimer(window);
596 /* if (!sequence.seqlen) {NextSeq(fin,lib,multi_lib,m,b,m_single,
597 b_single,bound,flog,fout_coils,fout);
600 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
601 if (report.xany.window == zoom) {XUnmapWindow (display, zoom); break;}
602 if (report.xany.window == window) {
603 /* Zoom in and show register */
605 tmp = XYtoPos(report.xbutton.x,report.xbutton.y);
607 ZoomPnt(); zoomptr=tmp; Zoomer(); ZoomPnt();
610 if (report.xany.window == helpb)
611 XMapWindow(display,help);
612 if (report.xany.window == next)
613 NextSeq(fin,lib,multi_lib,m,b,m_single, b_single,
614 mode,bound,flog,fout_coils,fout);
615 if (report.xany.window == prev)
616 PrevSeq(fin,lib,multi_lib,m,b,m_single,b_single,mode,bound, pfreqp);
617 if (report.xany.window == autob)
618 { if (autoadvance ^= 1) NextSeq(fin,lib,multi_lib,m,b,m_single,
619 b_single,mode,bound,flog,fout_coils,
623 if (report.xany.window == zoomb) {
624 if (zoomptr) ZoomPnt();
625 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
627 if (report.xany.window == printb) {
629 ps = sopen(print?print:defprint,"a");
630 PSDefs(ps, sequence.seqlen);
631 ShowScores(scores,&sequence,ps,bound,mode);
634 if (report.xany.window == tableb) {
635 if (method == MultiCoil)
636 table = (table +1) % NUMBER_CLASSES;
637 else if (method == HMMCoil)
638 table = (table +1) % NUMBER_TABLES;
640 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,0, bound);
641 if (zoomptr) Zoomer();
645 if (report.xany.window == methodb) {
646 method_cycle(&method,NUMBER_METHODS);
647 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,0, bound);
648 if (zoomptr) Zoomer();
651 if (report.xany.window == quit) {
652 if (flog || fout || ftotal_like || fout_coils)
655 finish_log(mode,log_bound,fin,flog,fout_coils,
656 fout,m,b,m_single,b_single,
657 lib,multi_lib,main_method,main_table,
658 main_preprocessor_method,
663 if (flog) sclose(flog); if (fout) sclose(fout);
664 if (fout_coils) sclose(fout_coils);
665 if (ftotal_like) sclose(ftotal_like);
672 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
673 switch (XLookupKeysym(&(report.xkey),0)) {
676 if (bound > 1) bound =1;
677 draw_bound_box(bound, paramb);
681 if (bound < 0) bound =0;
682 draw_bound_box(bound, paramb);
685 case XK_Return: /** Signal that bound has been set, should rescore **/
686 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
687 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,-1,-1, bound);
688 if (zoomptr) Zoomer();
692 case XK_C: /** Switch between register, avg, and max **/
693 case XK_c: /** coil score for MultiCoil. **/
694 Coil_Score = (Coil_Score +1) % 3;
696 ShowSeq(lib,multi_lib, /** Change to res score **/
697 m,b,m_single,b_single,mode,-2,-2, bound);
698 else ShowSeq(lib,multi_lib, /** Change to coil score with */
699 m,b,m_single,b_single,mode,-1,-1, bound); /* out recomputing */
700 if (zoomptr) Zoomer(); /* res score. */
703 case XK_O: /** Oligomerization ratio. ***/
705 mode ^= ONLY_COILED_CLASSES; /* To XOR to complement that bit. */
706 if (NUMBER_CLASSES >0) {
707 if (method == MultiCoil)
708 type_class_score_convert(sequence,all_scores,bound,
709 mode & ONLY_COILED_CLASSES);
710 if (preprocessor_method == MultiCoil)
711 type_class_score_convert(sequence,all_preprocess_like,bound,
712 mode & ONLY_COILED_CLASSES);
714 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,2, bound);
715 if (zoomptr) Zoomer();
720 ps = sopen(print?print:defprint,"w");
721 PSDefs(ps, sequence.seqlen);
722 ShowScores(scores,&sequence,ps,bound,mode);
728 {ZoomPnt(); --zoomptr; Zoomer(); ZoomPnt();}
731 if (zoomptr && zoomptr<sequence.seqlen)
732 {ZoomPnt(); ++zoomptr; Zoomer(); ZoomPnt();}
736 if (zoomptr) ZoomPnt();
737 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
740 #ifdef DEBUG_PRINTOUT
742 offset_to_use =-1; /* Do max scoring register */
743 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,1, bound);
744 if (zoomptr) Zoomer();
748 offset_to_use =7; /* Do MultiCoil max score as combo of max score */
749 ShowSeq(lib,multi_lib, /* in each dimension. */
750 m,b,m_single,b_single,mode,1,1, bound);
751 if (zoomptr) Zoomer();
758 XMapWindow(display,help);
763 if (report.xany.window == window)
764 {sclose(fin); if (flog) sclose(flog);
765 if (ftotal_like) sclose(ftotal_like);
769 if (fout) sclose(fout);
772 if (report.xany.window == zoom) {
774 XUnmapWindow(display, zoom);
780 method_cycle(&method,NUMBER_METHODS);
782 if (method == MultiCoil)
783 table = table % NUMBER_CLASSES;
784 else if (method == HMMCoil)
785 table = table % NUMBER_TABLES;
788 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,0, bound);
789 if (zoomptr) Zoomer();
794 method_cycle(&preprocessor_method,NUMBER_PREPROCESSORS);
796 if (preprocessor_method == MultiCoil)
797 preproc_table = preproc_table% NUMBER_CLASSES;
798 else if (preprocessor_method == HMMCoil)
799 preproc_table = preproc_table%NUMBER_TABLES;
800 else preproc_table =0;
802 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,1, bound);
803 if (zoomptr) Zoomer();
807 if (preprocessor_method == MultiCoil)
808 preproc_table = (preproc_table +1) % NUMBER_CLASSES;
809 else if (preprocessor_method == HMMCoil)
810 preproc_table = (preproc_table +1) % NUMBER_TABLES;
811 else preproc_table = 0;
812 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,2, bound);
813 if (zoomptr) Zoomer();
817 if (method == MultiCoil)
818 table = (table +1) % NUMBER_CLASSES;
819 else if (method == HMMCoil)
820 table = (table +1) % NUMBER_TABLES;
822 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,0, bound);
823 if (zoomptr) Zoomer();
828 if (report.xany.window == window) PrevSeq(fin, lib,multi_lib,m,b,
829 m_single,b_single,mode,bound,pfreqp);
837 if (report.xany.window == window)
838 NextSeq(fin,lib,multi_lib,m,b,m_single,b_single,mode,bound,flog,
851 if (XCheckTypedEvent(display,Expose,&report)) {
853 XDrawString(display,window,gcgray,500,17,"(initializing)",14);
859 struct itimerval value;
860 signal(SIGALRM,Alarm);
861 getitimer(ITIMER_REAL,&value);
862 value.it_value.tv_sec = 0;
863 value.it_value.tv_usec = 500;
864 value.it_interval.tv_sec = 0;
865 value.it_interval.tv_usec = 500;
866 setitimer(ITIMER_REAL,&value,NULL);
871 struct itimerval value;
872 getitimer(ITIMER_REAL,&value);
873 value.it_value.tv_sec = 0;
874 value.it_value.tv_usec = 0;
875 setitimer(ITIMER_REAL,&value,NULL);
876 XDrawImageString(display,window,gcgray,500,17," ",14);
877 XDrawImageString(display,window,gcgray,225,729,"press any key to continue",25);
883 if (!zoomptr) return;
884 PostoXY(zoomptr,&x,&y);
885 pts[0].x = x; pts[0].y = y;
886 pts[1].x = -5; pts[1].y = 5;
887 pts[2].x = 10; pts[2].y = 0;
888 XFillPolygon(display,window,gcdash,pts,3,Convex,CoordModePrevious);
889 XDrawLine(display,window,gcdash,x,y+5,x,y+8);
893 /* Find most likely residue, break ties arbitrarily */
898 for (i=0; i<sequence.seqlen; ++i)
899 if (scores[8*i+7] > like) {
900 like = scores[8*i+7];
911 XMapWindow(display, zoom);
912 /* Indicate where in the sequence we are */
913 sprintf(buf,"Position %d ",zoomptr);
914 XDrawImageString(display,zoom,gc,7,20,
916 for (i = -3; i<=3; ++i) {
917 if (zoomptr+i>0 && zoomptr+i<=sequence.seqlen)
918 buf[0] = numaa(sequence.seq[zoomptr-1+i]);
921 XDrawImageString(display,zoom,gc,159+9*i,20, buf,1);
923 XDrawLine(display,zoom,gc,159,22,168,22);
925 /* Indicate the probability of being in any given state */
926 sc = &(scores[8*(zoomptr-1)]);
929 if (sc[i] == sc[7]) {offset=i; break;} /* Used to print out maxscoring */
930 /* register as 'a' + offset. */
931 /* Offset 7 will be a flag that no particular offset scored max. */
933 /* Give likelihood */
934 sprintf(buf,"%5.1f%% chance",100*sc[7]);
935 XDrawImageString(display,zoom,gc,220,20, buf,strlen(buf));
938 fprintf(stderr,"\noffset= %d, reg = %c",all_offsets[table][zoomptr-1],
939 'a' + (zoomptr-1 + all_offsets[table][zoomptr-1])%POSNUM);
944 sprintf(buf,"register %c",'a'+offset);
946 sprintf(buf,"combo reg ");
947 else sprintf(buf," ");
948 XDrawImageString(display,zoom,gc,230,40, buf,strlen(buf));
950 for (i=0; i<7; ++i) {
951 sprintf(buf,"%c:%4.1f%% ",'a'+i,(100*sc[i]));
952 XDrawImageString(display,zoom,gc,15+92*(i&3),60+20*(i/4),buf,strlen(buf));
954 XDrawString(display,zoom,gc,10,110,"Click on viewgram or use arrow",30);
955 XDrawString(display,zoom,gc,10,127,"keys to view other positions.",29);
956 XDrawString(display,zoom,gc,10,144,"You may need to move this",25);
957 XDrawString(display,zoom,gc,10,161,"window to see the viewgram.",27);
967 /*********The following is for storing sequences. ***/
968 static char seqbuff[40000];
970 static Sequence seqs[30];
971 static int seqbufpos[30];
972 static int seqsizes[30];
975 { int size, pos, limsize;
977 size = seqsizes[sequp-seqnum];
978 pos = seqbufpos[sequp-seqnum];
980 sequence = seqs[sequp-seqnum]; /** Retrieve seq from memory. **/
983 limsize = (pos+size<=40000) ? size : 40000-pos;
984 bcopy(&(seqbuff[pos]),sequence.code,limsize);
985 bcopy(seqbuff,&(sequence.code[limsize]),size-limsize);
991 { int size, pos, limsize,i;
992 char *temp, *current;
994 bcopy(seqs,&(seqs[1]),29*sizeof(Sequence));
995 bcopy(seqbufpos,&(seqbufpos[1]),29*sizeof(int));
996 bcopy(seqsizes,&(seqsizes[1]),29*sizeof(int));
997 if (sequp-seqlow>=30) seqlow = sequp-29;
999 /* The 2 accounts for the end of string characters in .code and .title. */
1000 size = 2 + strlen(sequence.code) + strlen(sequence.title) +
1001 sequence.seqlen * (sequence.reg ? 2 : 1);
1003 pos = seqbufpos[1] + seqsizes[1]; /* Number of characters in array */
1005 /* if (pos>=40000) pos -= 40000;*/
1006 /* limsize = (pos+size<=40000) ? size : 40000-pos; */
1008 if (pos+size >=40000) pos=0; /* Array will be full so start at beginning. */
1009 limsize = (pos+size<=40000) ? size : 40000-pos;
1012 /* Save the "code+title+seq+reg" in something permanant, and wrap around
1015 /* David uses bcopy because wants to copy past end of string characters */
1017 bcopy(sequence.code,&(seqbuff[pos]),limsize);
1018 bcopy(&(sequence.code[limsize]),seqbuff,size-limsize);
1021 sequence.code=strcpy(&(seqbuff[pos]),sequence.code);
1023 strcpy(&(seqbuff[pos+strlen(sequence.code)+1]),sequence.title);
1027 &(seqbuff[pos+strlen(sequence.code) + strlen(sequence.title) + 2]);
1028 for(i=0; i<sequence.seqlen; i++)
1029 sequence.seq[i] = temp[i];
1034 sequence.reg=&(sequence.seq[sequence.seqlen]);
1035 for(i=0; i<sequence.seqlen; i++)
1036 sequence.reg[i] = temp[i];
1040 while ((seqlow<sequp) && (seqbufpos[sequp-seqlow] >= pos && seqbufpos[sequp-seqlow] < pos+limsize || seqbufpos[sequp-seqlow] < size-limsize))
1043 seqs[0].seqlen= sequence.seqlen;
1053 method_cycle(int *Method, int Number_of_Methods)
1056 *Method = (*Method+1)%Number_of_Methods;
1058 while ((*Method == MultiCoil && MultiCoil_Method == unavailable) ||
1059 (*Method == PairCoil && PairCoil_Method == unavailable) ||
1060 (*Method == NEWCOILS && NEWCOILS_Method == unavailable) ||
1061 (*Method == HMMCoil && HMMCoil_Method == unavailable) ||
1062 (*Method == ActualCoil && ActualCoil_Method == unavailable) ||
1063 (*Method == PairCoilDiff && PairCoilDiff_Method == unavailable)||
1064 (*Method == PairCoilDiffAvg && PairCoilDiffAvg_Method == unavailable) ||
1066 (*Method == ActualCoil) && (!sequence.reg) /* Don't do actual coil */
1067 ); /* if didn't input register */
1073 PrevSeq (FILE *fin, char lib[MAXFUNCTNUM],
1074 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1075 double *m, double *b,
1076 double *m_single, double *b_single,
1077 int mode, double bound, long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
1079 if (seqnum<=seqlow) {XBell(display,0); return;}
1082 NewSeq(mode,sequence, pfreqp);
1083 ShowSeq(lib, multi_lib,m, b, m_single, b_single, mode,1,1, bound);
1087 NextSeq (FILE *fin,char lib[MAXFUNCTNUM],
1088 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1089 double *m, double *b, double *m_single, double *b_single,
1090 int mode, double bound, FILE *flog, FILE *fout_coils,FILE *fout)
1092 char title[MAXLINE],code[MAXLINE];
1093 char seq[MAXSEQLEN], reg[MAXSEQLEN];
1095 int first_time_seq=0;
1098 if (seqnum < sequp) {
1104 sequence.seqlen = 0;
1107 sequence.title = title;
1108 sequence.code = code;
1111 if (!getseq2(fin,&sequence))
1112 {XBell(display,0); sequence= temp_seq; autoadvance=0; return;}
1120 NewSeq(mode,sequence);
1121 if (first_time_seq) {
1122 output_seq(lib,multi_lib,m,b,m_single, b_single,
1123 mode,log_bound,flog,fout_coils,fout,
1124 avg_max, main_method, main_preprocessor_method, main_table);
1130 ShowSeq(lib,multi_lib,m, b, m_single, b_single,mode,1,1, bound);
1135 /** This is stuff need to always do before rescore a sequence. **/
1136 NewSeq_nonX (int mode, Sequence sequence, long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
1138 if (mode & TST_MODE0)
1139 recalc_prob(sequence, 0, mode, pfreqp);
1140 if (mode & TST_MODE1)
1141 recalc_prob(sequence,1, mode, pfreqp);
1144 NewSeq (int mode, Sequence sequence, long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
1146 NewSeq_nonX (mode, sequence, pfreqp);
1148 XUnmapWindow(display, zoom);
1153 /* The likelihood line is mx +b */
1154 void output_seq(char lib[MAXFUNCTNUM],
1155 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1156 double *m,double *b,
1157 double *m_single, double *b_single,
1159 double log_bound,FILE *flog,
1163 int main_method, int main_preprocessor_method,
1167 double maxother=0; double maxscore;
1168 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1173 double junk_scores[MAXSEQLEN][POSNUM+1];
1176 fprintf(stderr,"\nIn outputseq");
1179 if (is_not_differentiator(main_method))
1183 all_scores = ScoreSeq (lib, multi_lib,m, b, m_single, b_single, mode,
1184 main_table, main_method, main_preprocessor_method,
1185 log_offset_to_use, &maxscore, 1, do_preproc, log_bound);
1189 add_to_total_likes(all_scores, total_gauss_like, sequence.seqlen,
1193 fprintf(stderr,"\nGot all_scores");
1195 scores = &all_scores[main_table* MAXSEQLEN*(POSNUM+1) ];
1197 /* This is all_scores[main_table + 0*number_tables][0][0] */
1198 /* so is score on main_table using library # 0 **/
1199 /* Similarly for preproces_like. **/
1201 preprocess_like = &all_preprocess_like[main_table][0][0];
1204 for (i=0; i<sequence.seqlen; i++)
1205 if (preprocess_like[i*(POSNUM+1) + 7] > maxother)
1206 maxother=preprocess_like[i*(POSNUM +1) + 7];
1207 /* This is preprocess_like[i][7] */
1210 switch (main_method) {
1212 if (mode & VER_MODE) {
1213 NEWCOILSScore(mode, sequence, &maxother,junk_scores,
1216 /* scores in junk_scores */
1217 log_output_ver(SC2SEQ | SCSTOCK, maxscore, maxother,
1218 maxscore,seqnum, sequence.title,
1219 sequence.code, flog);
1226 /*** For pos files coils, does a log file version of coil/seq score. **/
1227 output_pos_scores(fout_coils, sequence, all_scores, log_bound, mode,2,1);
1230 if (mode & PRN_MODE)
1233 /* If aren't dealing with coil scores, convert to weighted-avg **/
1234 /* coil_scores (by using average_score_over_coils3(....,0) **/
1235 /* so can output to log file. **/
1237 for (tablenum =0; tablenum <NUMBER_CLASSES; tablenum++) {
1238 average_score_over_coils3(sequence,
1239 &all_scores[(tablenum+NUMBER_CLASSES)* MAXSEQLEN*(POSNUM+1)],
1240 &all_scores[tablenum* MAXSEQLEN*(POSNUM+1)],
1241 &all_scores[(2+NUMBER_CLASSES)* MAXSEQLEN*(POSNUM+1)],
1242 /** Total Coil Like. **/
1244 all_offsets[tablenum],
1245 start_coil_at_reg_shift, log_bound,
1249 log_coil_output(all_scores, all_offsets,sequence, log_bound, flog,
1252 /** Now convert back for safety in case need to do an out file **/
1254 switch_back_to_res_score(log_bound, sequence,mode,all_scores);
1258 log_output2_prn(mode,maxscore, maxscore, scores, scores,
1260 seqnum, sequence, log_bound, log_bound, log_bound,
1268 if (mode & VER_MODE)
1269 log_output_ver(SC2SEQ , maxscore, 0,
1270 maxscore,seqnum, sequence.title,
1271 sequence.code, flog);
1274 if (mode & PRN_MODE)
1275 log_output2_prn(mode,maxscore, maxscore, scores, scores,sequence.seqlen,
1276 seqnum, sequence, log_bound, log_bound, log_bound,
1281 /* case PairCoilDiff: */
1282 case PairCoilDiffAvg:
1283 if (mode & PRN_MODE) {
1284 log_output2_prn(mode,maxscore, maxother, scores,preprocess_like,
1286 seqnum, sequence, .5, .5, log_bound,
1296 fprintf(stderr,"\n About to go into fout.");
1301 if ( (!(mode & RAW_OUT)) ||
1302 ( (main_method != MultiCoil) && (main_method != PairCoilDiff))) {
1303 if (main_method== MultiCoil)
1304 if (mode & WEB_OUT_MODE)
1305 web_output(sequence, fout, all_scores,all_offsets,2);
1306 /* Output table of score and seq and reg values */
1308 tuple_output2(sequence, mode, fout, all_scores,log_bound,2);
1309 /* Output dimer and trimer prob. in regions above bound and */
1310 /* in long enough windows only. **/
1312 else txt_output2(sequence, mode, fout, scores, log_bound);
1314 else { /* Want raw tuple_output for multi_coil and PairCoilDiff*/
1316 if (main_method == MultiCoil)
1317 if (log_offset_to_use==7) /* Combo offset just output max score */
1318 tuple_output2(sequence, mode, fout, /** for each dimension. */
1320 log_bound, number_score_dim);
1321 else /** Max reg for raw score can't really be determined. Instead */
1322 /* for coils, output scores for the known correct register, & */
1323 /* for pdb- output ALL possible 7 register scores (all negs). */
1324 tuple_output_for_max(sequence,mode,fout, all_scores, NUMBER_TABLES,
1325 log_bound, number_score_dim);
1332 int switch_back_to_res_score(double bound, Sequence sequence, int mode,
1333 double All_Scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1])
1338 for (tablenum =0; tablenum < NUMBER_CLASSES ; tablenum++)
1339 for (j =0; j <=POSNUM; j++)
1340 for (i=0; i<sequence.seqlen; i++) {
1341 All_Scores[tablenum][i][j] =
1342 All_Scores[tablenum+ NUMBER_CLASSES][i][j];
1343 if ((mode & ONLY_COILED_CLASSES) && (tablenum < (NUMBER_CLASSES-1))
1344 && (All_Scores[2+ NUMBER_CLASSES][i][j]>bound))
1345 All_Scores[tablenum][i][j]/= All_Scores[2+ NUMBER_CLASSES][i][j];
1350 int copy_offsets(char preproc_offsets[MAXSEQLEN],
1352 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1354 double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1359 for (tablenum =0; tablenum < number_dimens; tablenum++)
1360 for(i=0; i<sequence.seqlen; i++) {
1361 all_offsets[tablenum][i]=preproc_offsets[i];
1362 /*** Make sure the output score corresponds to this offset... **/
1363 /*** Also do it for backup residue scores in tablenum + NUMBER_CLASSES */
1364 if ((offset_to_use ==7) || (offset_to_use == -1)) {
1365 all_scores[tablenum][i][7] =
1366 all_scores[tablenum][i][(i+all_offsets[tablenum][i]) %POSNUM];
1367 all_scores[tablenum+ NUMBER_CLASSES][i][7] =
1368 all_scores[tablenum+NUMBER_CLASSES][i]
1369 [(i+all_offsets[tablenum][i]) %POSNUM]; }
1373 int get_offsets(double all_likes[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1374 Sequence sequence, char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1378 int got_offset, max_offset;
1379 double max_off_score;
1382 for (tablenum =0; tablenum < number_dimens; tablenum++)
1383 for(i=0; i<sequence.seqlen; i++) {
1384 offsets[tablenum][i]=-1;
1386 max_off_score=-HUGE_VAL;
1387 for (j=0; j<POSNUM; j++) {
1388 if (all_likes[tablenum][i][j] == all_likes[tablenum][i][7]) {
1389 offsets[tablenum][i] = (j-i) % POSNUM;
1390 if (offsets[tablenum][i]<0)
1391 offsets[tablenum][i] += POSNUM;
1394 if (all_likes[tablenum][i][j] > max_off_score) {
1395 max_off_score=all_likes[tablenum][i][j];
1396 max_offset=(j-i) % POSNUM;
1397 if (max_offset <0) max_offset +=POSNUM;
1400 if (!got_offset) offsets[tablenum][i] = max_offset;
1402 /* This is a new thing so that if an offset change is not real */
1403 /* we keep the old offset. **********/
1404 if ( (i>0) && (all_likes[tablenum][i][(i+offsets[tablenum][i-1]) %POSNUM]
1405 == all_likes[tablenum][i][(i+offsets[tablenum][i])%POSNUM]) )
1406 offsets[tablenum][i] = offsets[tablenum][i-1];
1418 void preprocessor_score(char lib[MAXFUNCTNUM],
1419 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1420 double *m, double *b,
1421 double *m_single, double *b_single, int mode,
1422 int Preprocessor_Method, int Offset_to_Use,
1423 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1425 char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1429 char offsets[MAXSEQLEN];
1433 int dimension, num_score_dim;
1436 for (tablenum =0; tablenum < NUMBER_TABLES; tablenum++) {
1437 switch_tables(tablenum); /* Set prob. to be for that table. */
1439 switch(Preprocessor_Method) {
1441 if (combine_dist[tablenum])
1442 functnum = num_dist[tablenum]; /* For use in scscore.c **/
1445 num_score_dim= num_table_dimensions(num_dist[tablenum],
1446 combine_dist[tablenum]);
1448 for (dist=0; dist<num_score_dim; dist++) {
1449 if (num_dist[tablenum]>1) { /* Hack: else don't access correctly */
1450 local_lib = &multi_lib[tablenum][dist];
1452 else local_lib = multi_lib[tablenum];
1453 dimension = compute_dimension2(tablenum, dist, num_dist,
1455 MultiCoilDimensionScore(mode, sequence, local_lib,
1456 &maxscore,tablenum, offsets,
1457 all_preprocess_like[dimension],
1463 if (mode & PAIRCOIL_PAIRS) {
1464 functnum = pair_functnum; /* For use in scscore.c **/
1465 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1466 &maxscore,tablenum, offsets,
1467 all_preprocess_like[tablenum],
1470 else /* Use singles probabilities... for likelihood use line 0 */
1471 SingleCoilScore(mode,sequence,m_single[tablenum],
1473 &maxscore, tablenum,
1474 all_preprocess_like[tablenum], Offset_to_Use,
1479 NEWCOILSScore(mode, sequence, &maxscore, all_preprocess_like[tablenum], Offset_to_Use);
1483 HMMScore(sequence.seq, sequence.seqlen, hmm, all_preprocess_like[tablenum],
1484 tablenum, Offset_to_Use);
1488 /* if (mode & POS_MODE) */ /* Score the coils in the posfile based on */
1489 /* their having 100% likelihood. */
1490 ActualCoils(sequence, all_preprocess_like[tablenum], Offset_to_Use,1);
1494 fprintf(stderr,"\nBad preprocessor, using PairCoil.\n");
1495 functnum = pair_functnum; /* For use in scscore.c **/
1496 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1497 &maxscore,tablenum, offsets,
1498 all_preprocess_like[tablenum],
1505 if (Preprocessor_Method !=MultiCoil)
1506 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1509 /***** Get copy of res scores in tablenum and tablenum+ NUMBER_CLASSES **/
1510 switch_gauss_param(Offset_to_Use, Preprocessor_Method);
1512 convert_raw_scores_to_gauss_prob_like2(all_preprocess_like,sequence.seqlen,
1513 NUMBER_TABLES,means_submatrix,inverse_covars,
1514 determinants, NUMBER_CLASSES, init_class_prob,
1516 mode & ONLY_COILED_CLASSES,
1517 bound,Offset_to_Use,1);
1519 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1522 if (Coil_Score) { /* Value of 1 means do average, 2 means do max. */
1523 /*** Need to save residue scores, since the coil scores **/
1524 /*** will replace them, but old residue scores might be needed **/
1525 /*** again in tablenum+ NUMBER_CLASSES. ****/
1526 for (tablenum =0; tablenum <3; tablenum++) {
1527 average_score_over_coils3(sequence,
1528 all_preprocess_like[tablenum+NUMBER_CLASSES],
1529 all_preprocess_like[tablenum],
1530 all_preprocess_like[2+NUMBER_CLASSES],
1531 /** Total Coil Like. **/
1533 all_preproc_offsets[tablenum],
1534 start_coil_at_reg_shift, bound,
1548 double *ScoreSeq (char lib[MAXFUNCTNUM],
1549 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1550 double *m_paircoil, double *b_paircoil,
1551 double *m_single, double *b_single,
1552 int mode, int Table, int Method,
1553 int Preprocessor_Method, int Offset_to_Use,
1555 int rescore_seq, int rescore_preproc, double bound)
1557 /** rescore_preproc == 1 means need to rescore. If it or rescore_seq are */
1558 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
1559 /** -2 means need to change from MultiCoil coil score back to res score. */
1560 /** If it is another non-zero */
1561 /** value, don't need to rescore it, but do need to rescore anything that */
1562 /** uses it as a preproc. */
1564 double max_paircoil_score;
1565 static char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1566 extern char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1567 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1568 static double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1575 int num_score_dim, dimension;
1576 int local_number_tables;
1578 double *preproc_res_like;
1579 int scored_already =0;
1581 /**************************************************************************/
1582 /********* PreprocLike is used in the 2-3 stranded differentiators. ****/
1585 if ( (rescore_preproc==1) &&
1586 (!ftotal_like)) { /* A hack to save computation time. */
1587 preprocessor_score(lib, multi_lib,
1588 m_paircoil, b_paircoil, m_single, b_single,
1589 mode, Preprocessor_Method,
1591 all_preprocess_like, sequence, all_preproc_offsets,
1594 /*****************************************************************************/
1595 if (Preprocessor_Method == MultiCoil) {
1596 if ( ( rescore_preproc==-1) && Coil_Score) {
1597 /* Coil_Score value of 1 means do average, 2 means do max. */
1598 /*** Need to save residue scores, since the coil scores **/
1599 /*** will replace them, but old residue scores might be needed **/
1600 /*** again in tablenum+ NUMBER_CLASSES. ****/
1601 for (tablenum =0; tablenum < NUMBER_CLASSES; tablenum++)
1602 average_score_over_coils3(sequence,
1603 all_preprocess_like[tablenum+ NUMBER_CLASSES],
1604 all_preprocess_like[tablenum],
1605 all_preprocess_like[2+NUMBER_CLASSES],
1606 /** Total Coil Like. **/
1608 all_preproc_offsets[tablenum],
1609 start_coil_at_reg_shift, bound,
1611 if (mode & ONLY_COILED_CLASSES) /** Convert to olig state. **/
1612 type_class_score_convert(sequence,all_preprocess_like,bound,
1613 mode & ONLY_COILED_CLASSES);
1615 else if (rescore_preproc == -2) /** Switch back to res score.*/
1616 switch_back_to_res_score(bound, sequence,mode,all_preprocess_like);
1619 /****************************************************************************/
1625 if (Method == PairCoilDiff)
1626 local_number_tables = 1;
1627 else local_number_tables= NUMBER_TABLES;
1630 /*** Clean up this conditional, since really the only thing that might ***/
1631 /*** change for paircoildiffer when change preproc is the offset to use **/
1632 /*** (i.e. all_scores[table][i][7]. ***/
1634 if ((rescore_seq>0) ||
1635 (!is_not_differentiator(Method) && rescore_preproc) ) {
1637 for (tablenum=0; tablenum<local_number_tables; tablenum++) {
1638 switch_tables(tablenum); /* Switch prob to right table. */
1641 maxsc=0; /* Initialization. */
1647 functnum = pair_functnum; /* For use in scscore.c **/
1648 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1649 b_paircoil[tablenum],
1650 &maxsc,tablenum, all_offsets[tablenum],all_scores[tablenum],
1656 if (combine_dist[tablenum])
1657 functnum = num_dist[tablenum]; /* For use in scscore.c **/
1660 num_score_dim= num_table_dimensions(num_dist[tablenum],
1661 combine_dist[tablenum]);
1663 for (dist=0; dist<num_score_dim; dist++) {
1664 if (num_dist[tablenum] >1)
1665 local_lib = &multi_lib[tablenum][dist];
1666 else local_lib = multi_lib[tablenum];
1668 dimension = compute_dimension2(tablenum, dist,
1669 num_dist, combine_dist);
1670 MultiCoilDimensionScore(mode, sequence, local_lib, &maxsc,
1671 tablenum,all_offsets[dimension],
1672 all_scores[dimension],
1679 if (mode & PAIRCOIL_PAIRS) {
1680 functnum = pair_functnum; /* For use in scscore.c **/
1681 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1682 b_paircoil[tablenum], &maxsc,
1683 tablenum,all_offsets[tablenum],
1684 all_scores[tablenum],
1685 Offset_to_Use, (mode & RAW_OUT));
1687 else /*Use singles probabilities... for likelihood use line number 0 */
1688 SingleCoilScore(mode,sequence,m_single[tablenum],
1691 all_scores[tablenum], Offset_to_Use,
1695 NEWCOILSScore(mode,sequence,&maxsc, all_scores[tablenum],
1700 HMMScore(sequence.seq, sequence.seqlen, hmm, all_scores[tablenum],
1701 tablenum, Offset_to_Use);
1706 ActualCoils(sequence, all_scores[tablenum], Offset_to_Use,0);
1709 if (tablenum == Table) *maxscore = maxsc;
1716 /*********************************************************************/
1717 /*********************************************************************/
1718 /********* This ends the scoring along all and all dimensions ********/
1719 /********* Now convert the raw scores that haven't been converted ****/
1720 /********* yet into likelihoods (i.e. the multidimensional scores.****/
1723 /***** Note: Can use same "init_class_prob" here, because since set number **/
1724 /*** of classes to 2, the 3rd class (Pdb-) will be ignored, so just get ****/
1725 /*** ratio of two to three stranded likelihoods, BOTH weighted by their ****/
1726 /*** init probaabiities. ********/
1728 if (Method != MultiCoil)
1729 get_offsets(all_scores, sequence,all_offsets,NUMBER_TABLES);
1730 else if (!(mode & RAW_OUT)) {
1732 switch_gauss_param(Offset_to_Use, Method);
1736 /**********Do conversions on MultiCoil dimension scores. ******/
1738 if (Method == MultiCoil) {
1739 convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
1740 NUMBER_TABLES, means_submatrix,inverse_covars,
1741 determinants, NUMBER_CLASSES,init_class_prob,
1743 mode & ONLY_COILED_CLASSES, bound,
1745 get_offsets(all_scores, sequence,all_offsets,3);
1748 if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
1749 for (tablenum=0; tablenum<NUMBER_CLASSES; tablenum++)
1750 average_score_over_coils3(sequence,
1751 all_scores[tablenum+NUMBER_CLASSES],
1752 all_scores[tablenum],
1753 all_scores[2+NUMBER_CLASSES],
1754 /** Total Coil Like. **/
1756 all_offsets[tablenum],
1757 start_coil_at_reg_shift, bound,
1763 } /* Ends if rescore_seq */
1766 /*****************************************************************************/
1767 if ( (Method == MultiCoil) || (Method == PairCoilDiff)){
1769 if (Method == MultiCoil)
1770 if ( (rescore_seq ==1 ) || /** Signals rescored MultiCoil **/
1771 (rescore_seq==-1)) /*Signals bound changed, need to redo seq_score */
1772 get_seq_scores(seq_scores, sequence, all_scores[0+ NUMBER_CLASSES],
1773 all_scores[1+ NUMBER_CLASSES],bound);
1776 if ( (rescore_seq==-1) && Coil_Score && !scored_already) {
1777 if (Method == MultiCoil) {
1778 for (tablenum=0; tablenum<NUMBER_CLASSES; tablenum++) {
1779 average_score_over_coils3(sequence,
1780 all_scores[tablenum+NUMBER_CLASSES],
1781 all_scores[tablenum],
1782 all_scores[2+NUMBER_CLASSES],
1783 /** Total Coil Like. **/
1785 all_offsets[tablenum],
1786 start_coil_at_reg_shift, bound,
1788 if (mode & ONLY_COILED_CLASSES) /** Convert to olig. state. **/
1789 type_class_score_convert(sequence,all_scores,bound,
1790 mode & ONLY_COILED_CLASSES);
1797 else if (rescore_seq == -2) /** Switch back to reg score.*/
1798 switch_back_to_res_score(bound, sequence,mode,all_scores);
1800 /****************************************************************************/
1803 return (double*)all_scores;
1810 ShowSeq(char lib[MAXFUNCTNUM],
1811 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1812 double *m, double *b,
1813 double *m_single, double *b_single,
1814 int mode, int rescore_sequence, int rescore_preproc, double bound)
1819 if ( (rescore_sequence) || (rescore_preproc))
1820 all_scores = ScoreSeq(lib, multi_lib, m, b,m_single, b_single,
1822 method, preprocessor_method,
1823 offset_to_use, &maxscore,
1824 rescore_sequence, rescore_preproc, bound);
1827 scores = &all_scores[table * MAXSEQLEN*(POSNUM+1)];
1828 /* This is all_scores[main_table][0][0] */
1829 preprocess_like = &all_preprocess_like[preproc_table][0][0];
1831 ShowScores(scores,&sequence,NULL, bound,mode);
1834 /*---------------------------------------------------------------------------*/
1839 OpenScreen(int mode, double log_bound, FILE *fin, FILE *flog,
1842 double *m, double *b, double *m_single, double *b_single,
1843 char lib[MAXFUNCTNUM],
1844 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1846 int main_preprocessor_method,
1847 int *seqcnt, long *pfreqp)
1851 XSetWindowAttributes attr;
1853 static unsigned char dashl[2] = {1,1};
1854 static unsigned char dashd[2] = {2,2};
1859 display = XOpenDisplay("");
1861 /* Check mode to see if want GUI. */
1862 if (!display || (mode & NO_GUI)) {
1863 /** COmplete logfile, etc. if can't open xdisplay. */
1864 if (!(mode & NO_GUI))
1865 fprintf(stderr,"Sorry, bad display.\nMaybe your DISPLAY environment variable isn't set properly,\nor you haven't xhosted this machine.\n");
1866 if ( fout || ( (mode & PRN_MODE) || (mode & VER_MODE)) ) {
1867 fprintf(stderr,"\nThe log file will still be written by:\n\n");
1868 for (i = 0; motd[i][0]; ++i) /* Print out title. */
1869 fprintf(stderr," %s\n",motd[i]);
1871 finish_log(mode,log_bound,fin,flog,fout_coils,
1873 m,b, m_single, b_single,
1874 lib,multi_lib, main_method, main_table,
1875 main_preprocessor_method, seqcnt, avg_max);
1879 if (flog) sclose(flog);
1880 if (fout) sclose(fout);
1881 if (ftotal_like) sclose(ftotal_like);
1888 font = XLoadFont(display,"9x15");
1889 window = XCreateSimpleWindow(display, RootWindow(display, 0),
1891 WIDTH-80, HEIGHT/4 +75,
1893 BlackPixel(display, 0),
1894 WhitePixel(display, 0));
1896 help = XCreateSimpleWindow(display, RootWindow(display, 0),
1898 WIDTH / 5 * 4 +60, HEIGHT / 5 * 3 + 400,
1900 BlackPixel(display, 0),
1901 WhitePixel(display, 0));
1902 helpb = XCreateSimpleWindow(display, window,
1904 BlackPixel(display, 0), WhitePixel(display, 0));
1905 next = XCreateSimpleWindow(display, window,
1907 BlackPixel(display, 0), WhitePixel(display, 0));
1908 prev = XCreateSimpleWindow(display, window,
1910 BlackPixel(display, 0), WhitePixel(display, 0));
1911 autob = XCreateSimpleWindow(display, window,
1913 BlackPixel(display, 0), WhitePixel(display, 0));
1914 zoomb = XCreateSimpleWindow(display, window,
1916 BlackPixel(display, 0), WhitePixel(display, 0));
1917 printb = XCreateSimpleWindow(display, window,
1919 BlackPixel(display, 0), WhitePixel(display, 0));
1920 tableb=XCreateSimpleWindow(display, window,
1922 BlackPixel(display, 0), WhitePixel(display, 0));
1924 methodb=XCreateSimpleWindow(display, window,
1926 BlackPixel(display, 0), WhitePixel(display, 0));
1928 paramb=XCreateSimpleWindow(display, window,
1930 BlackPixel(display, 0), WhitePixel(display, 0));
1931 quit = XCreateSimpleWindow(display, window,
1933 BlackPixel(display, 0), WhitePixel(display, 0));
1935 /** To get size: **/
1936 /** Printing " trimer x.xx" takes 12 char = 12*7 =84 **/
1937 /** Each line takes 15 of y coord. giving 3*15 =45 + 10 for blank line **/
1938 seqscoreb=XCreateSimpleWindow(display, window,
1939 XCOORD_SEQSCORES, YCOORD_SEQSCORES,
1941 BlackPixel(display, 0), WhitePixel(display, 0));
1943 zoom = XCreateSimpleWindow(display, RootWindow(display, 0),
1947 BlackPixel(display, 0),
1948 WhitePixel(display, 0));
1949 XSelectInput(display, helpb, ButtonPressMask);
1950 XSelectInput(display, next, ButtonPressMask);
1951 XSelectInput(display, prev, ButtonPressMask);
1952 XSelectInput(display, autob, ButtonPressMask);
1953 XSelectInput(display, zoomb, ButtonPressMask);
1954 XSelectInput(display, printb, ButtonPressMask);
1955 /**** XSelectInput(display, methodb, ButtonPressMask); ***/
1956 XSelectInput(display, tableb, ButtonPressMask);
1957 XSelectInput(display, paramb, ButtonPressMask);
1959 XSelectInput(display, quit, ButtonPressMask);
1960 XSelectInput(display, window, ExposureMask);
1961 XSelectInput(display, zoom, ExposureMask | KeyPressMask | ButtonPressMask);
1962 XSelectInput(display, help, ExposureMask | KeyPressMask | ButtonPressMask);
1963 attr.bit_gravity = NorthWestGravity;
1964 XChangeWindowAttributes(display,window,CWBitGravity,&attr);
1965 XMapWindow(display, window);
1968 values.foreground = BlackPixel(display, 0);
1969 values.background = WhitePixel(display, 0);
1971 gc = XCreateGC(display, window, GCForeground | GCBackground | GCFont, &values);
1972 values.line_style = LineDoubleDash;
1973 values.font = XLoadFont(display,"6x13");
1974 gcgray = XCreateGC(display, window, GCForeground | GCBackground
1975 | GCLineStyle | GCFont, &values);
1976 XSetDashes(display, gcgray, 1, dashl, 2);
1977 values.function = GXxor;
1978 values.foreground = 1;
1979 gcdash = XCreateGC(display, window, GCForeground
1980 | GCLineStyle | GCFunction, &values);
1981 XSetDashes(display, gcdash, 0, dashd, 2);
1985 icon = XCreateBitmapFromData(display,window,bart_bits,bart_width,bart_height);*/
1986 /* XSetStandardProperties(display,window,"Coiled Coil Finder","xcoil",icon,NULL,0,NULL);*/
1987 /* XSetStandardProperties(display,zoom,"XCoil Zoomer","xcoil",icon,NULL,0,NULL);*/
1988 /* XSetStandardProperties(display,help,"XCoil Help","xcoil",icon,NULL,0,NULL);*/
1992 PSDefs (FILE *ps, int seqlen)
1997 int nominallength= 300; /* Default value for number residues per line. */
1999 if ( (seqlen%nominallength <= 50) && (seqlen >= nominallength)) {
2000 nominallength = 350;
2003 /*************************/
2004 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
2005 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
2006 nominallength = seqlen;
2008 nominallength = ps_res_per_line; /* How many res per line. */
2009 /************************/
2011 if (seqlen ==0) number_lines=1;
2013 number_lines = (int)(seqlen-1)/nominallength +1;
2014 /* nominallength amino acids per line. */
2015 /* C rounds down and we want 0-nominallength to go to 1 */
2016 /* nominallength+1 to 2*noiminallength to go to 2.. */
2021 /*****Scale and translate to fit the page if want it to rotate page. */
2022 /* unscaled there are 9 lines per page. */
2023 /* .89 scaling if 300 res per line.
2024 xscaler = .89 * 300/nominallength;
2025 fputs("90 rotate\n",ps);
2026 fputs("78 -60 translate\n",ps);
2027 if (number_lines > 9) yscaler = (double) 9/number_lines;
2028 fprintf(ps,"%lf %lf scale\n", xscaler, yscaler);
2031 /*****Scale and translate to fit the page if don't want to rotate page. */
2032 /* unscaled there can fit 12 lines per page, so 3600. */
2033 xscaler = .75 * 300/nominallength; /* .75 scaling if 300 res per line. */
2034 fputs("25 740 translate\n",ps);
2035 if (number_lines > 12) yscaler = (double) 12/number_lines;
2036 fprintf(ps,"%lf %lf scale\n", xscaler,yscaler);
2040 fputs("/black {0 setgray} def\n",ps);
2041 fputs("/gray {0.6 setgray} def\n",ps);
2042 fputs("/tt /Courier findfont 15 scalefont def\n",ps);
2043 fputs("/rm /Times-Roman findfont 12 scalefont def\n",ps);
2044 fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
2045 fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
2046 fputs("/m {neg moveto} def\n",ps);
2047 fputs("/l {neg lineto stroke} def\n",ps);
2048 fputs("/rl {neg rlineto} def\n",ps);
2049 fputs("/rl_stroke {neg rlineto stroke} def\n",ps);
2050 fputs("/bar {dup 0.0019 lt {pop} {0 exch rlineto stroke} ifelse} def\n",ps);
2053 /*** static int xmin = 60; ***** Made it global variable *****/
2054 /*Originally was 29 */
2055 /** xmin determines where in window start seq. **/
2056 static int nominallength;
2059 static int ydelta=57;
2060 /* If click at x,y, determine which residue position it corresponds to.
2063 XYtoPos (int x, int y)
2066 if (x<xmin) return 0;
2068 if (i>= nominallength) return 0;
2069 for (y -= ybase+ydelta; y>0; y-=ydelta, i+=nominallength);
2070 if (y<-30) return 0;
2071 if (i>=length) return 0;
2074 PostoXY (int i, int *x, int *y)
2077 *x = xmin + 2*(i%nominallength);
2078 *y = ybase+ydelta + ydelta*(i/nominallength);
2080 psstr (char *str, int len, FILE *ps)
2083 for (i=0; i<len; ++i)
2095 /*** ybase determines where top of the seq score box should be so don't ***/
2096 /*** print over the sequenc title. **/
2097 draw_seqscore_box(double bound, double seq_scores[MAX_CLASS_NUMBER],
2098 Window *seqscoreb, FILE *ps, int ybase)
2101 /* XWindowChanges values; */
2105 XClearWindow(display,*seqscoreb);
2107 /** Relocate seqscoreb base on ybase**/
2108 /**This is complex way to move it.****
2109 values.y= YCOORD_SEQSCORES+ybase-30;
2110 XConfigureWindow(display,*seqscoreb,CWY,&values);
2113 XMoveWindow(display,*seqscoreb,XCOORD_SEQSCORES, YCOORD_SEQSCORES+ybase-30);
2117 /** 3 means print 3 into the box horiz, 10 means print 10 into box vert? */
2119 XMapWindow(display, *seqscoreb);
2120 XDrawString(display,*seqscoreb,gcgray,3,10,"Seq Prob:",9);
2122 if (ps) { /* Draw the box. **/
2123 fprintf(ps,"%d %d m\n",10, 50+ybase-30);
2124 fprintf(ps,"0 40 rl\n");
2125 fprintf(ps,"79 0 rl\n");
2126 fprintf(ps,"0 -40 rl\n");
2127 fprintf(ps,"-79 0 rl_stroke\n");
2130 if (seq_scores[0] + seq_scores[1] + seq_scores[2] > 0) { /* Signals above */
2131 sprintf(buf," Dimer %5.2lf", seq_scores[0]); /* bound. */
2132 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2133 sprintf(buf," Trimer %5.2lf", seq_scores[1]);
2134 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2137 fprintf(ps,"rm setfont ");
2138 fprintf(ps,"%d %d m (Seq Prob for) show\n", 15,60 +ybase-30);
2140 fprintf(ps,"%d %d m (bound %5.2lf:) show\n", 15,70 +ybase-30,bound);
2142 fprintf(ps,"%d %d m ( Dimer %5.2lf) show\n",15,80 +ybase-30,
2144 fprintf(ps,"%d %d m ( Trimer %5.2lf) show\n",15,90 +ybase-30,
2151 sprintf(buf," No Residue ");
2152 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2153 sprintf(buf," above %5.2lf ", bound);
2154 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2157 fprintf(ps,"rm setfont ");
2158 fprintf(ps,"%d %d m (Seq Prob:) show\n", 15,60 +ybase-30);
2160 fprintf(ps,"%d %d m ( No Residue) show\n", 15, 75 +ybase-30);
2161 fprintf(ps,"%d %d m ( above %5.2lf) show\n",15,85 +ybase-30, bound);
2169 draw_bound_box(double bound, Window paramb)
2173 XMapWindow(display, paramb);
2174 XDrawString(display,paramb,gcgray,3,10,"bound",5);
2175 sprintf(buf,"%5.2lf", bound);
2176 XDrawImageString(display,paramb,gcgray,3,25,buf,strlen(buf));
2179 ShowScores (double *scores, Sequence *sequence, FILE *ps, double bound,
2182 int i,q, x, y0, xmax;
2187 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2188 char *MethodTitle = MethodTitleArray;
2189 char PreprocTitleArray[500]; /* For printing the preproc.name to ps file. */
2190 char *PreprocTitle = PreprocTitleArray;
2193 MethodTitleArray[0] = '\0'; /* Initialization. */
2194 PreprocTitleArray[0]='\0';
2196 nominallength = 300;
2197 length = sequence->seqlen;
2199 /* This makes display look nice (so don't have a line with < 50 residues. */
2200 if ((length%nominallength <= 50) && (length>= nominallength)){
2201 nominallength = 350;
2202 /* if (ps) fputs("-50 0 translate\n",ps); */
2206 /*************************/
2207 if (ps) /* Do this stuff only for postscript printout. */
2208 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
2209 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
2210 nominallength = sequence->seqlen;
2212 nominallength = ps_res_per_line; /* How many res per line. */
2213 /************************/
2216 charperline = (xmin+2*nominallength-5)/9; /*Each char is "9" units witdth */
2217 rows = (length-1)/nominallength+1; /* Number of rows to display. */
2219 xmax = 2*length+xmin; /* The max x position in current row. */
2221 xmax = 2*nominallength+xmin;
2222 width = xmin+2*nominallength+ WIDTH_ADD_ON; /* Width of window to display. */
2223 /* The +WIDTH_ADD_ON used to be +11, but */
2224 /* was cutting of some numbers. */
2225 x = strlen(sequence->title);
2226 XClearWindow(display,window);
2228 /**********This is where window size is adjusted to fit the sequence. ***/
2229 /********** using width and number of rows. **************************/
2231 XResizeWindow(display,window,width,
2232 38+15*((x+charperline-1)/charperline) +ydelta*rows);
2234 /*************************************************************************/
2235 XMapWindow(display, helpb);
2236 XMapWindow(display, next);
2237 XMapWindow(display, prev);
2238 XMapWindow(display, autob);
2239 XMapWindow(display, zoomb);
2240 XMapWindow(display, printb);
2241 XMapWindow(display, tableb);
2243 XMapWindow(display, methodb);
2244 XMapWindow(display, paramb); Now done in draw_bound_box()
2246 XMapWindow(display, quit);
2249 while (XCheckTypedEvent(display,Expose,&report));
2251 XDrawString(display,helpb,gcgray,3,10,"help",4);
2252 XDrawString(display,next,gcgray,3,10,"next",4);
2253 XDrawString(display,prev,gcgray,3,10,"prev",4);
2254 XDrawString(display,autob,gcgray,3,10,"auto",4);
2255 XDrawString(display,zoomb,gcgray,3,10,"zoom",4);
2256 XDrawString(display,printb,gcgray,3,10,"print",5);
2257 /* XDrawString(display,methodb,gcgray,3,10,"method",6); */
2258 XDrawString(display,tableb,gcgray,3,10,"table",5);
2260 draw_bound_box(bound, paramb);
2261 /** draw_seqscore_box(bound, seq_scores,?? &seqscoreb, ps, ybase);
2262 MOVED TO LATER (AFTER ybase gets adjusted at end of this function). ****/
2265 XDrawString(display,quit,gcgray,3,10,"quit",4);
2268 /************This section prints out the method name..... ********************/
2269 /**** Add 9 to xcoord for each character. *****/
2270 XDrawString(display,window,gc, XCOORD_SCORED_BY,15,"Scored by ",10);
2272 make_methodtitle(MethodTitle, method, table, mode,0);
2273 XDrawString(display,window,gc,XCOORD_SCORED_BY +90,15,MethodTitle,
2274 strlen(MethodTitle));
2276 XDrawString(display,window,gc, XCOORD_SCORED_BY,30,"Filter by ",10);
2277 make_methodtitle(PreprocTitle, preprocessor_method, preproc_table,
2279 XDrawString(display,window,gc,XCOORD_SCORED_BY + 90,30,
2280 PreprocTitle,strlen(PreprocTitle));
2282 /***************************************************************************/
2284 XDrawString(display,window,gc, 5,30,sequence->code,strlen(sequence->code));
2285 for (i=0,ybase=30; x > i + charperline; i += charperline)
2286 XDrawString(display,window,gc, 5,ybase+=15,&(sequence->title[i]),charperline);
2287 XDrawString(display,window,gc, 5,ybase+15,&(sequence->title[i]),x-i);
2290 fprintf(ps,"tt setfont 5 15 m (%s) show\n",sequence->code);
2291 /* fprintf(ps,"175 15 m (Viewgram by XCoil) show\n"); */
2293 fprintf(ps,"430 15 m (Scored by %s) show\n",methodname[method]);
2295 fprintf(ps,"180 15 m (Scored by %s) show\n",MethodTitle);
2296 fprintf(ps,"180 26 m (Filter by %s) show\n",PreprocTitle);
2299 for (i=0,ybase=30; x > i + charperline; i += charperline) {
2300 fprintf(ps,"5 %d m (",ybase+=15);
2301 psstr(&(sequence->title[i]),charperline,ps);
2302 fprintf(ps,") show\n");
2304 fprintf(ps,"5 %d m (",ybase+15);
2305 psstr(&(sequence->title[i]),x-i,ps);
2306 fprintf(ps,") show\n");
2307 fputs("rm setfont 2 setlinewidth\n",ps);
2309 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2311 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2312 y = y0-30*scores[8*i+7];
2313 XDrawLine(display, window, gc, x, y0, x, (int)(y+0.5));
2314 XDrawLine(display, window, gc, x+1, y0, x+1, (int)(y+0.5));
2316 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2319 /********************* If in pos mode draws in the correct likelihoods every
2320 now and then? ******************
2321 if (sequence->reg && method==HMMCoil) {
2324 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2325 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2326 if (sequence->reg[i] != '.') {
2327 double y = y0-30*scores[8*i+sequence->reg[i]];
2328 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2329 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2331 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2334 if (ps) fputs("black\n",ps);
2336 **********************/
2338 /**********Now draw in the PreProcScore in gray every other line. ****/
2342 for (i=0,y0=ybase; i<length; i+=2,x+=4) { /* Normally is ++i,x+=4, but */
2343 double y; /* we do every other line here. */
2344 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2346 y = y0-30*preprocess_like[i*(POSNUM +1) + 7];
2347 /* This is preprocess_like[i][7]. */
2348 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2349 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2351 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2357 if (ps) fputs("black\n",ps);
2359 /*************************************************************************/
2361 if (ps) fputs("0 setlinewidth\n",ps);
2362 for (i=9,y0=ybase; i<length; i+=10) {
2363 if (!((i-9)%nominallength)) {x=xmin-1; y0 += ydelta;}
2365 XDrawLine(display, window, gc, x,y0,x,y0+2); /* Draw ticks every */
2366 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+2); /* 10 residues */
2369 for (i=49,y0=ybase; i<length; i+=50) {
2371 if (!((i-49)%nominallength)) {x=xmin-1; y0 += ydelta;}
2373 XDrawLine(display, window, gc, x,y0,x,y0+5);
2374 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+5);
2375 sprintf(buff,"%d",i+1);
2376 XDrawString(display, window, gcgray, x-7, y0+18, buff, strlen(buff));
2377 if (ps) fprintf(ps,"%d %d m (%d) cshow\n",x,y0+18,i+1);
2379 for (i=1,y0=ybase+ydelta; i<=rows; ++i,y0+= ydelta) {
2381 xmax = 2*(length-(rows-1)*nominallength)+xmin; /* Max x position in */
2383 XDrawString(display,window,gcgray,xmin-25,y0-26,"100%",4);
2384 XDrawString(display,window,gcgray,xmin-19,y0-11,"50%",3);
2385 XDrawString(display,window,gcgray,xmin-13,y0+4,"0%",2);
2386 XDrawLine(display,window,gcgray,xmin,y0-30,xmax-1,y0-30);
2387 XDrawLine(display,window,gcgray,xmin,y0-15,xmax-1,y0-15);
2389 fprintf(ps,"%d %d m (100%%) rshow\n",xmin-1,y0-26);
2390 fprintf(ps,"%d %d m (50%%) rshow\n",xmin-1,y0-11);
2391 fprintf(ps,"%d %d m (0%%) rshow\n",xmin-1,y0+4);
2392 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-30,xmax,y0-30);
2393 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-15,xmax,y0-15);
2394 fprintf(ps,"%d %d m %d %d l\n",xmin,y0,xmax,y0);
2399 draw_seqscore_box(bound, seq_scores, &seqscoreb, ps, ybase);
2401 if (zoomptr) ZoomPnt();
2404 fputs("showpage\n",ps);
2413 make_methodtitle(char *MethodTitle, int Method, int which_table,
2414 int mode, int for_logfile)
2421 if (which_table ==0) like2or3 =2;
2425 /************This section prints out the method name..... ********************/
2426 /**** Add 9 to xcoord for each character. *****/
2430 if ((Method == PairCoilDiff) || (Method == PairCoilDiffAvg)) {
2431 /* Does Differen. use dimer or trimer like. */
2433 /* XDrawString(display,window,gc,xcoord,15,"Dimer",5); */
2434 MethodTitle= strcat(MethodTitle, "Dimer");
2437 else if (like2or3==3) {
2438 /* XDrawString(display,window,gc,xcoord,15,"Trimer",6); */
2439 MethodTitle= strcat(MethodTitle, "Trimer");
2445 if (Method == PairCoil) {
2446 if (mode & PAIRCOIL_PAIRS) { /* Are we using pair or single prob */
2447 /* XDrawString(display,window,gc,xcoord,15,"Pair",4); */
2448 MethodTitle= strcat(MethodTitle, "Pair");
2452 /* XDrawString(display,window,gc,xcoord,15,"Single",6); */
2453 MethodTitle= strcat(MethodTitle, "Single");
2457 } /* Ends the stuff for dealing with PairCoil, TDRes, TDCoil. */
2460 /*****************************************************************************/
2461 /*************************** Now write the methodname. *********************/
2463 /* XDrawString(display,window,gc,xcoord,15,methodname[method],strlen(methodname[method])); */
2464 xcoord += 9*strlen(methodname[Method]);
2465 MethodTitle= strcat(MethodTitle, methodname[Method]);
2467 if (Method == MultiCoil) {
2468 if (mode & ONLY_COILED_CLASSES)
2469 strcat(MethodTitle,"Olig");
2472 strcat(MethodTitle,"WtdAvg");
2474 else if (Coil_Score ==2)
2475 strcat(MethodTitle, "MAX");
2477 if (start_coil_at_reg_shift)
2478 strcat(MethodTitle, "Shift");
2482 /* Now print Table are using. */
2483 /* And if are using a particular offset, print that. */
2484 if ( (Method==MultiCoil) || (Method==PairCoil) || (Method == HMMCoil)
2485 || (Method==PairCoilDiff) || (Method==PairCoilDiffAvg)) {
2488 if (NUMBER_TABLES >1) /** MultiCoil logfile uses all tables.... ***/
2489 if ( (Method != MultiCoil) || (!for_logfile)) {
2490 sprintf(buf,"%d",which_table+1);
2491 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2492 xcoord += 9*strlen(buf);
2493 MethodTitle= strcat(MethodTitle, buf);
2497 sprintf(buf,""); /* Clear out buf. */
2500 if ((offset_to_use != -1) && (offset_to_use != 7))
2501 sprintf(buf,"%c",'a' + offset_to_use);
2502 else if (Method == MultiCoil)
2503 if (offset_to_use == -1) {
2504 strcat(MethodTitle,"Max");
2508 else if (offset_to_use == 7) {
2509 strcat(MethodTitle,"Combo");
2514 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2515 xcoord += 9*strlen(buf);
2516 MethodTitle= strcat(MethodTitle, buf);
2520 /***************************************************************************/
2523 int is_not_differentiator(int Method)
2525 if ( (Method == MultiCoil) || (Method == PairCoil) || (Method==HMMCoil) ||
2526 (Method==NEWCOILS) || (Method == ActualCoil) )
2535 log_file_header(FILE *flog, int mode, int argc, char *argv[],
2536 int avg_max, double bound,
2537 int preprocessor_method, int preproc_table,
2538 int log_offset_to_use, int start_coil_at_reg_shift,
2539 int main_method, int main_table,
2540 char pair_lib[MAXFUNCTNUM],
2541 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM])
2544 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2545 char *MethodTitle = MethodTitleArray;
2547 MethodTitleArray[0] = '\0'; /* Initialization. */
2548 make_methodtitle(MethodTitle, method,main_table,mode,1);
2550 fprintf(flog,"\nInput file is %s", argc>1?argv[1]:"-");
2552 fprintf(flog, "\n\nScores by %s.", MethodTitle);
2554 "\nThe probability cutoff for the coiled-coil locater is %.2lf.",
2557 if (main_method == PairCoil) {
2558 if (mode & WEIGHTED_PROBS) {
2559 fprintf(flog,"\nStatistically weighted PairCoil score is ");
2562 fprintf(flog,"\nPairCoil score is ");
2565 fprintf(flog,"AVERAGE of residue scores.");
2568 fprintf(flog,"MAX of residue scores.");
2571 if (mode & USE_LIKE_LINE) {
2572 fprintf(flog,"\nPaircoil used the likelihood line.");
2575 fprintf(flog,"\nPaircoil used direct prob. formula instead of likelihood line.");
2578 fprintf(flog, "\nThe distances used by PairCoil are: ");
2580 for (i=0; i<pair_functnum; i++) {
2581 fprintf(flog,"%d ",pair_lib[i]);
2585 if ( (main_method==PairCoilDiff) || (main_method==PairCoilDiffAvg)) {
2586 char MethodTitleArray[500]; /* For printing the methodname. */
2587 char *MethodTitle = MethodTitleArray;
2588 int oligimerization;
2590 MethodTitleArray[0] = '\0'; /* Initialization. */
2591 make_methodtitle(MethodTitle, preprocessor_method, preproc_table,mode,0);
2594 oligimerization == 2;
2595 else oligimerization ==3;
2597 fprintf(flog, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2598 /* Whether logfile is for 2 or */
2599 /* 3 stranded predictions. */
2601 if (main_method==PairCoilDiffAvg) {
2602 if (start_coil_at_reg_shift) {
2603 fprintf(flog,"\n New coils were considered to start at register shifts.");
2606 fprintf(flog,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2612 fprintf(flog,"\nMain Table is Table %d\n\n",main_table);
2615 if (main_method == MultiCoil) {
2617 fprintf(flog, "\nThe scoring distances are: ");
2618 for (tablenum =0; tablenum<NUMBER_TABLES; tablenum++) {
2619 /* fprintf(flog,"\n For table %d: ",tablenum+1); */
2620 if (tablenum == 0) fprintf(flog,"\n For dimeric table: ");
2621 if (tablenum == 1) fprintf(flog,"\n For trimeric table: ");
2622 for (i=0; i<num_dist[tablenum]; i++) {
2623 fprintf(flog,"%d ",multi_lib[tablenum][i]);
2628 if (mode & POS_MODE)
2629 fprintf(flog,"\nIf the start or end of a coil is listed inside (), then that position does not occur in the pos file.\n");
2637 coil_tuple_output(int number_tables, char *lib, int mode, int avg_max,
2638 double bound, Sequence sequence, FILE *fout_coils,
2639 int number_multi_lib[MAX_TABLE_NUMBER],
2640 int number_score_dim)
2642 char offsets[MAXSEQLEN];
2646 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
2647 extern double *all_scores;
2650 for (tab=0; tab< number_tables; tab++) {
2651 for (dist=0; dist<number_multi_lib[tab]; dist++) {
2653 /*************************************************************************/
2654 /**********This computes offset of that table/dist method of scoring. ***/
2656 for (i=0; i<sequence.seqlen; i++) {
2657 for (j=0; j<POSNUM; j++)
2658 if (all_preprocess_like[tab + number_tables*(int)dist][i][j] ==
2659 all_preprocess_like[tab + number_tables*(int)dist][i][7]) {
2660 offsets[i] = (j-i) % POSNUM;
2661 if (offsets[i]<0) offsets[i] += POSNUM;
2663 /* This is a new thing so that if an offset change is not real */
2664 /* we keep the old offset. **********/
2665 if ( (i>0) && (offsets[i-1] != -1) &&
2666 (offsets[i] != offsets[i-1]) &&
2667 (all_preprocess_like[tab + number_tables*(int)dist]
2668 [i][(i+offsets[i-1]) %POSNUM]
2669 == all_preprocess_like[tab+number_tables*(int)dist][i][7]) )
2670 offsets[i] = offsets[i-1];
2673 /**************************************************************************/
2674 /****** Now compute coil scores for that method (table,dist) pair. ***/
2675 get_raw_coil_score(sequence,
2676 all_preprocess_like[tab+number_tables*(int)dist],
2677 all_scores[tab+number_tables*(int)dist], mode,
2678 avg_max,bound,offsets);
2682 /************* Now output the tuples of coil scores. ****/
2683 tuple_output2(sequence, mode, fout_coils,
2685 bound, number_score_dim);
2690 switch_gauss_param(int Offset_to_Use, int Method)
2692 if (Method == MultiCoil) {
2693 if (Offset_to_Use == 7) { /* Combo offset */
2694 determinants= &both_determinants[0][0];
2695 inverse_covars= &both_inverse_covars[0][0][0][0];
2696 means_submatrix= &both_means_submatrix[0][0][0];
2698 else { /* max offset or some fixed offset */
2699 determinants= &both_determinants[1][0];
2700 inverse_covars= &both_inverse_covars[1][0][0][0];
2701 means_submatrix= &both_means_submatrix[1][0][0];
2708 void usage(char *prog_command)
2713 /** Strip off directories to get just the program name. **/
2714 i=strlen(prog_command);
2715 while ((i>=0) && (prog_command[i] != '/') )
2717 prog_name= &prog_command[i+1];
2719 /************ This is for my extra stuff......
2720 fprintf(stderr,"\nUsage: %s <file-name>\n [Table to Remove seq from...\n could be dimer, trimer, dimer-1, trimer-1, negatives, config_file]\n [<distances with table1> <distances with table2>]\n <-config config_file>\n",prog_name);
2721 **********************/
2723 fprintf(stderr,"\n Usage: %s <file-name> [-config config_file_name]\n",prog_name);