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() */
5 #define XCOORD_SCORED_BY 430
6 #define XCOORD_SEQSCORES 0
7 #define YCOORD_SEQSCORES 55 /* Each char takes "15" so put us below printing */
8 #define WIDTH_ADD_ON 30 /** Empty space after seq in main window. **/
9 #define xmin 120 /* Want this far enough over so have room for seqscoreb */
10 /** seqscoreb is 108 wide. xmin used to be 60 **/
14 #include "sc2seq.h" /* This also includes scconst.h, sc.h, interface.h */
16 #include <X11/keysym.h>
22 Window window, zoom, printb, quit, next, prev, autob, zoomb, methodb, paramb, helpb, help, seqscoreb;
23 GC gc, gcgray, gcdash;
27 #include <sys/types.h>
36 /* FILE *debug_out; */
39 /*************************** MULTICOIL LIKELIHOOD STUFF. *********************/
40 double both_determinants[2][MAX_CLASS_NUMBER];
41 double *determinants= &both_determinants[0][0];
43 double both_inverse_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
44 double *inverse_covars= &both_inverse_covars[0][0][0][0];
46 double both_covars_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
48 double both_means_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
49 double *means_submatrix= &both_means_submatrix[0][0][0];
51 double init_class_prob[MAX_CLASS_NUMBER];
53 int multi_functnum[MAX_TABLE_NUMBER];
54 int number_score_dim=0;
56 int number_multi_lib[MAX_TABLE_NUMBER];
57 /* Used to detrmine the number of different library */
58 /* distance values to do gaussian param over. */
61 FILE *ftotal_like=NULL; /* Extern for now. ***/
62 double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER];
63 /* When ftotal_like is read from get_defaults(), the */
64 /* i,j entry will hold the total gaussian like when */
65 /* assuming dimers have init_prob = */
66 /* MIN_DIMER_PROB + i*DIMER_GRID_STEP */
67 /* and trimers have init_prob = */
68 /* MIN_TRIMER_PROB + j*TRIMER_GRID_STEP */
69 /* These constants are defined in scconst.h. */
70 /* The value is updated during "output_seq()". */
72 /*************************** PAIRDIFF LIKELIHOOD STUFF. *********************/
73 double differ_determinants[MAX_CLASS_NUMBER];
74 double differ_inverse_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
75 double differ_covars_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
76 double differ_means_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
78 int number_differ_score_dim =0;
79 int number_differ_lib=0;
80 int differ_window_length = 28;
81 int which_differentiator=1; /** Default is just the score for the window **/
82 /** 2 means to give each residue a pair of scores**/
83 /** corresponding to the max and min of the **/
84 /** windows containing it. **/
85 int differ_functnum=0;
88 /************************** General Scoring parameters. ********************/
89 int functnum=0, pair_functnum=0;
91 double log_bound = 0; /*In PRN_MODE only scores > bound output to logfile.*/
93 int WINDOW = PAIRWINDOW;
94 int window_length[MAX_TABLE_NUMBER];
96 double SCALE0 = DEFAULT_SCALE0;
97 double scale0s[MAX_TABLE_NUMBER]; /** Both used as bonnie's scale0 and **/
98 double scale0p[MAX_TABLE_NUMBER]; /** in computing posterior prob. from **/
99 /** prior as sigma^2(prior)=scale0*sigma^2 */
100 /** where these are the standard deviation *
101 /** of the n*(prior_prob) and the distib */
102 /** of number of hits when sampling */
104 int log_offset_to_use;
105 int offset_to_use=-1; /* Default offset in DimerTrimerScore is all. */
106 int avg_max=1; /* Default means use max coil scores in log file. */
107 /** Value of 0 means do average, 1 means do max. **/
108 /** 2 means do both. */
109 int number_tables =0;
111 int method = MultiCoil;
112 int main_method = MultiCoil;
113 int main_preprocessor_method = MultiCoil;
114 int preprocessor_method = MultiCoil;
119 /*********************** Parameters involved in Coil Scores ***********/
120 int by_coil_or_seq =0; /*** To determe what type of coil_out file to do. **/
121 int weighted_avg = 0; /* Determines which method for determining coil */
122 /* score will be used in DimerTrimerScore. */
123 int start_coil_at_reg_shift =0; /* For averaging scores over coil length. */
124 int Coil_Score =0; /* 0 means for MultiCoil do residues score, 1 means */
125 /* average scores over coils, 2 means take max coil */
131 char *methodname[] = {"MultiCoil", "Coil","NEWCOILS", "HMMCoil","ActualCoil",
133 "PairCoilDiffAvg"}; /* "Coil" is PairCoil or */
137 char table_diff_stat_filename[MAXLINE]="/dev/null";
138 char weightfile_outname[MAXLINE]="/dev/null";
148 double seq_scores[MAX_CLASS_NUMBER];
150 double *scores; /** The space these pointers points to is alocated by **/
151 /** using static variable arrays in ScoreSeq */
152 /** all_scores is scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM +1]*/
153 /** scores is scores[current_score_algor][MAXSEQLEN][POSNUM+1] **/
155 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
156 /** Since preprocess_likelihood is used everywhere, just */
157 /** make it an external variable. */
158 /** This is given a value in ScoreSeq and sometime in */
160 double *preprocess_like;
161 /* This will point to the current table version */
162 /* of all_preprocess_like */
163 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
167 /************************ XWindow stuff. ***********************************/
169 int ps_res_per_line=0; /* 0 is a flag that it wasn't set in get_defaults */
170 /* so should just print it using nominal_length */
171 /* as done on the screen. If it is set to -1 */
172 /* then should print it all on one line. */
177 " *********************************",
179 " ** MultiCoil (version 1.0) **",
182 " *********************************",
184 "MultiCoil is described in:",
185 " Ethan Wolf, Peter Kim, Bonnie Berger,",
186 " `MultiCoil: A Program for Predicting Two- and
187 " Three-stranded Coiled Coils'.",
189 "Email inquiries to: multicoil-help@theory.lcs.mit.edu.",
195 "help: Clicking on help displays this message.",
197 "next: Displays the the next sequence.",
198 " Same as pressing 'n' or down-arrow.",
200 "prev: Displays the previous sequence.",
201 " Same as pressing 'p' or up-arrow.",
203 "auto: Automatically steps through the sequences",
206 "zoom: Gives more information on the highest",
207 " scoring residue. Same as pressing 'z'.",
208 " Clicking on a specific region will give",
209 " more information on that region.",
211 "print: Prints the sequence.",
214 "table: Switches the table used by PairCoil and",
215 " by TableDiff methods, same as 't'. Note that",
216 " key 'i' (input_table) changes the preprocessor's table.",
219 "method: Switches the scoring method, same as 'm'.",
220 " 'w' swithces Which preprocessor is used. ",
222 "key 'Up' or 'Down' moves the bound up or down by .01",
223 "key 'Return' rescores the sequence after changing the bound",
225 "key 's' or 'l': switch between singles and pair",
226 " probabilities at lib distances.",
228 "key 'o': Calc paircoil class prob. with respect only to",
229 " positive oligimerization classes (ignore pdb- class)",
230 "key '7': MultiCoil max score is obtained by taking tuple",
231 " of each dimension's max score (combo register).",
232 "key 'x': MultiCoil computes scores for each offset, taking",
233 " max of those scores to be the max score (max reg).",
234 "key 'r': Cycle through offsets 0-7 for register (7=combo reg)",
236 "key 'c': Switch between MultiCoil register score, average coil",
237 " score, and max coil score",
238 "key 'b': For Coil score, toggle whether start new coils at",
239 " a register shift (break in structure).",
240 "key 'a': For Coil score, toggle between unweighted and weighted avg",
242 "quit: Quits the program and completes log file.",
243 " 'q' quits without completing the log.",
244 " Pressing 'q' in any window also closes ",
251 ShowTitle (Window win)
255 for (i = 0; motd[i][0]; ++i)
256 XDrawString(display, win, gc, 5, 17*i+20, motd[i], strlen(motd[i]));
257 XDrawString(display,win,gcgray,100,17*i+20,
258 "Press Any Key To Continue",25);
261 ShowHelp (Window win)
265 /* Places hmotd text in win */
266 for (i=0; hmotd[i][0]; ++i)
267 XDrawString(display,win,gc,5,17*i+20, hmotd[i],strlen(hmotd[i]));
272 static int seqnum=0, sequp=0, seqlow=1;
275 main (int argc, char *argv[])
277 struct itimerval interval;
281 char buff[MAXLINE+MAXSEQLEN];
285 /* Note that fout_raw[MAX_TABLE_NUMBER+1] and fout_coils[MAX_TABLE_NUMBER+1] */
286 /* will be files of tuples of each residue's scores according to all tables. */
288 char *command_line = NULL;
290 FILE *fgin=NULL, *fout=NULL,
291 *flog=NULL, *flog_coil_conflicts=NULL, *fout_coils=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 extern int number_classes;
302 char class_sc_filenames[2][MAX_CLASS_NUMBER][MAXLINE];
303 double both_class_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM]
304 double both_class_means[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
306 char gauss_param[2][MAXLINE];
307 /** 0 entry is for combo offset, 1 entry is for **/
308 /** max of other registers. ***/
309 extern int functnum; /* What is used in scscore.c. Takes value from */
310 /* pair_functnum and multi_functnum. */
312 char lib[MAXFUNCTNUM]; /* Use MAXFUNCTUM=7 here since are 7 distances */
313 /* in coil, so we know AT MOST 7 lib functions. */
314 extern int pair_functnum; /* The number of distances in lib[].*/
316 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM];
317 extern int multi_functnum[MAX_TABLE_NUMBER];
320 char differ_class_sc_filenames[MAX_CLASS_NUMBER][MAXLINE];
321 double differ_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
322 double differ_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
323 char differ_gauss_param[MAXLINE];
324 char differ_lib[MAXFUNCTNUM]; /* Distance values for PairCoilDiff */
325 extern int differ_functnum;
331 int mode = 0; /* If the i-th option is set, mode gets 1 in i-th bit. */
334 double upper_bound; /* In DEBUG_MODE (for printing out the predicted */
335 /* heptad repeat), scores in [bound,upper_bound] */
336 /* are printed in capital letters. */
338 double m[MAX_TABLE_NUMBER],
339 b[MAX_TABLE_NUMBER]; /*The likelihood line is mx+b for paircoil. */
340 double m_single[MAX_TABLE_NUMBER],
341 b_single[MAX_TABLE_NUMBER]; /* Like line for singlecoil. */
344 int table_to_remove_from = -1;
345 /* Used if want to remove ALL input seq */
346 /* from that table at outset. **/
348 int old_num_dim_table[2];
349 int new_num_dim_table[2];
350 int old_num_dim_differ =0;
351 /* The gauss param file for multiple distances is */
352 /* contains the parameters for all distances. */
356 double prior_freq_single[MAX_TABLE_NUMBER],
357 prior_freq_pair[MAX_TABLE_NUMBER];
358 int which_priors =0, which_priorp =0;
359 int good_turing_fixed_disc =0; /** If set to 1 do the good_turing_fixed_disc estimate of probs*/
360 int structural_pos[POSNUM+1]; /** Holds list of which positions are **/
361 /** structurally important. All others **/
362 /** get genbnk probabilities.... **/
363 /** Value -1 indicates no more entries. **/
367 /*** INITIALIZATION ***/
369 fprintf(stderr,"\nAbout to initialize");
370 initialize(&number_classes, &number_multi_lib[0], window_length, scale0s,
372 &likelihoods[0][0], &pir_name[0], print, gauss_param,
374 &mode, prior_freq_single, prior_freq_pair,
377 fprintf(stderr,"\nInitialized");
378 read_command_line(argc, argv, &command_line, multi_functnum,
379 &mode, number_multi_lib, multi_lib);
380 fprintf(stderr,"\n\nRead commandline, command_line=%s",command_line);
382 /* Get the default options in the .paircoil file!!! */
383 get_defaults(argv[1], &mode, &log_bound, &fgin, fpin,
385 &flog, &flog_coil_conflicts,
386 &fout_coils, &by_coil_or_seq, likelihoods,
387 pir_name, lib, differ_lib, multi_lib,
388 &pair_functnum, &differ_functnum, multi_functnum,
389 print, &main_method, &main_preprocessor_method,
390 &main_table, &number_tables, &offset_to_use,
391 &avg_max, &weighted_avg, &Coil_Score,
392 &start_coil_at_reg_shift,
393 table_diff_stat_filename, weightfile_outname,
394 &ps_res_per_line, &number_classes,
395 class_sc_filenames, differ_class_sc_filenames,
396 gauss_param, differ_gauss_param, number_multi_lib,
399 &table_to_remove_from,
400 command_line, window_length, scale0s, scale0p,
401 prior_freq_single, prior_freq_pair, &which_priors,
402 &which_priorp, &good_turing_fixed_disc,
403 structural_pos, &which_differentiator, &differ_window_length);
405 fprintf(stderr,"\nOut of get_def");
411 compute_multivariate_num_differ_dim(&number_differ_score_dim,
413 differ_functnum, &old_num_dim_differ,
414 which_differentiator);
418 compute_multivariate_num_dim(mode, &number_score_dim,
419 number_multi_lib,multi_functnum, number_tables,
420 new_num_dim_table, old_num_dim_table);
422 if (number_classes >0) { /*** Get gaussian param for the classes. */
423 if (ftotal_like) { /** Want to output total gauss like over classes*/
424 for (i=0; i<number_classes; i++) /* for various init_prob settings. */
425 init_class_prob[i]=1; /** All 1 values is flag to just output */
426 init_totals(total_gauss_like); /** gauss values, not like values in */
427 } /** convert_raw_scores_to_gauss_like. */
431 if (gauss_param[i][0] != ',') {
432 get_gauss_param_all_classes2(class_sc_filenames[i],
433 old_num_dim_table[0]+old_num_dim_table[1],
434 both_class_means[i], both_class_covars[i],
435 number_classes, gauss_param[i]);
438 /******Note we store the submatrices in the original matrices....****/
439 compute_submatrices2(both_class_means[i],both_means_submatrix[i],
440 both_class_covars[i],number_classes,
441 both_inverse_covars[i],
442 both_determinants[i], multi_lib,
443 old_num_dim_table[0], old_num_dim_table[1],
444 new_num_dim_table[0], new_num_dim_table[1],
445 mode & MULTI_TRIMER_PAIRS);
447 for (class=0; class<number_classes; class++) {
448 fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
449 print_matrix(both_means_submatrix[i][class],1,
450 new_num_dim_table[0] + new_num_dim_table[1],stderr);
451 fprintf(stderr,"\n\nInver Covar Submatrix %d:\n", class);
452 print_matrix(both_inverse_covars[i][class],
453 new_num_dim_table[0] + new_num_dim_table[1],
454 new_num_dim_table[0] + new_num_dim_table[1],stderr);
458 fprintf(stderr,"\nDet0 =%lf, det1= %lf, det2= %lf\n",
459 both_determinants[i][0],both_determinants[i][1],
460 both_determinants[i][2]);
462 fprintf(stderr,"\nnew_num_dim_tab0=%d, nuw_num_dim_tab1=%d, old_num_dim_tab0=%d, old_num_dim_tab1=%d, number_score_dim=%d", new_num_dim_table[0], new_num_dim_table[1], old_num_dim_table[0], old_num_dim_table[1], number_score_dim);
469 switch_gauss_param(offset_to_use, main_method);
471 fprintf(stderr, "\n Out of get_gauss_param");
476 /*************************************************************************/
477 /******************Now do gauss parameters for PairCoilDiff method. *****/
478 /*************************************************************************/
479 /**** Note the "2" parameter represents that are 2 classes (dimers, trimers)***/
480 /**** not three, since just differentiating so don't need pdb-. ************/
482 if (differ_gauss_param[0] != ',') {
483 char differ_lib_hack[2][MAXFUNCTNUM];
484 /** If are computing max and min window scores, **/
485 /** then need to extract both dimenion i and **/
486 /** i + old_num_dim_differ/2. Treat them like they **/
487 /** were from another table. ***/
490 get_gauss_param_all_classes2(differ_class_sc_filenames,
492 differ_means, differ_covars,
493 2, differ_gauss_param);
495 if (which_differentiator ==2) {
496 for (i=0; i<number_differ_score_dim/2; i++) {
497 differ_lib_hack[0][i] = differ_lib[i];
498 differ_lib_hack[1][i] = differ_lib[i];
500 compute_submatrices2(differ_means,differ_means_submatrix,
502 differ_inverse_covars,
503 differ_determinants, differ_lib_hack,
504 old_num_dim_differ/2, old_num_dim_differ/2,
505 number_differ_score_dim/2, number_differ_score_dim/2,
509 compute_submatrices2(differ_means,differ_means_submatrix,
511 differ_inverse_covars,
512 differ_determinants, differ_lib,
513 old_num_dim_differ, 0,
514 number_differ_score_dim, 0,
518 for (class=0; class<2; class++) {
519 fprintf(stderr,"\n\nDiffer Mean Submatrix %d:\n",class);
520 print_matrix(differ_means_submatrix[class],1,
521 number_differ_score_dim,stderr);
522 fprintf(stderr,"\n\nDiffer Inver Covar Submatrix %d:\n", class);
523 print_matrix(differ_inverse_covars[class],
524 number_differ_score_dim,number_differ_score_dim,
530 fprintf(stderr,"\nDet0 =%lf, det1= %lf\n",
531 differ_determinants[0],differ_determinants[1]);
533 fprintf(stderr,"\nnew_num_dim_differ=%d, old_num_dim_differ=%d,", number_differ_score_dim, old_num_dim_differ);
541 fprintf(stderr, "\n Out of get_gauss_param");
545 log_offset_to_use = offset_to_use;
546 method = main_method;
547 preprocessor_method = main_preprocessor_method;
551 if ( (mode & PRN_MODE) && (flog || flog_coil_conflicts)){
552 log_file_header(flog, flog_coil_conflicts, mode,
553 argc, argv, avg_max, log_bound,
554 preprocessor_method, preproc_table, log_offset_to_use,
555 start_coil_at_reg_shift, weighted_avg, main_method,
559 if (mode & POS_MODE) ActualCoil_Method = available;
562 /****** Now if didn't get locations from .paircoil, get from ENV variable */
563 if (!fgin) { /* Didn't get fgin in */
564 fprintf(stderr,"\n\n **********************");
565 fprintf(stderr,"\nDid not get genbnk location from config file.");
566 fprintf(stderr,"\n Will instead use uniform distrib. as underlying distrib.");
567 fprintf(stderr,"\n**********************\n\n");
569 /**** exit(-1); ****/ /* get_defaults. */
571 if (!fpin[0]) { /* Didn't get fpin in */
572 fprintf(stderr,"\nDid not get positive data set from config file.");
573 exit(-1); /* get_defaults. */
575 /******************************************/
578 PairCoilData(fgin,fpin,&number_tables, mode & PROLINE_FREE_MODE,
579 table_to_remove_from, argv[1], mode,
580 prior_freq_single, prior_freq_pair, which_priors,
581 which_priorp, good_turing_fixed_disc, structural_pos);
582 /* If fpin[1] != NULL this will call */
583 /* TrimerDimerDiff to set up table diff */
585 /* If it is NULL then table is average */
586 /* of probabilities from 1 and 2. */
590 hmm = CCHMM(fgin,fpin);
594 if (mode & USE_LIKE_LINE) {
595 /** Don't want to compute likelihood line if */
597 functnum = pair_functnum;
598 get_likelihood_line(likelihoods, m, b, m_single, b_single,
599 lib, pair_functnum,pir_name,mode,number_tables);
601 for (i=0; fpin[i]; i++)
602 fprintf(stderr, "\n m[%d]= %lf, b[%d]= %lf\n",i,m[i],i,b[i]);
611 if (fgin) sclose(fgin);
617 fin = sopen(argc>1?argv[1]:"-","r");
618 fprintf(stderr, "\nfin is opened to %s",argc>1?argv[1]:"-");
620 /*** Note: If can't open xdisplay then OpenScreen will write log and */
621 /*** out files if those options are set. */
622 fprintf(stderr,"\nInitializing\n");
623 OpenScreen(mode,log_bound,fin,flog,flog_coil_conflicts, fout_coils,
624 fout, m,b,m_single, b_single,
625 lib,differ_lib,multi_lib,
626 main_method,main_preprocessor_method, &seqnum, which_priors,
627 which_priorp, prior_freq_single, prior_freq_pair,
628 good_turing_fixed_disc, structural_pos);
630 /* Let user read message while initializing */
636 interval.it_value.tv_sec = 2;
637 interval.it_value.tv_usec = 0;
638 interval.it_interval.tv_sec = 2;
639 interval.it_interval.tv_usec = 0;
642 XSelectInput(display, window, KeyPressMask | ButtonPressMask | ExposureMask);
644 if (autoadvance) do {
645 FD_SET(ConnectionNumber(display),&fds);
646 if (!select(FD_SETSIZE,&fds,&fdnone,&fdnone,&interval))
647 NextSeq(fin,lib,differ_lib,multi_lib,
648 m,b,m_single,b_single,
649 mode,bound,flog,flog_coil_conflicts,fout_coils,
650 fout, which_priors, which_priorp, prior_freq_single,
651 prior_freq_pair, good_turing_fixed_disc, structural_pos);
653 while ( (!XPending(display)) && (autoadvance));
654 /* Autoadvance until button is hit or end of seq file. */
655 XNextEvent(display,&report);
656 switch (report.type) {
658 if (report.xany.window == zoom)
660 else if (report.xany.window == help)
662 else if (sequence.seqlen)
663 ShowScores(scores,&sequence,NULL,bound, mode);
668 /* if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
669 m,b,m_single, b_single,
670 bound,flog,flog_coil_conflicts,
671 fout_coils,fout, which_priors,
672 which_priorp, prior_freq_single,
673 prior_freq_pair, good_turing_fixed_disc,
677 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
678 if (report.xany.window == zoom) {XUnmapWindow (display, zoom); break;}
679 if (report.xany.window == window) {
680 /* Zoom in and show register */
682 tmp = XYtoPos(report.xbutton.x,report.xbutton.y);
684 ZoomPnt(); zoomptr=tmp; Zoomer(); ZoomPnt();
687 if (report.xany.window == helpb)
688 XMapWindow(display,help);
689 if (report.xany.window == next)
690 NextSeq(fin,lib,differ_lib,multi_lib,
691 m,b,m_single, b_single,
692 mode,bound,flog,flog_coil_conflicts,fout_coils,fout,
693 which_priors,which_priorp, prior_freq_single, prior_freq_pair,
694 good_turing_fixed_disc, structural_pos);
695 if (report.xany.window == prev)
696 PrevSeq(fin,lib,differ_lib,multi_lib,
697 m,b,m_single,b_single,
698 mode,bound, which_priors,which_priorp, prior_freq_single,
699 prior_freq_pair, good_turing_fixed_disc, structural_pos);
700 if (report.xany.window == autob)
701 /* toggle the automatic sequencer */
702 { if (autoadvance ^= 1) NextSeq(fin,lib,differ_lib,multi_lib,
703 m,b,m_single,b_single,
705 flog_coil_conflicts,fout_coils,
706 fout, which_priors,which_priorp,
707 prior_freq_single, prior_freq_pair,
708 good_turing_fixed_disc,
712 if (report.xany.window == zoomb) {
713 if (zoomptr) ZoomPnt();
714 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
716 if (report.xany.window == printb) {
718 ps = sopen(print?print:defprint,"a");
719 PSDefs(ps, sequence.seqlen);
720 ShowScores(scores,&sequence,ps,bound,mode);
724 if (report.xany.window == tableb) {
725 table = (table +1) % number_tables;
726 ShowSeq(lib,differ_lib,multi_lib,
727 m,b,m_single,b_single,mode,0,0, bound);
728 if (zoomptr) Zoomer();
731 if (report.xany.window == methodb) {
732 method_cycle(&method,NUMBER_METHODS);
733 ShowSeq(lib,differ_lib,multi_lib,
734 m,b,m_single,b_single,mode,1,0, bound);
735 if (zoomptr) Zoomer();
737 if (report.xany.window == quit) {
738 if (flog || flog_coil_conflicts || fout || ftotal_like || fout_coils)
741 finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
742 fout,m,b,m_single,b_single,
743 lib,differ_lib,multi_lib,main_method,main_table,
744 main_preprocessor_method,
745 &seqnum, avg_max, which_priors,which_priorp,
746 prior_freq_single, prior_freq_pair,
747 good_turing_fixed_disc, structural_pos);
751 if (flog) sclose(flog); if (fout) sclose(fout);
752 if (fout_coils) sclose(fout_coils);
753 if (flog_coil_conflicts) sclose(flog_coil_conflicts);
754 if (ftotal_like) sclose(ftotal_like);
760 /* if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
761 m,b,m_single, b_single,
762 mode,bound,flog,flog_coil_conflicts,
763 fout_coils,fout, which_priors,
764 which_priorp, prior_freq_single,
765 prior_freq_pair, good_turing_fixed_disc,
769 if (report.xany.window == help) {XUnmapWindow(display,help); break;}
770 switch (XLookupKeysym(&(report.xkey),0)) {
773 if (bound > 1) bound =1;
774 draw_bound_box(bound, paramb);
778 if (bound < 0) bound =0;
779 draw_bound_box(bound, paramb);
782 case XK_Return: /** Signal that bound has been set, should rescore **/
783 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
784 ShowSeq(lib,differ_lib,multi_lib,
785 m,b,m_single,b_single,mode,-1,-1, bound);
786 if (zoomptr) Zoomer();
790 case XK_C: /** Switch between register, avg, and max **/
791 case XK_c: /** coil score for MultiCoil. **/
792 Coil_Score = (Coil_Score +1) % 3;
794 ShowSeq(lib,differ_lib,multi_lib, /** Change to res score **/
795 m,b,m_single,b_single,mode,1,1, -2,-2, bound);
796 else ShowSeq(lib,differ_lib,multi_lib, /** Change to coil score with */
797 m,b,m_single,b_single,mode,-1,-1, bound); /* out recomputing */
798 if (zoomptr) Zoomer(); /* res score. */
801 case XK_O: /** Oligomerization ratio. ***/
803 mode ^= ONLY_COILED_CLASSES; /* To XOR to complement that bit. */
804 if (number_classes >0) {
805 if (method == MultiCoil)
806 type_class_score_convert(sequence,all_scores,bound,
807 mode & ONLY_COILED_CLASSES,number_tables,
809 if (preprocessor_method == MultiCoil)
810 type_class_score_convert(sequence,all_preprocess_like,bound,
811 mode & ONLY_COILED_CLASSES,number_tables,
814 ShowSeq(lib,differ_lib,multi_lib,
815 m,b,m_single,b_single,mode,0,2, bound);
816 if (zoomptr) Zoomer();
821 ps = sopen(print?print:defprint,"w");
822 PSDefs(ps, sequence.seqlen);
823 ShowScores(scores,&sequence,ps,bound,mode);
829 {ZoomPnt(); --zoomptr; Zoomer(); ZoomPnt();}
832 if (zoomptr && zoomptr<sequence.seqlen)
833 {ZoomPnt(); ++zoomptr; Zoomer(); ZoomPnt();}
837 if (zoomptr) ZoomPnt();
838 zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
842 case XK_A: /** Switch between unweighted avg and **/
843 case XK_a: /** weighted avg coil score. **/
844 weighted_avg = !weighted_avg;
845 if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
846 ShowSeq(lib,differ_lib,multi_lib,
847 m,b,m_single,b_single,mode,-1,-1, bound);
848 if (zoomptr) Zoomer();
853 case XK_B: /** For CoilScore, change how deal with reg shift. **/
854 start_coil_at_reg_shift = !start_coil_at_reg_shift;
855 if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
856 ShowSeq(lib,differ_lib,multi_lib,
857 m,b,m_single,b_single,mode,-1,-1, bound);
858 if (zoomptr) Zoomer();
865 offset_to_use = (offset_to_use +1) % 8; /* Cycle through offsets */
866 /*** switch_gauss_param(offset_to_use); **/
867 ShowSeq(lib,differ_lib,multi_lib,
868 m,b,m_single,b_single,mode,1,1, bound);
869 if (zoomptr) Zoomer();
872 offset_to_use =-1; /* Do max scoring register */
873 ShowSeq(lib,differ_lib,multi_lib,
874 m,b,m_single,b_single,mode,1,1, bound);
875 if (zoomptr) Zoomer();
878 offset_to_use =7; /* Do MultiCoil max score as combo of max score */
879 ShowSeq(lib,differ_lib,multi_lib, /* in each dimension. */
880 m,b,m_single,b_single,mode,1,1, bound);
881 if (zoomptr) Zoomer();
884 case XK_s: /* Use single probabilities. */
886 if (mode & PAIRCOIL_PAIRS) {
887 mode -= PAIRCOIL_PAIRS; /** Turn it off **/
888 ShowSeq(lib,differ_lib,multi_lib,
889 m,b,m_single,b_single,mode,1,1, bound);
890 if (zoomptr) Zoomer();
894 case XK_l: /* Use pair probabilities at lib distances */
896 if (!(mode & PAIRCOIL_PAIRS)) {
897 mode |= PAIRCOIL_PAIRS; /* Turn it on. */
898 ShowSeq(lib,differ_lib,multi_lib,
899 m,b,m_single,b_single,mode,1,1, bound);
900 if (zoomptr) Zoomer();
909 XMapWindow(display,help);
914 if (report.xany.window == window)
915 {sclose(fin); if (flog) sclose(flog);
916 if (flog_coil_conflicts) sclose(flog_coil_conflicts);
917 if (ftotal_like) sclose(ftotal_like);
921 if (fout) sclose(fout);
923 if (report.xany.window == zoom) {
925 XUnmapWindow(display, zoom);
931 method_cycle(&method,NUMBER_METHODS);
933 if (method == MultiCoil)
934 table = table % number_classes;
936 table = table % number_tables;
938 ShowSeq(lib,differ_lib,multi_lib,
939 m,b,m_single,b_single,mode,1,0, bound);
940 if (zoomptr) Zoomer();
945 method_cycle(&preprocessor_method,NUMBER_PREPROCESSORS);
947 if (preprocessor_method == MultiCoil)
948 preproc_table = preproc_table% number_classes;
950 preproc_table = preproc_table%number_tables;
952 ShowSeq(lib,differ_lib,multi_lib,
953 m,b,m_single,b_single,mode,0,1, bound);
954 if (zoomptr) Zoomer();
958 if (preprocessor_method == MultiCoil)
959 preproc_table = (preproc_table +1) % number_classes;
961 preproc_table = (preproc_table +1) %number_tables;
962 ShowSeq(lib,differ_lib,multi_lib,
963 m,b,m_single,b_single,mode,0,2, bound);
964 if (zoomptr) Zoomer();
968 if (method == MultiCoil)
969 table = (table +1) % number_classes;
971 table = (table +1) % number_tables;
972 ShowSeq(lib,differ_lib,multi_lib,
973 m,b,m_single,b_single,mode,0,0, bound);
974 if (zoomptr) Zoomer();
979 if (report.xany.window == window) PrevSeq(fin, lib,differ_lib,
980 multi_lib,m,b,m_single,b_single,mode,bound, which_priors,
981 which_priorp, prior_freq_single,prior_freq_pair,
982 good_turing_fixed_disc, structural_pos);
990 if (report.xany.window == window)
991 NextSeq(fin,lib,differ_lib,multi_lib,
992 m,b,m_single,b_single,
993 mode,bound,flog,flog_coil_conflicts,fout_coils,
994 fout, which_priors, which_priorp,
995 prior_freq_single,prior_freq_pair, good_turing_fixed_disc,
1004 static void Alarm ()
1007 if (XCheckTypedEvent(display,Expose,&report)) {
1009 XDrawString(display,window,gcgray,500,17,"(initializing)",14);
1015 struct itimerval value;
1016 signal(SIGALRM,Alarm);
1017 getitimer(ITIMER_REAL,&value);
1018 value.it_value.tv_sec = 0;
1019 value.it_value.tv_usec = 500;
1020 value.it_interval.tv_sec = 0;
1021 value.it_interval.tv_usec = 500;
1022 setitimer(ITIMER_REAL,&value,NULL);
1027 struct itimerval value;
1028 getitimer(ITIMER_REAL,&value);
1029 value.it_value.tv_sec = 0;
1030 value.it_value.tv_usec = 0;
1031 setitimer(ITIMER_REAL,&value,NULL);
1032 XDrawImageString(display,window,gcgray,500,17," ",14);
1033 XDrawImageString(display,window,gcgray,225,729,"press any key to continue",25);
1039 if (!zoomptr) return;
1040 PostoXY(zoomptr,&x,&y);
1041 pts[0].x = x; pts[0].y = y;
1042 pts[1].x = -5; pts[1].y = 5;
1043 pts[2].x = 10; pts[2].y = 0;
1044 XFillPolygon(display,window,gcdash,pts,3,Convex,CoordModePrevious);
1045 XDrawLine(display,window,gcdash,x,y+5,x,y+8);
1049 /* Find most likely residue, break ties arbitrarily */
1054 for (i=0; i<sequence.seqlen; ++i)
1055 if (scores[8*i+7] > like) {
1056 like = scores[8*i+7];
1067 XMapWindow(display, zoom);
1068 /* Indicate where in the sequence we are */
1069 sprintf(buf,"Position %d ",zoomptr);
1070 XDrawImageString(display,zoom,gc,7,20,
1072 for (i = -3; i<=3; ++i) {
1073 if (zoomptr+i>0 && zoomptr+i<=sequence.seqlen)
1074 buf[0] = numaa(sequence.seq[zoomptr-1+i]);
1077 XDrawImageString(display,zoom,gc,159+9*i,20, buf,1);
1079 XDrawLine(display,zoom,gc,159,22,168,22);
1081 /* Indicate the probability of being in any given state */
1082 sc = &(scores[8*(zoomptr-1)]);
1084 for (i=0; i<=7; ++i)
1085 if (sc[i] == sc[7]) {offset=i; break;} /* Used to print out maxscoring */
1086 /* register as 'a' + offset. */
1087 /* Offset 7 will be a flag that no particular offset scored max. */
1089 /* Give likelihood */
1090 sprintf(buf,"%5.1f%% chance",100*sc[7]);
1091 XDrawImageString(display,zoom,gc,220,20, buf,strlen(buf));
1094 fprintf(stderr,"\noffset= %d, reg = %c",all_offsets[table][zoomptr-1],
1095 'a' + (zoomptr-1 + all_offsets[table][zoomptr-1])%POSNUM);
1100 sprintf(buf,"register %c",'a'+offset);
1102 sprintf(buf,"combo reg ");
1103 else sprintf(buf," ");
1104 XDrawImageString(display,zoom,gc,230,40, buf,strlen(buf));
1106 for (i=0; i<7; ++i) {
1107 sprintf(buf,"%c:%4.1f%% ",'a'+i,(100*sc[i]));
1108 XDrawImageString(display,zoom,gc,15+92*(i&3),60+20*(i/4),buf,strlen(buf));
1110 XDrawString(display,zoom,gc,10,110,"Click on viewgram or use arrow",30);
1111 XDrawString(display,zoom,gc,10,127,"keys to view other positions.",29);
1112 XDrawString(display,zoom,gc,10,144,"You may need to move this",25);
1113 XDrawString(display,zoom,gc,10,161,"window to see the viewgram.",27);
1123 /*********The following is for storing sequences. ***/
1124 static char seqbuff[40000];
1126 static Sequence seqs[30];
1127 static int seqbufpos[30];
1128 static int seqsizes[30];
1131 { int size, pos, limsize;
1133 size = seqsizes[sequp-seqnum];
1134 pos = seqbufpos[sequp-seqnum];
1136 sequence = seqs[sequp-seqnum]; /** Retrieve seq from memory. **/
1139 limsize = (pos+size<=40000) ? size : 40000-pos;
1140 bcopy(&(seqbuff[pos]),sequence.code,limsize);
1141 bcopy(seqbuff,&(sequence.code[limsize]),size-limsize);
1147 { int size, pos, limsize,i;
1148 char *temp, *current;
1150 bcopy(seqs,&(seqs[1]),29*sizeof(Sequence));
1151 bcopy(seqbufpos,&(seqbufpos[1]),29*sizeof(int));
1152 bcopy(seqsizes,&(seqsizes[1]),29*sizeof(int));
1153 if (sequp-seqlow>=30) seqlow = sequp-29;
1155 /* The 2 accounts for the end of string characters in .code and .title. */
1156 size = 2 + strlen(sequence.code) + strlen(sequence.title) +
1157 sequence.seqlen * (sequence.reg ? 2 : 1);
1159 pos = seqbufpos[1] + seqsizes[1]; /* Number of characters in array */
1161 /* if (pos>=40000) pos -= 40000;*/
1162 /* limsize = (pos+size<=40000) ? size : 40000-pos; */
1164 if (pos+size >=40000) pos=0; /* Array will be full so start at beginning. */
1165 limsize = (pos+size<=40000) ? size : 40000-pos;
1168 /* Save the "code+title+seq+reg" in something permanant, and wrap around
1171 /* David uses bcopy because wants to copy past end of string characters */
1173 bcopy(sequence.code,&(seqbuff[pos]),limsize);
1174 bcopy(&(sequence.code[limsize]),seqbuff,size-limsize);
1177 sequence.code=strcpy(&(seqbuff[pos]),sequence.code);
1179 strcpy(&(seqbuff[pos+strlen(sequence.code)+1]),sequence.title);
1183 &(seqbuff[pos+strlen(sequence.code) + strlen(sequence.title) + 2]);
1184 for(i=0; i<sequence.seqlen; i++)
1185 sequence.seq[i] = temp[i];
1190 sequence.reg=&(sequence.seq[sequence.seqlen]);
1191 for(i=0; i<sequence.seqlen; i++)
1192 sequence.reg[i] = temp[i];
1196 while ((seqlow<sequp) && (seqbufpos[sequp-seqlow] >= pos && seqbufpos[sequp-seqlow] < pos+limsize || seqbufpos[sequp-seqlow] < size-limsize))
1199 seqs[0].seqlen= sequence.seqlen;
1209 method_cycle(int *Method, int Number_of_Methods)
1212 *Method = (*Method+1)%Number_of_Methods;
1214 while ((*Method == MultiCoil && MultiCoil_Method == unavailable) ||
1215 (*Method == PairCoil && PairCoil_Method == unavailable) ||
1216 (*Method == NEWCOILS && NEWCOILS_Method == unavailable) ||
1217 (*Method == HMMCoil && HMMCoil_Method == unavailable) ||
1218 /**** (*Method == TableDiffRes &&
1219 TableDiffRes_Method == unavailable) ||
1220 (*Method == TableDiffCoil &&
1221 TableDiffCoil_Method == unavailable) ||
1223 (*Method == ActualCoil && ActualCoil_Method == unavailable) ||
1224 (*Method == PairCoilDiff && PairCoilDiff_Method == unavailable)||
1225 (*Method == PairCoilDiffAvg && PairCoilDiffAvg_Method == unavailable) ||
1227 (*Method == ActualCoil) && (!sequence.reg) /* Don't do actual coil */
1228 ); /* if didn't input register */
1234 PrevSeq (FILE *fin, char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1235 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1236 double *m, double *b,
1237 double *m_single, double *b_single,
1238 int mode, double bound, int which_priors, int which_priorp,
1239 double prior_freq_single[MAX_TABLE_NUMBER],
1240 double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1241 int structural_pos[POSNUM+1])
1243 if (seqnum<=seqlow) {XBell(display,0); return;}
1246 NewSeq(mode,sequence, which_priors, which_priorp, prior_freq_single,
1247 prior_freq_pair, good_turing_fixed_disc, structural_pos);
1248 ShowSeq(lib, differ_lib, multi_lib,
1249 m, b, m_single, b_single, mode,1,1, bound);
1253 NextSeq (FILE *fin,char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1254 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1255 double *m, double *b, double *m_single, double *b_single,
1257 double bound, FILE *flog, FILE* flog_coil_conflicts,
1259 FILE *fout, int which_priors, int which_priorp,
1260 double prior_freq_single[MAX_TABLE_NUMBER],
1261 double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1262 int structural_pos[POSNUM+1])
1264 char title[MAXLINE],code[MAXLINE];
1265 char seq[MAXSEQLEN], reg[MAXSEQLEN];
1267 int first_time_seq=0;
1270 fprintf(stderr,"\nIn next seq");
1273 if (seqnum < sequp) {
1279 sequence.seqlen = 0;
1282 sequence.title = title;
1283 sequence.code = code;
1286 /* if ( ( (mode & POS_MODE) ? !getpos(fin,&sequence) :
1287 !getseq2(fin,&sequence) ) ) */
1288 /*** getseq2 now gets register info too... ***/
1290 if (!getseq2(fin,&sequence))
1291 {XBell(display,0); sequence= temp_seq; autoadvance=0; return;}
1299 NewSeq(mode,sequence, which_priors, which_priorp,
1300 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1302 if (first_time_seq) {
1304 output_seq(lib,differ_lib, multi_lib,
1305 m,b,m_single, b_single,
1306 mode,log_bound,flog,flog_coil_conflicts,fout_coils,fout,
1307 avg_max, main_method, main_preprocessor_method, main_table);
1312 ShowSeq(lib, differ_lib,multi_lib,
1313 m, b, m_single, b_single,mode,1,1, bound);
1318 /** This is stuff need to always do before rescore a sequence. **/
1319 NewSeq_nonX (int mode, Sequence sequence, int which_priors, int which_priorp,
1320 double prior_freq_single[MAX_TABLE_NUMBER],
1321 double prior_freq_pair[MAX_TABLE_NUMBER],
1322 int good_turing_fixed_disc, int structural_pos[POSNUM+1])
1324 if (mode & TST_MODE0)
1325 recalc_prob(sequence, 0, mode, which_priors, which_priorp,
1326 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1328 if (mode & TST_MODE1)
1329 recalc_prob(sequence,1, mode, which_priors, which_priorp,
1330 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1334 NewSeq (int mode, Sequence sequence, int which_priors, int which_priorp,
1335 double prior_freq_single[MAX_TABLE_NUMBER],
1336 double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1337 int structural_pos[POSNUM+1])
1339 NewSeq_nonX (mode, sequence, which_priors, which_priorp,
1340 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1343 XUnmapWindow(display, zoom);
1348 /* The likelihood line is mx +b */
1349 void output_seq(char lib[MAXFUNCTNUM],char differ_lib[MAXFUNCTNUM],
1350 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1351 double *m,double *b,
1352 double *m_single, double *b_single,
1354 double log_bound,FILE *flog, FILE *flog_coil_conflicts,
1358 int main_method, int main_preprocessor_method,
1362 double maxother=0; double maxscore;
1363 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1368 double junk_scores[MAXSEQLEN][POSNUM+1];
1371 fprintf(stderr,"\nIn outputseq");
1374 if (is_not_differentiator(main_method))
1378 all_scores = ScoreSeq (lib, differ_lib, multi_lib,
1379 m, b, m_single, b_single, mode,
1381 main_method, main_preprocessor_method,
1382 log_offset_to_use, &maxscore,
1383 1, do_preproc, log_bound);
1387 add_to_total_likes(all_scores, total_gauss_like, sequence.seqlen,
1391 fprintf(stderr,"\nGot all_scores");
1393 scores = &all_scores[main_table* MAXSEQLEN*(POSNUM+1) ];
1395 /* This is all_scores[main_table + 0*number_tables][0][0] */
1396 /* so is score on main_table using library # 0 **/
1397 /* Similarly for preproces_like. **/
1399 preprocess_like = &all_preprocess_like[main_table][0][0];
1402 for (i=0; i<sequence.seqlen; i++)
1403 if (preprocess_like[i*(POSNUM+1) + 7] > maxother)
1404 maxother=preprocess_like[i*(POSNUM +1) + 7];
1405 /* This is preprocess_like[i][7] */
1408 switch (main_method) {
1410 if (mode & VER_MODE) {
1411 NEWCOILSScore(mode, sequence, &maxother,junk_scores,
1414 /* scores in junk_scores */
1415 log_output_ver(SC2SEQ | SCSTOCK, maxscore, maxother,
1416 maxscore,seqnum, sequence.title,
1417 sequence.code, flog);
1424 /*** For pos files coils, does a log file version of coil/seq score. **/
1425 output_pos_scores(fout_coils, sequence, all_scores, log_bound, mode,
1426 by_coil_or_seq, weighted_avg);
1429 if (mode & PRN_MODE)
1431 if (Coil_Score && flog) {
1432 /*** Were having trouble with this in printing out coils, since offset
1433 *** would sometime change at ends of coils between residue offset and
1434 *** coil score offset. To make log file print out nice, use res. offset.
1435 char all_offs[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1436 get_offsets(all_scores, sequence,all_offs,3);
1438 log_coil_output(all_scores, all_offsets,sequence, log_bound, flog,
1442 log_output2_prn(mode,maxscore, maxscore, scores, scores,
1444 seqnum, sequence, log_bound, log_bound, log_bound,
1452 if (mode & VER_MODE)
1453 log_output_ver(SC2SEQ , maxscore, 0,
1454 maxscore,seqnum, sequence.title,
1455 sequence.code, flog);
1458 if (mode & PRN_MODE)
1459 log_output2_prn(mode,maxscore, maxscore, scores, scores,sequence.seqlen,
1460 seqnum, sequence, log_bound, log_bound, log_bound,
1464 if (flog_coil_conflicts)
1465 coil_conflicts(mode,all_scores,all_scores,sequence, seqnum,log_bound,
1466 log_bound,number_tables, flog_coil_conflicts,avg_max);
1469 /* case PairCoilDiff: */
1470 case PairCoilDiffAvg:
1471 if (mode & PRN_MODE) {
1472 log_output2_prn(mode,maxscore, maxother, scores,preprocess_like,
1474 seqnum, sequence, .5, .5, log_bound,
1478 if (flog_coil_conflicts) {
1480 coil_conflicts(mode,all_scores,all_preprocess_like,sequence, seqnum,
1481 .5, log_bound, number_tables, flog_coil_conflicts,
1493 fprintf(stderr,"\n About to go into fout.");
1498 if ( (!(mode & RAW_OUT)) ||
1499 ( (main_method != MultiCoil) && (main_method != PairCoilDiff))) {
1500 if (main_method== MultiCoil)
1501 tuple_output2(sequence, mode, fout, all_scores,3, 1,log_bound,3);
1502 /* Output based on 3 tables, meaning */
1503 /* output the dimer and trimer likes */
1504 /* and total like. */
1505 else txt_output2(sequence, mode, fout, scores, log_bound);
1507 else { /* Want raw tuple_output for multi_coil and PairCoilDiff*/
1509 if (main_method == MultiCoil)
1510 if (log_offset_to_use==7) /* Combo offset just output max score */
1511 tuple_output2(sequence, mode, fout, /** for each dimension. */
1512 all_scores, number_tables,
1513 number_multi_lib,log_bound, number_score_dim);
1514 else /** Max reg for raw score can't really be determined. So instead */
1515 /* for coils, output scores for the known correct register, and */
1516 /* for pdb- output ALL possible 7 register scores (all negs). */
1517 tuple_output_for_max(sequence,mode,fout, all_scores, number_tables,
1518 number_multi_lib, log_bound, number_score_dim);
1520 else if (main_method == PairCoilDiff)
1521 /* Output max score for each dimension. */
1522 /** Had to do a hack since tuple_out2 expects number_lib */
1523 /** to be array of size num_tab, so pass in 1 for num_tab */
1524 tuple_output2(sequence, mode, fout, all_scores, 1,
1525 &number_differ_lib,log_bound, number_differ_score_dim);
1532 int copy_offsets(char preproc_offsets[MAXSEQLEN],
1534 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1536 double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1541 for (tablenum =0; tablenum < number_dimens; tablenum++)
1542 for(i=0; i<sequence.seqlen; i++) {
1543 all_offsets[tablenum][i]=preproc_offsets[i];
1544 /*** Make sure the output score corresponds to this offset... **/
1545 /*** Also do it for backup residue scores in tablenum +3 ***/
1546 if ((offset_to_use ==7) || (offset_to_use == -1)) {
1547 all_scores[tablenum][i][7] =
1548 all_scores[tablenum][i][(i+all_offsets[tablenum][i]) %POSNUM];
1549 all_scores[tablenum+3][i][7] =
1550 all_scores[tablenum+3][i][(i+all_offsets[tablenum][i]) %POSNUM]; }
1554 int get_offsets(double all_likes[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1555 Sequence sequence, char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1559 int got_offset, max_offset;
1560 double max_off_score;
1563 for (tablenum =0; tablenum < number_dimens; tablenum++)
1564 for(i=0; i<sequence.seqlen; i++) {
1565 offsets[tablenum][i]=-1;
1567 max_off_score=-HUGE_VAL;
1568 for (j=0; j<POSNUM; j++) {
1569 if (all_likes[tablenum][i][j] == all_likes[tablenum][i][7]) {
1570 offsets[tablenum][i] = (j-i) % POSNUM;
1571 if (offsets[tablenum][i]<0)
1572 offsets[tablenum][i] += POSNUM;
1575 if (all_likes[tablenum][i][j] > max_off_score) {
1576 max_off_score=all_likes[tablenum][i][j];
1577 max_offset=(j-i) % POSNUM;
1578 if (max_offset <0) max_offset +=POSNUM;
1581 if (!got_offset) offsets[tablenum][i] = max_offset;
1583 /* This is a new thing so that if an offset change is not real */
1584 /* we keep the old offset. **********/
1585 if ( (i>0) && (all_likes[tablenum][i][(i+offsets[tablenum][i-1]) %POSNUM]
1586 == all_likes[tablenum][i][(i+offsets[tablenum][i])%POSNUM]) )
1587 offsets[tablenum][i] = offsets[tablenum][i-1];
1599 void preprocessor_score(char lib[MAXFUNCTNUM],
1600 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1601 double *m, double *b,
1602 double *m_single, double *b_single, int mode,
1603 int Preprocessor_Method, int Offset_to_Use,
1604 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1606 char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1610 char offsets[MAXSEQLEN];
1615 int dimension, num_multilib;
1619 for (tablenum =0; tablenum < number_tables; tablenum++) {
1620 switch_tables(tablenum); /* Set prob. to be for that table. */
1621 if (tablenum ==0) use_pairs =1; /* Use pairs for cctb */
1622 else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
1623 /* Do what config file said to do for trimers. **/
1625 switch(Preprocessor_Method) {
1627 functnum = multi_functnum[tablenum]; /* For use in scscore.c **/
1628 num_multilib= compute_num_multilib(tablenum, number_multi_lib,
1629 mode & MULTI_TRIMER_PAIRS);
1630 for (dist=0; dist<num_multilib; dist++) {
1631 if (number_multi_lib[tablenum]>1) {
1632 local_lib = &multi_lib[tablenum][dist];
1634 else local_lib = multi_lib[tablenum];
1635 dimension = compute_dimension(tablenum, number_tables, dist,
1636 number_score_dim, number_multi_lib,
1637 mode & MULTI_TRIMER_PAIRS);
1639 MultiCoilDimensionScore(mode, sequence, local_lib,
1640 &maxscore,tablenum, offsets,
1641 all_preprocess_like[dimension],
1642 Offset_to_Use,use_pairs);
1647 if (mode & PAIRCOIL_PAIRS) {
1648 functnum = pair_functnum; /* For use in scscore.c **/
1649 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1650 &maxscore,tablenum, offsets,
1651 all_preprocess_like[tablenum],
1654 else /* Use singles probabilities... for likelihood use line 0 */
1655 SingleCoilScore(mode,sequence,m_single[tablenum],
1657 &maxscore, tablenum,
1658 all_preprocess_like[tablenum], Offset_to_Use,
1663 NEWCOILSScore(mode, sequence, &maxscore, all_preprocess_like[tablenum], Offset_to_Use);
1667 HMMScore(sequence.seq, sequence.seqlen, hmm, all_preprocess_like[tablenum],
1668 tablenum, Offset_to_Use);
1672 /* if (mode & POS_MODE) */ /* Score the coils in the posfile based on */
1673 /* their having 100% likelihood. */
1674 ActualCoils(sequence, all_preprocess_like[tablenum], Offset_to_Use,1);
1678 fprintf(stderr,"\nBad preprocessor, using PairCoil.\n");
1679 functnum = pair_functnum; /* For use in scscore.c **/
1680 PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1681 &maxscore,tablenum, offsets,
1682 all_preprocess_like[tablenum],
1688 if (Preprocessor_Method !=MultiCoil)
1689 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1692 /***** Get copy of res scores in tablenum and tablenum+3 **/
1693 switch_gauss_param(Offset_to_Use, Preprocessor_Method);
1695 convert_raw_scores_to_gauss_prob_like2(all_preprocess_like,sequence.seqlen,
1696 number_tables,means_submatrix,inverse_covars,
1697 determinants, number_classes, init_class_prob,
1699 mode & ONLY_COILED_CLASSES,
1700 bound,Offset_to_Use, 1);
1702 get_offsets(all_preprocess_like, sequence,all_preproc_offsets,3);
1704 if (Coil_Score) { /* Value of 1 means do average, 2 means do max. */
1705 /*** Need to save residue scores, since the coil scores **/
1706 /*** will replace them, but old residue scores might be needed **/
1707 /*** again in tablenum+3. ****/
1708 for (tablenum =0; tablenum <3; tablenum++) {
1709 average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
1710 all_preprocess_like[tablenum],
1711 all_preprocess_like[2+3],/** Total Coil Like. **/
1712 weighted_avg, offset_to_use,
1713 all_preproc_offsets[tablenum],
1714 start_coil_at_reg_shift, bound,
1728 double *ScoreSeq (char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1729 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1730 double *m_paircoil, double *b_paircoil,
1731 double *m_single, double *b_single,
1732 int mode, int Table, int Method,
1733 int Preprocessor_Method, int Offset_to_Use,
1735 int rescore_seq, int rescore_preproc, double bound)
1737 /** rescore_preproc == 1 means need to rescore. If it or rescore_seq are */
1738 /** -1 then means just need to recompute MultiCoil coil and seq scores. */
1739 /** -2 means need to change from MultiCoil coil score back to res score. */
1740 /** If it is another non-zero */
1741 /** value, don't need to rescore it, but do need to rescore anything that */
1742 /** uses it as a preproc. */
1744 double max_paircoil_score;
1745 static char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1746 extern char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1747 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1748 static double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1756 int num_multilib, dimension;
1757 int local_number_tables;
1759 double *preproc_res_like;
1760 int scored_already =0;
1762 /**************************************************************************/
1763 /********* PreprocLike is used in the 2-3 stranded differentiators. ****/
1766 if ( (rescore_preproc==1) &&
1767 (!ftotal_like)) { /* A hack to save computation time. */
1768 preprocessor_score(lib, multi_lib,
1769 m_paircoil, b_paircoil, m_single, b_single,
1770 mode, Preprocessor_Method,
1772 all_preprocess_like, sequence, all_preproc_offsets,
1776 /*****************************************************************************/
1777 if (Preprocessor_Method == MultiCoil) {
1778 if ( ( rescore_preproc==-1) && Coil_Score) {
1779 /* Coil_Score value of 1 means do average, 2 means do max. */
1780 /*** Need to save residue scores, since the coil scores **/
1781 /*** will replace them, but old residue scores might be needed **/
1782 /*** again in tablenum+3. ****/
1783 for (tablenum =0; tablenum <3; tablenum++)
1784 average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
1785 all_preprocess_like[tablenum],
1786 all_preprocess_like[2+3],/** Total Coil Like. **/
1787 weighted_avg, offset_to_use,
1788 all_preproc_offsets[tablenum],
1789 start_coil_at_reg_shift, bound,
1791 if (mode & ONLY_COILED_CLASSES) /** Convert to olig state. **/
1792 type_class_score_convert(sequence,all_preprocess_like,bound,
1793 mode & ONLY_COILED_CLASSES,number_tables,
1796 else if (rescore_preproc == -2) /** Switch back to res score.*/
1797 for (tablenum =0; tablenum <3 ; tablenum++)
1798 for (j =0; j <=POSNUM; j++)
1799 for (i=0; i<sequence.seqlen; i++) {
1800 all_preprocess_like[tablenum][i][j] =
1801 all_preprocess_like[tablenum+3][i][j];
1802 if ((mode & ONLY_COILED_CLASSES) && (tablenum <2)&&
1803 (all_preprocess_like[2+3][i][j]>bound))
1804 all_preprocess_like[tablenum][i][j]/=
1805 all_preprocess_like[2+3][i][j];
1809 /****************************************************************************/
1815 if (Method == PairCoilDiff)
1816 local_number_tables = 1;
1817 else local_number_tables= number_tables;
1820 /*** Clean up this conditional, since really the only thing that might ***/
1821 /*** change for paircoildiffer when change preproc is the offset to use **/
1822 /*** (i.e. all_scores[table][i][7]. ***/
1824 if ((rescore_seq>0) ||
1825 (!is_not_differentiator(Method) && rescore_preproc) ) {
1827 for (tablenum=0; tablenum<local_number_tables; tablenum++) {
1828 switch_tables(tablenum); /* Switch prob to right table. */
1830 if (tablenum ==0) use_pairs =1; /* Use pairs for cctb */
1831 else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
1832 /* Do what config file said to do for trimers. **/
1834 maxsc=0; /* Initialization. */
1840 functnum = pair_functnum; /* For use in scscore.c **/
1841 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1842 b_paircoil[tablenum],
1843 &maxsc,tablenum, all_offsets[tablenum],all_scores[tablenum],
1849 functnum = multi_functnum[tablenum]; /* For use in scscore.c **/
1850 num_multilib = compute_num_multilib(tablenum,number_multi_lib,
1851 mode & MULTI_TRIMER_PAIRS);
1852 for (dist=0; dist<num_multilib; dist++) {
1853 if (number_multi_lib[tablenum] >1)
1854 local_lib = &multi_lib[tablenum][dist];
1855 else local_lib = multi_lib[tablenum];
1856 dimension = compute_dimension(tablenum, number_tables, dist,
1857 number_score_dim, number_multi_lib,
1858 mode & MULTI_TRIMER_PAIRS);
1860 MultiCoilDimensionScore(mode, sequence, local_lib, &maxsc,
1861 tablenum,all_offsets[dimension],
1862 all_scores[dimension],
1863 Offset_to_Use,use_pairs);
1869 if (mode & PAIRCOIL_PAIRS) {
1870 functnum = pair_functnum; /* For use in scscore.c **/
1871 PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1872 b_paircoil[tablenum], &maxsc,
1873 tablenum,all_offsets[tablenum],
1874 all_scores[tablenum],
1875 Offset_to_Use, (mode & RAW_OUT));
1877 else /*Use singles probabilities... for likelihood use line number 0 */
1878 SingleCoilScore(mode,sequence,m_single[tablenum],
1881 all_scores[tablenum], Offset_to_Use,
1885 NEWCOILSScore(mode,sequence,&maxsc, all_scores[tablenum],
1890 HMMScore(sequence.seq, sequence.seqlen, hmm, all_scores[tablenum],
1891 tablenum, Offset_to_Use);
1895 functnum = differ_functnum; /* For use in scscore.c **/
1896 for (dimension=0; dimension<number_differ_lib; dimension++) {
1897 if (number_differ_lib >1) {
1898 local_lib = &differ_lib[dimension];
1900 else local_lib = differ_lib;
1903 PairCoilDifferDimension(mode, sequence, local_lib, &maxsc,
1904 all_preproc_offsets[preproc_table],
1905 all_preprocess_like[preproc_table],
1907 all_offsets, dimension,
1908 Offset_to_Use, number_differ_score_dim,
1912 /***********************
1913 PairCoilDiffer(mode, sequence, local_lib, &maxsc,
1914 all_preproc_offsets[tablenum],
1915 all_preprocess_like[preproc_table],
1916 Offset_to_Use, tablenum,
1917 all_scores[tablenum],0,
1918 weighted_avg,start_coil_at_reg_shift);
1920 *************************/
1924 case PairCoilDiffAvg:
1925 functnum = differ_functnum;
1926 PairCoilDiffer(mode, sequence, differ_lib, &maxsc,
1927 all_preproc_offsets[preproc_table],
1928 all_preprocess_like[preproc_table],
1929 Offset_to_Use, tablenum,
1930 all_scores[tablenum],
1931 1, weighted_avg,start_coil_at_reg_shift);
1936 ActualCoils(sequence, all_scores[tablenum], Offset_to_Use,0);
1939 if (tablenum == Table) *maxscore = maxsc;
1947 /*********************************************************************/
1948 /*********************************************************************/
1949 /********* This ends the scoring along all and all dimensions ********/
1950 /********* Now convert the raw scores that haven't been converted ****/
1951 /********* yet into likelihoods (i.e. the multidimensional scores.****/
1954 /***** Note: Can use same "init_class_prob" here, because since set number **/
1955 /*** of classes to 2, the 3rd class (Pdb-) will be ignored, so just get ****/
1956 /*** ratio of two to three stranded likelihoods, BOTH weighted by their ****/
1957 /*** init probaabiities. ********/
1959 if ( (Method != MultiCoil) && (Method !=PairCoilDiff))
1960 get_offsets(all_scores, sequence,all_offsets,number_tables);
1961 else if (!(mode & RAW_OUT)) {
1963 switch_gauss_param(Offset_to_Use, Method);
1966 /***********Do conversions on PairCoilDiff scores, using preprocessor. **/
1967 /*********** Note for differentiator there are only 2 possible classes **/
1968 /*********** (dimer and trimer). **********/
1969 if (Method == PairCoilDiff) {
1970 convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
1971 number_tables, differ_means_submatrix,
1972 differ_inverse_covars,differ_determinants, 2,
1973 init_class_prob, number_differ_score_dim,
1974 0, bound, Offset_to_Use,0);
1975 if (Preprocessor_Method == MultiCoil)
1976 preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
1977 else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
1980 /* Only use window scores that lie completely in coil */
1981 if ( (which_differentiator == 3) || (which_differentiator ==4))
1982 zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
1983 sequence.seqlen, bound);
1985 if (which_differentiator == 1)
1986 zero_out_bad_windows(differ_window_length,
1987 all_scores,all_preprocess_like[preproc_table],
1988 sequence.seqlen,bound);
1990 copy_offsets(all_preproc_offsets[preproc_table],sequence,
1991 all_offsets,2, all_scores, Offset_to_Use);
1994 if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
1995 for (tablenum=0; tablenum<2; tablenum++) {
1996 average_score_over_coils3(sequence, all_scores[tablenum+3],
1997 all_scores[tablenum],
1999 /** Total Coil Like. **/
2000 weighted_avg, offset_to_use,
2001 all_preproc_offsets[preproc_table],
2002 start_coil_at_reg_shift, bound,
2007 /**********Do conversions on MultiCoil dimension scores. ******/
2009 else if (Method == MultiCoil) {
2010 convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
2011 number_tables, means_submatrix,inverse_covars,
2012 determinants, number_classes,init_class_prob,
2014 mode & ONLY_COILED_CLASSES, bound,
2016 get_offsets(all_scores, sequence,all_offsets,3);
2020 if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
2021 for (tablenum=0; tablenum<3; tablenum++)
2022 average_score_over_coils3(sequence, all_scores[tablenum+3],
2023 all_scores[tablenum],
2024 all_scores[2+3], /** Total Coil Like. **/
2025 weighted_avg, offset_to_use,
2026 all_offsets[tablenum],
2027 start_coil_at_reg_shift, bound,
2033 } /* Ends if rescore_seq */
2036 /*****************************************************************************/
2037 if ( (Method == MultiCoil) || (Method == PairCoilDiff)){
2039 if (Method == MultiCoil)
2040 if ( (rescore_seq ==1 ) || /** Signals rescored MultiCoil **/
2041 (rescore_seq==-1)) /*Signals bound changed, need to redo seq_score */
2042 get_seq_scores(seq_scores, sequence, all_scores[0+3],
2043 all_scores[1+3],bound, weighted_avg);
2045 if (rescore_seq==-1 && !scored_already && (Method==PairCoilDiff)) {
2046 if ( (which_differentiator == 3) || (which_differentiator ==4))
2047 zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
2048 sequence.seqlen, bound);
2049 if (which_differentiator == 1) /* Only use window scores that lie */
2050 zero_out_bad_windows(differ_window_length, /* completely in coil */
2051 all_scores, all_preprocess_like[preproc_table],
2052 sequence.seqlen,bound);
2055 if ( (rescore_seq==-1) && Coil_Score && !scored_already) {
2056 if (Method == MultiCoil) {
2057 for (tablenum=0; tablenum<3; tablenum++) {
2058 average_score_over_coils3(sequence, all_scores[tablenum+3],
2059 all_scores[tablenum],
2060 all_scores[2+3], /** Total Coil Like. **/
2061 weighted_avg, offset_to_use,
2062 all_offsets[tablenum],
2063 start_coil_at_reg_shift, bound,
2065 if (mode & ONLY_COILED_CLASSES) /** Convert to olig. state. **/
2066 type_class_score_convert(sequence,all_scores,bound,
2067 mode & ONLY_COILED_CLASSES, number_tables,
2071 else if (Method == PairCoilDiff) {
2072 if (Preprocessor_Method == MultiCoil)
2073 preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
2074 else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
2077 for (tablenum=0; tablenum<2; tablenum++)
2078 average_score_over_coils3(sequence, all_scores[tablenum+3],
2079 all_scores[tablenum],
2081 /** Total Coil Like. **/
2082 weighted_avg, offset_to_use,
2083 all_offsets[tablenum],
2084 start_coil_at_reg_shift, bound,
2090 else if (rescore_seq == -2) /** Switch back to reg score.*/
2091 for (tablenum =0; tablenum <3 ; tablenum++)
2092 for (j =0; j <=POSNUM; j++)
2093 for (i=0; i<sequence.seqlen; i++) {
2094 all_scores[tablenum][i][j] =
2095 all_scores[tablenum+3][i][j];
2096 if ((Method == MultiCoil) &&
2097 (mode & ONLY_COILED_CLASSES) && (tablenum <2) &&
2098 (all_scores[2+3][i][j]>bound))
2099 all_scores[tablenum][i][j]/= all_scores[2+3][i][j];
2102 /****************************************************************************/
2107 return (double*)all_scores;
2114 ShowSeq(char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
2115 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
2116 double *m, double *b,
2117 double *m_single, double *b_single,
2118 int mode, int rescore_sequence, int rescore_preproc, double bound)
2122 /* fprintf(stderr,"\nIn show seq"); */
2125 if ( (rescore_sequence) || (rescore_preproc))
2126 all_scores = ScoreSeq(lib, differ_lib, multi_lib, m, b,m_single, b_single,
2128 method, preprocessor_method,
2129 offset_to_use, &maxscore,
2130 rescore_sequence, rescore_preproc, bound);
2133 /* fprintf(stderr,"\n Got all_scores in showseq"); */
2135 scores = &all_scores[table * MAXSEQLEN*(POSNUM+1)];
2136 /* This is all_scores[main_table][0][0] */
2137 preprocess_like = &all_preprocess_like[preproc_table][0][0];
2139 ShowScores(scores,&sequence,NULL, bound,mode);
2142 /*---------------------------------------------------------------------------*/
2147 OpenScreen(int mode, double log_bound, FILE *fin, FILE *flog,
2148 FILE* flog_coil_conflicts,FILE *fout_coils,
2150 double *m, double *b, double *m_single, double *b_single,
2151 char lib[MAXFUNCTNUM],
2152 char differ_lib[MAXFUNCTNUM],
2153 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
2155 int main_preprocessor_method,
2156 int *seqcnt, int which_priors, int which_priorp,
2157 double prior_freq_single[MAX_TABLE_NUMBER],
2158 double prior_freq_pair[MAX_TABLE_NUMBER],
2159 int good_turing_fixed_disc, int structural_pos[POSNUM+1])
2163 XSetWindowAttributes attr;
2165 static unsigned char dashl[2] = {1,1};
2166 static unsigned char dashd[2] = {2,2};
2171 display = XOpenDisplay("");
2173 /* Check mode to see if want GUI. */
2174 if (!display || (mode & NO_GUI)) {
2175 /** COmplete logfile, etc. if can't open xdisplay. */
2176 if (!(mode & NO_GUI))
2177 fprintf(stderr,"Sorry, bad display.\nMaybe your DISPLAY environment variable isn't set properly,\nor you haven't xhosted this machine.\n");
2178 if ( fout || ( (mode & PRN_MODE) || (mode & VER_MODE)) ) {
2179 fprintf(stderr,"\nThe log file will still be written by:\n\n");
2180 for (i = 0; motd[i][0]; ++i) /* Print out title. */
2181 fprintf(stderr," %s\n",motd[i]);
2183 finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
2185 m,b, m_single, b_single,
2186 lib,differ_lib,multi_lib, main_method, main_table,
2187 main_preprocessor_method, seqcnt, avg_max, which_priors,
2188 which_priorp, prior_freq_single, prior_freq_pair,
2189 good_turing_fixed_disc, structural_pos);
2193 if (flog) sclose(flog);
2194 if (flog_coil_conflicts) sclose(flog_coil_conflicts);
2195 if (fout) sclose(fout);
2196 if (ftotal_like) sclose(ftotal_like);
2202 font = XLoadFont(display,"9x15");
2203 window = XCreateSimpleWindow(display, RootWindow(display, 0),
2205 WIDTH, HEIGHT/4 - 20,
2207 BlackPixel(display, 0),
2208 WhitePixel(display, 0));
2210 help = XCreateSimpleWindow(display, RootWindow(display, 0),
2212 WIDTH / 5 * 4 +60, HEIGHT / 5 * 3 + 200,
2214 BlackPixel(display, 0),
2215 WhitePixel(display, 0));
2216 helpb = XCreateSimpleWindow(display, window,
2218 BlackPixel(display, 0), WhitePixel(display, 0));
2219 next = XCreateSimpleWindow(display, window,
2221 BlackPixel(display, 0), WhitePixel(display, 0));
2222 prev = XCreateSimpleWindow(display, window,
2224 BlackPixel(display, 0), WhitePixel(display, 0));
2225 autob = XCreateSimpleWindow(display, window,
2227 BlackPixel(display, 0), WhitePixel(display, 0));
2228 zoomb = XCreateSimpleWindow(display, window,
2230 BlackPixel(display, 0), WhitePixel(display, 0));
2231 printb = XCreateSimpleWindow(display, window,
2233 BlackPixel(display, 0), WhitePixel(display, 0));
2234 methodb=XCreateSimpleWindow(display, window,
2236 BlackPixel(display, 0), WhitePixel(display, 0));
2238 tableb=XCreateSimpleWindow(display, window,
2240 BlackPixel(display, 0), WhitePixel(display, 0));
2242 paramb=XCreateSimpleWindow(display, window,
2244 BlackPixel(display, 0), WhitePixel(display, 0));
2245 quit = XCreateSimpleWindow(display, window,
2247 BlackPixel(display, 0), WhitePixel(display, 0));
2249 /** To get size: **/
2250 /** Printing " trimer x.xx" takes 12 char = 12*7 =84 **/
2251 /** Each line takes 15 of y coord. giving 3*15 =45 + 10 for blank line **/
2252 seqscoreb=XCreateSimpleWindow(display, window,
2253 XCOORD_SEQSCORES, YCOORD_SEQSCORES,
2255 BlackPixel(display, 0), WhitePixel(display, 0));
2257 zoom = XCreateSimpleWindow(display, RootWindow(display, 0),
2261 BlackPixel(display, 0),
2262 WhitePixel(display, 0));
2263 XSelectInput(display, helpb, ButtonPressMask);
2264 XSelectInput(display, next, ButtonPressMask);
2265 XSelectInput(display, prev, ButtonPressMask);
2266 XSelectInput(display, autob, ButtonPressMask);
2267 XSelectInput(display, zoomb, ButtonPressMask);
2268 XSelectInput(display, printb, ButtonPressMask);
2269 XSelectInput(display, methodb, ButtonPressMask);
2271 XSelectInput(display, tableb, ButtonPressMask);
2273 XSelectInput(display, paramb, ButtonPressMask);
2275 XSelectInput(display, quit, ButtonPressMask);
2276 XSelectInput(display, window, ExposureMask);
2277 XSelectInput(display, zoom, ExposureMask | KeyPressMask | ButtonPressMask);
2278 XSelectInput(display, help, ExposureMask | KeyPressMask | ButtonPressMask);
2279 attr.bit_gravity = NorthWestGravity;
2280 XChangeWindowAttributes(display,window,CWBitGravity,&attr);
2281 XMapWindow(display, window);
2284 values.foreground = BlackPixel(display, 0);
2285 values.background = WhitePixel(display, 0);
2287 gc = XCreateGC(display, window, GCForeground | GCBackground | GCFont, &values);
2288 values.line_style = LineDoubleDash;
2289 values.font = XLoadFont(display,"6x13");
2290 gcgray = XCreateGC(display, window, GCForeground | GCBackground
2291 | GCLineStyle | GCFont, &values);
2292 XSetDashes(display, gcgray, 1, dashl, 2);
2293 values.function = GXxor;
2294 values.foreground = 1;
2295 gcdash = XCreateGC(display, window, GCForeground
2296 | GCLineStyle | GCFunction, &values);
2297 XSetDashes(display, gcdash, 0, dashd, 2);
2301 icon = XCreateBitmapFromData(display,window,bart_bits,bart_width,bart_height);*/
2302 /* XSetStandardProperties(display,window,"Coiled Coil Finder","xcoil",icon,NULL,0,NULL);*/
2303 /* XSetStandardProperties(display,zoom,"XCoil Zoomer","xcoil",icon,NULL,0,NULL);*/
2304 /* XSetStandardProperties(display,help,"XCoil Help","xcoil",icon,NULL,0,NULL);*/
2308 PSDefs (FILE *ps, int seqlen)
2313 int nominallength= 300; /* Default value for number residues per line. */
2315 if ( (seqlen%nominallength <= 50) && (seqlen >= nominallength)) {
2316 nominallength = 350;
2319 /*************************/
2320 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
2321 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
2322 nominallength = seqlen;
2324 nominallength = ps_res_per_line; /* How many res per line. */
2325 /************************/
2327 if (seqlen ==0) number_lines=1;
2329 number_lines = (int)(seqlen-1)/nominallength +1;
2330 /* nominallength amino acids per line. */
2331 /* C rounds down and we want 0-nominallength to go to 1 */
2332 /* nominallength+1 to 2*noiminallength to go to 2.. */
2337 /*****Scale and translate to fit the page if want it to rotate page. */
2338 /* unscaled there are 9 lines per page. */
2339 /* .89 scaling if 300 res per line.
2340 xscaler = .89 * 300/nominallength;
2341 fputs("90 rotate\n",ps);
2342 fputs("78 -60 translate\n",ps);
2343 if (number_lines > 9) yscaler = (double) 9/number_lines;
2344 fprintf(ps,"%lf %lf scale\n", xscaler, yscaler);
2347 /*****Scale and translate to fit the page if don't want to rotate page. */
2348 /* unscaled there can fit 12 lines per page, so 3600. */
2349 xscaler = .75 * 300/nominallength; /* .75 scaling if 300 res per line. */
2350 fputs("25 740 translate\n",ps);
2351 if (number_lines > 12) yscaler = (double) 12/number_lines;
2352 fprintf(ps,"%lf %lf scale\n", xscaler,yscaler);
2356 fputs("/black {0 setgray} def\n",ps);
2357 fputs("/gray {0.6 setgray} def\n",ps);
2358 fputs("/tt /Courier findfont 15 scalefont def\n",ps);
2359 fputs("/rm /Times-Roman findfont 12 scalefont def\n",ps);
2360 fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
2361 fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
2362 fputs("/m {neg moveto} def\n",ps);
2363 fputs("/l {neg lineto stroke} def\n",ps);
2364 fputs("/rl {neg rlineto} def\n",ps);
2365 fputs("/rl_stroke {neg rlineto stroke} def\n",ps);
2366 fputs("/bar {dup 0.0019 lt {pop} {0 exch rlineto stroke} ifelse} def\n",ps);
2369 /*** static int xmin = 60; ***** Made it global variable *****/
2370 /*Originally was 29 */
2371 /** xmin determines where in window start seq. **/
2372 static int nominallength;
2375 static int ydelta=57;
2376 /* If click at x,y, determine which residue position it corresponds to.
2379 XYtoPos (int x, int y)
2382 if (x<xmin) return 0;
2384 if (i>= nominallength) return 0;
2385 for (y -= ybase+ydelta; y>0; y-=ydelta, i+=nominallength);
2386 if (y<-30) return 0;
2387 if (i>=length) return 0;
2390 PostoXY (int i, int *x, int *y)
2393 *x = xmin + 2*(i%nominallength);
2394 *y = ybase+ydelta + ydelta*(i/nominallength);
2396 psstr (char *str, int len, FILE *ps)
2399 for (i=0; i<len; ++i)
2411 /*** ybase determines where top of the seq score box should be so don't ***/
2412 /*** print over the sequenc title. **/
2413 draw_seqscore_box(double bound, double seq_scores[MAX_CLASS_NUMBER],
2414 Window *seqscoreb, FILE *ps, int ybase)
2417 /* XWindowChanges values; */
2421 XClearWindow(display,*seqscoreb);
2423 /** Relocate seqscoreb base on ybase**/
2424 /**This is complex way to move it.****
2425 values.y= YCOORD_SEQSCORES+ybase-30;
2426 XConfigureWindow(display,*seqscoreb,CWY,&values);
2429 XMoveWindow(display,*seqscoreb,XCOORD_SEQSCORES, YCOORD_SEQSCORES+ybase-30);
2433 /** 3 means print 3 into the box horiz, 10 means print 10 into box vert? */
2435 XMapWindow(display, *seqscoreb);
2436 XDrawString(display,*seqscoreb,gcgray,3,10,"Seq Prob:",9);
2438 if (ps) { /* Draw the box. **/
2439 fprintf(ps,"%d %d m\n",10, 50+ybase-30);
2440 fprintf(ps,"0 40 rl\n");
2441 fprintf(ps,"79 0 rl\n");
2442 fprintf(ps,"0 -40 rl\n");
2443 fprintf(ps,"-79 0 rl_stroke\n");
2446 if (seq_scores[0] + seq_scores[1] + seq_scores[2] > 0) { /* Signals above */
2447 sprintf(buf," Dimer %5.2lf", seq_scores[0]); /* bound. */
2448 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2449 sprintf(buf," Trimer %5.2lf", seq_scores[1]);
2450 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2453 fprintf(ps,"rm setfont ");
2454 fprintf(ps,"%d %d m (Seq Prob for) show\n", 15,60 +ybase-30);
2456 fprintf(ps,"%d %d m (bound %5.2lf:) show\n", 15,70 +ybase-30,bound);
2458 fprintf(ps,"%d %d m ( Dimer %5.2lf) show\n",15,80 +ybase-30,
2460 fprintf(ps,"%d %d m ( Trimer %5.2lf) show\n",15,90 +ybase-30,
2467 sprintf(buf," No Residue ");
2468 XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2469 sprintf(buf," above %5.2lf ", bound);
2470 XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2473 fprintf(ps,"rm setfont ");
2474 fprintf(ps,"%d %d m (Seq Prob:) show\n", 15,60 +ybase-30);
2476 fprintf(ps,"%d %d m ( No Residue) show\n", 15, 75 +ybase-30);
2477 fprintf(ps,"%d %d m ( above %5.2lf) show\n",15,85 +ybase-30, bound);
2485 draw_bound_box(double bound, Window paramb)
2489 XMapWindow(display, paramb);
2490 XDrawString(display,paramb,gcgray,3,10,"bound",5);
2491 sprintf(buf,"%5.2lf", bound);
2492 XDrawImageString(display,paramb,gcgray,3,25,buf,strlen(buf));
2495 ShowScores (double *scores, Sequence *sequence, FILE *ps, double bound,
2498 int i,q, x, y0, xmax;
2503 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2504 char *MethodTitle = MethodTitleArray;
2505 char PreprocTitleArray[500]; /* For printing the preproc.name to ps file. */
2506 char *PreprocTitle = PreprocTitleArray;
2509 MethodTitleArray[0] = '\0'; /* Initialization. */
2510 PreprocTitleArray[0]='\0';
2512 nominallength = 300;
2513 length = sequence->seqlen;
2515 /* This makes display look nice (so don't have a line with < 50 residues. */
2516 if ((length%nominallength <= 50) && (length>= nominallength)){
2517 nominallength = 350;
2518 /* if (ps) fputs("-50 0 translate\n",ps); */
2522 /*************************/
2523 if (ps) /* Do this stuff only for postscript printout. */
2524 if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
2525 if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
2526 nominallength = sequence->seqlen;
2528 nominallength = ps_res_per_line; /* How many res per line. */
2529 /************************/
2532 charperline = (xmin+2*nominallength-5)/9; /*Each char is "9" units witdth */
2533 rows = (length-1)/nominallength+1; /* Number of rows to display. */
2535 xmax = 2*length+xmin; /* The max x position in current row. */
2537 xmax = 2*nominallength+xmin;
2538 width = xmin+2*nominallength+ WIDTH_ADD_ON; /* Width of window to display. */
2539 /* The +WIDTH_ADD_ON used to be +11, but */
2540 /* was cutting of some numbers. */
2541 x = strlen(sequence->title);
2542 XClearWindow(display,window);
2544 /**********This is where window size is adjusted to fit the sequence. ***/
2545 /********** using width and number of rows. **************************/
2547 XResizeWindow(display,window,width,
2548 38+15*((x+charperline-1)/charperline) +ydelta*rows);
2550 /*************************************************************************/
2551 XMapWindow(display, helpb);
2552 XMapWindow(display, next);
2553 XMapWindow(display, prev);
2554 XMapWindow(display, autob);
2555 XMapWindow(display, zoomb);
2556 XMapWindow(display, printb);
2557 XMapWindow(display, methodb);
2559 XMapWindow(display, tableb);
2560 XMapWindow(display, paramb); Now done in draw_bound_box()
2562 XMapWindow(display, quit);
2565 while (XCheckTypedEvent(display,Expose,&report));
2567 XDrawString(display,helpb,gcgray,3,10,"help",4);
2568 XDrawString(display,next,gcgray,3,10,"next",4);
2569 XDrawString(display,prev,gcgray,3,10,"prev",4);
2570 XDrawString(display,autob,gcgray,3,10,"auto",4);
2571 XDrawString(display,zoomb,gcgray,3,10,"zoom",4);
2572 XDrawString(display,printb,gcgray,3,10,"print",5);
2573 XDrawString(display,methodb,gcgray,3,10,"method",6);
2575 /* XDrawString(display,tableb,gcgray,3,10,"table",5); */
2577 draw_bound_box(bound, paramb);
2578 /** draw_seqscore_box(bound, seq_scores,?? &seqscoreb, ps, ybase);
2579 MOVED TO LATER (AFTER ybase gets adjusted at end of this function). ****/
2582 XDrawString(display,quit,gcgray,3,10,"quit",4);
2585 /************This section prints out the method name..... ********************/
2586 /**** Add 9 to xcoord for each character. *****/
2587 XDrawString(display,window,gc, XCOORD_SCORED_BY,15,"Scored by ",10);
2589 make_methodtitle(MethodTitle, method, table, mode);
2590 XDrawString(display,window,gc,XCOORD_SCORED_BY +90,15,MethodTitle,
2591 strlen(MethodTitle));
2593 XDrawString(display,window,gc, XCOORD_SCORED_BY,30,"Filter by ",10);
2594 make_methodtitle(PreprocTitle, preprocessor_method, preproc_table,
2596 XDrawString(display,window,gc,XCOORD_SCORED_BY + 90,30,
2597 PreprocTitle,strlen(PreprocTitle));
2599 /***************************************************************************/
2601 XDrawString(display,window,gc, 5,30,sequence->code,strlen(sequence->code));
2602 for (i=0,ybase=30; x > i + charperline; i += charperline)
2603 XDrawString(display,window,gc, 5,ybase+=15,&(sequence->title[i]),charperline);
2604 XDrawString(display,window,gc, 5,ybase+15,&(sequence->title[i]),x-i);
2607 fprintf(ps,"tt setfont 5 15 m (%s) show\n",sequence->code);
2608 /* fprintf(ps,"175 15 m (Viewgram by XCoil) show\n"); */
2610 fprintf(ps,"430 15 m (Scored by %s) show\n",methodname[method]);
2612 fprintf(ps,"180 15 m (Scored by %s) show\n",MethodTitle);
2613 fprintf(ps,"180 26 m (Filter by %s) show\n",PreprocTitle);
2616 for (i=0,ybase=30; x > i + charperline; i += charperline) {
2617 fprintf(ps,"5 %d m (",ybase+=15);
2618 psstr(&(sequence->title[i]),charperline,ps);
2619 fprintf(ps,") show\n");
2621 fprintf(ps,"5 %d m (",ybase+15);
2622 psstr(&(sequence->title[i]),x-i,ps);
2623 fprintf(ps,") show\n");
2624 fputs("rm setfont 2 setlinewidth\n",ps);
2626 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2628 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2629 y = y0-30*scores[8*i+7];
2630 XDrawLine(display, window, gc, x, y0, x, (int)(y+0.5));
2631 XDrawLine(display, window, gc, x+1, y0, x+1, (int)(y+0.5));
2633 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2636 /********************* If in pos mode draws in the correct likelihoods every
2637 now and then? ******************
2638 if (sequence->reg && method==HMMCoil) {
2641 for (i=0,y0=ybase; i<length; ++i,x+=2) {
2642 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2643 if (sequence->reg[i] != '.') {
2644 double y = y0-30*scores[8*i+sequence->reg[i]];
2645 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2646 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2648 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2651 if (ps) fputs("black\n",ps);
2653 **********************/
2655 /**********Now draw in the PreProcScore in gray every other line. ****/
2659 for (i=0,y0=ybase; i<length; i+=2,x+=4) { /* Normally is ++i,x+=4, but */
2660 double y; /* we do every other line here. */
2661 if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2663 y = y0-30*preprocess_like[i*(POSNUM +1) + 7];
2664 /* This is preprocess_like[i][7]. */
2665 XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2666 XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2668 fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2671 if (ps) fputs("black\n",ps);
2673 /*************************************************************************/
2675 if (ps) fputs("0 setlinewidth\n",ps);
2676 for (i=9,y0=ybase; i<length; i+=10) {
2677 if (!((i-9)%nominallength)) {x=xmin-1; y0 += ydelta;}
2679 XDrawLine(display, window, gc, x,y0,x,y0+2); /* Draw ticks every */
2680 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+2); /* 10 residues */
2683 for (i=49,y0=ybase; i<length; i+=50) {
2685 if (!((i-49)%nominallength)) {x=xmin-1; y0 += ydelta;}
2687 XDrawLine(display, window, gc, x,y0,x,y0+5);
2688 if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+5);
2689 sprintf(buff,"%d",i+1);
2690 XDrawString(display, window, gcgray, x-7, y0+18, buff, strlen(buff));
2691 if (ps) fprintf(ps,"%d %d m (%d) cshow\n",x,y0+18,i+1);
2693 for (i=1,y0=ybase+ydelta; i<=rows; ++i,y0+= ydelta) {
2695 xmax = 2*(length-(rows-1)*nominallength)+xmin; /* Max x position in */
2697 XDrawString(display,window,gcgray,xmin-25,y0-26,"100%",4);
2698 XDrawString(display,window,gcgray,xmin-19,y0-11,"50%",3);
2699 XDrawString(display,window,gcgray,xmin-13,y0+4,"0%",2);
2700 XDrawLine(display,window,gcgray,xmin,y0-30,xmax-1,y0-30);
2701 XDrawLine(display,window,gcgray,xmin,y0-15,xmax-1,y0-15);
2703 fprintf(ps,"%d %d m (100%%) rshow\n",xmin-1,y0-26);
2704 fprintf(ps,"%d %d m (50%%) rshow\n",xmin-1,y0-11);
2705 fprintf(ps,"%d %d m (0%%) rshow\n",xmin-1,y0+4);
2706 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-30,xmax,y0-30);
2707 fprintf(ps,"%d %d m %d %d l\n",xmin,y0-15,xmax,y0-15);
2708 fprintf(ps,"%d %d m %d %d l\n",xmin,y0,xmax,y0);
2713 draw_seqscore_box(bound, seq_scores, &seqscoreb, ps, ybase);
2715 if (zoomptr) ZoomPnt();
2718 fputs("showpage\n",ps);
2727 make_methodtitle(char *MethodTitle, int Method, int which_table,
2735 if (which_table ==0) like2or3 =2;
2739 /************This section prints out the method name..... ********************/
2740 /**** Add 9 to xcoord for each character. *****/
2744 if ((Method == PairCoilDiff) || (Method == PairCoilDiffAvg)) {
2745 /* Does Differen. use dimer or trimer like. */
2747 /* XDrawString(display,window,gc,xcoord,15,"Dimer",5); */
2748 MethodTitle= strcat(MethodTitle, "Dimer");
2751 else if (like2or3==3) {
2752 /* XDrawString(display,window,gc,xcoord,15,"Trimer",6); */
2753 MethodTitle= strcat(MethodTitle, "Trimer");
2759 if (Method == PairCoil) {
2760 if (mode & PAIRCOIL_PAIRS) { /* Are we using pair or single prob */
2761 /* XDrawString(display,window,gc,xcoord,15,"Pair",4); */
2762 MethodTitle= strcat(MethodTitle, "Pair");
2766 /* XDrawString(display,window,gc,xcoord,15,"Single",6); */
2767 MethodTitle= strcat(MethodTitle, "Single");
2771 } /* Ends the stuff for dealing with PairCoil, TDRes, TDCoil. */
2774 /*****************************************************************************/
2775 /*************************** Now write the methodname. *********************/
2777 /* XDrawString(display,window,gc,xcoord,15,methodname[method],strlen(methodname[method])); */
2778 xcoord += 9*strlen(methodname[Method]);
2779 MethodTitle= strcat(MethodTitle, methodname[Method]);
2781 if (Method == MultiCoil) {
2782 if (mode & ONLY_COILED_CLASSES)
2783 strcat(MethodTitle,"Olig");
2787 strcat(MethodTitle,"WtdAvg");
2788 else strcat(MethodTitle,"Avg");
2790 else if (Coil_Score ==2)
2791 strcat(MethodTitle, "MAX");
2793 if (start_coil_at_reg_shift)
2794 strcat(MethodTitle, "Shift");
2796 strcat(MethodTitle, "NoShifts");
2800 /* Now print Table are using. */
2801 /* And if are using a particular offset, print that. */
2802 if ( (Method==MultiCoil) || (Method==PairCoil) || (Method == HMMCoil)
2803 || (Method==PairCoilDiff) || (Method==PairCoilDiffAvg)) {
2806 if (number_tables >1) {
2807 sprintf(buf,"%d",which_table+1);
2808 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2809 xcoord += 9*strlen(buf);
2810 MethodTitle= strcat(MethodTitle, buf);
2814 sprintf(buf,""); /* Clear out buf. */
2817 if ((offset_to_use != -1) && (offset_to_use != 7))
2818 sprintf(buf,"%c",'a' + offset_to_use);
2819 else if (Method == MultiCoil)
2820 if (offset_to_use == -1) {
2821 strcat(MethodTitle,"Max");
2824 else if (offset_to_use == 7) {
2825 strcat(MethodTitle,"Combo");
2829 /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2830 xcoord += 9*strlen(buf);
2831 MethodTitle= strcat(MethodTitle, buf);
2835 /***************************************************************************/
2838 int is_not_differentiator(int Method)
2840 if ( (Method == MultiCoil) || (Method == PairCoil) || (Method==HMMCoil) ||
2841 (Method==NEWCOILS) || (Method == ActualCoil) )
2850 log_file_header(FILE *flog, FILE *flog_coil_conflicts,
2851 int mode, int argc, char *argv[],
2852 int avg_max, double bound,
2853 int preprocessor_method, int preproc_table,
2854 int log_offset_to_use, int start_coil_at_reg_shift,
2855 int weighted_avg, int main_method, int main_table)
2858 char MethodTitleArray[500]; /* For printing the methodname to ps file. */
2859 char *MethodTitle = MethodTitleArray;
2861 MethodTitleArray[0] = '\0'; /* Initialization. */
2862 make_methodtitle(MethodTitle, method,main_table,mode);
2865 fprintf(flog,"\nInput file is %s", argc>1?argv[1]:"-");
2866 fprintf(flog, "\n\nMethod Used is %s", MethodTitle);
2867 fprintf(flog, "\nThe bound for the coil locater is %.2lf",bound);
2869 if (flog_coil_conflicts) {
2870 fprintf(flog_coil_conflicts,"\nInput file is %s", argc>1?argv[1]:"-");
2871 if (mode & TST_MODE0)
2872 fprintf(flog_coil_conflicts,"\nSubtracting sequence off from cctb before scoring it.");
2873 if (mode & TST_MODE1)
2874 fprintf(flog_coil_conflicts,"\nSubtracting sequence off from trimer table before scoring it.");
2878 if (mode & WEIGHTED_PROBS) {
2879 if (flog_coil_conflicts)
2880 fprintf(flog_coil_conflicts,"\nStatistically weighted PairCoil score is ");
2882 fprintf(flog,"\nStatistically weighted PairCoil score is ");
2885 if (flog_coil_conflicts)
2886 fprintf(flog_coil_conflicts,"\nPairCoil score is ");
2888 fprintf(flog,"\nPairCoil score is ");
2891 if (flog_coil_conflicts)
2892 fprintf(flog_coil_conflicts,"AVERAGE of residue scores.");
2894 fprintf(flog,"AVERAGE of residue scores.");
2897 if (flog_coil_conflicts)
2898 fprintf(flog_coil_conflicts,"MAX of residue scores.");
2900 fprintf(flog,"MAX of residue scores.");
2903 if (mode & USE_LIKE_LINE) {
2904 if (flog_coil_conflicts)
2905 fprintf(flog_coil_conflicts,"\nPaircoil used the likelihood line.");
2907 fprintf(flog,"\nPaircoil used the likelihood line.");
2910 if (flog_coil_conflicts)
2911 fprintf(flog_coil_conflicts,"\nPaircoil used direct prob. formula instead of likelihood line.");
2913 fprintf(flog,"\nPaircoil used direct prob. formula instead of likelihood line.");
2916 if (flog_coil_conflicts)
2917 fprintf(flog_coil_conflicts, "\nThe distances used by PairCoil are: ");
2919 fprintf(flog, "\nThe distances used by PairCoil are: ");
2921 for (i=0; i<pair_functnum; i++) {
2922 if (flog_coil_conflicts)
2923 fprintf(flog_coil_conflicts,"%d ",i);
2925 fprintf(flog,"%d ",i);
2928 if (flog_coil_conflicts) {
2929 fprintf(flog_coil_conflicts, "\n\nMethod Used is %s", MethodTitle);
2930 fprintf(flog_coil_conflicts, "\nThe bound for the coil locater is %.2lf",bound);
2933 if ( (main_method==PairCoilDiff) || (main_method==PairCoilDiffAvg)) {
2934 char MethodTitleArray[500]; /* For printing the methodname. */
2935 char *MethodTitle = MethodTitleArray;
2936 int oligimerization;
2938 MethodTitleArray[0] = '\0'; /* Initialization. */
2939 make_methodtitle(MethodTitle, preprocessor_method, preproc_table,mode);
2942 oligimerization == 2;
2943 else oligimerization ==3;
2946 fprintf(flog, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2947 if (flog_coil_conflicts)
2948 fprintf(flog_coil_conflicts, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2949 /* Whether logfile is for 2 or */
2950 /* 3 stranded predictions. */
2952 if (main_method==PairCoilDiffAvg) {
2953 if (weighted_avg == 0) {
2955 fprintf(flog, "\n Averaging over coils was done by coil length.");
2956 if (flog_coil_conflicts)
2957 fprintf(flog_coil_conflicts, "\n Averaging over coils was done by coil length.");
2959 else { /* weighted_avg ==1 */
2961 fprintf(flog, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
2962 if (flog_coil_conflicts)
2963 fprintf(flog_coil_conflicts, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
2966 if (start_coil_at_reg_shift) {
2968 fprintf(flog,"\n New coils were considered to start at register shifts.");
2969 if (flog_coil_conflicts)
2970 fprintf(flog_coil_conflicts,"\n New coils were considered to start at register shifts.");
2974 fprintf(flog,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2975 if (flog_coil_conflicts)
2976 fprintf(flog_coil_conflicts,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2982 fprintf(flog,"\nMain Table is Table %d\n\n",main_table);
2984 if (mode & POS_MODE)
2985 fprintf(flog,"\nIf the start or end of a coil is listed inside (), then that position does not occur in the pos file.\n");
2993 coil_tuple_output(int number_tables, char *lib, int mode, int avg_max,
2994 double bound, Sequence sequence, FILE *fout_coils,
2995 int number_multi_lib[MAX_TABLE_NUMBER],
2996 int number_score_dim)
2998 char offsets[MAXSEQLEN];
3002 extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
3003 extern double *all_scores;
3006 for (tab=0; tab< number_tables; tab++) {
3007 for (dist=0; dist<number_multi_lib[tab]; dist++) {
3009 /*************************************************************************/
3010 /**********This computes offset of that table/dist method of scoring. ***/
3012 for (i=0; i<sequence.seqlen; i++) {
3013 for (j=0; j<POSNUM; j++)
3014 if (all_preprocess_like[tab + number_tables*(int)dist][i][j] ==
3015 all_preprocess_like[tab + number_tables*(int)dist][i][7]) {
3016 offsets[i] = (j-i) % POSNUM;
3017 if (offsets[i]<0) offsets[i] += POSNUM;
3019 /* This is a new thing so that if an offset change is not real */
3020 /* we keep the old offset. **********/
3021 if ( (i>0) && (offsets[i-1] != -1) &&
3022 (offsets[i] != offsets[i-1]) &&
3023 (all_preprocess_like[tab + number_tables*(int)dist]
3024 [i][(i+offsets[i-1]) %POSNUM]
3025 == all_preprocess_like[tab+number_tables*(int)dist][i][7]) )
3026 offsets[i] = offsets[i-1];
3029 /**************************************************************************/
3030 /****** Now compute coil scores for that method (table,dist) pair. ***/
3031 get_raw_coil_score(sequence,
3032 all_preprocess_like[tab+number_tables*(int)dist],
3033 all_scores[tab+number_tables*(int)dist], mode,
3034 avg_max,bound,offsets);
3038 /************* Now output the tuples of coil scores. ****/
3039 tuple_output2(sequence, mode, fout_coils,
3040 all_scores, number_tables,
3041 number_multi_lib,bound, number_score_dim);
3046 switch_gauss_param(int Offset_to_Use, int Method)
3048 if (Method == MultiCoil) {
3049 if (Offset_to_Use == 7) { /* Combo offset */
3050 determinants= &both_determinants[0][0];
3051 inverse_covars= &both_inverse_covars[0][0][0][0];
3052 means_submatrix= &both_means_submatrix[0][0][0];
3054 else { /* max offset or some fixed offset */
3055 determinants= &both_determinants[1][0];
3056 inverse_covars= &both_inverse_covars[1][0][0][0];
3057 means_submatrix= &both_means_submatrix[1][0][0];
3060 else if (Method == PairCoilDiff) {
3061 determinants= &differ_determinants[0];
3062 inverse_covars= &differ_inverse_covars[0][0][0];
3063 means_submatrix= &differ_means_submatrix[0][0];
3070 void usage(char *prog_command)
3075 /** Strip off directories to get just the program name. **/
3076 i=strlen(prog_command);
3077 while ((i>=0) && (prog_command[i] != '/') )
3079 prog_name= &prog_command[i+1];
3081 fprintf(stderr,"\nUsage: %s <file-name> [dimer, trimer, dimer-1, trimer-1, negatives, config_file] [multi_lib_dim1] [multi_lib_dim2] \n",prog_name);