--- /dev/null
+/* Code by Ethan Wolf 1995. Graphics code by David Wilson. */
+/*** See dimer_trimer_like.c, function calculate_likelihood to see how */
+/*** to compute likelihood line number before call DimerTrimerScore() */
+
+#define XCOORD_SCORED_BY 430
+#define XCOORD_SEQSCORES 0
+#define YCOORD_SEQSCORES 55 /* Each char takes "15" so put us below printing */
+#define WIDTH_ADD_ON 30 /** Empty space after seq in main window. **/
+#define xmin 120 /* Want this far enough over so have room for seqscoreb */
+ /** seqscoreb is 108 wide. xmin used to be 60 **/
+
+#include "sc.h"
+#include <stdio.h>
+#include "sc2seq.h" /* This also includes scconst.h, sc.h, interface.h */
+#include <X11/Xlib.h>
+#include <X11/keysym.h>
+#include <math.h>
+#include "scio.h"
+#include "options.h"
+
+Display *display;
+Window window, zoom, printb, quit, next, prev, autob, zoomb, methodb, paramb, helpb, help, seqscoreb;
+GC gc, gcgray, gcdash;
+
+
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/time.h>
+#include <string.h>
+#include "config.h"
+#include "switches.h"
+
+#include "hmm.h"
+
+
+/* FILE *debug_out; */
+
+
+/*************************** MULTICOIL LIKELIHOOD STUFF. *********************/
+double both_determinants[2][MAX_CLASS_NUMBER];
+double *determinants= &both_determinants[0][0];
+
+double both_inverse_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+double *inverse_covars= &both_inverse_covars[0][0][0][0];
+
+double both_covars_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+
+double both_means_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
+double *means_submatrix= &both_means_submatrix[0][0][0];
+
+double init_class_prob[MAX_CLASS_NUMBER];
+
+int multi_functnum[MAX_TABLE_NUMBER];
+int number_score_dim=0;
+int number_classes=0;
+int number_multi_lib[MAX_TABLE_NUMBER];
+ /* Used to detrmine the number of different library */
+ /* distance values to do gaussian param over. */
+
+#ifdef TEST_VERSION
+ FILE *ftotal_like=NULL; /* Extern for now. ***/
+ double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER];
+ /* When ftotal_like is read from get_defaults(), the */
+ /* i,j entry will hold the total gaussian like when */
+ /* assuming dimers have init_prob = */
+ /* MIN_DIMER_PROB + i*DIMER_GRID_STEP */
+ /* and trimers have init_prob = */
+ /* MIN_TRIMER_PROB + j*TRIMER_GRID_STEP */
+ /* These constants are defined in scconst.h. */
+ /* The value is updated during "output_seq()". */
+
+ /*************************** PAIRDIFF LIKELIHOOD STUFF. *********************/
+ double differ_determinants[MAX_CLASS_NUMBER];
+ double differ_inverse_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ double differ_covars_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ double differ_means_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
+
+ int number_differ_score_dim =0;
+ int number_differ_lib=0;
+ int differ_window_length = 28;
+ int which_differentiator=1; /** Default is just the score for the window **/
+ /** 2 means to give each residue a pair of scores**/
+ /** corresponding to the max and min of the **/
+ /** windows containing it. **/
+ int differ_functnum=0;
+#endif
+
+/************************** General Scoring parameters. ********************/
+int functnum=0, pair_functnum=0;
+int preproc_table=0;
+double log_bound = 0; /*In PRN_MODE only scores > bound output to logfile.*/
+
+int WINDOW = PAIRWINDOW;
+int window_length[MAX_TABLE_NUMBER];
+
+double SCALE0 = DEFAULT_SCALE0;
+double scale0s[MAX_TABLE_NUMBER]; /** Both used as bonnie's scale0 and **/
+double scale0p[MAX_TABLE_NUMBER]; /** in computing posterior prob. from **/
+ /** prior as sigma^2(prior)=scale0*sigma^2 */
+ /** where these are the standard deviation *
+ /** of the n*(prior_prob) and the distib */
+ /** of number of hits when sampling */
+
+int log_offset_to_use;
+int offset_to_use=-1; /* Default offset in DimerTrimerScore is all. */
+int avg_max=1; /* Default means use max coil scores in log file. */
+ /** Value of 0 means do average, 1 means do max. **/
+ /** 2 means do both. */
+int number_tables =0;
+
+int method = MultiCoil;
+int main_method = MultiCoil;
+int main_preprocessor_method = MultiCoil;
+int preprocessor_method = MultiCoil;
+
+int main_table = 0;
+int table = 0;
+
+/*********************** Parameters involved in Coil Scores ***********/
+int by_coil_or_seq =0; /*** To determe what type of coil_out file to do. **/
+int weighted_avg = 0; /* Determines which method for determining coil */
+ /* score will be used in DimerTrimerScore. */
+int start_coil_at_reg_shift =0; /* For averaging scores over coil length. */
+int Coil_Score =0; /* 0 means for MultiCoil do residues score, 1 means */
+ /* average scores over coils, 2 means take max coil */
+ /* score. */
+
+
+Sequence sequence;
+
+char *methodname[] = {"MultiCoil", "Coil","NEWCOILS", "HMMCoil","ActualCoil",
+ "PairCoilDiff",
+ "PairCoilDiffAvg"}; /* "Coil" is PairCoil or */
+
+
+
+char table_diff_stat_filename[MAXLINE]="/dev/null";
+char weightfile_outname[MAXLINE]="/dev/null";
+
+
+
+
+
+/***/
+ void *hmm;
+/***/
+
+double seq_scores[MAX_CLASS_NUMBER];
+double *all_scores;
+double *scores; /** The space these pointers points to is alocated by **/
+ /** using static variable arrays in ScoreSeq */
+ /** all_scores is scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM +1]*/
+ /** scores is scores[current_score_algor][MAXSEQLEN][POSNUM+1] **/
+
+double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
+ /** Since preprocess_likelihood is used everywhere, just */
+ /** make it an external variable. */
+ /** This is given a value in ScoreSeq and sometime in */
+ /** ActualCoil(). */
+double *preprocess_like;
+ /* This will point to the current table version */
+ /* of all_preprocess_like */
+char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
+
+
+
+/************************ XWindow stuff. ***********************************/
+int zoomptr=0;
+int ps_res_per_line=0; /* 0 is a flag that it wasn't set in get_defaults */
+ /* so should just print it using nominal_length */
+ /* as done on the screen. If it is set to -1 */
+ /* then should print it all on one line. */
+
+
+char *motd[] = {
+
+ " *********************************",
+ " ** **",
+ " ** MultiCoil (version 1.0) **",
+ " ** **",
+ " ** **",
+ " *********************************",
+ "",
+ "MultiCoil is described in:",
+ " Ethan Wolf, Peter Kim, Bonnie Berger,",
+ " `MultiCoil: A Program for Predicting Two- and
+ " Three-stranded Coiled Coils'.",
+ " ",
+ "Email inquiries to: multicoil-help@theory.lcs.mit.edu.",
+ ""
+
+};
+
+char *hmotd[] = {
+ "help: Clicking on help displays this message.",
+ " ",
+ "next: Displays the the next sequence.",
+ " Same as pressing 'n' or down-arrow.",
+ " ",
+ "prev: Displays the previous sequence.",
+ " Same as pressing 'p' or up-arrow.",
+ " ",
+ "auto: Automatically steps through the sequences",
+ " of a file.",
+ " ",
+ "zoom: Gives more information on the highest",
+ " scoring residue. Same as pressing 'z'.",
+ " Clicking on a specific region will give",
+ " more information on that region.",
+ " ",
+ "print: Prints the sequence.",
+ " ",
+/**********
+ "table: Switches the table used by PairCoil and",
+ " by TableDiff methods, same as 't'. Note that",
+ " key 'i' (input_table) changes the preprocessor's table.",
+ " ",
+**********/
+ "method: Switches the scoring method, same as 'm'.",
+ " 'w' swithces Which preprocessor is used. ",
+ " ",
+ "key 'Up' or 'Down' moves the bound up or down by .01",
+ "key 'Return' rescores the sequence after changing the bound",
+ " ",
+ "key 's' or 'l': switch between singles and pair",
+ " probabilities at lib distances.",
+ " ",
+ "key 'o': Calc paircoil class prob. with respect only to",
+ " positive oligimerization classes (ignore pdb- class)",
+ "key '7': MultiCoil max score is obtained by taking tuple",
+ " of each dimension's max score (combo register).",
+ "key 'x': MultiCoil computes scores for each offset, taking",
+ " max of those scores to be the max score (max reg).",
+ "key 'r': Cycle through offsets 0-7 for register (7=combo reg)",
+ "",
+ "key 'c': Switch between MultiCoil register score, average coil",
+ " score, and max coil score",
+ "key 'b': For Coil score, toggle whether start new coils at",
+ " a register shift (break in structure).",
+ "key 'a': For Coil score, toggle between unweighted and weighted avg",
+ "",
+ "quit: Quits the program and completes log file.",
+ " 'q' quits without completing the log.",
+ " Pressing 'q' in any window also closes ",
+ " that window.",
+ ""
+
+};
+
+
+ShowTitle (Window win)
+{
+ int i;
+
+ for (i = 0; motd[i][0]; ++i)
+ XDrawString(display, win, gc, 5, 17*i+20, motd[i], strlen(motd[i]));
+ XDrawString(display,win,gcgray,100,17*i+20,
+ "Press Any Key To Continue",25);
+}
+
+ShowHelp (Window win)
+{
+ int i;
+
+ /* Places hmotd text in win */
+ for (i=0; hmotd[i][0]; ++i)
+ XDrawString(display,win,gc,5,17*i+20, hmotd[i],strlen(hmotd[i]));
+
+}
+
+
+static int seqnum=0, sequp=0, seqlow=1;
+
+char autoadvance=0;
+main (int argc, char *argv[])
+{
+ struct itimerval interval;
+ fd_set fds, fdnone;
+ FILE *fin, *ps;
+ XEvent report;
+ char buff[MAXLINE+MAXSEQLEN];
+
+
+
+/* Note that fout_raw[MAX_TABLE_NUMBER+1] and fout_coils[MAX_TABLE_NUMBER+1] */
+/* will be files of tuples of each residue's scores according to all tables. */
+/* Never mind... **/
+ char *command_line = NULL;
+
+ FILE *fgin=NULL, *fout=NULL,
+ *flog=NULL, *flog_coil_conflicts=NULL, *fout_coils=NULL;
+ FILE *fpin[MAX_TABLE_NUMBER];
+ FILE *pir=NULL;
+ char pir_name[MAXLINE];
+ char likelihoods[MAX_TABLE_NUMBER][MAXLINE];
+ /* Files of likelihood lines. */
+ char printfile[500]; /* The filename to output postscript print to. */
+ char *print= printfile;
+
+ extern int number_classes;
+
+ char class_sc_filenames[2][MAX_CLASS_NUMBER][MAXLINE];
+ double both_class_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM]
+ double both_class_means[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
+ [MAX_NUM_SCORE_DIM];
+ char gauss_param[2][MAXLINE];
+ /** 0 entry is for combo offset, 1 entry is for **/
+ /** max of other registers. ***/
+ extern int functnum; /* What is used in scscore.c. Takes value from */
+ /* pair_functnum and multi_functnum. */
+
+ char lib[MAXFUNCTNUM]; /* Use MAXFUNCTUM=7 here since are 7 distances */
+ /* in coil, so we know AT MOST 7 lib functions. */
+ extern int pair_functnum; /* The number of distances in lib[].*/
+
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM];
+ extern int multi_functnum[MAX_TABLE_NUMBER];
+
+#ifdef TEST_VERSION
+ char differ_class_sc_filenames[MAX_CLASS_NUMBER][MAXLINE];
+ double differ_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
+ double differ_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ char differ_gauss_param[MAXLINE];
+ char differ_lib[MAXFUNCTNUM]; /* Distance values for PairCoilDiff */
+ extern int differ_functnum;
+#ifdef TEST_VERSION
+
+
+
+ int acnt = 1;
+ int mode = 0; /* If the i-th option is set, mode gets 1 in i-th bit. */
+ double bound =0;
+
+ double upper_bound; /* In DEBUG_MODE (for printing out the predicted */
+ /* heptad repeat), scores in [bound,upper_bound] */
+ /* are printed in capital letters. */
+
+ double m[MAX_TABLE_NUMBER],
+ b[MAX_TABLE_NUMBER]; /*The likelihood line is mx+b for paircoil. */
+ double m_single[MAX_TABLE_NUMBER],
+ b_single[MAX_TABLE_NUMBER]; /* Like line for singlecoil. */
+ int i, class;
+
+ int table_to_remove_from = -1;
+ /* Used if want to remove ALL input seq */
+ /* from that table at outset. **/
+
+ int old_num_dim_table[2];
+ int new_num_dim_table[2];
+ int old_num_dim_differ =0;
+ /* The gauss param file for multiple distances is */
+ /* contains the parameters for all distances. */
+
+
+#ifdef TEST_VERSION
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ prior_freq_pair[MAX_TABLE_NUMBER];
+ int which_priors =0, which_priorp =0;
+ int good_turing_fixed_disc =0; /** If set to 1 do the good_turing_fixed_disc estimate of probs*/
+ int structural_pos[POSNUM+1]; /** Holds list of which positions are **/
+ /** structurally important. All others **/
+ /** get genbnk probabilities.... **/
+ /** Value -1 indicates no more entries. **/
+#endif
+
+
+/*** INITIALIZATION ***/
+
+ fprintf(stderr,"\nAbout to initialize");
+ initialize(&number_classes, &number_multi_lib[0], window_length, scale0s,
+ scale0p,fpin[0],
+ &likelihoods[0][0], &pir_name[0], print, gauss_param,
+ differ_gauss_param,
+ &mode, prior_freq_single, prior_freq_pair,
+ structural_pos);
+
+ fprintf(stderr,"\nInitialized");
+ read_command_line(argc, argv, &command_line, multi_functnum,
+ &mode, number_multi_lib, multi_lib);
+ fprintf(stderr,"\n\nRead commandline, command_line=%s",command_line);
+
+/* Get the default options in the .paircoil file!!! */
+ get_defaults(argv[1], &mode, &log_bound, &fgin, fpin,
+ &ftotal_like, &fout,
+ &flog, &flog_coil_conflicts,
+ &fout_coils, &by_coil_or_seq, likelihoods,
+ pir_name, lib, differ_lib, multi_lib,
+ &pair_functnum, &differ_functnum, multi_functnum,
+ print, &main_method, &main_preprocessor_method,
+ &main_table, &number_tables, &offset_to_use,
+ &avg_max, &weighted_avg, &Coil_Score,
+ &start_coil_at_reg_shift,
+ table_diff_stat_filename, weightfile_outname,
+ &ps_res_per_line, &number_classes,
+ class_sc_filenames, differ_class_sc_filenames,
+ gauss_param, differ_gauss_param, number_multi_lib,
+ &number_differ_lib,
+ init_class_prob,
+ &table_to_remove_from,
+ command_line, window_length, scale0s, scale0p,
+ prior_freq_single, prior_freq_pair, &which_priors,
+ &which_priorp, &good_turing_fixed_disc,
+ structural_pos, &which_differentiator, &differ_window_length);
+
+ fprintf(stderr,"\nOut of get_def");
+
+
+ bound = log_bound;
+
+#ifdef TEST_VERSION
+ compute_multivariate_num_differ_dim(&number_differ_score_dim,
+ number_differ_lib,
+ differ_functnum, &old_num_dim_differ,
+ which_differentiator);
+#endif
+
+
+ compute_multivariate_num_dim(mode, &number_score_dim,
+ number_multi_lib,multi_functnum, number_tables,
+ new_num_dim_table, old_num_dim_table);
+
+ if (number_classes >0) { /*** Get gaussian param for the classes. */
+ if (ftotal_like) { /** Want to output total gauss like over classes*/
+ for (i=0; i<number_classes; i++) /* for various init_prob settings. */
+ init_class_prob[i]=1; /** All 1 values is flag to just output */
+ init_totals(total_gauss_like); /** gauss values, not like values in */
+ } /** convert_raw_scores_to_gauss_like. */
+
+
+ for (i=0; i<2; i++)
+ if (gauss_param[i][0] != ',') {
+ get_gauss_param_all_classes2(class_sc_filenames[i],
+ old_num_dim_table[0]+old_num_dim_table[1],
+ both_class_means[i], both_class_covars[i],
+ number_classes, gauss_param[i]);
+
+
+/******Note we store the submatrices in the original matrices....****/
+ compute_submatrices2(both_class_means[i],both_means_submatrix[i],
+ both_class_covars[i],number_classes,
+ both_inverse_covars[i],
+ both_determinants[i], multi_lib,
+ old_num_dim_table[0], old_num_dim_table[1],
+ new_num_dim_table[0], new_num_dim_table[1],
+ mode & MULTI_TRIMER_PAIRS);
+/************/
+ for (class=0; class<number_classes; class++) {
+ fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
+ print_matrix(both_means_submatrix[i][class],1,
+ new_num_dim_table[0] + new_num_dim_table[1],stderr);
+ fprintf(stderr,"\n\nInver Covar Submatrix %d:\n", class);
+ print_matrix(both_inverse_covars[i][class],
+ new_num_dim_table[0] + new_num_dim_table[1],
+ new_num_dim_table[0] + new_num_dim_table[1],stderr);
+
+ }
+/**********/
+ fprintf(stderr,"\nDet0 =%lf, det1= %lf, det2= %lf\n",
+ both_determinants[i][0],both_determinants[i][1],
+ both_determinants[i][2]);
+
+ 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);
+
+
+
+ }
+
+ }
+ switch_gauss_param(offset_to_use, main_method);
+
+ fprintf(stderr, "\n Out of get_gauss_param");
+
+
+
+#ifdef TEST_VERSION
+/*************************************************************************/
+/******************Now do gauss parameters for PairCoilDiff method. *****/
+/*************************************************************************/
+/**** Note the "2" parameter represents that are 2 classes (dimers, trimers)***/
+/**** not three, since just differentiating so don't need pdb-. ************/
+
+ if (differ_gauss_param[0] != ',') {
+ char differ_lib_hack[2][MAXFUNCTNUM];
+ /** If are computing max and min window scores, **/
+ /** then need to extract both dimenion i and **/
+ /** i + old_num_dim_differ/2. Treat them like they **/
+ /** were from another table. ***/
+
+
+ get_gauss_param_all_classes2(differ_class_sc_filenames,
+ old_num_dim_differ,
+ differ_means, differ_covars,
+ 2, differ_gauss_param);
+
+ if (which_differentiator ==2) {
+ for (i=0; i<number_differ_score_dim/2; i++) {
+ differ_lib_hack[0][i] = differ_lib[i];
+ differ_lib_hack[1][i] = differ_lib[i];
+ }
+ compute_submatrices2(differ_means,differ_means_submatrix,
+ differ_covars,2,
+ differ_inverse_covars,
+ differ_determinants, differ_lib_hack,
+ old_num_dim_differ/2, old_num_dim_differ/2,
+ number_differ_score_dim/2, number_differ_score_dim/2,
+ 1);
+ }
+ else
+ compute_submatrices2(differ_means,differ_means_submatrix,
+ differ_covars,2,
+ differ_inverse_covars,
+ differ_determinants, differ_lib,
+ old_num_dim_differ, 0,
+ number_differ_score_dim, 0,
+ 1);
+
+/************/
+ for (class=0; class<2; class++) {
+ fprintf(stderr,"\n\nDiffer Mean Submatrix %d:\n",class);
+ print_matrix(differ_means_submatrix[class],1,
+ number_differ_score_dim,stderr);
+ fprintf(stderr,"\n\nDiffer Inver Covar Submatrix %d:\n", class);
+ print_matrix(differ_inverse_covars[class],
+ number_differ_score_dim,number_differ_score_dim,
+ stderr);
+
+ }
+/**********/
+
+ fprintf(stderr,"\nDet0 =%lf, det1= %lf\n",
+ differ_determinants[0],differ_determinants[1]);
+
+ fprintf(stderr,"\nnew_num_dim_differ=%d, old_num_dim_differ=%d,", number_differ_score_dim, old_num_dim_differ);
+
+
+
+ }
+#endif
+
+
+ fprintf(stderr, "\n Out of get_gauss_param");
+
+
+
+ log_offset_to_use = offset_to_use;
+ method = main_method;
+ preprocessor_method = main_preprocessor_method;
+ table = main_table;
+
+
+ if ( (mode & PRN_MODE) && (flog || flog_coil_conflicts)){
+ log_file_header(flog, flog_coil_conflicts, mode,
+ argc, argv, avg_max, log_bound,
+ preprocessor_method, preproc_table, log_offset_to_use,
+ start_coil_at_reg_shift, weighted_avg, main_method,
+ main_table);
+ }
+
+ if (mode & POS_MODE) ActualCoil_Method = available;
+
+
+/****** Now if didn't get locations from .paircoil, get from ENV variable */
+ if (!fgin) { /* Didn't get fgin in */
+ fprintf(stderr,"\n\n **********************");
+ fprintf(stderr,"\nDid not get genbnk location from config file.");
+ fprintf(stderr,"\n Will instead use uniform distrib. as underlying distrib.");
+ fprintf(stderr,"\n**********************\n\n");
+
+/**** exit(-1); ****/ /* get_defaults. */
+ }
+ if (!fpin[0]) { /* Didn't get fpin in */
+ fprintf(stderr,"\nDid not get positive data set from config file.");
+ exit(-1); /* get_defaults. */
+ }
+/******************************************/
+
+
+ PairCoilData(fgin,fpin,&number_tables, mode & PROLINE_FREE_MODE,
+ table_to_remove_from, argv[1], mode,
+ prior_freq_single, prior_freq_pair, which_priors,
+ which_priorp, good_turing_fixed_disc, structural_pos);
+ /* If fpin[1] != NULL this will call */
+ /* TrimerDimerDiff to set up table diff */
+ /* statistics. */
+ /* If it is NULL then table is average */
+ /* of probabilities from 1 and 2. */
+
+
+/***/
+ hmm = CCHMM(fgin,fpin);
+/***/
+
+
+ if (mode & USE_LIKE_LINE) {
+ /** Don't want to compute likelihood line if */
+ /** not using. */
+ functnum = pair_functnum;
+ get_likelihood_line(likelihoods, m, b, m_single, b_single,
+ lib, pair_functnum,pir_name,mode,number_tables);
+ }
+ for (i=0; fpin[i]; i++)
+ fprintf(stderr, "\n m[%d]= %lf, b[%d]= %lf\n",i,m[i],i,b[i]);
+
+
+
+ i=0;
+ while (fpin[i]) {
+ sclose(fpin[i]);
+ i++;
+ }
+ if (fgin) sclose(fgin);
+ NEWCOILSinit();
+
+ sequence.seqlen = 0;
+
+
+ fin = sopen(argc>1?argv[1]:"-","r");
+ fprintf(stderr, "\nfin is opened to %s",argc>1?argv[1]:"-");
+
+ /*** Note: If can't open xdisplay then OpenScreen will write log and */
+ /*** out files if those options are set. */
+ fprintf(stderr,"\nInitializing\n");
+ OpenScreen(mode,log_bound,fin,flog,flog_coil_conflicts, fout_coils,
+ fout, m,b,m_single, b_single,
+ lib,differ_lib,multi_lib,
+ main_method,main_preprocessor_method, &seqnum, which_priors,
+ which_priorp, prior_freq_single, prior_freq_pair,
+ good_turing_fixed_disc, structural_pos);
+ painter();
+ /* Let user read message while initializing */
+
+
+
+ FD_ZERO(&fds);
+ FD_ZERO(&fdnone);
+ interval.it_value.tv_sec = 2;
+ interval.it_value.tv_usec = 0;
+ interval.it_interval.tv_sec = 2;
+ interval.it_interval.tv_usec = 0;
+
+ nopainter();
+ XSelectInput(display, window, KeyPressMask | ButtonPressMask | ExposureMask);
+ while (1) {
+ if (autoadvance) do {
+ FD_SET(ConnectionNumber(display),&fds);
+ if (!select(FD_SETSIZE,&fds,&fdnone,&fdnone,&interval))
+ NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,
+ mode,bound,flog,flog_coil_conflicts,fout_coils,
+ fout, which_priors, which_priorp, prior_freq_single,
+ prior_freq_pair, good_turing_fixed_disc, structural_pos);
+ }
+ while ( (!XPending(display)) && (autoadvance));
+ /* Autoadvance until button is hit or end of seq file. */
+ XNextEvent(display,&report);
+ switch (report.type) {
+ case Expose:
+ if (report.xany.window == zoom)
+ Zoomer();
+ else if (report.xany.window == help)
+ ShowHelp(help);
+ else if (sequence.seqlen)
+ ShowScores(scores,&sequence,NULL,bound, mode);
+ else
+ ShowTitle(window);
+ break;
+ case ButtonPress:
+/* if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single, b_single,
+ bound,flog,flog_coil_conflicts,
+ fout_coils,fout, which_priors,
+ which_priorp, prior_freq_single,
+ prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+ break;}
+*/
+ if (report.xany.window == help) {XUnmapWindow(display,help); break;}
+ if (report.xany.window == zoom) {XUnmapWindow (display, zoom); break;}
+ if (report.xany.window == window) {
+ /* Zoom in and show register */
+ int tmp;
+ tmp = XYtoPos(report.xbutton.x,report.xbutton.y);
+ if (tmp) {
+ ZoomPnt(); zoomptr=tmp; Zoomer(); ZoomPnt();
+ }
+ }
+ if (report.xany.window == helpb)
+ XMapWindow(display,help);
+ if (report.xany.window == next)
+ NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single, b_single,
+ mode,bound,flog,flog_coil_conflicts,fout_coils,fout,
+ which_priors,which_priorp, prior_freq_single, prior_freq_pair,
+ good_turing_fixed_disc, structural_pos);
+ if (report.xany.window == prev)
+ PrevSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,
+ mode,bound, which_priors,which_priorp, prior_freq_single,
+ prior_freq_pair, good_turing_fixed_disc, structural_pos);
+ if (report.xany.window == autob)
+ /* toggle the automatic sequencer */
+ { if (autoadvance ^= 1) NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,
+ mode,bound,flog,
+ flog_coil_conflicts,fout_coils,
+ fout, which_priors,which_priorp,
+ prior_freq_single, prior_freq_pair,
+ good_turing_fixed_disc,
+ structural_pos); }
+ else
+ autoadvance = 0;
+ if (report.xany.window == zoomb) {
+ if (zoomptr) ZoomPnt();
+ zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
+ }
+ if (report.xany.window == printb) {
+ XBell(display,-10);
+ ps = sopen(print?print:defprint,"a");
+ PSDefs(ps, sequence.seqlen);
+ ShowScores(scores,&sequence,ps,bound,mode);
+ sclose(ps);
+ }
+/***
+ if (report.xany.window == tableb) {
+ table = (table +1) % number_tables;
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,0,0, bound);
+ if (zoomptr) Zoomer();
+ }
+*****/
+ if (report.xany.window == methodb) {
+ method_cycle(&method,NUMBER_METHODS);
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,0, bound);
+ if (zoomptr) Zoomer();
+ }
+ if (report.xany.window == quit) {
+ if (flog || flog_coil_conflicts || fout || ftotal_like || fout_coils)
+ {
+ seqnum=sequp;
+ finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
+ fout,m,b,m_single,b_single,
+ lib,differ_lib,multi_lib,main_method,main_table,
+ main_preprocessor_method,
+ &seqnum, avg_max, which_priors,which_priorp,
+ prior_freq_single, prior_freq_pair,
+ good_turing_fixed_disc, structural_pos);
+ }
+
+ sclose(fin);
+ if (flog) sclose(flog); if (fout) sclose(fout);
+ if (fout_coils) sclose(fout_coils);
+ if (flog_coil_conflicts) sclose(flog_coil_conflicts);
+ if (ftotal_like) sclose(ftotal_like);
+ exit(0); }
+
+ break;
+ case KeyPress:
+ autoadvance = 0;
+/* if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single, b_single,
+ mode,bound,flog,flog_coil_conflicts,
+ fout_coils,fout, which_priors,
+ which_priorp, prior_freq_single,
+ prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+ break;}
+*/
+ if (report.xany.window == help) {XUnmapWindow(display,help); break;}
+ switch (XLookupKeysym(&(report.xkey),0)) {
+ case XK_Up:
+ bound = bound + .01;
+ if (bound > 1) bound =1;
+ draw_bound_box(bound, paramb);
+ break;
+ case XK_Down:
+ bound = bound - .01;
+ if (bound < 0) bound =0;
+ draw_bound_box(bound, paramb);
+ break;
+
+ case XK_Return: /** Signal that bound has been set, should rescore **/
+ /** -1 then means just need to recompute MultiCoil coil and seq scores. */
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,-1,-1, bound);
+ if (zoomptr) Zoomer();
+ break;
+
+
+ case XK_C: /** Switch between register, avg, and max **/
+ case XK_c: /** coil score for MultiCoil. **/
+ Coil_Score = (Coil_Score +1) % 3;
+ if (Coil_Score == 0)
+ ShowSeq(lib,differ_lib,multi_lib, /** Change to res score **/
+ m,b,m_single,b_single,mode,1,1, -2,-2, bound);
+ else ShowSeq(lib,differ_lib,multi_lib, /** Change to coil score with */
+ m,b,m_single,b_single,mode,-1,-1, bound); /* out recomputing */
+ if (zoomptr) Zoomer(); /* res score. */
+ break;
+
+ case XK_O: /** Oligomerization ratio. ***/
+ case XK_o:
+ mode ^= ONLY_COILED_CLASSES; /* To XOR to complement that bit. */
+ if (number_classes >0) {
+ if (method == MultiCoil)
+ type_class_score_convert(sequence,all_scores,bound,
+ mode & ONLY_COILED_CLASSES,number_tables,
+ number_classes);
+ if (preprocessor_method == MultiCoil)
+ type_class_score_convert(sequence,all_preprocess_like,bound,
+ mode & ONLY_COILED_CLASSES,number_tables,
+ number_classes);
+ }
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,0,2, bound);
+ if (zoomptr) Zoomer();
+ break;
+
+ case XK_Print:
+ XBell(display,-10);
+ ps = sopen(print?print:defprint,"w");
+ PSDefs(ps, sequence.seqlen);
+ ShowScores(scores,&sequence,ps,bound,mode);
+ sclose(ps);
+ XBell(display,-10);
+ break;
+ case XK_Left:
+ if (zoomptr>1)
+ {ZoomPnt(); --zoomptr; Zoomer(); ZoomPnt();}
+ break;
+ case XK_Right:
+ if (zoomptr && zoomptr<sequence.seqlen)
+ {ZoomPnt(); ++zoomptr; Zoomer(); ZoomPnt();}
+ break;
+ case XK_Z:
+ case XK_z:
+ if (zoomptr) ZoomPnt();
+ zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
+ break;
+
+#ifdef TEST_VERSION
+ case XK_A: /** Switch between unweighted avg and **/
+ case XK_a: /** weighted avg coil score. **/
+ weighted_avg = !weighted_avg;
+ if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,-1,-1, bound);
+ if (zoomptr) Zoomer();
+ }
+ break;
+
+ case XK_b:
+ case XK_B: /** For CoilScore, change how deal with reg shift. **/
+ start_coil_at_reg_shift = !start_coil_at_reg_shift;
+ if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,-1,-1, bound);
+ if (zoomptr) Zoomer();
+ }
+ break;
+
+
+ case XK_R:
+ case XK_r:
+ offset_to_use = (offset_to_use +1) % 8; /* Cycle through offsets */
+/*** switch_gauss_param(offset_to_use); **/
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,1, bound);
+ if (zoomptr) Zoomer();
+ break;
+ case XK_x:
+ offset_to_use =-1; /* Do max scoring register */
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,1, bound);
+ if (zoomptr) Zoomer();
+ break;
+ case XK_7:
+ offset_to_use =7; /* Do MultiCoil max score as combo of max score */
+ ShowSeq(lib,differ_lib,multi_lib, /* in each dimension. */
+ m,b,m_single,b_single,mode,1,1, bound);
+ if (zoomptr) Zoomer();
+ break;
+
+ case XK_s: /* Use single probabilities. */
+ case XK_S:
+ if (mode & PAIRCOIL_PAIRS) {
+ mode -= PAIRCOIL_PAIRS; /** Turn it off **/
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,1, bound);
+ if (zoomptr) Zoomer();
+ }
+ break;
+
+ case XK_l: /* Use pair probabilities at lib distances */
+ case XK_L:
+ if (!(mode & PAIRCOIL_PAIRS)) {
+ mode |= PAIRCOIL_PAIRS; /* Turn it on. */
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,1, bound);
+ if (zoomptr) Zoomer();
+ }
+ break;
+#endif
+
+ case XK_h:
+ case XK_H:
+ case XK_question:
+ case XK_Help:
+ XMapWindow(display,help);
+ break;
+ case XK_Q:
+ case XK_q:
+ case XK_Cancel:
+ if (report.xany.window == window)
+ {sclose(fin); if (flog) sclose(flog);
+ if (flog_coil_conflicts) sclose(flog_coil_conflicts);
+ if (ftotal_like) sclose(ftotal_like);
+ if (fout_coils)
+ sclose(fout_coils);
+
+ if (fout) sclose(fout);
+ exit(0);}
+ if (report.xany.window == zoom) {
+ ZoomPnt();
+ XUnmapWindow(display, zoom);
+ zoomptr=0;
+ }
+ break;
+ case XK_M:
+ case XK_m:
+ method_cycle(&method,NUMBER_METHODS);
+
+ if (method == MultiCoil)
+ table = table % number_classes;
+ else
+ table = table % number_tables;
+
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,1,0, bound);
+ if (zoomptr) Zoomer();
+ break;
+
+ case XK_W:
+ case XK_w:
+ method_cycle(&preprocessor_method,NUMBER_PREPROCESSORS);
+
+ if (preprocessor_method == MultiCoil)
+ preproc_table = preproc_table% number_classes;
+ else
+ preproc_table = preproc_table%number_tables;
+
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,0,1, bound);
+ if (zoomptr) Zoomer();
+ break;
+
+ case XK_i:
+ if (preprocessor_method == MultiCoil)
+ preproc_table = (preproc_table +1) % number_classes;
+ else
+ preproc_table = (preproc_table +1) %number_tables;
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,0,2, bound);
+ if (zoomptr) Zoomer();
+ break;
+ case XK_T:
+ case XK_t:
+ if (method == MultiCoil)
+ table = (table +1) % number_classes;
+ else
+ table = (table +1) % number_tables;
+ ShowSeq(lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,mode,0,0, bound);
+ if (zoomptr) Zoomer();
+ break;
+ case XK_P:
+ case XK_p:
+ case XK_Prior:
+ if (report.xany.window == window) PrevSeq(fin, lib,differ_lib,
+ multi_lib,m,b,m_single,b_single,mode,bound, which_priors,
+ which_priorp, prior_freq_single,prior_freq_pair,
+ good_turing_fixed_disc, structural_pos);
+ break;
+
+ case XK_N:
+ case XK_n:
+ case XK_Next:
+ case XK_space:
+ case XK_Linefeed:
+ if (report.xany.window == window)
+ NextSeq(fin,lib,differ_lib,multi_lib,
+ m,b,m_single,b_single,
+ mode,bound,flog,flog_coil_conflicts,fout_coils,
+ fout, which_priors, which_priorp,
+ prior_freq_single,prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+ break;
+ }
+ break;
+ }
+ }
+}
+
+static void Alarm ()
+{
+ XEvent report;
+ if (XCheckTypedEvent(display,Expose,&report)) {
+ ShowTitle(window);
+ XDrawString(display,window,gcgray,500,17,"(initializing)",14);
+ }
+}
+#include <signal.h>
+painter ()
+{
+ struct itimerval value;
+ signal(SIGALRM,Alarm);
+ getitimer(ITIMER_REAL,&value);
+ value.it_value.tv_sec = 0;
+ value.it_value.tv_usec = 500;
+ value.it_interval.tv_sec = 0;
+ value.it_interval.tv_usec = 500;
+ setitimer(ITIMER_REAL,&value,NULL);
+ Alarm();
+}
+nopainter ()
+{
+ struct itimerval value;
+ getitimer(ITIMER_REAL,&value);
+ value.it_value.tv_sec = 0;
+ value.it_value.tv_usec = 0;
+ setitimer(ITIMER_REAL,&value,NULL);
+ XDrawImageString(display,window,gcgray,500,17," ",14);
+ XDrawImageString(display,window,gcgray,225,729,"press any key to continue",25);
+}
+ZoomPnt ()
+{
+ int x,y;
+ XPoint pts[3];
+ if (!zoomptr) return;
+ PostoXY(zoomptr,&x,&y);
+ pts[0].x = x; pts[0].y = y;
+ pts[1].x = -5; pts[1].y = 5;
+ pts[2].x = 10; pts[2].y = 0;
+ XFillPolygon(display,window,gcdash,pts,3,Convex,CoordModePrevious);
+ XDrawLine(display,window,gcdash,x,y+5,x,y+8);
+}
+
+
+/* Find most likely residue, break ties arbitrarily */
+GoodZoom ()
+{
+ int i, zoom = -1;
+ double like=-1;
+ for (i=0; i<sequence.seqlen; ++i)
+ if (scores[8*i+7] > like) {
+ like = scores[8*i+7];
+ zoom = i;
+ }
+ return zoom+1;
+}
+
+
+Zoomer ()
+{ int i, offset;
+ double *sc;
+ char buf[32];
+ XMapWindow(display, zoom);
+ /* Indicate where in the sequence we are */
+ sprintf(buf,"Position %d ",zoomptr);
+ XDrawImageString(display,zoom,gc,7,20,
+ buf,strlen(buf));
+ for (i = -3; i<=3; ++i) {
+ if (zoomptr+i>0 && zoomptr+i<=sequence.seqlen)
+ buf[0] = numaa(sequence.seq[zoomptr-1+i]);
+ else
+ buf[0] = ' ';
+ XDrawImageString(display,zoom,gc,159+9*i,20, buf,1);
+ }
+ XDrawLine(display,zoom,gc,159,22,168,22);
+
+ /* Indicate the probability of being in any given state */
+ sc = &(scores[8*(zoomptr-1)]);
+
+ for (i=0; i<=7; ++i)
+ if (sc[i] == sc[7]) {offset=i; break;} /* Used to print out maxscoring */
+ /* register as 'a' + offset. */
+ /* Offset 7 will be a flag that no particular offset scored max. */
+
+ /* Give likelihood */
+ sprintf(buf,"%5.1f%% chance",100*sc[7]);
+ XDrawImageString(display,zoom,gc,220,20, buf,strlen(buf));
+
+/*****
+ fprintf(stderr,"\noffset= %d, reg = %c",all_offsets[table][zoomptr-1],
+ 'a' + (zoomptr-1 + all_offsets[table][zoomptr-1])%POSNUM);
+*******/
+
+ if ( sc[7] != 0)
+ if (offset != 7)
+ sprintf(buf,"register %c",'a'+offset);
+ else
+ sprintf(buf,"combo reg ");
+ else sprintf(buf," ");
+ XDrawImageString(display,zoom,gc,230,40, buf,strlen(buf));
+
+ for (i=0; i<7; ++i) {
+ sprintf(buf,"%c:%4.1f%% ",'a'+i,(100*sc[i]));
+ XDrawImageString(display,zoom,gc,15+92*(i&3),60+20*(i/4),buf,strlen(buf));
+ }
+ XDrawString(display,zoom,gc,10,110,"Click on viewgram or use arrow",30);
+ XDrawString(display,zoom,gc,10,127,"keys to view other positions.",29);
+ XDrawString(display,zoom,gc,10,144,"You may need to move this",25);
+ XDrawString(display,zoom,gc,10,161,"window to see the viewgram.",27);
+}
+
+
+
+
+
+
+
+
+/*********The following is for storing sequences. ***/
+static char seqbuff[40000];
+
+static Sequence seqs[30];
+static int seqbufpos[30];
+static int seqsizes[30];
+
+ReadSeq()
+{ int size, pos, limsize;
+/*
+ size = seqsizes[sequp-seqnum];
+ pos = seqbufpos[sequp-seqnum];
+*/
+ sequence = seqs[sequp-seqnum]; /** Retrieve seq from memory. **/
+
+/*
+ limsize = (pos+size<=40000) ? size : 40000-pos;
+ bcopy(&(seqbuff[pos]),sequence.code,limsize);
+ bcopy(seqbuff,&(sequence.code[limsize]),size-limsize);
+*/
+}
+
+
+StoreSeq()
+{ int size, pos, limsize,i;
+ char *temp, *current;
+
+ bcopy(seqs,&(seqs[1]),29*sizeof(Sequence));
+ bcopy(seqbufpos,&(seqbufpos[1]),29*sizeof(int));
+ bcopy(seqsizes,&(seqsizes[1]),29*sizeof(int));
+ if (sequp-seqlow>=30) seqlow = sequp-29;
+
+/* The 2 accounts for the end of string characters in .code and .title. */
+ size = 2 + strlen(sequence.code) + strlen(sequence.title) +
+ sequence.seqlen * (sequence.reg ? 2 : 1);
+
+ pos = seqbufpos[1] + seqsizes[1]; /* Number of characters in array */
+
+/* if (pos>=40000) pos -= 40000;*/
+/* limsize = (pos+size<=40000) ? size : 40000-pos; */
+
+ if (pos+size >=40000) pos=0; /* Array will be full so start at beginning. */
+ limsize = (pos+size<=40000) ? size : 40000-pos;
+
+
+/* Save the "code+title+seq+reg" in something permanant, and wrap around
+at end of array. */
+
+/* David uses bcopy because wants to copy past end of string characters */
+/*
+ bcopy(sequence.code,&(seqbuff[pos]),limsize);
+ bcopy(&(sequence.code[limsize]),seqbuff,size-limsize);
+*/
+
+ sequence.code=strcpy(&(seqbuff[pos]),sequence.code);
+ sequence.title=
+ strcpy(&(seqbuff[pos+strlen(sequence.code)+1]),sequence.title);
+
+ temp=sequence.seq;
+ sequence.seq=
+ &(seqbuff[pos+strlen(sequence.code) + strlen(sequence.title) + 2]);
+ for(i=0; i<sequence.seqlen; i++)
+ sequence.seq[i] = temp[i];
+
+
+ if (sequence.reg) {
+ temp=sequence.reg;
+ sequence.reg=&(sequence.seq[sequence.seqlen]);
+ for(i=0; i<sequence.seqlen; i++)
+ sequence.reg[i] = temp[i];
+ }
+
+
+ while ((seqlow<sequp) && (seqbufpos[sequp-seqlow] >= pos && seqbufpos[sequp-seqlow] < pos+limsize || seqbufpos[sequp-seqlow] < size-limsize))
+ ++seqlow;
+ seqs[0] = sequence;
+ seqs[0].seqlen= sequence.seqlen;
+ seqsizes[0] = size;
+ seqbufpos[0] = pos;
+}
+
+
+
+
+
+
+method_cycle(int *Method, int Number_of_Methods)
+{
+ do {
+ *Method = (*Method+1)%Number_of_Methods;
+ }
+ while ((*Method == MultiCoil && MultiCoil_Method == unavailable) ||
+ (*Method == PairCoil && PairCoil_Method == unavailable) ||
+ (*Method == NEWCOILS && NEWCOILS_Method == unavailable) ||
+ (*Method == HMMCoil && HMMCoil_Method == unavailable) ||
+/**** (*Method == TableDiffRes &&
+ TableDiffRes_Method == unavailable) ||
+ (*Method == TableDiffCoil &&
+ TableDiffCoil_Method == unavailable) ||
+****/
+ (*Method == ActualCoil && ActualCoil_Method == unavailable) ||
+ (*Method == PairCoilDiff && PairCoilDiff_Method == unavailable)||
+ (*Method == PairCoilDiffAvg && PairCoilDiffAvg_Method == unavailable) ||
+
+ (*Method == ActualCoil) && (!sequence.reg) /* Don't do actual coil */
+ ); /* if didn't input register */
+
+}
+
+
+
+PrevSeq (FILE *fin, char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m, double *b,
+ double *m_single, double *b_single,
+ int mode, double bound, int which_priors, int which_priorp,
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
+ int structural_pos[POSNUM+1])
+{
+ if (seqnum<=seqlow) {XBell(display,0); return;}
+ --seqnum;
+ ReadSeq();
+ NewSeq(mode,sequence, which_priors, which_priorp, prior_freq_single,
+ prior_freq_pair, good_turing_fixed_disc, structural_pos);
+ ShowSeq(lib, differ_lib, multi_lib,
+ m, b, m_single, b_single, mode,1,1, bound);
+}
+
+
+NextSeq (FILE *fin,char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m, double *b, double *m_single, double *b_single,
+ int mode,
+ double bound, FILE *flog, FILE* flog_coil_conflicts,
+ FILE *fout_coils,
+ FILE *fout, int which_priors, int which_priorp,
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
+ int structural_pos[POSNUM+1])
+{
+ char title[MAXLINE],code[MAXLINE];
+ char seq[MAXSEQLEN], reg[MAXSEQLEN];
+ Sequence temp_seq;
+ int first_time_seq=0;
+
+/********
+ fprintf(stderr,"\nIn next seq");
+*********/
+
+ if (seqnum < sequp) {
+ ++seqnum;
+ ReadSeq();
+ }
+ else {
+ temp_seq= sequence;
+ sequence.seqlen = 0;
+ sequence.seq = seq;
+ sequence.reg = reg;
+ sequence.title = title;
+ sequence.code = code;
+
+
+/* if ( ( (mode & POS_MODE) ? !getpos(fin,&sequence) :
+ !getseq2(fin,&sequence) ) ) */
+/*** getseq2 now gets register info too... ***/
+
+ if (!getseq2(fin,&sequence))
+ {XBell(display,0); sequence= temp_seq; autoadvance=0; return;}
+
+
+ ++seqnum; ++sequp;
+ StoreSeq();
+ first_time_seq =1;
+ }
+
+ NewSeq(mode,sequence, which_priors, which_priorp,
+ prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+ if (first_time_seq) {
+
+ output_seq(lib,differ_lib, multi_lib,
+ m,b,m_single, b_single,
+ mode,log_bound,flog,flog_coil_conflicts,fout_coils,fout,
+ avg_max, main_method, main_preprocessor_method, main_table);
+
+ }
+
+
+ ShowSeq(lib, differ_lib,multi_lib,
+ m, b, m_single, b_single,mode,1,1, bound);
+
+}
+
+
+/** This is stuff need to always do before rescore a sequence. **/
+NewSeq_nonX (int mode, Sequence sequence, int which_priors, int which_priorp,
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ double prior_freq_pair[MAX_TABLE_NUMBER],
+ int good_turing_fixed_disc, int structural_pos[POSNUM+1])
+{
+ if (mode & TST_MODE0)
+ recalc_prob(sequence, 0, mode, which_priors, which_priorp,
+ prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+ if (mode & TST_MODE1)
+ recalc_prob(sequence,1, mode, which_priors, which_priorp,
+ prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+}
+
+NewSeq (int mode, Sequence sequence, int which_priors, int which_priorp,
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
+ int structural_pos[POSNUM+1])
+{
+ NewSeq_nonX (mode, sequence, which_priors, which_priorp,
+ prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
+ structural_pos);
+
+ XUnmapWindow(display, zoom);
+ zoomptr=0;
+}
+
+
+/* The likelihood line is mx +b */
+void output_seq(char lib[MAXFUNCTNUM],char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m,double *b,
+ double *m_single, double *b_single,
+ int mode,
+ double log_bound,FILE *flog, FILE *flog_coil_conflicts,
+ FILE *fout_coils,
+ FILE *fout,
+ int avg_max,
+ int main_method, int main_preprocessor_method,
+ int main_table)
+{
+
+ double maxother=0; double maxscore;
+ extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
+
+ int i,j;
+ int do_preproc;
+
+ double junk_scores[MAXSEQLEN][POSNUM+1];
+
+/***
+ fprintf(stderr,"\nIn outputseq");
+****/
+
+ if (is_not_differentiator(main_method))
+ do_preproc = 0;
+ else do_preproc =1;
+
+ all_scores = ScoreSeq (lib, differ_lib, multi_lib,
+ m, b, m_single, b_single, mode,
+ main_table,
+ main_method, main_preprocessor_method,
+ log_offset_to_use, &maxscore,
+ 1, do_preproc, log_bound);
+
+
+ if (ftotal_like) {
+ add_to_total_likes(all_scores, total_gauss_like, sequence.seqlen,
+ ftotal_like);
+ }
+/******
+ fprintf(stderr,"\nGot all_scores");
+*******/
+ scores = &all_scores[main_table* MAXSEQLEN*(POSNUM+1) ];
+
+ /* This is all_scores[main_table + 0*number_tables][0][0] */
+ /* so is score on main_table using library # 0 **/
+ /* Similarly for preproces_like. **/
+
+ preprocess_like = &all_preprocess_like[main_table][0][0];
+
+
+ for (i=0; i<sequence.seqlen; i++)
+ if (preprocess_like[i*(POSNUM+1) + 7] > maxother)
+ maxother=preprocess_like[i*(POSNUM +1) + 7];
+ /* This is preprocess_like[i][7] */
+
+
+ switch (main_method) {
+ case STOCK_PAIR:
+ if (mode & VER_MODE) {
+ NEWCOILSScore(mode, sequence, &maxother,junk_scores,
+ log_offset_to_use);
+ /* just dump */
+ /* scores in junk_scores */
+ log_output_ver(SC2SEQ | SCSTOCK, maxscore, maxother,
+ maxscore,seqnum, sequence.title,
+ sequence.code, flog);
+ }
+ break;
+
+ case MultiCoil:
+/****/
+ if (fout_coils)
+ /*** For pos files coils, does a log file version of coil/seq score. **/
+ output_pos_scores(fout_coils, sequence, all_scores, log_bound, mode,
+ by_coil_or_seq, weighted_avg);
+
+
+ if (mode & PRN_MODE)
+
+ if (Coil_Score && flog) {
+/*** Were having trouble with this in printing out coils, since offset
+ *** would sometime change at ends of coils between residue offset and
+ *** coil score offset. To make log file print out nice, use res. offset.
+ char all_offs[MAX_NUM_SCORE_DIM][MAXSEQLEN];
+ get_offsets(all_scores, sequence,all_offs,3);
+***********/
+ log_coil_output(all_scores, all_offsets,sequence, log_bound, flog,
+ seqnum,mode);
+ }
+ else if (flog)
+ log_output2_prn(mode,maxscore, maxscore, scores, scores,
+ sequence.seqlen,
+ seqnum, sequence, log_bound, log_bound, log_bound,
+ flog, avg_max);
+ /****/
+ break;
+
+ case PairCoil:
+ case NEWCOILS:
+ /****
+ if (mode & VER_MODE)
+ log_output_ver(SC2SEQ , maxscore, 0,
+ maxscore,seqnum, sequence.title,
+ sequence.code, flog);
+***/
+
+ if (mode & PRN_MODE)
+ log_output2_prn(mode,maxscore, maxscore, scores, scores,sequence.seqlen,
+ seqnum, sequence, log_bound, log_bound, log_bound,
+ flog, avg_max);
+
+
+ if (flog_coil_conflicts)
+ coil_conflicts(mode,all_scores,all_scores,sequence, seqnum,log_bound,
+ log_bound,number_tables, flog_coil_conflicts,avg_max);
+ break;
+
+/* case PairCoilDiff: */
+ case PairCoilDiffAvg:
+ if (mode & PRN_MODE) {
+ log_output2_prn(mode,maxscore, maxother, scores,preprocess_like,
+ sequence.seqlen,
+ seqnum, sequence, .5, .5, log_bound,
+ flog, avg_max);
+ }
+
+ if (flog_coil_conflicts) {
+
+ coil_conflicts(mode,all_scores,all_preprocess_like,sequence, seqnum,
+ .5, log_bound, number_tables, flog_coil_conflicts,
+ avg_max);
+
+ }
+ break;
+
+ case ActualCoil:
+ break;
+ }
+
+
+/***************
+ fprintf(stderr,"\n About to go into fout.");
+***************/
+
+
+ if ( fout)
+ if ( (!(mode & RAW_OUT)) ||
+ ( (main_method != MultiCoil) && (main_method != PairCoilDiff))) {
+ if (main_method== MultiCoil)
+ tuple_output2(sequence, mode, fout, all_scores,3, 1,log_bound,3);
+ /* Output based on 3 tables, meaning */
+ /* output the dimer and trimer likes */
+ /* and total like. */
+ else txt_output2(sequence, mode, fout, scores, log_bound);
+ }
+ else { /* Want raw tuple_output for multi_coil and PairCoilDiff*/
+
+ if (main_method == MultiCoil)
+ if (log_offset_to_use==7) /* Combo offset just output max score */
+ tuple_output2(sequence, mode, fout, /** for each dimension. */
+ all_scores, number_tables,
+ number_multi_lib,log_bound, number_score_dim);
+ else /** Max reg for raw score can't really be determined. So instead */
+ /* for coils, output scores for the known correct register, and */
+ /* for pdb- output ALL possible 7 register scores (all negs). */
+ tuple_output_for_max(sequence,mode,fout, all_scores, number_tables,
+ number_multi_lib, log_bound, number_score_dim);
+
+ else if (main_method == PairCoilDiff)
+ /* Output max score for each dimension. */
+ /** Had to do a hack since tuple_out2 expects number_lib */
+ /** to be array of size num_tab, so pass in 1 for num_tab */
+ tuple_output2(sequence, mode, fout, all_scores, 1,
+ &number_differ_lib,log_bound, number_differ_score_dim);
+
+ }
+
+}
+
+
+int copy_offsets(char preproc_offsets[MAXSEQLEN],
+ Sequence sequence,
+ char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
+ int number_dimens,
+ double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+ int offset_to_use)
+{
+ int tablenum,i;
+
+ for (tablenum =0; tablenum < number_dimens; tablenum++)
+ for(i=0; i<sequence.seqlen; i++) {
+ all_offsets[tablenum][i]=preproc_offsets[i];
+ /*** Make sure the output score corresponds to this offset... **/
+ /*** Also do it for backup residue scores in tablenum +3 ***/
+ if ((offset_to_use ==7) || (offset_to_use == -1)) {
+ all_scores[tablenum][i][7] =
+ all_scores[tablenum][i][(i+all_offsets[tablenum][i]) %POSNUM];
+ all_scores[tablenum+3][i][7] =
+ all_scores[tablenum+3][i][(i+all_offsets[tablenum][i]) %POSNUM]; }
+ }
+}
+
+int get_offsets(double all_likes[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+ Sequence sequence, char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
+ int number_dimens)
+{
+ int tablenum;
+ int got_offset, max_offset;
+ double max_off_score;
+ int i,j;
+
+ for (tablenum =0; tablenum < number_dimens; tablenum++)
+ for(i=0; i<sequence.seqlen; i++) {
+ offsets[tablenum][i]=-1;
+ got_offset=0;
+ max_off_score=-HUGE_VAL;
+ for (j=0; j<POSNUM; j++) {
+ if (all_likes[tablenum][i][j] == all_likes[tablenum][i][7]) {
+ offsets[tablenum][i] = (j-i) % POSNUM;
+ if (offsets[tablenum][i]<0)
+ offsets[tablenum][i] += POSNUM;
+ got_offset=1;
+ }
+ if (all_likes[tablenum][i][j] > max_off_score) {
+ max_off_score=all_likes[tablenum][i][j];
+ max_offset=(j-i) % POSNUM;
+ if (max_offset <0) max_offset +=POSNUM;
+ }
+ }
+ if (!got_offset) offsets[tablenum][i] = max_offset;
+
+ /* This is a new thing so that if an offset change is not real */
+ /* we keep the old offset. **********/
+ if ( (i>0) && (all_likes[tablenum][i][(i+offsets[tablenum][i-1]) %POSNUM]
+ == all_likes[tablenum][i][(i+offsets[tablenum][i])%POSNUM]) )
+ offsets[tablenum][i] = offsets[tablenum][i-1];
+
+ }
+}
+
+
+
+
+
+
+
+
+void preprocessor_score(char lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m, double *b,
+ double *m_single, double *b_single, int mode,
+ int Preprocessor_Method, int Offset_to_Use,
+ double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+ Sequence sequence,
+ char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
+ double bound)
+{
+ double maxscore;
+ char offsets[MAXSEQLEN];
+ int tablenum;
+ char *local_lib;
+ int dist;
+ int use_pairs=1;
+ int dimension, num_multilib;
+
+
+
+ for (tablenum =0; tablenum < number_tables; tablenum++) {
+ switch_tables(tablenum); /* Set prob. to be for that table. */
+ if (tablenum ==0) use_pairs =1; /* Use pairs for cctb */
+ else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
+ /* Do what config file said to do for trimers. **/
+
+ switch(Preprocessor_Method) {
+ case MultiCoil:
+ functnum = multi_functnum[tablenum]; /* For use in scscore.c **/
+ num_multilib= compute_num_multilib(tablenum, number_multi_lib,
+ mode & MULTI_TRIMER_PAIRS);
+ for (dist=0; dist<num_multilib; dist++) {
+ if (number_multi_lib[tablenum]>1) {
+ local_lib = &multi_lib[tablenum][dist];
+ }
+ else local_lib = multi_lib[tablenum];
+ dimension = compute_dimension(tablenum, number_tables, dist,
+ number_score_dim, number_multi_lib,
+ mode & MULTI_TRIMER_PAIRS);
+
+ MultiCoilDimensionScore(mode, sequence, local_lib,
+ &maxscore,tablenum, offsets,
+ all_preprocess_like[dimension],
+ Offset_to_Use,use_pairs);
+ }
+ break;
+
+ case PairCoil:
+ if (mode & PAIRCOIL_PAIRS) {
+ functnum = pair_functnum; /* For use in scscore.c **/
+ PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
+ &maxscore,tablenum, offsets,
+ all_preprocess_like[tablenum],
+ Offset_to_Use, 0);
+ }
+ else /* Use singles probabilities... for likelihood use line 0 */
+ SingleCoilScore(mode,sequence,m_single[tablenum],
+ b_single[tablenum],
+ &maxscore, tablenum,
+ all_preprocess_like[tablenum], Offset_to_Use,
+ mode & RAW_OUT);
+ break;
+
+ case NEWCOILS:
+ NEWCOILSScore(mode, sequence, &maxscore, all_preprocess_like[tablenum], Offset_to_Use);
+ break;
+
+ case HMMCoil:
+ HMMScore(sequence.seq, sequence.seqlen, hmm, all_preprocess_like[tablenum],
+ tablenum, Offset_to_Use);
+ break;
+
+ case ActualCoil:
+ /* if (mode & POS_MODE) */ /* Score the coils in the posfile based on */
+ /* their having 100% likelihood. */
+ ActualCoils(sequence, all_preprocess_like[tablenum], Offset_to_Use,1);
+ break;
+
+ default:
+ fprintf(stderr,"\nBad preprocessor, using PairCoil.\n");
+ functnum = pair_functnum; /* For use in scscore.c **/
+ PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
+ &maxscore,tablenum, offsets,
+ all_preprocess_like[tablenum],
+ Offset_to_Use, 0);
+ break;
+ }
+ }
+
+ if (Preprocessor_Method !=MultiCoil)
+ get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
+ number_tables);
+ else {
+/***** Get copy of res scores in tablenum and tablenum+3 **/
+ switch_gauss_param(Offset_to_Use, Preprocessor_Method);
+
+ convert_raw_scores_to_gauss_prob_like2(all_preprocess_like,sequence.seqlen,
+ number_tables,means_submatrix,inverse_covars,
+ determinants, number_classes, init_class_prob,
+ number_score_dim,
+ mode & ONLY_COILED_CLASSES,
+ bound,Offset_to_Use, 1);
+
+ get_offsets(all_preprocess_like, sequence,all_preproc_offsets,3);
+
+ if (Coil_Score) { /* Value of 1 means do average, 2 means do max. */
+ /*** Need to save residue scores, since the coil scores **/
+ /*** will replace them, but old residue scores might be needed **/
+ /*** again in tablenum+3. ****/
+ for (tablenum =0; tablenum <3; tablenum++) {
+ average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
+ all_preprocess_like[tablenum],
+ all_preprocess_like[2+3],/** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_preproc_offsets[tablenum],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ }
+ }
+ }
+
+
+
+
+}
+
+
+
+
+double *ScoreSeq (char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m_paircoil, double *b_paircoil,
+ double *m_single, double *b_single,
+ int mode, int Table, int Method,
+ int Preprocessor_Method, int Offset_to_Use,
+ double *maxscore,
+ int rescore_seq, int rescore_preproc, double bound)
+{
+ /** rescore_preproc == 1 means need to rescore. If it or rescore_seq are */
+ /** -1 then means just need to recompute MultiCoil coil and seq scores. */
+ /** -2 means need to change from MultiCoil coil score back to res score. */
+ /** If it is another non-zero */
+ /** value, don't need to rescore it, but do need to rescore anything that */
+ /** uses it as a preproc. */
+
+ double max_paircoil_score;
+ static char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
+ extern char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
+ extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
+ static double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
+ int i,j;
+ int tablenum;
+ double maxsc;
+ char *local_lib;
+ int dist;
+
+ int use_pairs=1;
+ int num_multilib, dimension;
+ int local_number_tables;
+
+ double *preproc_res_like;
+ int scored_already =0;
+
+/**************************************************************************/
+/********* PreprocLike is used in the 2-3 stranded differentiators. ****/
+
+
+ if ( (rescore_preproc==1) &&
+ (!ftotal_like)) { /* A hack to save computation time. */
+ preprocessor_score(lib, multi_lib,
+ m_paircoil, b_paircoil, m_single, b_single,
+ mode, Preprocessor_Method,
+ Offset_to_Use,
+ all_preprocess_like, sequence, all_preproc_offsets,
+ bound);
+ }
+
+/*****************************************************************************/
+ if (Preprocessor_Method == MultiCoil) {
+ if ( ( rescore_preproc==-1) && Coil_Score) {
+ /* Coil_Score value of 1 means do average, 2 means do max. */
+ /*** Need to save residue scores, since the coil scores **/
+ /*** will replace them, but old residue scores might be needed **/
+ /*** again in tablenum+3. ****/
+ for (tablenum =0; tablenum <3; tablenum++)
+ average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
+ all_preprocess_like[tablenum],
+ all_preprocess_like[2+3],/** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_preproc_offsets[tablenum],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ if (mode & ONLY_COILED_CLASSES) /** Convert to olig state. **/
+ type_class_score_convert(sequence,all_preprocess_like,bound,
+ mode & ONLY_COILED_CLASSES,number_tables,
+ number_classes);
+ }
+ else if (rescore_preproc == -2) /** Switch back to res score.*/
+ for (tablenum =0; tablenum <3 ; tablenum++)
+ for (j =0; j <=POSNUM; j++)
+ for (i=0; i<sequence.seqlen; i++) {
+ all_preprocess_like[tablenum][i][j] =
+ all_preprocess_like[tablenum+3][i][j];
+ if ((mode & ONLY_COILED_CLASSES) && (tablenum <2)&&
+ (all_preprocess_like[2+3][i][j]>bound))
+ all_preprocess_like[tablenum][i][j]/=
+ all_preprocess_like[2+3][i][j];
+ }
+
+ }
+/****************************************************************************/
+
+
+
+
+
+ if (Method == PairCoilDiff)
+ local_number_tables = 1;
+ else local_number_tables= number_tables;
+
+
+/*** Clean up this conditional, since really the only thing that might ***/
+/*** change for paircoildiffer when change preproc is the offset to use **/
+/*** (i.e. all_scores[table][i][7]. ***/
+
+ if ((rescore_seq>0) ||
+ (!is_not_differentiator(Method) && rescore_preproc) ) {
+ scored_already=1;
+ for (tablenum=0; tablenum<local_number_tables; tablenum++) {
+ switch_tables(tablenum); /* Switch prob to right table. */
+
+ if (tablenum ==0) use_pairs =1; /* Use pairs for cctb */
+ else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
+ /* Do what config file said to do for trimers. **/
+
+ maxsc=0; /* Initialization. */
+
+
+ switch (Method) {
+
+ case STOCK_PAIR:
+ functnum = pair_functnum; /* For use in scscore.c **/
+ PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
+ b_paircoil[tablenum],
+ &maxsc,tablenum, all_offsets[tablenum],all_scores[tablenum],
+ Offset_to_Use, 0);
+ break;
+
+
+ case MultiCoil:
+ functnum = multi_functnum[tablenum]; /* For use in scscore.c **/
+ num_multilib = compute_num_multilib(tablenum,number_multi_lib,
+ mode & MULTI_TRIMER_PAIRS);
+ for (dist=0; dist<num_multilib; dist++) {
+ if (number_multi_lib[tablenum] >1)
+ local_lib = &multi_lib[tablenum][dist];
+ else local_lib = multi_lib[tablenum];
+ dimension = compute_dimension(tablenum, number_tables, dist,
+ number_score_dim, number_multi_lib,
+ mode & MULTI_TRIMER_PAIRS);
+
+ MultiCoilDimensionScore(mode, sequence, local_lib, &maxsc,
+ tablenum,all_offsets[dimension],
+ all_scores[dimension],
+ Offset_to_Use,use_pairs);
+ }
+
+ break;
+
+ case PairCoil:
+ if (mode & PAIRCOIL_PAIRS) {
+ functnum = pair_functnum; /* For use in scscore.c **/
+ PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
+ b_paircoil[tablenum], &maxsc,
+ tablenum,all_offsets[tablenum],
+ all_scores[tablenum],
+ Offset_to_Use, (mode & RAW_OUT));
+ }
+ else /*Use singles probabilities... for likelihood use line number 0 */
+ SingleCoilScore(mode,sequence,m_single[tablenum],
+ b_single[tablenum],
+ &maxsc, tablenum,
+ all_scores[tablenum], Offset_to_Use,
+ (mode & RAW_OUT));
+ break;
+ case NEWCOILS:
+ NEWCOILSScore(mode,sequence,&maxsc, all_scores[tablenum],
+ Offset_to_Use);
+ break;
+
+ case HMMCoil:
+ HMMScore(sequence.seq, sequence.seqlen, hmm, all_scores[tablenum],
+ tablenum, Offset_to_Use);
+ break;
+
+ case PairCoilDiff:
+ functnum = differ_functnum; /* For use in scscore.c **/
+ for (dimension=0; dimension<number_differ_lib; dimension++) {
+ if (number_differ_lib >1) {
+ local_lib = &differ_lib[dimension];
+ }
+ else local_lib = differ_lib;
+
+
+ PairCoilDifferDimension(mode, sequence, local_lib, &maxsc,
+ all_preproc_offsets[preproc_table],
+ all_preprocess_like[preproc_table],
+ all_scores,
+ all_offsets, dimension,
+ Offset_to_Use, number_differ_score_dim,
+ bound);
+
+ }
+/***********************
+ PairCoilDiffer(mode, sequence, local_lib, &maxsc,
+ all_preproc_offsets[tablenum],
+ all_preprocess_like[preproc_table],
+ Offset_to_Use, tablenum,
+ all_scores[tablenum],0,
+ weighted_avg,start_coil_at_reg_shift);
+
+*************************/
+ break;
+
+/*************
+ case PairCoilDiffAvg:
+ functnum = differ_functnum;
+ PairCoilDiffer(mode, sequence, differ_lib, &maxsc,
+ all_preproc_offsets[preproc_table],
+ all_preprocess_like[preproc_table],
+ Offset_to_Use, tablenum,
+ all_scores[tablenum],
+ 1, weighted_avg,start_coil_at_reg_shift);
+ break;
+****************/
+
+ case ActualCoil:
+ ActualCoils(sequence, all_scores[tablenum], Offset_to_Use,0);
+ break;
+ }
+ if (tablenum == Table) *maxscore = maxsc;
+
+ }
+
+
+
+
+
+/*********************************************************************/
+/*********************************************************************/
+/********* This ends the scoring along all and all dimensions ********/
+/********* Now convert the raw scores that haven't been converted ****/
+/********* yet into likelihoods (i.e. the multidimensional scores.****/
+
+
+/***** Note: Can use same "init_class_prob" here, because since set number **/
+/*** of classes to 2, the 3rd class (Pdb-) will be ignored, so just get ****/
+/*** ratio of two to three stranded likelihoods, BOTH weighted by their ****/
+/*** init probaabiities. ********/
+
+ if ( (Method != MultiCoil) && (Method !=PairCoilDiff))
+ get_offsets(all_scores, sequence,all_offsets,number_tables);
+ else if (!(mode & RAW_OUT)) {
+
+ switch_gauss_param(Offset_to_Use, Method);
+
+
+/***********Do conversions on PairCoilDiff scores, using preprocessor. **/
+/*********** Note for differentiator there are only 2 possible classes **/
+/*********** (dimer and trimer). **********/
+ if (Method == PairCoilDiff) {
+ convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
+ number_tables, differ_means_submatrix,
+ differ_inverse_covars,differ_determinants, 2,
+ init_class_prob, number_differ_score_dim,
+ 0, bound, Offset_to_Use,0);
+ if (Preprocessor_Method == MultiCoil)
+ preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
+ else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
+
+
+/* Only use window scores that lie completely in coil */
+ if ( (which_differentiator == 3) || (which_differentiator ==4))
+ zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
+ sequence.seqlen, bound);
+
+ if (which_differentiator == 1)
+ zero_out_bad_windows(differ_window_length,
+ all_scores,all_preprocess_like[preproc_table],
+ sequence.seqlen,bound);
+
+ copy_offsets(all_preproc_offsets[preproc_table],sequence,
+ all_offsets,2, all_scores, Offset_to_Use);
+
+
+ if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
+ for (tablenum=0; tablenum<2; tablenum++) {
+ average_score_over_coils3(sequence, all_scores[tablenum+3],
+ all_scores[tablenum],
+ preproc_res_like,
+ /** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_preproc_offsets[preproc_table],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ }
+ }
+
+/**********Do conversions on MultiCoil dimension scores. ******/
+
+ else if (Method == MultiCoil) {
+ convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
+ number_tables, means_submatrix,inverse_covars,
+ determinants, number_classes,init_class_prob,
+ number_score_dim,
+ mode & ONLY_COILED_CLASSES, bound,
+ Offset_to_Use,1);
+ get_offsets(all_scores, sequence,all_offsets,3);
+
+
+
+ if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
+ for (tablenum=0; tablenum<3; tablenum++)
+ average_score_over_coils3(sequence, all_scores[tablenum+3],
+ all_scores[tablenum],
+ all_scores[2+3], /** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_offsets[tablenum],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ }
+ }
+
+
+ } /* Ends if rescore_seq */
+
+
+/*****************************************************************************/
+ if ( (Method == MultiCoil) || (Method == PairCoilDiff)){
+
+ if (Method == MultiCoil)
+ if ( (rescore_seq ==1 ) || /** Signals rescored MultiCoil **/
+ (rescore_seq==-1)) /*Signals bound changed, need to redo seq_score */
+ get_seq_scores(seq_scores, sequence, all_scores[0+3],
+ all_scores[1+3],bound, weighted_avg);
+
+ if (rescore_seq==-1 && !scored_already && (Method==PairCoilDiff)) {
+ if ( (which_differentiator == 3) || (which_differentiator ==4))
+ zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
+ sequence.seqlen, bound);
+ if (which_differentiator == 1) /* Only use window scores that lie */
+ zero_out_bad_windows(differ_window_length, /* completely in coil */
+ all_scores, all_preprocess_like[preproc_table],
+ sequence.seqlen,bound);
+ }
+
+ if ( (rescore_seq==-1) && Coil_Score && !scored_already) {
+ if (Method == MultiCoil) {
+ for (tablenum=0; tablenum<3; tablenum++) {
+ average_score_over_coils3(sequence, all_scores[tablenum+3],
+ all_scores[tablenum],
+ all_scores[2+3], /** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_offsets[tablenum],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ if (mode & ONLY_COILED_CLASSES) /** Convert to olig. state. **/
+ type_class_score_convert(sequence,all_scores,bound,
+ mode & ONLY_COILED_CLASSES, number_tables,
+ number_classes);
+ }
+ }
+ else if (Method == PairCoilDiff) {
+ if (Preprocessor_Method == MultiCoil)
+ preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
+ else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
+
+
+ for (tablenum=0; tablenum<2; tablenum++)
+ average_score_over_coils3(sequence, all_scores[tablenum+3],
+ all_scores[tablenum],
+ preproc_res_like,
+ /** Total Coil Like. **/
+ weighted_avg, offset_to_use,
+ all_offsets[tablenum],
+ start_coil_at_reg_shift, bound,
+ Coil_Score -1);
+ }
+ }
+
+
+ else if (rescore_seq == -2) /** Switch back to reg score.*/
+ for (tablenum =0; tablenum <3 ; tablenum++)
+ for (j =0; j <=POSNUM; j++)
+ for (i=0; i<sequence.seqlen; i++) {
+ all_scores[tablenum][i][j] =
+ all_scores[tablenum+3][i][j];
+ if ((Method == MultiCoil) &&
+ (mode & ONLY_COILED_CLASSES) && (tablenum <2) &&
+ (all_scores[2+3][i][j]>bound))
+ all_scores[tablenum][i][j]/= all_scores[2+3][i][j];
+ }
+ }
+/****************************************************************************/
+
+
+
+
+ return (double*)all_scores;
+
+}
+
+
+
+
+ShowSeq(char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ double *m, double *b,
+ double *m_single, double *b_single,
+ int mode, int rescore_sequence, int rescore_preproc, double bound)
+{
+ double maxscore;
+
+/* fprintf(stderr,"\nIn show seq"); */
+
+
+ if ( (rescore_sequence) || (rescore_preproc))
+ all_scores = ScoreSeq(lib, differ_lib, multi_lib, m, b,m_single, b_single,
+ mode, table,
+ method, preprocessor_method,
+ offset_to_use, &maxscore,
+ rescore_sequence, rescore_preproc, bound);
+
+
+/* fprintf(stderr,"\n Got all_scores in showseq"); */
+
+ scores = &all_scores[table * MAXSEQLEN*(POSNUM+1)];
+ /* This is all_scores[main_table][0][0] */
+ preprocess_like = &all_preprocess_like[preproc_table][0][0];
+
+ ShowScores(scores,&sequence,NULL, bound,mode);
+
+}
+/*---------------------------------------------------------------------------*/
+
+
+Font font;
+
+OpenScreen(int mode, double log_bound, FILE *fin, FILE *flog,
+ FILE* flog_coil_conflicts,FILE *fout_coils,
+ FILE *fout,
+ double *m, double *b, double *m_single, double *b_single,
+ char lib[MAXFUNCTNUM],
+ char differ_lib[MAXFUNCTNUM],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ int main_method,
+ int main_preprocessor_method,
+ int *seqcnt, int which_priors, int which_priorp,
+ double prior_freq_single[MAX_TABLE_NUMBER],
+ double prior_freq_pair[MAX_TABLE_NUMBER],
+ int good_turing_fixed_disc, int structural_pos[POSNUM+1])
+{
+ int WIDTH, HEIGHT;
+ XGCValues values;
+ XSetWindowAttributes attr;
+ Pixmap icon;
+ static unsigned char dashl[2] = {1,1};
+ static unsigned char dashd[2] = {2,2};
+ int i;
+
+ WIDTH = 597;
+ HEIGHT = 729;
+ display = XOpenDisplay("");
+
+ /* Check mode to see if want GUI. */
+ if (!display || (mode & NO_GUI)) {
+ /** COmplete logfile, etc. if can't open xdisplay. */
+ if (!(mode & NO_GUI))
+ fprintf(stderr,"Sorry, bad display.\nMaybe your DISPLAY environment variable isn't set properly,\nor you haven't xhosted this machine.\n");
+ if ( fout || ( (mode & PRN_MODE) || (mode & VER_MODE)) ) {
+ fprintf(stderr,"\nThe log file will still be written by:\n\n");
+ for (i = 0; motd[i][0]; ++i) /* Print out title. */
+ fprintf(stderr," %s\n",motd[i]);
+ seqnum=sequp;
+ finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
+ fout,
+ m,b, m_single, b_single,
+ lib,differ_lib,multi_lib, main_method, main_table,
+ main_preprocessor_method, seqcnt, avg_max, which_priors,
+ which_priorp, prior_freq_single, prior_freq_pair,
+ good_turing_fixed_disc, structural_pos);
+ }
+
+ sclose(fin);
+ if (flog) sclose(flog);
+ if (flog_coil_conflicts) sclose(flog_coil_conflicts);
+ if (fout) sclose(fout);
+ if (ftotal_like) sclose(ftotal_like);
+ if (fout_coils)
+ sclose(fout_coils);
+ exit(0);
+ }
+
+ font = XLoadFont(display,"9x15");
+ window = XCreateSimpleWindow(display, RootWindow(display, 0),
+ 0, 0,
+ WIDTH, HEIGHT/4 - 20,
+ 0,
+ BlackPixel(display, 0),
+ WhitePixel(display, 0));
+
+ help = XCreateSimpleWindow(display, RootWindow(display, 0),
+ 0, 0,
+ WIDTH / 5 * 4 +60, HEIGHT / 5 * 3 + 200,
+ 0,
+ BlackPixel(display, 0),
+ WhitePixel(display, 0));
+ helpb = XCreateSimpleWindow(display, window,
+ 90, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ next = XCreateSimpleWindow(display, window,
+ 125, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ prev = XCreateSimpleWindow(display, window,
+ 160, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ autob = XCreateSimpleWindow(display, window,
+ 195, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ zoomb = XCreateSimpleWindow(display, window,
+ 230, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ printb = XCreateSimpleWindow(display, window,
+ 265, 0, 36, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ methodb=XCreateSimpleWindow(display, window,
+ 306, 0, 42, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+/***
+ tableb=XCreateSimpleWindow(display, window,
+ 353, 0, 36, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+***/
+ paramb=XCreateSimpleWindow(display, window,
+ 353, 0, 36, 30, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+ quit = XCreateSimpleWindow(display, window,
+ 393, 0, 30, 15, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+
+/** To get size: **/
+/** Printing " trimer x.xx" takes 12 char = 12*7 =84 **/
+/** Each line takes 15 of y coord. giving 3*15 =45 + 10 for blank line **/
+ seqscoreb=XCreateSimpleWindow(display, window,
+ XCOORD_SEQSCORES, YCOORD_SEQSCORES,
+ 84, 55, 1,
+ BlackPixel(display, 0), WhitePixel(display, 0));
+
+ zoom = XCreateSimpleWindow(display, RootWindow(display, 0),
+ 0, 0,
+ 400, 169,
+ 0,
+ BlackPixel(display, 0),
+ WhitePixel(display, 0));
+ XSelectInput(display, helpb, ButtonPressMask);
+ XSelectInput(display, next, ButtonPressMask);
+ XSelectInput(display, prev, ButtonPressMask);
+ XSelectInput(display, autob, ButtonPressMask);
+ XSelectInput(display, zoomb, ButtonPressMask);
+ XSelectInput(display, printb, ButtonPressMask);
+ XSelectInput(display, methodb, ButtonPressMask);
+/***
+ XSelectInput(display, tableb, ButtonPressMask);
+****/
+ XSelectInput(display, paramb, ButtonPressMask);
+
+ XSelectInput(display, quit, ButtonPressMask);
+ XSelectInput(display, window, ExposureMask);
+ XSelectInput(display, zoom, ExposureMask | KeyPressMask | ButtonPressMask);
+ XSelectInput(display, help, ExposureMask | KeyPressMask | ButtonPressMask);
+ attr.bit_gravity = NorthWestGravity;
+ XChangeWindowAttributes(display,window,CWBitGravity,&attr);
+ XMapWindow(display, window);
+ XFlush(display);
+
+ values.foreground = BlackPixel(display, 0);
+ values.background = WhitePixel(display, 0);
+ values.font = font;
+ gc = XCreateGC(display, window, GCForeground | GCBackground | GCFont, &values);
+ values.line_style = LineDoubleDash;
+ values.font = XLoadFont(display,"6x13");
+ gcgray = XCreateGC(display, window, GCForeground | GCBackground
+ | GCLineStyle | GCFont, &values);
+ XSetDashes(display, gcgray, 1, dashl, 2);
+ values.function = GXxor;
+ values.foreground = 1;
+ gcdash = XCreateGC(display, window, GCForeground
+ | GCLineStyle | GCFunction, &values);
+ XSetDashes(display, gcdash, 0, dashd, 2);
+
+
+/*
+ icon = XCreateBitmapFromData(display,window,bart_bits,bart_width,bart_height);*/
+ /* XSetStandardProperties(display,window,"Coiled Coil Finder","xcoil",icon,NULL,0,NULL);*/
+/* XSetStandardProperties(display,zoom,"XCoil Zoomer","xcoil",icon,NULL,0,NULL);*/
+/* XSetStandardProperties(display,help,"XCoil Help","xcoil",icon,NULL,0,NULL);*/
+ XFlush(display);
+}
+
+PSDefs (FILE *ps, int seqlen)
+{
+ double yscaler=1;
+ double xscaler;
+ int number_lines;
+ int nominallength= 300; /* Default value for number residues per line. */
+
+ if ( (seqlen%nominallength <= 50) && (seqlen >= nominallength)) {
+ nominallength = 350;
+ }
+
+/*************************/
+ if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
+ if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
+ nominallength = seqlen;
+ else
+ nominallength = ps_res_per_line; /* How many res per line. */
+/************************/
+
+ if (seqlen ==0) number_lines=1;
+ else
+ number_lines = (int)(seqlen-1)/nominallength +1;
+ /* nominallength amino acids per line. */
+ /* C rounds down and we want 0-nominallength to go to 1 */
+ /* nominallength+1 to 2*noiminallength to go to 2.. */
+
+
+ fputs("%!\n",ps);
+
+/*****Scale and translate to fit the page if want it to rotate page. */
+/* unscaled there are 9 lines per page. */
+/* .89 scaling if 300 res per line.
+ xscaler = .89 * 300/nominallength;
+ fputs("90 rotate\n",ps);
+ fputs("78 -60 translate\n",ps);
+ if (number_lines > 9) yscaler = (double) 9/number_lines;
+ fprintf(ps,"%lf %lf scale\n", xscaler, yscaler);
+*/
+
+/*****Scale and translate to fit the page if don't want to rotate page. */
+/* unscaled there can fit 12 lines per page, so 3600. */
+ xscaler = .75 * 300/nominallength; /* .75 scaling if 300 res per line. */
+ fputs("25 740 translate\n",ps);
+ if (number_lines > 12) yscaler = (double) 12/number_lines;
+ fprintf(ps,"%lf %lf scale\n", xscaler,yscaler);
+
+
+
+ fputs("/black {0 setgray} def\n",ps);
+ fputs("/gray {0.6 setgray} def\n",ps);
+ fputs("/tt /Courier findfont 15 scalefont def\n",ps);
+ fputs("/rm /Times-Roman findfont 12 scalefont def\n",ps);
+ fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
+ fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
+ fputs("/m {neg moveto} def\n",ps);
+ fputs("/l {neg lineto stroke} def\n",ps);
+ fputs("/rl {neg rlineto} def\n",ps);
+ fputs("/rl_stroke {neg rlineto stroke} def\n",ps);
+ fputs("/bar {dup 0.0019 lt {pop} {0 exch rlineto stroke} ifelse} def\n",ps);
+}
+
+/*** static int xmin = 60; ***** Made it global variable *****/
+ /*Originally was 29 */
+ /** xmin determines where in window start seq. **/
+static int nominallength;
+static int length;
+static int ybase;
+static int ydelta=57;
+/* If click at x,y, determine which residue position it corresponds to.
+ Return 0 if none.
+ */
+XYtoPos (int x, int y)
+{
+ int i;
+ if (x<xmin) return 0;
+ i = (x-xmin)/2;
+ if (i>= nominallength) return 0;
+ for (y -= ybase+ydelta; y>0; y-=ydelta, i+=nominallength);
+ if (y<-30) return 0;
+ if (i>=length) return 0;
+ return i+1;
+}
+PostoXY (int i, int *x, int *y)
+{
+ --i;
+ *x = xmin + 2*(i%nominallength);
+ *y = ybase+ydelta + ydelta*(i/nominallength);
+}
+psstr (char *str, int len, FILE *ps)
+{
+ int i;
+ for (i=0; i<len; ++i)
+ switch (str[i]) {
+ case '(':
+ case ')':
+ case '\\':
+ fputc('\\',ps);
+ default:
+ fputc(str[i],ps);
+ }
+}
+
+
+/*** ybase determines where top of the seq score box should be so don't ***/
+/*** print over the sequenc title. **/
+draw_seqscore_box(double bound, double seq_scores[MAX_CLASS_NUMBER],
+ Window *seqscoreb, FILE *ps, int ybase)
+{
+ char buf[14];
+/* XWindowChanges values; */
+
+
+
+ XClearWindow(display,*seqscoreb);
+
+ /** Relocate seqscoreb base on ybase**/
+/**This is complex way to move it.****
+ values.y= YCOORD_SEQSCORES+ybase-30;
+ XConfigureWindow(display,*seqscoreb,CWY,&values);
+***/
+
+ XMoveWindow(display,*seqscoreb,XCOORD_SEQSCORES, YCOORD_SEQSCORES+ybase-30);
+
+
+
+/** 3 means print 3 into the box horiz, 10 means print 10 into box vert? */
+
+ XMapWindow(display, *seqscoreb);
+ XDrawString(display,*seqscoreb,gcgray,3,10,"Seq Prob:",9);
+
+ if (ps) { /* Draw the box. **/
+ fprintf(ps,"%d %d m\n",10, 50+ybase-30);
+ fprintf(ps,"0 40 rl\n");
+ fprintf(ps,"79 0 rl\n");
+ fprintf(ps,"0 -40 rl\n");
+ fprintf(ps,"-79 0 rl_stroke\n");
+ }
+
+ if (seq_scores[0] + seq_scores[1] + seq_scores[2] > 0) { /* Signals above */
+ sprintf(buf," Dimer %5.2lf", seq_scores[0]); /* bound. */
+ XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
+ sprintf(buf," Trimer %5.2lf", seq_scores[1]);
+ XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
+
+ if (ps) {
+ fprintf(ps,"rm setfont ");
+ fprintf(ps,"%d %d m (Seq Prob for) show\n", 15,60 +ybase-30);
+
+ fprintf(ps,"%d %d m (bound %5.2lf:) show\n", 15,70 +ybase-30,bound);
+
+ fprintf(ps,"%d %d m ( Dimer %5.2lf) show\n",15,80 +ybase-30,
+ seq_scores[0]);
+ fprintf(ps,"%d %d m ( Trimer %5.2lf) show\n",15,90 +ybase-30,
+ seq_scores[1]);
+ }
+
+
+ }
+ else {
+ sprintf(buf," No Residue ");
+ XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
+ sprintf(buf," above %5.2lf ", bound);
+ XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
+
+ if (ps) {
+ fprintf(ps,"rm setfont ");
+ fprintf(ps,"%d %d m (Seq Prob:) show\n", 15,60 +ybase-30);
+
+ fprintf(ps,"%d %d m ( No Residue) show\n", 15, 75 +ybase-30);
+ fprintf(ps,"%d %d m ( above %5.2lf) show\n",15,85 +ybase-30, bound);
+ }
+
+ }
+
+
+}
+
+draw_bound_box(double bound, Window paramb)
+{
+ char buf[6];
+
+ XMapWindow(display, paramb);
+ XDrawString(display,paramb,gcgray,3,10,"bound",5);
+ sprintf(buf,"%5.2lf", bound);
+ XDrawImageString(display,paramb,gcgray,3,25,buf,strlen(buf));
+}
+
+ShowScores (double *scores, Sequence *sequence, FILE *ps, double bound,
+ int mode)
+{
+ int i,q, x, y0, xmax;
+ int width;
+ int rows;
+ int charperline;
+ int xcoord;
+ char MethodTitleArray[500]; /* For printing the methodname to ps file. */
+ char *MethodTitle = MethodTitleArray;
+ char PreprocTitleArray[500]; /* For printing the preproc.name to ps file. */
+ char *PreprocTitle = PreprocTitleArray;
+
+
+ MethodTitleArray[0] = '\0'; /* Initialization. */
+ PreprocTitleArray[0]='\0';
+
+ nominallength = 300;
+ length = sequence->seqlen;
+
+/* This makes display look nice (so don't have a line with < 50 residues. */
+ if ((length%nominallength <= 50) && (length>= nominallength)){
+ nominallength = 350;
+/* if (ps) fputs("-50 0 translate\n",ps); */
+ }
+
+
+/*************************/
+ if (ps) /* Do this stuff only for postscript printout. */
+ if (ps_res_per_line != 0) /* 0 signifies do the default as above. */
+ if (ps_res_per_line == -1) /* -1 signifies put the seq all on 1 line. */
+ nominallength = sequence->seqlen;
+ else
+ nominallength = ps_res_per_line; /* How many res per line. */
+/************************/
+
+
+ charperline = (xmin+2*nominallength-5)/9; /*Each char is "9" units witdth */
+ rows = (length-1)/nominallength+1; /* Number of rows to display. */
+ if (rows==1)
+ xmax = 2*length+xmin; /* The max x position in current row. */
+ else
+ xmax = 2*nominallength+xmin;
+ width = xmin+2*nominallength+ WIDTH_ADD_ON; /* Width of window to display. */
+ /* The +WIDTH_ADD_ON used to be +11, but */
+ /* was cutting of some numbers. */
+ x = strlen(sequence->title);
+ XClearWindow(display,window);
+
+/**********This is where window size is adjusted to fit the sequence. ***/
+/********** using width and number of rows. **************************/
+
+ XResizeWindow(display,window,width,
+ 38+15*((x+charperline-1)/charperline) +ydelta*rows);
+
+/*************************************************************************/
+ XMapWindow(display, helpb);
+ XMapWindow(display, next);
+ XMapWindow(display, prev);
+ XMapWindow(display, autob);
+ XMapWindow(display, zoomb);
+ XMapWindow(display, printb);
+ XMapWindow(display, methodb);
+/****
+ XMapWindow(display, tableb);
+ XMapWindow(display, paramb); Now done in draw_bound_box()
+******/
+ XMapWindow(display, quit);
+ XSync(display,0);
+ { XEvent report;
+ while (XCheckTypedEvent(display,Expose,&report));
+ }
+ XDrawString(display,helpb,gcgray,3,10,"help",4);
+ XDrawString(display,next,gcgray,3,10,"next",4);
+ XDrawString(display,prev,gcgray,3,10,"prev",4);
+ XDrawString(display,autob,gcgray,3,10,"auto",4);
+ XDrawString(display,zoomb,gcgray,3,10,"zoom",4);
+ XDrawString(display,printb,gcgray,3,10,"print",5);
+ XDrawString(display,methodb,gcgray,3,10,"method",6);
+
+/* XDrawString(display,tableb,gcgray,3,10,"table",5); */
+
+ draw_bound_box(bound, paramb);
+/** draw_seqscore_box(bound, seq_scores,?? &seqscoreb, ps, ybase);
+ MOVED TO LATER (AFTER ybase gets adjusted at end of this function). ****/
+
+
+ XDrawString(display,quit,gcgray,3,10,"quit",4);
+
+
+/************This section prints out the method name..... ********************/
+/**** Add 9 to xcoord for each character. *****/
+ XDrawString(display,window,gc, XCOORD_SCORED_BY,15,"Scored by ",10);
+
+ make_methodtitle(MethodTitle, method, table, mode);
+ XDrawString(display,window,gc,XCOORD_SCORED_BY +90,15,MethodTitle,
+ strlen(MethodTitle));
+
+ XDrawString(display,window,gc, XCOORD_SCORED_BY,30,"Filter by ",10);
+ make_methodtitle(PreprocTitle, preprocessor_method, preproc_table,
+ mode);
+ XDrawString(display,window,gc,XCOORD_SCORED_BY + 90,30,
+ PreprocTitle,strlen(PreprocTitle));
+
+/***************************************************************************/
+
+ XDrawString(display,window,gc, 5,30,sequence->code,strlen(sequence->code));
+ for (i=0,ybase=30; x > i + charperline; i += charperline)
+ XDrawString(display,window,gc, 5,ybase+=15,&(sequence->title[i]),charperline);
+ XDrawString(display,window,gc, 5,ybase+15,&(sequence->title[i]),x-i);
+
+ if (ps) {
+ fprintf(ps,"tt setfont 5 15 m (%s) show\n",sequence->code);
+/* fprintf(ps,"175 15 m (Viewgram by XCoil) show\n"); */
+/*
+ fprintf(ps,"430 15 m (Scored by %s) show\n",methodname[method]);
+*/
+ fprintf(ps,"180 15 m (Scored by %s) show\n",MethodTitle);
+ fprintf(ps,"180 26 m (Filter by %s) show\n",PreprocTitle);
+
+
+ for (i=0,ybase=30; x > i + charperline; i += charperline) {
+ fprintf(ps,"5 %d m (",ybase+=15);
+ psstr(&(sequence->title[i]),charperline,ps);
+ fprintf(ps,") show\n");
+ }
+ fprintf(ps,"5 %d m (",ybase+15);
+ psstr(&(sequence->title[i]),x-i,ps);
+ fprintf(ps,") show\n");
+ fputs("rm setfont 2 setlinewidth\n",ps);
+ }
+ for (i=0,y0=ybase; i<length; ++i,x+=2) {
+ double y;
+ if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
+ y = y0-30*scores[8*i+7];
+ XDrawLine(display, window, gc, x, y0, x, (int)(y+0.5));
+ XDrawLine(display, window, gc, x+1, y0, x+1, (int)(y+0.5));
+ if (ps) {
+ fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
+ }
+ }
+/********************* If in pos mode draws in the correct likelihoods every
+now and then? ******************
+ if (sequence->reg && method==HMMCoil) {
+ if (ps)
+ fputs("gray\n",ps);
+ for (i=0,y0=ybase; i<length; ++i,x+=2) {
+ if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
+ if (sequence->reg[i] != '.') {
+ double y = y0-30*scores[8*i+sequence->reg[i]];
+ XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
+ XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
+ if (ps)
+ fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
+ }
+ }
+ if (ps) fputs("black\n",ps);
+ }
+**********************/
+
+/**********Now draw in the PreProcScore in gray every other line. ****/
+
+ if (ps)
+ fputs("gray\n",ps);
+ for (i=0,y0=ybase; i<length; i+=2,x+=4) { /* Normally is ++i,x+=4, but */
+ double y; /* we do every other line here. */
+ if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
+
+ y = y0-30*preprocess_like[i*(POSNUM +1) + 7];
+ /* This is preprocess_like[i][7]. */
+ XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
+ XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
+ if (ps)
+ fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
+ }
+
+ if (ps) fputs("black\n",ps);
+
+/*************************************************************************/
+
+ if (ps) fputs("0 setlinewidth\n",ps);
+ for (i=9,y0=ybase; i<length; i+=10) {
+ if (!((i-9)%nominallength)) {x=xmin-1; y0 += ydelta;}
+ x += 20;
+ XDrawLine(display, window, gc, x,y0,x,y0+2); /* Draw ticks every */
+ if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+2); /* 10 residues */
+
+ }
+ for (i=49,y0=ybase; i<length; i+=50) {
+ char buff[6];
+ if (!((i-49)%nominallength)) {x=xmin-1; y0 += ydelta;}
+ x += 100;
+ XDrawLine(display, window, gc, x,y0,x,y0+5);
+ if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+5);
+ sprintf(buff,"%d",i+1);
+ XDrawString(display, window, gcgray, x-7, y0+18, buff, strlen(buff));
+ if (ps) fprintf(ps,"%d %d m (%d) cshow\n",x,y0+18,i+1);
+ }
+ for (i=1,y0=ybase+ydelta; i<=rows; ++i,y0+= ydelta) {
+ if (i==rows)
+ xmax = 2*(length-(rows-1)*nominallength)+xmin; /* Max x position in */
+ /* last row. */
+ XDrawString(display,window,gcgray,xmin-25,y0-26,"100%",4);
+ XDrawString(display,window,gcgray,xmin-19,y0-11,"50%",3);
+ XDrawString(display,window,gcgray,xmin-13,y0+4,"0%",2);
+ XDrawLine(display,window,gcgray,xmin,y0-30,xmax-1,y0-30);
+ XDrawLine(display,window,gcgray,xmin,y0-15,xmax-1,y0-15);
+ if (ps) {
+ fprintf(ps,"%d %d m (100%%) rshow\n",xmin-1,y0-26);
+ fprintf(ps,"%d %d m (50%%) rshow\n",xmin-1,y0-11);
+ fprintf(ps,"%d %d m (0%%) rshow\n",xmin-1,y0+4);
+ fprintf(ps,"%d %d m %d %d l\n",xmin,y0-30,xmax,y0-30);
+ fprintf(ps,"%d %d m %d %d l\n",xmin,y0-15,xmax,y0-15);
+ fprintf(ps,"%d %d m %d %d l\n",xmin,y0,xmax,y0);
+ }
+ }
+
+
+ draw_seqscore_box(bound, seq_scores, &seqscoreb, ps, ybase);
+
+ if (zoomptr) ZoomPnt();
+ XFlush(display);
+ if (ps) {
+ fputs("showpage\n",ps);
+ fflush(ps);
+ }
+}
+
+
+
+
+
+make_methodtitle(char *MethodTitle, int Method, int which_table,
+ int mode)
+{
+ char buf[4];
+ int xcoord;
+ int like2or3;
+
+
+ if (which_table ==0) like2or3 =2;
+ else like2or3=3;
+
+
+/************This section prints out the method name..... ********************/
+/**** Add 9 to xcoord for each character. *****/
+
+
+ xcoord =520;
+ if ((Method == PairCoilDiff) || (Method == PairCoilDiffAvg)) {
+ /* Does Differen. use dimer or trimer like. */
+ if (like2or3==2) {
+/* XDrawString(display,window,gc,xcoord,15,"Dimer",5); */
+ MethodTitle= strcat(MethodTitle, "Dimer");
+ xcoord += 45;
+ }
+ else if (like2or3==3) {
+/* XDrawString(display,window,gc,xcoord,15,"Trimer",6); */
+ MethodTitle= strcat(MethodTitle, "Trimer");
+ xcoord += 54;
+ }
+ }
+
+
+ if (Method == PairCoil) {
+ if (mode & PAIRCOIL_PAIRS) { /* Are we using pair or single prob */
+/* XDrawString(display,window,gc,xcoord,15,"Pair",4); */
+ MethodTitle= strcat(MethodTitle, "Pair");
+ xcoord += 36;
+ }
+ else {
+/* XDrawString(display,window,gc,xcoord,15,"Single",6); */
+ MethodTitle= strcat(MethodTitle, "Single");
+ xcoord += 54;
+ }
+
+ } /* Ends the stuff for dealing with PairCoil, TDRes, TDCoil. */
+
+
+/*****************************************************************************/
+/*************************** Now write the methodname. *********************/
+
+/* XDrawString(display,window,gc,xcoord,15,methodname[method],strlen(methodname[method])); */
+ xcoord += 9*strlen(methodname[Method]);
+ MethodTitle= strcat(MethodTitle, methodname[Method]);
+
+ if (Method == MultiCoil) {
+ if (mode & ONLY_COILED_CLASSES)
+ strcat(MethodTitle,"Olig");
+ if (Coil_Score) {
+ if (Coil_Score ==1)
+ if (weighted_avg)
+ strcat(MethodTitle,"WtdAvg");
+ else strcat(MethodTitle,"Avg");
+
+ else if (Coil_Score ==2)
+ strcat(MethodTitle, "MAX");
+
+ if (start_coil_at_reg_shift)
+ strcat(MethodTitle, "Shift");
+ else
+ strcat(MethodTitle, "NoShifts");
+ }
+ }
+
+ /* Now print Table are using. */
+ /* And if are using a particular offset, print that. */
+ if ( (Method==MultiCoil) || (Method==PairCoil) || (Method == HMMCoil)
+ || (Method==PairCoilDiff) || (Method==PairCoilDiffAvg)) {
+
+
+ if (number_tables >1) {
+ sprintf(buf,"%d",which_table+1);
+/* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
+ xcoord += 9*strlen(buf);
+ MethodTitle= strcat(MethodTitle, buf);
+ }
+ }
+
+ sprintf(buf,""); /* Clear out buf. */
+
+
+ if ((offset_to_use != -1) && (offset_to_use != 7))
+ sprintf(buf,"%c",'a' + offset_to_use);
+ else if (Method == MultiCoil)
+ if (offset_to_use == -1) {
+ strcat(MethodTitle,"Max");
+ xcoord += 9*3;
+ }
+ else if (offset_to_use == 7) {
+ strcat(MethodTitle,"Combo");
+ xcoord += 9*5;
+ }
+
+ /* XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
+ xcoord += 9*strlen(buf);
+ MethodTitle= strcat(MethodTitle, buf);
+
+}
+
+/***************************************************************************/
+
+
+int is_not_differentiator(int Method)
+{
+ if ( (Method == MultiCoil) || (Method == PairCoil) || (Method==HMMCoil) ||
+ (Method==NEWCOILS) || (Method == ActualCoil) )
+ return (1);
+ else return(0);
+}
+
+
+
+
+
+log_file_header(FILE *flog, FILE *flog_coil_conflicts,
+ int mode, int argc, char *argv[],
+ int avg_max, double bound,
+ int preprocessor_method, int preproc_table,
+ int log_offset_to_use, int start_coil_at_reg_shift,
+ int weighted_avg, int main_method, int main_table)
+{
+ int i;
+ char MethodTitleArray[500]; /* For printing the methodname to ps file. */
+ char *MethodTitle = MethodTitleArray;
+
+ MethodTitleArray[0] = '\0'; /* Initialization. */
+ make_methodtitle(MethodTitle, method,main_table,mode);
+
+ if (flog) {
+ fprintf(flog,"\nInput file is %s", argc>1?argv[1]:"-");
+ fprintf(flog, "\n\nMethod Used is %s", MethodTitle);
+ fprintf(flog, "\nThe bound for the coil locater is %.2lf",bound);
+ }
+ if (flog_coil_conflicts) {
+ fprintf(flog_coil_conflicts,"\nInput file is %s", argc>1?argv[1]:"-");
+ if (mode & TST_MODE0)
+ fprintf(flog_coil_conflicts,"\nSubtracting sequence off from cctb before scoring it.");
+ if (mode & TST_MODE1)
+ fprintf(flog_coil_conflicts,"\nSubtracting sequence off from trimer table before scoring it.");
+ }
+
+
+ if (mode & WEIGHTED_PROBS) {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\nStatistically weighted PairCoil score is ");
+ if (flog)
+ fprintf(flog,"\nStatistically weighted PairCoil score is ");
+ }
+ else {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\nPairCoil score is ");
+ if (flog)
+ fprintf(flog,"\nPairCoil score is ");
+ }
+ if (avg_max ==0) {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"AVERAGE of residue scores.");
+ if (flog)
+ fprintf(flog,"AVERAGE of residue scores.");
+ }
+ else {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"MAX of residue scores.");
+ if (flog)
+ fprintf(flog,"MAX of residue scores.");
+ }
+
+ if (mode & USE_LIKE_LINE) {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\nPaircoil used the likelihood line.");
+ if (flog)
+ fprintf(flog,"\nPaircoil used the likelihood line.");
+ }
+ else {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\nPaircoil used direct prob. formula instead of likelihood line.");
+ if (flog)
+ fprintf(flog,"\nPaircoil used direct prob. formula instead of likelihood line.");
+ }
+
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts, "\nThe distances used by PairCoil are: ");
+ if (flog)
+ fprintf(flog, "\nThe distances used by PairCoil are: ");
+
+ for (i=0; i<pair_functnum; i++) {
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"%d ",i);
+ if (flog)
+ fprintf(flog,"%d ",i);
+ }
+
+ if (flog_coil_conflicts) {
+ fprintf(flog_coil_conflicts, "\n\nMethod Used is %s", MethodTitle);
+ fprintf(flog_coil_conflicts, "\nThe bound for the coil locater is %.2lf",bound);
+ }
+
+ if ( (main_method==PairCoilDiff) || (main_method==PairCoilDiffAvg)) {
+ char MethodTitleArray[500]; /* For printing the methodname. */
+ char *MethodTitle = MethodTitleArray;
+ int oligimerization;
+
+ MethodTitleArray[0] = '\0'; /* Initialization. */
+ make_methodtitle(MethodTitle, preprocessor_method, preproc_table,mode);
+
+ if (main_table ==0)
+ oligimerization == 2;
+ else oligimerization ==3;
+
+ if (flog)
+ fprintf(flog, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
+ /* Whether logfile is for 2 or */
+ /* 3 stranded predictions. */
+
+ if (main_method==PairCoilDiffAvg) {
+ if (weighted_avg == 0) {
+ if (flog)
+ fprintf(flog, "\n Averaging over coils was done by coil length.");
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts, "\n Averaging over coils was done by coil length.");
+ }
+ else { /* weighted_avg ==1 */
+ if (flog)
+ fprintf(flog, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
+ }
+
+ if (start_coil_at_reg_shift) {
+ if (flog)
+ fprintf(flog,"\n New coils were considered to start at register shifts.");
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\n New coils were considered to start at register shifts.");
+ }
+ else {
+ if (flog)
+ fprintf(flog,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
+ if (flog_coil_conflicts)
+ fprintf(flog_coil_conflicts,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
+ }
+ }
+ }
+
+ if (flog) {
+ fprintf(flog,"\nMain Table is Table %d\n\n",main_table);
+
+ if (mode & POS_MODE)
+ fprintf(flog,"\nIf the start or end of a coil is listed inside (), then that position does not occur in the pos file.\n");
+ }
+}
+
+
+
+
+
+coil_tuple_output(int number_tables, char *lib, int mode, int avg_max,
+ double bound, Sequence sequence, FILE *fout_coils,
+ int number_multi_lib[MAX_TABLE_NUMBER],
+ int number_score_dim)
+{
+ char offsets[MAXSEQLEN];
+ int dist;
+ int tab;
+ int i,j;
+ extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
+ extern double *all_scores;
+
+
+ for (tab=0; tab< number_tables; tab++) {
+ for (dist=0; dist<number_multi_lib[tab]; dist++) {
+
+/*************************************************************************/
+/**********This computes offset of that table/dist method of scoring. ***/
+
+ for (i=0; i<sequence.seqlen; i++) {
+ for (j=0; j<POSNUM; j++)
+ if (all_preprocess_like[tab + number_tables*(int)dist][i][j] ==
+ all_preprocess_like[tab + number_tables*(int)dist][i][7]) {
+ offsets[i] = (j-i) % POSNUM;
+ if (offsets[i]<0) offsets[i] += POSNUM;
+ }
+ /* This is a new thing so that if an offset change is not real */
+ /* we keep the old offset. **********/
+ if ( (i>0) && (offsets[i-1] != -1) &&
+ (offsets[i] != offsets[i-1]) &&
+ (all_preprocess_like[tab + number_tables*(int)dist]
+ [i][(i+offsets[i-1]) %POSNUM]
+ == all_preprocess_like[tab+number_tables*(int)dist][i][7]) )
+ offsets[i] = offsets[i-1];
+ }
+
+/**************************************************************************/
+/****** Now compute coil scores for that method (table,dist) pair. ***/
+ get_raw_coil_score(sequence,
+ all_preprocess_like[tab+number_tables*(int)dist],
+ all_scores[tab+number_tables*(int)dist], mode,
+ avg_max,bound,offsets);
+ }
+ }
+
+/************* Now output the tuples of coil scores. ****/
+ tuple_output2(sequence, mode, fout_coils,
+ all_scores, number_tables,
+ number_multi_lib,bound, number_score_dim);
+}
+
+
+
+switch_gauss_param(int Offset_to_Use, int Method)
+{
+ if (Method == MultiCoil) {
+ if (Offset_to_Use == 7) { /* Combo offset */
+ determinants= &both_determinants[0][0];
+ inverse_covars= &both_inverse_covars[0][0][0][0];
+ means_submatrix= &both_means_submatrix[0][0][0];
+ }
+ else { /* max offset or some fixed offset */
+ determinants= &both_determinants[1][0];
+ inverse_covars= &both_inverse_covars[1][0][0][0];
+ means_submatrix= &both_means_submatrix[1][0][0];
+ }
+ }
+ else if (Method == PairCoilDiff) {
+ determinants= &differ_determinants[0];
+ inverse_covars= &differ_inverse_covars[0][0][0];
+ means_submatrix= &differ_means_submatrix[0][0];
+ }
+
+
+}
+
+
+void usage(char *prog_command)
+{
+ int i;
+ char *prog_name;
+
+/** Strip off directories to get just the program name. **/
+ i=strlen(prog_command);
+ while ((i>=0) && (prog_command[i] != '/') )
+ i--;
+ prog_name= &prog_command[i+1];
+
+ fprintf(stderr,"\nUsage: %s <file-name> [dimer, trimer, dimer-1, trimer-1, negatives, config_file] [multi_lib_dim1] [multi_lib_dim2] \n",prog_name);
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
git://source.jalview.org / jpred.git / blobdiff