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[])
280 struct itimerval interval;
284 char buff[MAXLINE+MAXSEQLEN];
288 char *command_line = NULL;
290 FILE *fgin=NULL, *fout=NULL, *fout_coils=NULL, *flog=NULL;
291 FILE *fpin[MAX_TABLE_NUMBER];
293 char pir_name[MAXLINE];
294 char likelihoods[MAX_TABLE_NUMBER][MAXLINE];
295 /* Files of likelihood lines. */
296 char printfile[500]; /* The filename to output postscript print to. */
297 char *print= printfile;
299 char class_sc_filenames[2][MAX_CLASS_NUMBER][MAXLINE];
300 double both_class_covars[2][MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX]
301 [NUM_DIM_IN_ORIG_MATRIX];
302 double both_class_means[2][MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX];
304 char gauss_param[2][MAXLINE];
305 /** 0 entry is for combo offset, 1 entry is for **/
306 /** max of other registers. ***/
307 extern int functnum; /* What is used in scscore.c. Takes value from */
308 /* pair_functnum and dimension stuff for multicoil. */
310 char lib[MAXFUNCTNUM]; /* Use MAXFUNCTUM=7 here since are 7 distances */
311 /* in coil, so we know AT MOST 7 lib functions. */
312 extern int pair_functnum; /* The number of distances in lib[].*/
314 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM];
315 extern int num_dist[MAX_TABLE_NUMBER];
319 int mode = 0; /* If the i-th option is set, mode gets 1 in i-th bit. */
322 double m[MAX_TABLE_NUMBER], m_single[MAX_TABLE_NUMBER],
323 b[MAX_TABLE_NUMBER], b_single[MAX_TABLE_NUMBER];
324 /*The likelihood line is mx+b for paircoil. */
328 int number_tables =0; /** Currently not used (replaced by NUMBER_TABLES */
330 int table_to_remove_from = -1;
331 /* Used if want to remove ALL input seq */
332 /* from that table at outset. **/
333 char *command_line_config = NULL;
335 /*** INITIALIZATION ***/
338 printf("writing environment\n");
345 printf("%s\n",equality);
348 initialize(&num_dist[0], window_length, scale0s,
350 likelihoods, &pir_name[0], print, gauss_param,
351 &mode, init_class_prob);
353 printf("Printing arguments\n");
354 for (i=0;i<argc;i++){
355 printf("%s\n",argv[i]);
357 read_command_line(argc, argv, &command_line, combine_dist,
358 &mode, num_dist, multi_lib, &command_line_config);
360 fprintf(stderr,"\n\nRead commandline, command_line=%s",command_line);
362 /* Get the default options in the .paircoil file!!! */
363 get_defaults(command_line_config, argv[1], &mode, &log_bound, &fgin, fpin,
365 &flog, &fout_coils, &by_coil_or_seq, likelihoods,
366 pir_name, lib, multi_lib,
367 &pair_functnum, combine_dist,
368 print, &main_method, &main_preprocessor_method,
369 &main_table, &offset_to_use,
370 &avg_max, &Coil_Score,
373 gauss_param, num_dist,
375 &table_to_remove_from,
376 command_line, window_length, scale0s, scale0p);
378 if ((mode & RAW_OUT) && flog) {
379 fprintf(stderr,"\n Log file will not be written when doing raw scores.");
383 if ((mode & RAW_OUT) && fout_coils) {
385 "\n Sequence scores will not be written when doing raw scores.");
392 fprintf(stderr,"\nTable %d, num_dist = %d, combine_dist =%d", i,
393 num_dist[i], combine_dist[i]);
395 fprintf(stderr,"\nOut of get_def");
401 if (ftotal_like) { /** Want to output total gauss like over classes*/
402 for (i=0; i<NUMBER_CLASSES; i++) /* for various init_prob settings. */
403 init_class_prob[i]=1; /** All 1 values is flag to just output */
404 init_totals(total_gauss_like); /** gauss values, not like values in */
405 } /** convert_raw_scores_to_gauss_like. */
407 /** 7*(1-combine_dist) + 1*combine_dist gives the number of **/
408 /** dimensions started with, and **/
409 /** num_dist*(1-combine_dist) + 1*combine_dist gives the new number. **/
410 number_score_dim = num_table_dimensions(num_dist[0],combine_dist[0]) +
411 num_table_dimensions(num_dist[1],combine_dist[1]);
414 if (gauss_param[i][0] != ',') {
415 get_gauss_param_all_classes2(class_sc_filenames[i],
416 num_table_dimensions(7,combine_dist[0]) +
417 num_table_dimensions(7,combine_dist[1]),
418 both_class_means[i], both_class_covars[i],
421 /******Note we store the submatrices in the original matrices....****/
422 compute_submatrices2(both_class_means[i],both_means_submatrix[i],
423 both_class_covars[i],
424 both_inverse_covars[i],
425 both_determinants[i], multi_lib,
426 num_table_dimensions(7,combine_dist[0]),
427 num_table_dimensions(7,combine_dist[1]),
428 num_table_dimensions(num_dist[0],combine_dist[0]),
429 num_table_dimensions(num_dist[1],combine_dist[1]));
430 #ifdef DEBUG_PRINTOUT
431 for (class=0; class<NUMBER_CLASSES; class++) {
432 fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
433 print_matrix(both_means_submatrix[i][class],1,
434 number_score_dim,stderr);
435 fprintf(stderr,"\n\nInver Covar Submatrix %d:\n", class);
436 print_matrix(both_inverse_covars[i][class],
437 number_score_dim, number_score_dim,stderr);
441 fprintf(stderr,"\nDet0 =%lf, det1= %lf, det2= %lf\n",
442 both_determinants[i][0],both_determinants[i][1],
443 both_determinants[i][2]);
445 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",
446 num_table_dimensions(num_dist[0],combine_dist[0]),
447 num_table_dimensions(num_dist[1],combine_dist[1]),
448 num_table_dimensions(7,combine_dist[0]),
449 num_table_dimensions(7,combine_dist[1]),
456 switch_gauss_param(offset_to_use, main_method);
461 log_offset_to_use = offset_to_use;
462 method = main_method;
463 preprocessor_method = main_preprocessor_method;
467 if ( (mode & PRN_MODE) && flog ){
468 log_file_header(flog, mode,
469 argc, argv, avg_max, log_bound,
470 preprocessor_method, preproc_table, log_offset_to_use,
471 start_coil_at_reg_shift, main_method,
472 main_table, lib, multi_lib);
475 if (mode & POS_MODE) ActualCoil_Method = available;
477 /****** Now if didn't get locations from .paircoil, get from ENV variable */
478 if (!fgin) { /* Didn't get fgin in */
479 fprintf(stderr,"\n\n **********************");
480 fprintf(stderr,"\nDid not get genbnk location from config file.");
481 fprintf(stderr,"\n Will instead use uniform distrib. as underlying distrib.");
482 fprintf(stderr,"\n**********************\n\n");
484 /**** exit(-1); ****/ /* get_defaults. */
486 if (!fpin[0]) { /* Didn't get fpin in */
487 fprintf(stderr,"\nDid not get positive data set from config file.");
488 exit(-1); /* get_defaults. */
490 /******************************************/
493 PairCoilData(fgin,fpin, mode & PROLINE_FREE_MODE,
494 table_to_remove_from, argv[1], mode);
495 /* If fpin[1] != NULL this will call */
496 /* TrimerDimerDiff to set up table diff */
498 /* If it is NULL then table is average */
499 /* of probabilities from 1 and 2. */
503 hmm = CCHMM(fgin,fpin);
507 if (mode & USE_LIKE_LINE) {
508 /** Don't want to compute likelihood line if */
510 functnum = pair_functnum;
511 get_likelihood_line(likelihoods, m, b, m_single, b_single,
512 lib, pair_functnum,pir_name,mode,NUMBER_TABLES);
515 #ifdef DEBUG_PRINTOUT
516 for (i=0; i<NUMBER_TABLES; i++)
517 fprintf(stderr, "\n m[%d]= %lf, b[%d]= %lf\n",i,m[i],i,b[i]);
522 while (i< NUMBER_TABLES) {
526 if (fgin) sclose(fgin);
532 fin = sopen(argc>1?argv[1]:"-","r");
534 fprintf(stderr, "\nfin is opened to %s",argc>1?argv[1]:"-");
537 /*** Note: If can't open xdisplay then OpenScreen will write log and */
538 /*** out files if those options are set. */
540 OpenScreen(mode,log_bound,fin,flog, fout_coils,
541 fout, m,b,m_single, b_single,
543 main_method,main_preprocessor_method, &seqnum);
546 /* Let user read message while initializing */
552 interval.it_value.tv_sec = 2;
553 interval.it_value.tv_usec = 0;
554 interval.it_interval.tv_sec = 2;
555 interval.it_interval.tv_usec = 0;
558 XSelectInput(display, window, KeyPressMask | ButtonPressMask | ExposureMask);
560 if (autoadvance) do {
561 FD_SET(ConnectionNumber(display),&fds);
562 if (!select(FD_SETSIZE,&fds,&fdnone,&fdnone,&interval))
563 NextSeq(fin,lib,multi_lib,m,b,m_single,b_single,
564 mode,bound,flog,fout_coils,fout);
566 while ( (!XPending(display)) && (autoadvance));
567 /* Autoadvance until button is hit or end of seq file. */
568 XNextEvent(display,&report);
569 switch (report.type) {
571 if (report.xany.window == zoom)
573 else if (report.xany.window == help)
575 else if (sequence.seqlen)
576 ShowScores(scores,&sequence,NULL,bound, mode);
579 #ifdef SHOW_DISCLAIMER
580 XNextEvent(display, &report);
581 ShowDisclaimer(window);
586 /* if (!sequence.seqlen) {NextSeq(fin,lib,multi_lib,m,b,m_single,
587 b_single,bound,flog,fout_coils,fout);
590 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
591 if (report.xany.window == zoom) {XUnmapWindow (display, zoom); break;}
592 if (report.xany.window == window) {
593 /* Zoom in and show register */
595 tmp = XYtoPos(report.xbutton.x,report.xbutton.y);
597 ZoomPnt(); zoomptr=tmp; Zoomer(); ZoomPnt();
600 if (report.xany.window == helpb)
601 XMapWindow(display,help);
602 if (report.xany.window == next)
603 NextSeq(fin,lib,multi_lib,m,b,m_single, b_single,
604 mode,bound,flog,fout_coils,fout);
605 if (report.xany.window == prev)
606 PrevSeq(fin,lib,multi_lib,m,b,m_single,b_single,mode,bound);
607 if (report.xany.window == autob)
608 { if (autoadvance ^= 1) NextSeq(fin,lib,multi_lib,m,b,m_single,
609 b_single,mode,bound,flog,fout_coils,
613 if (report.xany.window == zoomb) {
614 if (zoomptr) ZoomPnt();
615 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
617 if (report.xany.window == printb) {
619 ps = sopen(print?print:defprint,"a");
620 PSDefs(ps, sequence.seqlen);
621 ShowScores(scores,&sequence,ps,bound,mode);
624 if (report.xany.window == tableb) {
625 if (method == MultiCoil)
626 table = (table +1) % NUMBER_CLASSES;
627 else if (method == HMMCoil)
628 table = (table +1) % NUMBER_TABLES;
630 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,0, bound);
631 if (zoomptr) Zoomer();
635 if (report.xany.window == methodb) {
636 method_cycle(&method,NUMBER_METHODS);
637 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,0, bound);
638 if (zoomptr) Zoomer();
641 if (report.xany.window == quit) {
642 if (flog || fout || ftotal_like || fout_coils)
645 finish_log(mode,log_bound,fin,flog,fout_coils,
646 fout,m,b,m_single,b_single,
647 lib,multi_lib,main_method,main_table,
648 main_preprocessor_method,
653 if (flog) sclose(flog); if (fout) sclose(fout);
654 if (fout_coils) sclose(fout_coils);
655 if (ftotal_like) sclose(ftotal_like);
661 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
662 switch (XLookupKeysym(&(report.xkey),0)) {
665 if (bound > 1) bound =1;
666 draw_bound_box(bound, paramb);
670 if (bound < 0) bound =0;
671 draw_bound_box(bound, paramb);
674 case XK_Return: /** Signal that bound has been set, should rescore **/
675 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
676 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,-1,-1, bound);
677 if (zoomptr) Zoomer();
681 case XK_C: /** Switch between register, avg, and max **/
682 case XK_c: /** coil score for MultiCoil. **/
683 Coil_Score = (Coil_Score +1) % 3;
685 ShowSeq(lib,multi_lib, /** Change to res score **/
686 m,b,m_single,b_single,mode,-2,-2, bound);
687 else ShowSeq(lib,multi_lib, /** Change to coil score with */
688 m,b,m_single,b_single,mode,-1,-1, bound); /* out recomputing */
689 if (zoomptr) Zoomer(); /* res score. */
692 case XK_O: /** Oligomerization ratio. ***/
694 mode ^= ONLY_COILED_CLASSES; /* To XOR to complement that bit. */
695 if (NUMBER_CLASSES >0) {
696 if (method == MultiCoil)
697 type_class_score_convert(sequence,all_scores,bound,
698 mode & ONLY_COILED_CLASSES);
699 if (preprocessor_method == MultiCoil)
700 type_class_score_convert(sequence,all_preprocess_like,bound,
701 mode & ONLY_COILED_CLASSES);
703 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,2, bound);
704 if (zoomptr) Zoomer();
709 ps = sopen(print?print:defprint,"w");
710 PSDefs(ps, sequence.seqlen);
711 ShowScores(scores,&sequence,ps,bound,mode);
717 {ZoomPnt(); --zoomptr; Zoomer(); ZoomPnt();}
720 if (zoomptr && zoomptr<sequence.seqlen)
721 {ZoomPnt(); ++zoomptr; Zoomer(); ZoomPnt();}
725 if (zoomptr) ZoomPnt();
726 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
729 #ifdef DEBUG_PRINTOUT
731 offset_to_use =-1; /* Do max scoring register */
732 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,1, bound);
733 if (zoomptr) Zoomer();
737 offset_to_use =7; /* Do MultiCoil max score as combo of max score */
738 ShowSeq(lib,multi_lib, /* in each dimension. */
739 m,b,m_single,b_single,mode,1,1, bound);
740 if (zoomptr) Zoomer();
747 XMapWindow(display,help);
752 if (report.xany.window == window)
753 {sclose(fin); if (flog) sclose(flog);
754 if (ftotal_like) sclose(ftotal_like);
758 if (fout) sclose(fout);
760 if (report.xany.window == zoom) {
762 XUnmapWindow(display, zoom);
768 method_cycle(&method,NUMBER_METHODS);
770 if (method == MultiCoil)
771 table = table % NUMBER_CLASSES;
772 else if (method == HMMCoil)
773 table = table % NUMBER_TABLES;
776 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,1,0, bound);
777 if (zoomptr) Zoomer();
782 method_cycle(&preprocessor_method,NUMBER_PREPROCESSORS);
784 if (preprocessor_method == MultiCoil)
785 preproc_table = preproc_table% NUMBER_CLASSES;
786 else if (preprocessor_method == HMMCoil)
787 preproc_table = preproc_table%NUMBER_TABLES;
788 else preproc_table =0;
790 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,1, bound);
791 if (zoomptr) Zoomer();
795 if (preprocessor_method == MultiCoil)
796 preproc_table = (preproc_table +1) % NUMBER_CLASSES;
797 else if (preprocessor_method == HMMCoil)
798 preproc_table = (preproc_table +1) % NUMBER_TABLES;
799 else preproc_table = 0;
800 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,2, bound);
801 if (zoomptr) Zoomer();
805 if (method == MultiCoil)
806 table = (table +1) % NUMBER_CLASSES;
807 else if (method == HMMCoil)
808 table = (table +1) % NUMBER_TABLES;
810 ShowSeq(lib,multi_lib,m,b,m_single,b_single,mode,0,0, bound);
811 if (zoomptr) Zoomer();
816 if (report.xany.window == window) PrevSeq(fin, lib,multi_lib,m,b,
817 m_single,b_single,mode,bound);
825 if (report.xany.window == window)
826 NextSeq(fin,lib,multi_lib,m,b,m_single,b_single,mode,bound,flog,
838 if (XCheckTypedEvent(display,Expose,&report)) {
840 XDrawString(display,window,gcgray,500,17,"(initializing)",14);
846 struct itimerval value;
847 signal(SIGALRM,Alarm);
848 getitimer(ITIMER_REAL,&value);
849 value.it_value.tv_sec = 0;
850 value.it_value.tv_usec = 500;
851 value.it_interval.tv_sec = 0;
852 value.it_interval.tv_usec = 500;
853 setitimer(ITIMER_REAL,&value,NULL);
858 struct itimerval value;
859 getitimer(ITIMER_REAL,&value);
860 value.it_value.tv_sec = 0;
861 value.it_value.tv_usec = 0;
862 setitimer(ITIMER_REAL,&value,NULL);
863 XDrawImageString(display,window,gcgray,500,17," ",14);
864 XDrawImageString(display,window,gcgray,225,729,"press any key to continue",25);
870 if (!zoomptr) return;
871 PostoXY(zoomptr,&x,&y);
872 pts[0].x = x; pts[0].y = y;
873 pts[1].x = -5; pts[1].y = 5;
874 pts[2].x = 10; pts[2].y = 0;
875 XFillPolygon(display,window,gcdash,pts,3,Convex,CoordModePrevious);
876 XDrawLine(display,window,gcdash,x,y+5,x,y+8);
880 /* Find most likely residue, break ties arbitrarily */
885 for (i=0; i<sequence.seqlen; ++i)
886 if (scores[8*i+7] > like) {
887 like = scores[8*i+7];
898 XMapWindow(display, zoom);
899 /* Indicate where in the sequence we are */
900 sprintf(buf,"Position %d ",zoomptr);
901 XDrawImageString(display,zoom,gc,7,20,
903 for (i = -3; i<=3; ++i) {
904 if (zoomptr+i>0 && zoomptr+i<=sequence.seqlen)
905 buf[0] = numaa(sequence.seq[zoomptr-1+i]);
908 XDrawImageString(display,zoom,gc,159+9*i,20, buf,1);
910 XDrawLine(display,zoom,gc,159,22,168,22);
912 /* Indicate the probability of being in any given state */
913 sc = &(scores[8*(zoomptr-1)]);
916 if (sc[i] == sc[7]) {offset=i; break;} /* Used to print out maxscoring */
917 /* register as 'a' + offset. */
918 /* Offset 7 will be a flag that no particular offset scored max. */
920 /* Give likelihood */
921 sprintf(buf,"%5.1f%% chance",100*sc[7]);
922 XDrawImageString(display,zoom,gc,220,20, buf,strlen(buf));
925 fprintf(stderr,"\noffset= %d, reg = %c",all_offsets[table][zoomptr-1],
926 'a' + (zoomptr-1 + all_offsets[table][zoomptr-1])%POSNUM);
931 sprintf(buf,"register %c",'a'+offset);
933 sprintf(buf,"combo reg ");
934 else sprintf(buf," ");
935 XDrawImageString(display,zoom,gc,230,40, buf,strlen(buf));
937 for (i=0; i<7; ++i) {
938 sprintf(buf,"%c:%4.1f%% ",'a'+i,(100*sc[i]));
939 XDrawImageString(display,zoom,gc,15+92*(i&3),60+20*(i/4),buf,strlen(buf));
941 XDrawString(display,zoom,gc,10,110,"Click on viewgram or use arrow",30);
942 XDrawString(display,zoom,gc,10,127,"keys to view other positions.",29);
943 XDrawString(display,zoom,gc,10,144,"You may need to move this",25);
944 XDrawString(display,zoom,gc,10,161,"window to see the viewgram.",27);
954 /*********The following is for storing sequences. ***/
955 static char seqbuff[40000];
957 static Sequence seqs[30];
958 static int seqbufpos[30];
959 static int seqsizes[30];
962 { int size, pos, limsize;
964 size = seqsizes[sequp-seqnum];
965 pos = seqbufpos[sequp-seqnum];
967 sequence = seqs[sequp-seqnum]; /** Retrieve seq from memory. **/
970 limsize = (pos+size<=40000) ? size : 40000-pos;
971 bcopy(&(seqbuff[pos]),sequence.code,limsize);
972 bcopy(seqbuff,&(sequence.code[limsize]),size-limsize);
978 { int size, pos, limsize,i;
979 char *temp, *current;
981 bcopy(seqs,&(seqs[1]),29*sizeof(Sequence));
982 bcopy(seqbufpos,&(seqbufpos[1]),29*sizeof(int));
983 bcopy(seqsizes,&(seqsizes[1]),29*sizeof(int));
984 if (sequp-seqlow>=30) seqlow = sequp-29;
986 /* The 2 accounts for the end of string characters in .code and .title. */
987 size = 2 + strlen(sequence.code) + strlen(sequence.title) +
988 sequence.seqlen * (sequence.reg ? 2 : 1);
990 pos = seqbufpos[1] + seqsizes[1]; /* Number of characters in array */
992 /* if (pos>=40000) pos -= 40000;*/
993 /* limsize = (pos+size<=40000) ? size : 40000-pos; */
995 if (pos+size >=40000) pos=0; /* Array will be full so start at beginning. */
996 limsize = (pos+size<=40000) ? size : 40000-pos;
999 /* Save the "code+title+seq+reg" in something permanant, and wrap around
1002 /* David uses bcopy because wants to copy past end of string characters */
1004 bcopy(sequence.code,&(seqbuff[pos]),limsize);
1005 bcopy(&(sequence.code[limsize]),seqbuff,size-limsize);
1008 sequence.code=strcpy(&(seqbuff[pos]),sequence.code);
1010 strcpy(&(seqbuff[pos+strlen(sequence.code)+1]),sequence.title);
1014 &(seqbuff[pos+strlen(sequence.code) + strlen(sequence.title) + 2]);
1015 for(i=0; i<sequence.seqlen; i++)
1016 sequence.seq[i] = temp[i];
1021 sequence.reg=&(sequence.seq[sequence.seqlen]);
1022 for(i=0; i<sequence.seqlen; i++)
1023 sequence.reg[i] = temp[i];
1027 while ((seqlow<sequp) && (seqbufpos[sequp-seqlow] >= pos && seqbufpos[sequp-seqlow] < pos+limsize || seqbufpos[sequp-seqlow] < size-limsize))
1030 seqs[0].seqlen= sequence.seqlen;
1040 method_cycle(int *Method, int Number_of_Methods)
1043 *Method = (*Method+1)%Number_of_Methods;
1045 while ((*Method == MultiCoil && MultiCoil_Method == unavailable) ||
1046 (*Method == PairCoil && PairCoil_Method == unavailable) ||
1047 (*Method == NEWCOILS && NEWCOILS_Method == unavailable) ||
1048 (*Method == HMMCoil && HMMCoil_Method == unavailable) ||
1049 (*Method == ActualCoil && ActualCoil_Method == unavailable) ||
1050 (*Method == PairCoilDiff && PairCoilDiff_Method == unavailable)||
1051 (*Method == PairCoilDiffAvg && PairCoilDiffAvg_Method == unavailable) ||
1053 (*Method == ActualCoil) && (!sequence.reg) /* Don't do actual coil */
1054 ); /* if didn't input register */
1060 PrevSeq (FILE *fin, char lib[MAXFUNCTNUM],
1061 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1062 double *m, double *b,
1063 double *m_single, double *b_single,
1064 int mode, double bound)
1066 if (seqnum<=seqlow) {XBell(display,0); return;}
1069 NewSeq(mode,sequence);
1070 ShowSeq(lib, multi_lib,m, b, m_single, b_single, mode,1,1, bound);
1074 NextSeq (FILE *fin,char lib[MAXFUNCTNUM],
1075 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1076 double *m, double *b, double *m_single, double *b_single,
1077 int mode, double bound, FILE *flog, FILE *fout_coils,FILE *fout)
1079 char title[MAXLINE],code[MAXLINE];
1080 char seq[MAXSEQLEN], reg[MAXSEQLEN];
1082 int first_time_seq=0;
1085 if (seqnum < sequp) {
1091 sequence.seqlen = 0;
1094 sequence.title = title;
1095 sequence.code = code;
1098 if (!getseq2(fin,&sequence))
1099 {XBell(display,0); sequence= temp_seq; autoadvance=0; return;}
1107 NewSeq(mode,sequence);
1108 if (first_time_seq) {
1109 output_seq(lib,multi_lib,m,b,m_single, b_single,
1110 mode,log_bound,flog,fout_coils,fout,
1111 avg_max, main_method, main_preprocessor_method, main_table);
1117 ShowSeq(lib,multi_lib,m, b, m_single, b_single,mode,1,1, bound);
1122 /** This is stuff need to always do before rescore a sequence. **/
1123 NewSeq_nonX (int mode, Sequence sequence)
1125 if (mode & TST_MODE0)
1126 recalc_prob(sequence, 0, mode);
1127 if (mode & TST_MODE1)
1128 recalc_prob(sequence,1, mode);
1131 NewSeq (int mode, Sequence sequence)
1133 NewSeq_nonX (mode, sequence);
1135 XUnmapWindow(display, zoom);
1140 /* The likelihood line is mx +b */
1141 void output_seq(char lib[MAXFUNCTNUM],
1142 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1143 double *m,double *b,
1144 double *m_single, double *b_single,
1146 double log_bound,FILE *flog,
1150 int main_method, int main_preprocessor_method,
1154 double maxother=0; double maxscore;
1155 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1160 double junk_scores[MAXSEQLEN][POSNUM+1];
1163 fprintf(stderr,"\nIn outputseq");
1166 if (is_not_differentiator(main_method))
1170 all_scores = ScoreSeq (lib, multi_lib,m, b, m_single, b_single, mode,
1171 main_table, main_method, main_preprocessor_method,
1172 log_offset_to_use, &maxscore, 1, do_preproc, log_bound);
1176 add_to_total_likes(all_scores, total_gauss_like, sequence.seqlen,
1180 fprintf(stderr,"\nGot all_scores");
1182 scores = &all_scores[main_table* MAXSEQLEN*(POSNUM+1) ];
1184 /* This is all_scores[main_table + 0*number_tables][0][0] */
1185 /* so is score on main_table using library # 0 **/
1186 /* Similarly for preproces_like. **/
1188 preprocess_like = &all_preprocess_like[main_table][0][0];
1191 for (i=0; i<sequence.seqlen; i++)
1192 if (preprocess_like[i*(POSNUM+1) + 7] > maxother)
1193 maxother=preprocess_like[i*(POSNUM +1) + 7];
1194 /* This is preprocess_like[i][7] */
1197 switch (main_method) {
1199 if (mode & VER_MODE) {
1200 NEWCOILSScore(mode, sequence, &maxother,junk_scores,
1203 /* scores in junk_scores */
1204 log_output_ver(SC2SEQ | SCSTOCK, maxscore, maxother,
1205 maxscore,seqnum, sequence.title,
1206 sequence.code, flog);
1213 /*** For pos files coils, does a log file version of coil/seq score. **/
1214 output_pos_scores(fout_coils, sequence, all_scores, log_bound, mode,2,1);
1217 if (mode & PRN_MODE)
1220 /* If aren't dealing with coil scores, convert to weighted-avg **/
1221 /* coil_scores (by using average_score_over_coils3(....,0) **/
1222 /* so can output to log file. **/
1224 for (tablenum =0; tablenum <NUMBER_CLASSES; tablenum++) {
1225 average_score_over_coils3(sequence,
1226 &all_scores[(tablenum+NUMBER_CLASSES)* MAXSEQLEN*(POSNUM+1)],
1227 &all_scores[tablenum* MAXSEQLEN*(POSNUM+1)],
1228 &all_scores[(2+NUMBER_CLASSES)* MAXSEQLEN*(POSNUM+1)],
1229 /** Total Coil Like. **/
1231 all_offsets[tablenum],
1232 start_coil_at_reg_shift, log_bound,
1236 log_coil_output(all_scores, all_offsets,sequence, log_bound, flog,
1239 /** Now convert back for safety in case need to do an out file **/
1241 switch_back_to_res_score(log_bound, sequence,mode,all_scores);
1245 log_output2_prn(mode,maxscore, maxscore, scores, scores,
1247 seqnum, sequence, log_bound, log_bound, log_bound,
1255 if (mode & VER_MODE)
1256 log_output_ver(SC2SEQ , maxscore, 0,
1257 maxscore,seqnum, sequence.title,
1258 sequence.code, flog);
1261 if (mode & PRN_MODE)
1262 log_output2_prn(mode,maxscore, maxscore, scores, scores,sequence.seqlen,
1263 seqnum, sequence, log_bound, log_bound, log_bound,
1268 /* case PairCoilDiff: */
1269 case PairCoilDiffAvg:
1270 if (mode & PRN_MODE) {
1271 log_output2_prn(mode,maxscore, maxother, scores,preprocess_like,
1273 seqnum, sequence, .5, .5, log_bound,
1283 fprintf(stderr,"\n About to go into fout.");
1288 if ( (!(mode & RAW_OUT)) ||
1289 ( (main_method != MultiCoil) && (main_method != PairCoilDiff))) {
1290 if (main_method== MultiCoil)
1291 if (mode & WEB_OUT_MODE)
1292 web_output(sequence, fout, all_scores,all_offsets,2);
1293 /* Output table of score and seq and reg values */
1295 tuple_output2(sequence, mode, fout, all_scores,log_bound,2);
1296 /* Output dimer and trimer prob. in regions above bound and */
1297 /* in long enough windows only. **/
1299 else txt_output2(sequence, mode, fout, scores, log_bound);
1301 else { /* Want raw tuple_output for multi_coil and PairCoilDiff*/
1303 if (main_method == MultiCoil)
1304 if (log_offset_to_use==7) /* Combo offset just output max score */
1305 tuple_output2(sequence, mode, fout, /** for each dimension. */
1307 log_bound, number_score_dim);
1308 else /** Max reg for raw score can't really be determined. Instead */
1309 /* for coils, output scores for the known correct register, & */
1310 /* for pdb- output ALL possible 7 register scores (all negs). */
1311 tuple_output_for_max(sequence,mode,fout, all_scores, NUMBER_TABLES,
1312 log_bound, number_score_dim);
1319 int switch_back_to_res_score(double bound, Sequence sequence, int mode,
1320 double All_Scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1])
1325 for (tablenum =0; tablenum < NUMBER_CLASSES ; tablenum++)
1326 for (j =0; j <=POSNUM; j++)
1327 for (i=0; i<sequence.seqlen; i++) {
1328 All_Scores[tablenum][i][j] =
1329 All_Scores[tablenum+ NUMBER_CLASSES][i][j];
1330 if ((mode & ONLY_COILED_CLASSES) && (tablenum < (NUMBER_CLASSES-1))
1331 && (All_Scores[2+ NUMBER_CLASSES][i][j]>bound))
1332 All_Scores[tablenum][i][j]/= All_Scores[2+ NUMBER_CLASSES][i][j];
1337 int copy_offsets(char preproc_offsets[MAXSEQLEN],
1339 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1341 double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1346 for (tablenum =0; tablenum < number_dimens; tablenum++)
1347 for(i=0; i<sequence.seqlen; i++) {
1348 all_offsets[tablenum][i]=preproc_offsets[i];
1349 /*** Make sure the output score corresponds to this offset... **/
1350 /*** Also do it for backup residue scores in tablenum + NUMBER_CLASSES */
1351 if ((offset_to_use ==7) || (offset_to_use == -1)) {
1352 all_scores[tablenum][i][7] =
1353 all_scores[tablenum][i][(i+all_offsets[tablenum][i]) %POSNUM];
1354 all_scores[tablenum+ NUMBER_CLASSES][i][7] =
1355 all_scores[tablenum+NUMBER_CLASSES][i]
1356 [(i+all_offsets[tablenum][i]) %POSNUM]; }
1360 int get_offsets(double all_likes[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1361 Sequence sequence, char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1365 int got_offset, max_offset;
1366 double max_off_score;
1369 for (tablenum =0; tablenum < number_dimens; tablenum++)
1370 for(i=0; i<sequence.seqlen; i++) {
1371 offsets[tablenum][i]=-1;
1373 max_off_score=-HUGE_VAL;
1374 for (j=0; j<POSNUM; j++) {
1375 if (all_likes[tablenum][i][j] == all_likes[tablenum][i][7]) {
1376 offsets[tablenum][i] = (j-i) % POSNUM;
1377 if (offsets[tablenum][i]<0)
1378 offsets[tablenum][i] += POSNUM;
1381 if (all_likes[tablenum][i][j] > max_off_score) {
1382 max_off_score=all_likes[tablenum][i][j];
1383 max_offset=(j-i) % POSNUM;
1384 if (max_offset <0) max_offset +=POSNUM;
1387 if (!got_offset) offsets[tablenum][i] = max_offset;
1389 /* This is a new thing so that if an offset change is not real */
1390 /* we keep the old offset. **********/
1391 if ( (i>0) && (all_likes[tablenum][i][(i+offsets[tablenum][i-1]) %POSNUM]
1392 == all_likes[tablenum][i][(i+offsets[tablenum][i])%POSNUM]) )
1393 offsets[tablenum][i] = offsets[tablenum][i-1];
1405 void preprocessor_score(char lib[MAXFUNCTNUM],
1406 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1407 double *m, double *b,
1408 double *m_single, double *b_single, int mode,
1409 int Preprocessor_Method, int Offset_to_Use,
1410 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1412 char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1416 char offsets[MAXSEQLEN];
1420 int dimension, num_score_dim;
1423 for (tablenum =0; tablenum < NUMBER_TABLES; tablenum++) {
1424 switch_tables(tablenum); /* Set prob. to be for that table. */
1426 switch(Preprocessor_Method) {
1428 if (combine_dist[tablenum])
1429 functnum = num_dist[tablenum]; /* For use in scscore.c **/
1432 num_score_dim= num_table_dimensions(num_dist[tablenum],
1433 combine_dist[tablenum]);
1435 for (dist=0; dist<num_score_dim; dist++) {
1436 if (num_dist[tablenum]>1) { /* Hack: else don't access correctly */
1437 local_lib = &multi_lib[tablenum][dist];
1439 else local_lib = multi_lib[tablenum];
1440 dimension = compute_dimension2(tablenum, dist, num_dist,
1442 MultiCoilDimensionScore(mode, sequence, local_lib,
1443 &maxscore,tablenum, offsets,
1444 all_preprocess_like[dimension],
1450 if (mode & PAIRCOIL_PAIRS) {
1451 functnum = pair_functnum; /* For use in scscore.c **/
1452 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1453 &maxscore,tablenum, offsets,
1454 all_preprocess_like[tablenum],
1457 else /* Use singles probabilities... for likelihood use line 0 */
1458 SingleCoilScore(mode,sequence,m_single[tablenum],
1460 &maxscore, tablenum,
1461 all_preprocess_like[tablenum], Offset_to_Use,
1466 NEWCOILSScore(mode, sequence, &maxscore, all_preprocess_like[tablenum], Offset_to_Use);
1470 HMMScore(sequence.seq, sequence.seqlen, hmm, all_preprocess_like[tablenum],
1471 tablenum, Offset_to_Use);
1475 /* if (mode & POS_MODE) */ /* Score the coils in the posfile based on */
1476 /* their having 100% likelihood. */
1477 ActualCoils(sequence, all_preprocess_like[tablenum], Offset_to_Use,1);
1481 fprintf(stderr,"\nBad preprocessor, using PairCoil.\n");
1482 functnum = pair_functnum; /* For use in scscore.c **/
1483 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1484 &maxscore,tablenum, offsets,
1485 all_preprocess_like[tablenum],
1492 if (Preprocessor_Method !=MultiCoil)
1493 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1496 /***** Get copy of res scores in tablenum and tablenum+ NUMBER_CLASSES **/
1497 switch_gauss_param(Offset_to_Use, Preprocessor_Method);
1499 convert_raw_scores_to_gauss_prob_like2(all_preprocess_like,sequence.seqlen,
1500 NUMBER_TABLES,means_submatrix,inverse_covars,
1501 determinants, NUMBER_CLASSES, init_class_prob,
1503 mode & ONLY_COILED_CLASSES,
1504 bound,Offset_to_Use,1);
1506 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1509 if (Coil_Score) { /* Value of 1 means do average, 2 means do max. */
1510 /*** Need to save residue scores, since the coil scores **/
1511 /*** will replace them, but old residue scores might be needed **/
1512 /*** again in tablenum+ NUMBER_CLASSES. ****/
1513 for (tablenum =0; tablenum <3; tablenum++) {
1514 average_score_over_coils3(sequence,
1515 all_preprocess_like[tablenum+NUMBER_CLASSES],
1516 all_preprocess_like[tablenum],
1517 all_preprocess_like[2+NUMBER_CLASSES],
1518 /** Total Coil Like. **/
1520 all_preproc_offsets[tablenum],
1521 start_coil_at_reg_shift, bound,
1535 double *ScoreSeq (char lib[MAXFUNCTNUM],
1536 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1537 double *m_paircoil, double *b_paircoil,
1538 double *m_single, double *b_single,
1539 int mode, int Table, int Method,
1540 int Preprocessor_Method, int Offset_to_Use,
1542 int rescore_seq, int rescore_preproc, double bound)
1544 /** rescore_preproc == 1 means need to rescore. If it or rescore_seq are */
1545 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
1546 /** -2 means need to change from MultiCoil coil score back to res score. */
1547 /** If it is another non-zero */
1548 /** value, don't need to rescore it, but do need to rescore anything that */
1549 /** uses it as a preproc. */
1551 double max_paircoil_score;
1552 static char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1553 extern char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1554 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1555 static double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1562 int num_score_dim, dimension;
1563 int local_number_tables;
1565 double *preproc_res_like;
1566 int scored_already =0;
1568 /**************************************************************************/
1569 /********* PreprocLike is used in the 2-3 stranded differentiators. ****/
1572 if ( (rescore_preproc==1) &&
1573 (!ftotal_like)) { /* A hack to save computation time. */
1574 preprocessor_score(lib, multi_lib,
1575 m_paircoil, b_paircoil, m_single, b_single,
1576 mode, Preprocessor_Method,
1578 all_preprocess_like, sequence, all_preproc_offsets,
1581 /*****************************************************************************/
1582 if (Preprocessor_Method == MultiCoil) {
1583 if ( ( rescore_preproc==-1) && Coil_Score) {
1584 /* Coil_Score value of 1 means do average, 2 means do max. */
1585 /*** Need to save residue scores, since the coil scores **/
1586 /*** will replace them, but old residue scores might be needed **/
1587 /*** again in tablenum+ NUMBER_CLASSES. ****/
1588 for (tablenum =0; tablenum < NUMBER_CLASSES; tablenum++)
1589 average_score_over_coils3(sequence,
1590 all_preprocess_like[tablenum+ NUMBER_CLASSES],
1591 all_preprocess_like[tablenum],
1592 all_preprocess_like[2+NUMBER_CLASSES],
1593 /** Total Coil Like. **/
1595 all_preproc_offsets[tablenum],
1596 start_coil_at_reg_shift, bound,
1598 if (mode & ONLY_COILED_CLASSES) /** Convert to olig state. **/
1599 type_class_score_convert(sequence,all_preprocess_like,bound,
1600 mode & ONLY_COILED_CLASSES);
1602 else if (rescore_preproc == -2) /** Switch back to res score.*/
1603 switch_back_to_res_score(bound, sequence,mode,all_preprocess_like);
1606 /****************************************************************************/
1612 if (Method == PairCoilDiff)
1613 local_number_tables = 1;
1614 else local_number_tables= NUMBER_TABLES;
1617 /*** Clean up this conditional, since really the only thing that might ***/
1618 /*** change for paircoildiffer when change preproc is the offset to use **/
1619 /*** (i.e. all_scores[table][i][7]. ***/
1621 if ((rescore_seq>0) ||
1622 (!is_not_differentiator(Method) && rescore_preproc) ) {
1624 for (tablenum=0; tablenum<local_number_tables; tablenum++) {
1625 switch_tables(tablenum); /* Switch prob to right table. */
1628 maxsc=0; /* Initialization. */
1634 functnum = pair_functnum; /* For use in scscore.c **/
1635 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1636 b_paircoil[tablenum],
1637 &maxsc,tablenum, all_offsets[tablenum],all_scores[tablenum],
1643 if (combine_dist[tablenum])
1644 functnum = num_dist[tablenum]; /* For use in scscore.c **/
1647 num_score_dim= num_table_dimensions(num_dist[tablenum],
1648 combine_dist[tablenum]);
1650 for (dist=0; dist<num_score_dim; dist++) {
1651 if (num_dist[tablenum] >1)
1652 local_lib = &multi_lib[tablenum][dist];
1653 else local_lib = multi_lib[tablenum];
1655 dimension = compute_dimension2(tablenum, dist,
1656 num_dist, combine_dist);
1657 MultiCoilDimensionScore(mode, sequence, local_lib, &maxsc,
1658 tablenum,all_offsets[dimension],
1659 all_scores[dimension],
1666 if (mode & PAIRCOIL_PAIRS) {
1667 functnum = pair_functnum; /* For use in scscore.c **/
1668 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1669 b_paircoil[tablenum], &maxsc,
1670 tablenum,all_offsets[tablenum],
1671 all_scores[tablenum],
1672 Offset_to_Use, (mode & RAW_OUT));
1674 else /*Use singles probabilities... for likelihood use line number 0 */
1675 SingleCoilScore(mode,sequence,m_single[tablenum],
1678 all_scores[tablenum], Offset_to_Use,
1682 NEWCOILSScore(mode,sequence,&maxsc, all_scores[tablenum],
1687 HMMScore(sequence.seq, sequence.seqlen, hmm, all_scores[tablenum],
1688 tablenum, Offset_to_Use);
1693 ActualCoils(sequence, all_scores[tablenum], Offset_to_Use,0);
1696 if (tablenum == Table) *maxscore = maxsc;
1703 /*********************************************************************/
1704 /*********************************************************************/
1705 /********* This ends the scoring along all and all dimensions ********/
1706 /********* Now convert the raw scores that haven't been converted ****/
1707 /********* yet into likelihoods (i.e. the multidimensional scores.****/
1710 /***** Note: Can use same "init_class_prob" here, because since set number **/
1711 /*** of classes to 2, the 3rd class (Pdb-) will be ignored, so just get ****/
1712 /*** ratio of two to three stranded likelihoods, BOTH weighted by their ****/
1713 /*** init probaabiities. ********/
1715 if (Method != MultiCoil)
1716 get_offsets(all_scores, sequence,all_offsets,NUMBER_TABLES);
1717 else if (!(mode & RAW_OUT)) {
1719 switch_gauss_param(Offset_to_Use, Method);
1723 /**********Do conversions on MultiCoil dimension scores. ******/
1725 if (Method == MultiCoil) {
1726 convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
1727 NUMBER_TABLES, means_submatrix,inverse_covars,
1728 determinants, NUMBER_CLASSES,init_class_prob,
1730 mode & ONLY_COILED_CLASSES, bound,
1732 get_offsets(all_scores, sequence,all_offsets,3);
1735 if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
1736 for (tablenum=0; tablenum<NUMBER_CLASSES; tablenum++)
1737 average_score_over_coils3(sequence,
1738 all_scores[tablenum+NUMBER_CLASSES],
1739 all_scores[tablenum],
1740 all_scores[2+NUMBER_CLASSES],
1741 /** Total Coil Like. **/
1743 all_offsets[tablenum],
1744 start_coil_at_reg_shift, bound,
1750 } /* Ends if rescore_seq */
1753 /*****************************************************************************/
1754 if ( (Method == MultiCoil) || (Method == PairCoilDiff)){
1756 if (Method == MultiCoil)
1757 if ( (rescore_seq ==1 ) || /** Signals rescored MultiCoil **/
1758 (rescore_seq==-1)) /*Signals bound changed, need to redo seq_score */
1759 get_seq_scores(seq_scores, sequence, all_scores[0+ NUMBER_CLASSES],
1760 all_scores[1+ NUMBER_CLASSES],bound);
1763 if ( (rescore_seq==-1) && Coil_Score && !scored_already) {
1764 if (Method == MultiCoil) {
1765 for (tablenum=0; tablenum<NUMBER_CLASSES; tablenum++) {
1766 average_score_over_coils3(sequence,
1767 all_scores[tablenum+NUMBER_CLASSES],
1768 all_scores[tablenum],
1769 all_scores[2+NUMBER_CLASSES],
1770 /** Total Coil Like. **/
1772 all_offsets[tablenum],
1773 start_coil_at_reg_shift, bound,
1775 if (mode & ONLY_COILED_CLASSES) /** Convert to olig. state. **/
1776 type_class_score_convert(sequence,all_scores,bound,
1777 mode & ONLY_COILED_CLASSES);
1784 else if (rescore_seq == -2) /** Switch back to reg score.*/
1785 switch_back_to_res_score(bound, sequence,mode,all_scores);
1787 /****************************************************************************/
1790 return (double*)all_scores;
1797 ShowSeq(char lib[MAXFUNCTNUM],
1798 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1799 double *m, double *b,
1800 double *m_single, double *b_single,
1801 int mode, int rescore_sequence, int rescore_preproc, double bound)
1806 if ( (rescore_sequence) || (rescore_preproc))
1807 all_scores = ScoreSeq(lib, multi_lib, m, b,m_single, b_single,
1809 method, preprocessor_method,
1810 offset_to_use, &maxscore,
1811 rescore_sequence, rescore_preproc, bound);
1814 scores = &all_scores[table * MAXSEQLEN*(POSNUM+1)];
1815 /* This is all_scores[main_table][0][0] */
1816 preprocess_like = &all_preprocess_like[preproc_table][0][0];
1818 ShowScores(scores,&sequence,NULL, bound,mode);
1821 /*---------------------------------------------------------------------------*/
1826 OpenScreen(int mode, double log_bound, FILE *fin, FILE *flog,
1829 double *m, double *b, double *m_single, double *b_single,
1830 char lib[MAXFUNCTNUM],
1831 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1833 int main_preprocessor_method,
1838 XSetWindowAttributes attr;
1840 static unsigned char dashl[2] = {1,1};
1841 static unsigned char dashd[2] = {2,2};
1846 display = XOpenDisplay("");
1848 /* Check mode to see if want GUI. */
1849 if (!display || (mode & NO_GUI)) {
1850 /** COmplete logfile, etc. if can't open xdisplay. */
1851 if (!(mode & NO_GUI))
1852 fprintf(stderr,"Sorry, bad display.\nMaybe your DISPLAY environment variable isn't set properly,\nor you haven't xhosted this machine.\n");
1853 if ( fout || ( (mode & PRN_MODE) || (mode & VER_MODE)) ) {
1854 fprintf(stderr,"\nThe log file will still be written by:\n\n");
1855 for (i = 0; motd[i][0]; ++i) /* Print out title. */
1856 fprintf(stderr," %s\n",motd[i]);
1858 finish_log(mode,log_bound,fin,flog,fout_coils,
1860 m,b, m_single, b_single,
1861 lib,multi_lib, main_method, main_table,
1862 main_preprocessor_method, seqcnt, avg_max);
1866 if (flog) sclose(flog);
1867 if (fout) sclose(fout);
1868 if (ftotal_like) sclose(ftotal_like);
1874 font = XLoadFont(display,"9x15");
1875 window = XCreateSimpleWindow(display, RootWindow(display, 0),
1877 WIDTH-80, HEIGHT/4 +75,
1879 BlackPixel(display, 0),
1880 WhitePixel(display, 0));
1882 help = XCreateSimpleWindow(display, RootWindow(display, 0),
1884 WIDTH / 5 * 4 +60, HEIGHT / 5 * 3 + 400,
1886 BlackPixel(display, 0),
1887 WhitePixel(display, 0));
1888 helpb = XCreateSimpleWindow(display, window,
1890 BlackPixel(display, 0), WhitePixel(display, 0));
1891 next = XCreateSimpleWindow(display, window,
1893 BlackPixel(display, 0), WhitePixel(display, 0));
1894 prev = XCreateSimpleWindow(display, window,
1896 BlackPixel(display, 0), WhitePixel(display, 0));
1897 autob = XCreateSimpleWindow(display, window,
1899 BlackPixel(display, 0), WhitePixel(display, 0));
1900 zoomb = XCreateSimpleWindow(display, window,
1902 BlackPixel(display, 0), WhitePixel(display, 0));
1903 printb = XCreateSimpleWindow(display, window,
1905 BlackPixel(display, 0), WhitePixel(display, 0));
1906 tableb=XCreateSimpleWindow(display, window,
1908 BlackPixel(display, 0), WhitePixel(display, 0));
1910 methodb=XCreateSimpleWindow(display, window,
1912 BlackPixel(display, 0), WhitePixel(display, 0));
1914 paramb=XCreateSimpleWindow(display, window,
1916 BlackPixel(display, 0), WhitePixel(display, 0));
1917 quit = XCreateSimpleWindow(display, window,
1919 BlackPixel(display, 0), WhitePixel(display, 0));
1921 /** To get size: **/
1922 /** Printing " trimer x.xx" takes 12 char = 12*7 =84 **/
1923 /** Each line takes 15 of y coord. giving 3*15 =45 + 10 for blank line **/
1924 seqscoreb=XCreateSimpleWindow(display, window,
1925 XCOORD_SEQSCORES, YCOORD_SEQSCORES,
1927 BlackPixel(display, 0), WhitePixel(display, 0));
1929 zoom = XCreateSimpleWindow(display, RootWindow(display, 0),
1933 BlackPixel(display, 0),
1934 WhitePixel(display, 0));
1935 XSelectInput(display, helpb, ButtonPressMask);
1936 XSelectInput(display, next, ButtonPressMask);
1937 XSelectInput(display, prev, ButtonPressMask);
1938 XSelectInput(display, autob, ButtonPressMask);
1939 XSelectInput(display, zoomb, ButtonPressMask);
1940 XSelectInput(display, printb, ButtonPressMask);
1941 /**** XSelectInput(display, methodb, ButtonPressMask); ***/
1942 XSelectInput(display, tableb, ButtonPressMask);
1943 XSelectInput(display, paramb, ButtonPressMask);
1945 XSelectInput(display, quit, ButtonPressMask);
1946 XSelectInput(display, window, ExposureMask);
1947 XSelectInput(display, zoom, ExposureMask | KeyPressMask | ButtonPressMask);
1948 XSelectInput(display, help, ExposureMask | KeyPressMask | ButtonPressMask);
1949 attr.bit_gravity = NorthWestGravity;
1950 XChangeWindowAttributes(display,window,CWBitGravity,&attr);
1951 XMapWindow(display, window);
1954 values.foreground = BlackPixel(display, 0);
1955 values.background = WhitePixel(display, 0);
1957 gc = XCreateGC(display, window, GCForeground | GCBackground | GCFont, &values);
1958 values.line_style = LineDoubleDash;
1959 values.font = XLoadFont(display,"6x13");
1960 gcgray = XCreateGC(display, window, GCForeground | GCBackground
1961 | GCLineStyle | GCFont, &values);
1962 XSetDashes(display, gcgray, 1, dashl, 2);
1963 values.function = GXxor;
1964 values.foreground = 1;
1965 gcdash = XCreateGC(display, window, GCForeground
1966 | GCLineStyle | GCFunction, &values);
1967 XSetDashes(display, gcdash, 0, dashd, 2);
1971 icon = XCreateBitmapFromData(display,window,bart_bits,bart_width,bart_height);*/
1972 /* XSetStandardProperties(display,window,"Coiled Coil Finder","xcoil",icon,NULL,0,NULL);*/
1973 /* XSetStandardProperties(display,zoom,"XCoil Zoomer","xcoil",icon,NULL,0,NULL);*/
1974 /* XSetStandardProperties(display,help,"XCoil Help","xcoil",icon,NULL,0,NULL);*/
1978 PSDefs (FILE *ps, int seqlen)
1983 int nominallength= 300; /* Default value for number residues per line. */
1985 if ( (seqlen%nominallength <= 50) && (seqlen >= nominallength)) {
1986 nominallength = 350;
1989 /*************************/
1990 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
1991 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
1992 nominallength = seqlen;
1994 nominallength = ps_res_per_line; /* How many res per line. */
1995 /************************/
1997 if (seqlen ==0) number_lines=1;
1999 number_lines = (int)(seqlen-1)/nominallength +1;
2000 /* nominallength amino acids per line. */
2001 /* C rounds down and we want 0-nominallength to go to 1 */
2002 /* nominallength+1 to 2*noiminallength to go to 2.. */
2007 /*****Scale and translate to fit the page if want it to rotate page. */
2008 /* unscaled there are 9 lines per page. */
2009 /* .89 scaling if 300 res per line.
2010 xscaler = .89 * 300/nominallength;
2011 fputs("90 rotate\n",ps);
2012 fputs("78 -60 translate\n",ps);
2013 if (number_lines > 9) yscaler = (double) 9/number_lines;
2014 fprintf(ps,"%lf %lf scale\n", xscaler, yscaler);
2017 /*****Scale and translate to fit the page if don't want to rotate page. */
2018 /* unscaled there can fit 12 lines per page, so 3600. */
2019 xscaler = .75 * 300/nominallength; /* .75 scaling if 300 res per line. */
2020 fputs("25 740 translate\n",ps);
2021 if (number_lines > 12) yscaler = (double) 12/number_lines;
2022 fprintf(ps,"%lf %lf scale\n", xscaler,yscaler);
2026 fputs("/black {0 setgray} def\n",ps);
2027 fputs("/gray {0.6 setgray} def\n",ps);
2028 fputs("/tt /Courier findfont 15 scalefont def\n",ps);
2029 fputs("/rm /Times-Roman findfont 12 scalefont def\n",ps);
2030 fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
2031 fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
2032 fputs("/m {neg moveto} def\n",ps);
2033 fputs("/l {neg lineto stroke} def\n",ps);
2034 fputs("/rl {neg rlineto} def\n",ps);
2035 fputs("/rl_stroke {neg rlineto stroke} def\n",ps);
2036 fputs("/bar {dup 0.0019 lt {pop} {0 exch rlineto stroke} ifelse} def\n",ps);
2039 /*** static int xmin = 60; ***** Made it global variable *****/
2040 /*Originally was 29 */
2041 /** xmin determines where in window start seq. **/
2042 static int nominallength;
2045 static int ydelta=57;
2046 /* If click at x,y, determine which residue position it corresponds to.
2049 XYtoPos (int x, int y)
2052 if (x<xmin) return 0;
2054 if (i>= nominallength) return 0;
2055 for (y -= ybase+ydelta; y>0; y-=ydelta, i+=nominallength);
2056 if (y<-30) return 0;
2057 if (i>=length) return 0;
2060 PostoXY (int i, int *x, int *y)
2063 *x = xmin + 2*(i%nominallength);
2064 *y = ybase+ydelta + ydelta*(i/nominallength);
2066 psstr (char *str, int len, FILE *ps)
2069 for (i=0; i<len; ++i)
2081 /*** ybase determines where top of the seq score box should be so don't ***/
2082 /*** print over the sequenc title. **/
2083 draw_seqscore_box(double bound, double seq_scores[MAX_CLASS_NUMBER],
2084 Window *seqscoreb, FILE *ps, int ybase)
2087 /* XWindowChanges values; */
2091 XClearWindow(display,*seqscoreb);
2093 /** Relocate seqscoreb base on ybase**/
2094 /**This is complex way to move it.****
2095 values.y= YCOORD_SEQSCORES+ybase-30;
2096 XConfigureWindow(display,*seqscoreb,CWY,&values);
2099 XMoveWindow(display,*seqscoreb,XCOORD_SEQSCORES, YCOORD_SEQSCORES+ybase-30);
2103 /** 3 means print 3 into the box horiz, 10 means print 10 into box vert? */
2105 XMapWindow(display, *seqscoreb);
2106 XDrawString(display,*seqscoreb,gcgray,3,10,"Seq Prob:",9);
2108 if (ps) { /* Draw the box. **/
2109 fprintf(ps,"%d %d m\n",10, 50+ybase-30);
2110 fprintf(ps,"0 40 rl\n");
2111 fprintf(ps,"79 0 rl\n");
2112 fprintf(ps,"0 -40 rl\n");
2113 fprintf(ps,"-79 0 rl_stroke\n");
2116 if (seq_scores[0] + seq_scores[1] + seq_scores[2] > 0) { /* Signals above */
2117 sprintf(buf," Dimer %5.2lf", seq_scores[0]); /* bound. */
2118 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2119 sprintf(buf," Trimer %5.2lf", seq_scores[1]);
2120 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2123 fprintf(ps,"rm setfont ");
2124 fprintf(ps,"%d %d m (Seq Prob for) show\n", 15,60 +ybase-30);
2126 fprintf(ps,"%d %d m (bound %5.2lf:) show\n", 15,70 +ybase-30,bound);
2128 fprintf(ps,"%d %d m ( Dimer %5.2lf) show\n",15,80 +ybase-30,
2130 fprintf(ps,"%d %d m ( Trimer %5.2lf) show\n",15,90 +ybase-30,
2137 sprintf(buf," No Residue ");
2138 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2139 sprintf(buf," above %5.2lf ", bound);
2140 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2143 fprintf(ps,"rm setfont ");
2144 fprintf(ps,"%d %d m (Seq Prob:) show\n", 15,60 +ybase-30);
2146 fprintf(ps,"%d %d m ( No Residue) show\n", 15, 75 +ybase-30);
2147 fprintf(ps,"%d %d m ( above %5.2lf) show\n",15,85 +ybase-30, bound);
2155 draw_bound_box(double bound, Window paramb)
2159 XMapWindow(display, paramb);
2160 XDrawString(display,paramb,gcgray,3,10,"bound",5);
2161 sprintf(buf,"%5.2lf", bound);
2162 XDrawImageString(display,paramb,gcgray,3,25,buf,strlen(buf));
2165 ShowScores (double *scores, Sequence *sequence, FILE *ps, double bound,
2168 int i,q, x, y0, xmax;
2173 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2174 char *MethodTitle = MethodTitleArray;
2175 char PreprocTitleArray[500]; /* For printing the preproc.name to ps file. */
2176 char *PreprocTitle = PreprocTitleArray;
2179 MethodTitleArray[0] = '\0'; /* Initialization. */
2180 PreprocTitleArray[0]='\0';
2182 nominallength = 300;
2183 length = sequence->seqlen;
2185 /* This makes display look nice (so don't have a line with < 50 residues. */
2186 if ((length%nominallength <= 50) && (length>= nominallength)){
2187 nominallength = 350;
2188 /* if (ps) fputs("-50 0 translate\n",ps); */
2192 /*************************/
2193 if (ps) /* Do this stuff only for postscript printout. */
2194 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
2195 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
2196 nominallength = sequence->seqlen;
2198 nominallength = ps_res_per_line; /* How many res per line. */
2199 /************************/
2202 charperline = (xmin+2*nominallength-5)/9; /*Each char is "9" units witdth */
2203 rows = (length-1)/nominallength+1; /* Number of rows to display. */
2205 xmax = 2*length+xmin; /* The max x position in current row. */
2207 xmax = 2*nominallength+xmin;
2208 width = xmin+2*nominallength+ WIDTH_ADD_ON; /* Width of window to display. */
2209 /* The +WIDTH_ADD_ON used to be +11, but */
2210 /* was cutting of some numbers. */
2211 x = strlen(sequence->title);
2212 XClearWindow(display,window);
2214 /**********This is where window size is adjusted to fit the sequence. ***/
2215 /********** using width and number of rows. **************************/
2217 XResizeWindow(display,window,width,
2218 38+15*((x+charperline-1)/charperline) +ydelta*rows);
2220 /*************************************************************************/
2221 XMapWindow(display, helpb);
2222 XMapWindow(display, next);
2223 XMapWindow(display, prev);
2224 XMapWindow(display, autob);
2225 XMapWindow(display, zoomb);
2226 XMapWindow(display, printb);
2227 XMapWindow(display, tableb);
2229 XMapWindow(display, methodb);
2230 XMapWindow(display, paramb); Now done in draw_bound_box()
2232 XMapWindow(display, quit);
2235 while (XCheckTypedEvent(display,Expose,&report));
2237 XDrawString(display,helpb,gcgray,3,10,"help",4);
2238 XDrawString(display,next,gcgray,3,10,"next",4);
2239 XDrawString(display,prev,gcgray,3,10,"prev",4);
2240 XDrawString(display,autob,gcgray,3,10,"auto",4);
2241 XDrawString(display,zoomb,gcgray,3,10,"zoom",4);
2242 XDrawString(display,printb,gcgray,3,10,"print",5);
2243 /* XDrawString(display,methodb,gcgray,3,10,"method",6); */
2244 XDrawString(display,tableb,gcgray,3,10,"table",5);
2246 draw_bound_box(bound, paramb);
2247 /** draw_seqscore_box(bound, seq_scores,?? &seqscoreb, ps, ybase);
2248 MOVED TO LATER (AFTER ybase gets adjusted at end of this function). ****/
2251 XDrawString(display,quit,gcgray,3,10,"quit",4);
2254 /************This section prints out the method name..... ********************/
2255 /**** Add 9 to xcoord for each character. *****/
2256 XDrawString(display,window,gc, XCOORD_SCORED_BY,15,"Scored by ",10);
2258 make_methodtitle(MethodTitle, method, table, mode,0);
2259 XDrawString(display,window,gc,XCOORD_SCORED_BY +90,15,MethodTitle,
2260 strlen(MethodTitle));
2262 XDrawString(display,window,gc, XCOORD_SCORED_BY,30,"Filter by ",10);
2263 make_methodtitle(PreprocTitle, preprocessor_method, preproc_table,
2265 XDrawString(display,window,gc,XCOORD_SCORED_BY + 90,30,
2266 PreprocTitle,strlen(PreprocTitle));
2268 /***************************************************************************/
2270 XDrawString(display,window,gc, 5,30,sequence->code,strlen(sequence->code));
2271 for (i=0,ybase=30; x > i + charperline; i += charperline)
2272 XDrawString(display,window,gc, 5,ybase+=15,&(sequence->title[i]),charperline);
2273 XDrawString(display,window,gc, 5,ybase+15,&(sequence->title[i]),x-i);
2276 fprintf(ps,"tt setfont 5 15 m (%s) show\n",sequence->code);
2277 /* fprintf(ps,"175 15 m (Viewgram by XCoil) show\n"); */
2279 fprintf(ps,"430 15 m (Scored by %s) show\n",methodname[method]);
2281 fprintf(ps,"180 15 m (Scored by %s) show\n",MethodTitle);
2282 fprintf(ps,"180 26 m (Filter by %s) show\n",PreprocTitle);
2285 for (i=0,ybase=30; x > i + charperline; i += charperline) {
2286 fprintf(ps,"5 %d m (",ybase+=15);
2287 psstr(&(sequence->title[i]),charperline,ps);
2288 fprintf(ps,") show\n");
2290 fprintf(ps,"5 %d m (",ybase+15);
2291 psstr(&(sequence->title[i]),x-i,ps);
2292 fprintf(ps,") show\n");
2293 fputs("rm setfont 2 setlinewidth\n",ps);
2295 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2297 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2298 y = y0-30*scores[8*i+7];
2299 XDrawLine(display, window, gc, x, y0, x, (int)(y+0.5));
2300 XDrawLine(display, window, gc, x+1, y0, x+1, (int)(y+0.5));
2302 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2305 /********************* If in pos mode draws in the correct likelihoods every
2306 now and then? ******************
2307 if (sequence->reg && method==HMMCoil) {
2310 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2311 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2312 if (sequence->reg[i] != '.') {
2313 double y = y0-30*scores[8*i+sequence->reg[i]];
2314 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2315 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2317 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2320 if (ps) fputs("black\n",ps);
2322 **********************/
2324 /**********Now draw in the PreProcScore in gray every other line. ****/
2328 for (i=0,y0=ybase; i<length; i+=2,x+=4) { /* Normally is ++i,x+=4, but */
2329 double y; /* we do every other line here. */
2330 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2332 y = y0-30*preprocess_like[i*(POSNUM +1) + 7];
2333 /* This is preprocess_like[i][7]. */
2334 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2335 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2337 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2343 if (ps) fputs("black\n",ps);
2345 /*************************************************************************/
2347 if (ps) fputs("0 setlinewidth\n",ps);
2348 for (i=9,y0=ybase; i<length; i+=10) {
2349 if (!((i-9)%nominallength)) {x=xmin-1; y0 += ydelta;}
2351 XDrawLine(display, window, gc, x,y0,x,y0+2); /* Draw ticks every */
2352 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+2); /* 10 residues */
2355 for (i=49,y0=ybase; i<length; i+=50) {
2357 if (!((i-49)%nominallength)) {x=xmin-1; y0 += ydelta;}
2359 XDrawLine(display, window, gc, x,y0,x,y0+5);
2360 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+5);
2361 sprintf(buff,"%d",i+1);
2362 XDrawString(display, window, gcgray, x-7, y0+18, buff, strlen(buff));
2363 if (ps) fprintf(ps,"%d %d m (%d) cshow\n",x,y0+18,i+1);
2365 for (i=1,y0=ybase+ydelta; i<=rows; ++i,y0+= ydelta) {
2367 xmax = 2*(length-(rows-1)*nominallength)+xmin; /* Max x position in */
2369 XDrawString(display,window,gcgray,xmin-25,y0-26,"100%",4);
2370 XDrawString(display,window,gcgray,xmin-19,y0-11,"50%",3);
2371 XDrawString(display,window,gcgray,xmin-13,y0+4,"0%",2);
2372 XDrawLine(display,window,gcgray,xmin,y0-30,xmax-1,y0-30);
2373 XDrawLine(display,window,gcgray,xmin,y0-15,xmax-1,y0-15);
2375 fprintf(ps,"%d %d m (100%%) rshow\n",xmin-1,y0-26);
2376 fprintf(ps,"%d %d m (50%%) rshow\n",xmin-1,y0-11);
2377 fprintf(ps,"%d %d m (0%%) rshow\n",xmin-1,y0+4);
2378 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-30,xmax,y0-30);
2379 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-15,xmax,y0-15);
2380 fprintf(ps,"%d %d m %d %d l\n",xmin,y0,xmax,y0);
2385 draw_seqscore_box(bound, seq_scores, &seqscoreb, ps, ybase);
2387 if (zoomptr) ZoomPnt();
2390 fputs("showpage\n",ps);
2399 make_methodtitle(char *MethodTitle, int Method, int which_table,
2400 int mode, int for_logfile)
2407 if (which_table ==0) like2or3 =2;
2411 /************This section prints out the method name..... ********************/
2412 /**** Add 9 to xcoord for each character. *****/
2416 if ((Method == PairCoilDiff) || (Method == PairCoilDiffAvg)) {
2417 /* Does Differen. use dimer or trimer like. */
2419 /* XDrawString(display,window,gc,xcoord,15,"Dimer",5); */
2420 MethodTitle= strcat(MethodTitle, "Dimer");
2423 else if (like2or3==3) {
2424 /* XDrawString(display,window,gc,xcoord,15,"Trimer",6); */
2425 MethodTitle= strcat(MethodTitle, "Trimer");
2431 if (Method == PairCoil) {
2432 if (mode & PAIRCOIL_PAIRS) { /* Are we using pair or single prob */
2433 /* XDrawString(display,window,gc,xcoord,15,"Pair",4); */
2434 MethodTitle= strcat(MethodTitle, "Pair");
2438 /* XDrawString(display,window,gc,xcoord,15,"Single",6); */
2439 MethodTitle= strcat(MethodTitle, "Single");
2443 } /* Ends the stuff for dealing with PairCoil, TDRes, TDCoil. */
2446 /*****************************************************************************/
2447 /*************************** Now write the methodname. *********************/
2449 /* XDrawString(display,window,gc,xcoord,15,methodname[method],strlen(methodname[method])); */
2450 xcoord += 9*strlen(methodname[Method]);
2451 MethodTitle= strcat(MethodTitle, methodname[Method]);
2453 if (Method == MultiCoil) {
2454 if (mode & ONLY_COILED_CLASSES)
2455 strcat(MethodTitle,"Olig");
2458 strcat(MethodTitle,"WtdAvg");
2460 else if (Coil_Score ==2)
2461 strcat(MethodTitle, "MAX");
2463 if (start_coil_at_reg_shift)
2464 strcat(MethodTitle, "Shift");
2468 /* Now print Table are using. */
2469 /* And if are using a particular offset, print that. */
2470 if ( (Method==MultiCoil) || (Method==PairCoil) || (Method == HMMCoil)
2471 || (Method==PairCoilDiff) || (Method==PairCoilDiffAvg)) {
2474 if (NUMBER_TABLES >1) /** MultiCoil logfile uses all tables.... ***/
2475 if ( (Method != MultiCoil) || (!for_logfile)) {
2476 sprintf(buf,"%d",which_table+1);
2477 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2478 xcoord += 9*strlen(buf);
2479 MethodTitle= strcat(MethodTitle, buf);
2483 sprintf(buf,""); /* Clear out buf. */
2486 if ((offset_to_use != -1) && (offset_to_use != 7))
2487 sprintf(buf,"%c",'a' + offset_to_use);
2488 else if (Method == MultiCoil)
2489 if (offset_to_use == -1) {
2490 strcat(MethodTitle,"Max");
2494 else if (offset_to_use == 7) {
2495 strcat(MethodTitle,"Combo");
2500 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2501 xcoord += 9*strlen(buf);
2502 MethodTitle= strcat(MethodTitle, buf);
2506 /***************************************************************************/
2509 int is_not_differentiator(int Method)
2511 if ( (Method == MultiCoil) || (Method == PairCoil) || (Method==HMMCoil) ||
2512 (Method==NEWCOILS) || (Method == ActualCoil) )
2521 log_file_header(FILE *flog, int mode, int argc, char *argv[],
2522 int avg_max, double bound,
2523 int preprocessor_method, int preproc_table,
2524 int log_offset_to_use, int start_coil_at_reg_shift,
2525 int main_method, int main_table,
2526 char pair_lib[MAXFUNCTNUM],
2527 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM])
2530 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2531 char *MethodTitle = MethodTitleArray;
2533 MethodTitleArray[0] = '\0'; /* Initialization. */
2534 make_methodtitle(MethodTitle, method,main_table,mode,1);
2536 fprintf(flog,"\nInput file is %s", argc>1?argv[1]:"-");
2538 fprintf(flog, "\n\nScores by %s.", MethodTitle);
2540 "\nThe probability cutoff for the coiled-coil locater is %.2lf.",
2543 if (main_method == PairCoil) {
2544 if (mode & WEIGHTED_PROBS) {
2545 fprintf(flog,"\nStatistically weighted PairCoil score is ");
2548 fprintf(flog,"\nPairCoil score is ");
2551 fprintf(flog,"AVERAGE of residue scores.");
2554 fprintf(flog,"MAX of residue scores.");
2557 if (mode & USE_LIKE_LINE) {
2558 fprintf(flog,"\nPaircoil used the likelihood line.");
2561 fprintf(flog,"\nPaircoil used direct prob. formula instead of likelihood line.");
2564 fprintf(flog, "\nThe distances used by PairCoil are: ");
2566 for (i=0; i<pair_functnum; i++) {
2567 fprintf(flog,"%d ",pair_lib[i]);
2571 if ( (main_method==PairCoilDiff) || (main_method==PairCoilDiffAvg)) {
2572 char MethodTitleArray[500]; /* For printing the methodname. */
2573 char *MethodTitle = MethodTitleArray;
2574 int oligimerization;
2576 MethodTitleArray[0] = '\0'; /* Initialization. */
2577 make_methodtitle(MethodTitle, preprocessor_method, preproc_table,mode,0);
2580 oligimerization == 2;
2581 else oligimerization ==3;
2583 fprintf(flog, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2584 /* Whether logfile is for 2 or */
2585 /* 3 stranded predictions. */
2587 if (main_method==PairCoilDiffAvg) {
2588 if (start_coil_at_reg_shift) {
2589 fprintf(flog,"\n New coils were considered to start at register shifts.");
2592 fprintf(flog,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2598 fprintf(flog,"\nMain Table is Table %d\n\n",main_table);
2601 if (main_method == MultiCoil) {
2603 fprintf(flog, "\nThe scoring distances are: ");
2604 for (tablenum =0; tablenum<NUMBER_TABLES; tablenum++) {
2605 /* fprintf(flog,"\n For table %d: ",tablenum+1); */
2606 if (tablenum == 0) fprintf(flog,"\n For dimeric table: ");
2607 if (tablenum == 1) fprintf(flog,"\n For trimeric table: ");
2608 for (i=0; i<num_dist[tablenum]; i++) {
2609 fprintf(flog,"%d ",multi_lib[tablenum][i]);
2614 if (mode & POS_MODE)
2615 fprintf(flog,"\nIf the start or end of a coil is listed inside (), then that position does not occur in the pos file.\n");
2623 coil_tuple_output(int number_tables, char *lib, int mode, int avg_max,
2624 double bound, Sequence sequence, FILE *fout_coils,
2625 int number_multi_lib[MAX_TABLE_NUMBER],
2626 int number_score_dim)
2628 char offsets[MAXSEQLEN];
2632 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
2633 extern double *all_scores;
2636 for (tab=0; tab< number_tables; tab++) {
2637 for (dist=0; dist<number_multi_lib[tab]; dist++) {
2639 /*************************************************************************/
2640 /**********This computes offset of that table/dist method of scoring. ***/
2642 for (i=0; i<sequence.seqlen; i++) {
2643 for (j=0; j<POSNUM; j++)
2644 if (all_preprocess_like[tab + number_tables*(int)dist][i][j] ==
2645 all_preprocess_like[tab + number_tables*(int)dist][i][7]) {
2646 offsets[i] = (j-i) % POSNUM;
2647 if (offsets[i]<0) offsets[i] += POSNUM;
2649 /* This is a new thing so that if an offset change is not real */
2650 /* we keep the old offset. **********/
2651 if ( (i>0) && (offsets[i-1] != -1) &&
2652 (offsets[i] != offsets[i-1]) &&
2653 (all_preprocess_like[tab + number_tables*(int)dist]
2654 [i][(i+offsets[i-1]) %POSNUM]
2655 == all_preprocess_like[tab+number_tables*(int)dist][i][7]) )
2656 offsets[i] = offsets[i-1];
2659 /**************************************************************************/
2660 /****** Now compute coil scores for that method (table,dist) pair. ***/
2661 get_raw_coil_score(sequence,
2662 all_preprocess_like[tab+number_tables*(int)dist],
2663 all_scores[tab+number_tables*(int)dist], mode,
2664 avg_max,bound,offsets);
2668 /************* Now output the tuples of coil scores. ****/
2669 tuple_output2(sequence, mode, fout_coils,
2671 bound, number_score_dim);
2676 switch_gauss_param(int Offset_to_Use, int Method)
2678 if (Method == MultiCoil) {
2679 if (Offset_to_Use == 7) { /* Combo offset */
2680 determinants= &both_determinants[0][0];
2681 inverse_covars= &both_inverse_covars[0][0][0][0];
2682 means_submatrix= &both_means_submatrix[0][0][0];
2684 else { /* max offset or some fixed offset */
2685 determinants= &both_determinants[1][0];
2686 inverse_covars= &both_inverse_covars[1][0][0][0];
2687 means_submatrix= &both_means_submatrix[1][0][0];
2694 void usage(char *prog_command)
2699 /** Strip off directories to get just the program name. **/
2700 i=strlen(prog_command);
2701 while ((i>=0) && (prog_command[i] != '/') )
2703 prog_name= &prog_command[i+1];
2705 /************ This is for my extra stuff......
2706 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);
2707 **********************/
2709 fprintf(stderr,"\n Usage: %s <file-name> [-config config_file_name]\n",prog_name);