1 /* Interface to sc2seq created by David Wilson and modified by Ethan Wolf */
2 /* Other interfaces by Ethan Wolf 1995. */
10 #include "likelihood.h"
14 #include "compute_like.h"
15 #include "interface.h"
18 #include "debug_printout_flag.h"
20 #define PIR_LIKES "/tmp/pir.likes"
21 #define PS_OF_LIKELIHOOD "~/likelihood.ps"
22 #define PAIR_PIR_TXT "/tmp/pir1.txt"
23 #define SINGLE_PIR_TXT "/tmp/pir2.txt"
25 double many_pprobs[MAX_TABLE_NUMBER][NUM_RES_TYPE][POSNUM],
26 many_pprobp[MAX_TABLE_NUMBER][NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM];
29 double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];
31 double (*pprobs)[POSNUM];
32 double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];
35 long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
36 long ptotals[MAX_TABLE_NUMBER][POSNUM];
37 long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
38 long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
42 #define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
50 ptr = (void*)malloc(n);
52 fprintf(stderr,"Memory shortage!!\n");
59 void switch_tables(int table)
61 extern int WINDOW, window_length[MAX_TABLE_NUMBER];
63 pprobs= many_pprobs[table];
64 pprobp= many_pprobp[table];
65 WINDOW = window_length[table];
70 /******** Used to subtract off a sequence from the table. ***/
72 void recalc_prob(Sequence sequence, int input_table, int mode, long pfreqp2[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
74 int ProlineFreeWindow = mode & PROLINE_FREE_MODE;
76 /* Remove the sequence from the positive table before computing on it. */
77 add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
79 sequence.reg, sequence.seq, sequence.seqlen, -1);
81 if (!(mode & WEIGHTED_PROBS))
82 estimate_database_probs2(input_table, pfreqs[input_table],
83 ptotals[input_table], many_pprobs,
84 pfreqp[input_table], ptotalp[input_table],
85 many_pprobp, ProlineFreeWindow);
86 /** Add seq back into the frequence counts. **/
87 add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
89 sequence.reg, sequence.seq, sequence.seqlen, 1);
91 #ifdef AVG_WEIRD_PROBS
92 calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
94 calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
101 void normalize_gen(double gprobs[NUM_RES_TYPE],
102 double gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM])
104 int res1, res2, dist;
106 for (res1 =0; res1<Proline; res1++) {
107 gprobs[res1] = getdlog(exp(gprobs[res1])/(1-exp(gprobs[Proline])));
108 for (res2= 0; res2<Proline; res2++)
109 for (dist=0;dist<7;dist++) {
110 /** Subtract off for Proline with anything in other position once **/
111 /** for each position, then add back on the extra bit subtracted. **/
112 gprobp[res1][res2][dist] = getdlog(exp(gprobp[res1][res2][dist]) /
113 ((1-exp(gprobs[Proline])-exp(gprobs[Proline]) +
114 exp(gprobp[Proline][Proline][dist]))) );
123 void PairCoilData (FILE *fgin, FILE **fpin,
124 int ProlineFree, int table_to_remove_from,
125 char *input_filename, int mode,
126 long pfreqp2[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
128 long gfreqs[AANUM]; long gtotals;
129 long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM];
131 extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
132 extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
133 //extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
134 extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
135 extern double SCALE0;
141 readgenfreq(fgin, >otals, gfreqs, gtotalp, gfreqp, stderr);
142 calcgprobs(gfreqs, gtotals, gprobs);
143 calcgprobp(gfreqp, gtotalp, gprobp);
145 normalize_gen(gprobs, gprobp);
148 else { setgprob_to_uniform(gprobs, gprobp); }
150 #ifdef AVG_WEIRD_PROBS
151 calc_weird_gprob_avg(gprobp, gprobs);
153 calc_weird_gprob_sum(gprobp, gprobs);
157 for (i=0; i<NUMBER_TABLES; i++) { /* Go through all input tables. */
158 if (fpin[i] != NULL) { /* Calculate probabilities for that table. */
159 readposfreq(fpin[i],ptotals[i],pfreqs[i], ptotalp[i], pfreqp[i], stderr);
161 if (i == table_to_remove_from)
162 remove_all_input_seq_from_table(pfreqs[table_to_remove_from],
163 ptotals[table_to_remove_from],
164 pfreqp[table_to_remove_from],
165 ptotalp[table_to_remove_from], input_filename, -1);
167 estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
168 pfreqp[i], ptotalp[i], many_pprobp,
172 fprintf(stderr,"\nCould not make input table number %d.",i);
177 #ifdef AVG_WEIRD_PROBS /* make table. */
178 calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);
180 calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);
192 initialize_relative();
195 #define NEGCCRATIO 30
200 #define SQRT2PI 2.5066282746310002
204 /** To compute the paircoil likelihood line. */
207 void likelihood_from_txt(char *txt_filename, double *m, double *b1)
211 double percent_coil=PERCENT_COILED;
214 long n[2]={0,0}, n_halfgauss[2]={0,0};
215 double sum[2]={0,0}, sum2[2]={0,0};
216 double a[2],b[2]; /* ax+b is the likelihood line approximation */
217 double mean[2], sd[2], down[2];
218 long _counts[2][-4*MINSCORE +1];
219 double _means[2][-4*MINSCORE+1];
222 int i, maxplaces[2]= {2*MINSCORE, 2*MINSCORE}, minplaces[2]={-2*MINSCORE,-2*MINSCORE};
223 int minplace, maxplace;
224 int number_strange=0, number_high=0;
225 FILE *devnull, *outfile, *outpsfile;
228 counts[0] = &(_counts[0][-2*MINSCORE]);
229 counts[1] = &(_counts[1][-2*MINSCORE]);
230 means[0] = &(_means[0][-2*MINSCORE]);
231 means[1] = &(_means[1][-2*MINSCORE]);
234 mode2 |= RATIO_TRAP_METHOD;
235 mode2 |= PERCENT_COIL_MODE;
239 for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
245 txt_file = sopen (txt_filename, "r");
246 while (fread(&score1,sizeof(double),1,txt_file)) {
247 if (score1<MINSCORE || score1>-(MINSCORE)) {
249 fprintf(stderr,"Strange score: %lf\n",score1);
254 place = ( (int)((score1-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
255 if (place > -2*MINSCORE) {
256 fprintf(stderr, "Out of bounds: %d.\n", place);
257 fprintf(stderr, "Score: %lf\n",score1);
259 counts[0][-2*MINSCORE]++;
260 means[0][-2*MINSCORE]++;
262 maxplaces[0]=-2*MINSCORE;
264 else if (place < 2*MINSCORE) {
265 fprintf(stderr, "Out of bounds: %d.\n", place);
266 fprintf(stderr, "Score: %lf\n",score1);
267 counts[0][2*MINSCORE]++;
268 means[0][2*MINSCORE]++;
269 minplaces[0]=2*MINSCORE;
275 means[0][place] += score1;
276 if (place>maxplaces[0]) maxplaces[0]=place;
277 if (place<minplaces[0]) minplaces[0]=place;
282 fprintf(stderr,"\n%d out of bounds scores > %d were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE));
283 fprintf(stderr,"\n%d strange scores < %d were scored as %d\n",number_strange-number_high,MINSCORE, MINSCORE);
286 /** Compute the means of the scores in each box. **/
287 for (i= minplaces[0]; i<=maxplaces[0]; ++i)
288 if (counts[0][i] != 0)
289 means[0][i] /= counts[0][i]; /** Each bucket gets average score*/
292 /*** We need a mean***/
294 mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
298 /*** Calculate gaussian fit to scores below the mean. **/
299 txt_file=sopen(txt_filename, "r");
300 while (fread(&score1,sizeof(double),1,txt_file))
301 if ( (score1 > MINSCORE) && (score1 < mean[0])) {
303 sum2[0] += (score1-mean[0])*(score1-mean[0]);
305 sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
306 down[0]= (double) 2*n_halfgauss[0]/n[0];
308 fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
311 minplace=minplaces[0]; maxplace=maxplaces[0];
313 devnull=sopen("/dev/null","r+");
314 outfile=sopen(PIR_LIKES,"w");
315 outpsfile=sopen(PS_OF_LIKELIHOOD,"w");
316 /* outfile=sopen("~/pir.out","w"); */
317 likelihood(n,mode2, filenum, 0.0 , minplace, maxplace, means,
318 counts, mean, sd, down, outfile,outpsfile,devnull,devnull,
323 char clean_command[200];
325 strcpy(clean_command, "rm -f ");
327 /*** The following removes the pir.txt output file. */
328 system(strcat(clean_command,txt_filename));
330 } /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
331 /** WRITE TO SAME FILESNAMES. **/
333 char clean_command[200];
335 strcpy(clean_command, "rm -f ");
337 system(strcat(clean_command,PIR_LIKES));
348 /* Look in file .likelihood for likelihood line. If not there, compute */
350 void get_likelihood_line(char likelihoods[MAX_TABLE_NUMBER][MAXLINE], double *m, double *b, double *m_single, double *b_single, char lib[MAXFUNCTNUM], int functnum, char *pir_name, int mode, int number_tables)
360 FILE *txt_file1, *txt_file2;
361 double score, scs[MAXSEQLEN][POSNUM];
363 char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
364 char title[500], code[500];
367 int found_like_lib=0; int found_single=0;
371 sequence.title=title;
375 for (i=0; i<functnum; i++)
376 libnumber |= (int) pow(2,lib[i]);
377 #ifdef DEBUG_PRINTOUT
378 fprintf(stderr,"\nlibnumber= %d",libnumber);
381 /* What directory to look for paircoil defauts in. Set environment */
382 /* variable LIKELINE_DIR to be a path to this file. */
383 for (likefile=0; likefile <number_tables; likefile++) {
387 #ifdef DEBUG_PRINTOUT
388 fprintf(stderr,"\nLikelihood %d=%s",likefile,likelihoods[likefile]);
390 if (likelihoods[likefile][0] == ',') { /* Don't bother loading likefile */
391 #ifdef DEBUG_PRINTOUT
393 "\n Not doing likelihood conversion for paircoil with table %d\n",
398 if (linefile=sopen(likelihoods[likefile],"r")) {
399 while ((fgets(buf,80,linefile)) && !(found_single && found_like_lib)) {
401 sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
404 else if (!found_like_lib) {
405 sscanf(buf, "%d %lf %lf",¤t_lib, &m[likefile], &b[likefile]);
406 if (current_lib == libnumber) found_like_lib=1;;
411 if (!(found_like_lib && found_single)) {
413 fprintf(stderr,"\nSingles likelihood function %d unknown.",
416 if (current_lib != libnumber ) /* Wasn't in file so should create it
417 by running on pir. Add this later.*/
418 fprintf(stderr, "\nLikelihood function %d unknown...", likefile);
421 if (pir_name[0]==',') {
422 fprintf(stderr,"\nUnable to continue. Please input the location of the pir in the .paircoil file.");
426 fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
428 pir = sopen(pir_name, "r");
430 if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
431 if (!found_single) txt_file2 = sopen(SINGLE_PIR_TXT,"w");
433 for (cnt=1;getseq2(pir,&sequence); cnt++)
434 if (sequence.seqlen >= MINWINDOW) {
435 seqlen= sequence.seqlen;
442 if (!found_like_lib) {
443 switch_tables(likefile);
444 scseqadj(lib,seq,seqlen,scs,offsets, &score,
446 txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
449 switch_tables(likefile);
451 scorestock(seq,seqlen,scs,offsets, &score,
452 mode & PROLINE_FREE_MODE, 0);
453 /* Score stock = singles method using our table. */
454 txt_output(sequence,0,txt_file2, scs, NULL,offsets,0);
457 } /* Ends if sequence > MINWINDOW, and hence the for */
458 /* loop that reads from pir. */
462 if (!found_like_lib) { /* Now compute the lib likelihood line. */
465 likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
467 linefile=sopen(likelihoods[likefile],"a");
468 fprintf(linefile,"%d %lf %lf\n",
469 libnumber,m[likefile], b[likefile]);
471 fprintf(linefile,"\t (lib =");
472 /* Now print human readable lib number. */
473 for (i=0; i<functnum; i++)
474 fprintf(linefile," %d",lib[i]);
475 fprintf(linefile," )\n");
480 if (!found_single) { /* Now compute the singles method like line */
483 likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
484 &b_single[likefile]);
486 linefile=sopen(likelihoods[likefile],"a");
487 fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile],
490 fprintf(linefile,"\t (lib =");
491 /* Now print human readable lib number. */
492 fprintf(linefile," Singles Method");
493 fprintf(linefile," )\n");
497 } /* Ends else from if (pir_name == ',') */
499 } /* Ends section for file not having lib or singles like line */
501 } /* Ends section where search likelihoods[likefile] file. */
503 fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
506 } /** Ends else the likelihood file exists to load a likelihood line... **/
508 } /* Ends count through tables=linefiles. */
513 double stocksc2likelihood (double sc)
518 gg = (sc - NEGMEAN)/NEGSD;
519 gcc = (sc - CCMEAN)/CCSD;
520 gg = exp((-gg * gg) / 2) / NEGSD / SQRT2PI;
521 gcc = exp((-gcc * gcc) / 2) / CCSD / SQRT2PI;
522 return(gcc / (NEGCCRATIO * gg + gcc));
527 /* Likelihood line is mx +b */
528 double sc2likelihood (double sc, double m, double b)
532 likelihood = m*sc + b;
534 if (likelihood>1) return(1);
535 else if (likelihood<0) return(0);
536 else return(likelihood);
540 static int current_seq_number=0;
542 void PairCoilScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double m, double b, double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use, int raw_score)
545 double sc2scores[MAXSEQLEN][POSNUM];
547 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
549 double frac_coiled_this_table;
550 extern double init_class_prob[MAX_CLASS_NUMBER];
552 frac_coiled_this_table = init_class_prob[table];
555 switch_tables(table);
557 /* score the sequence using sc2seq */
559 sequence.seq, sequence.seqlen,
566 if (mode & USE_LIKE_LINE)
567 *maxscore= sc2likelihood(*maxscore,m,b);
569 *maxscore= frac_coiled_this_table * exp(*maxscore);
572 for (i=0; i<sequence.seqlen; ++i) {
573 /* Do the most probabable register if offset_to_use is -1 or 7. */
574 if ((offset_to_use == -1) || (offset_to_use ==7))
575 if (raw_score) /*** Return rawscore. **/
576 scores[i][7] = sc2scores[i][offsets[i]];
578 if (mode & USE_LIKE_LINE)
579 scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
581 scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
582 if (scores[i][7] >1) scores[i][7]=1;
585 else /* Offset to use is 0-6. */
587 scores[i][7] = sc2scores[i][offset_to_use];
589 if (mode & USE_LIKE_LINE)
590 scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
592 scores[i][7] = frac_coiled_this_table *
593 exp(sc2scores[i][offset_to_use]);
594 if (scores[i][7] >1) scores[i][7]=1;
597 /* do other registers */
598 for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
601 scores[i][(i+j)%POSNUM] = sc2scores[i][j];
603 if (mode & USE_LIKE_LINE)
604 scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
606 scores[i][(i+j)%POSNUM] = frac_coiled_this_table *
607 exp(sc2scores[i][j]);
608 if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
619 /***************************************************************************/
621 /** Score using the singles method on our tables. **/
622 double *SingleCoilScore(int mode, Sequence sequence, double m, double b,
623 double *maxscore, int table,
624 double scores[MAXSEQLEN][POSNUM+1], int offset_to_use,
628 double singlescores[MAXSEQLEN][POSNUM];
629 char offsets[MAXSEQLEN];
631 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
633 double frac_coiled_this_table;
636 switch_tables(table);
638 /* score the thing */
640 sequence.seq, sequence.seqlen,
643 maxscore,ProlineFreeWindow, 0); /** Score using our table. */
647 if (mode & USE_LIKE_LINE)
648 *maxscore= sc2likelihood(*maxscore,m,b);
650 /* fprintf(stderr,"\nmaxscore = %lf", *maxscore);
652 *maxscore= frac_coiled_this_table * exp(*maxscore);
653 /* fprintf(stderr," maxlike = %lf", *maxscore);
657 for (i=0; i<sequence.seqlen; ++i) {
658 /* Do the most likely register. */
659 if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
661 if (mode & USE_LIKE_LINE)
662 scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b);
664 scores[i][7] = frac_coiled_this_table*
665 exp(singlescores[i][offsets[i]]);
666 if (scores[i][7] >1) scores[i][7]=1;
668 else scores[i][7] = singlescores[i][offsets[i]];
670 else /* Use offset_to_use. */
672 if (mode & USE_LIKE_LINE)
673 scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
675 scores[i][7] = frac_coiled_this_table*
676 exp(singlescores[i][offsets[i]]);
677 if (scores[i][7] >1) scores[i][7]=1;
679 else scores[i][7] = singlescores[i][offset_to_use];
681 /* do the other registers */
684 if (mode & USE_LIKE_LINE)
685 scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b);
687 scores[i][(i+j)%POSNUM] = frac_coiled_this_table*
688 exp(singlescores[i][j]);
689 if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
691 else scores[i][(i+j)%POSNUM] = singlescores[i][j]; /*Raw_score*/
697 /***************************************************************************/
701 void MultiCoilDimensionScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use)
704 double sc2scores[MAXSEQLEN][POSNUM];
706 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
709 switch_tables(table);
711 /* score the sequence using sc2seq */
713 sequence.seq, sequence.seqlen,
718 /*** Do a SingleCoil type score. ****/
720 scorestock(sequence.seq, sequence.seqlen,
722 maxscore,ProlineFreeWindow, 0); **/ /** Score using our table. */
725 for (i=0; i<sequence.seqlen; ++i) {
726 /* Do the most probabable register if offset_to_use is -1 or 7. */
727 if ((offset_to_use == -1) || (offset_to_use ==7))
728 scores[i][7] = sc2scores[i][offsets[i]];
729 else /* Offset to use is 0-6. */
730 scores[i][7] = sc2scores[i][offset_to_use];
732 /* do other registers */
733 for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
735 scores[i][(i+j)%POSNUM] = sc2scores[i][j];
745 /***************************************************************************/
747 /*** Score using the singles method on stock table and likelihood. **/
748 void NEWCOILSScore(int mode, Sequence sequence, double *maxscore,
749 double scores[MAXSEQLEN][POSNUM+1],
753 double stockscores[MAXSEQLEN][POSNUM];
754 char offsets[MAXSEQLEN];
757 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
759 /* score the thing */
761 sequence.seq, sequence.seqlen,
764 maxscore,ProlineFreeWindow, 1); /** Score using STOCK table. */
767 *maxscore= stocksc2likelihood(*maxscore);
770 for (i=0; i<sequence.seqlen; ++i) {
771 /* Do the most likely register. */
772 if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
773 scores[i][7] = stocksc2likelihood(stockscores[i][offsets[i]]);
774 else /* Use offset_to_use. */
775 scores[i][7] = stocksc2likelihood(stockscores[i][offset_to_use]);
776 /* do the other registers */
778 scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]);
785 /***********************************************************/
792 void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1],
793 int offset_to_use, int preprocessor)
798 fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");
801 for (i=0; i<sequence.seqlen; i++) {
802 for (j=0; j<POSNUM; j++)
803 if (!sequence.reg) scores[i][j] =0; /* No coils in input file. */
805 else if ( (offset_to_use ==7) || /* max offset is pos file reg */
806 ((offset_to_use == -1) && !preprocessor) )
807 if (j== sequence.reg[i])
811 else if ( (offset_to_use == -1) && preprocessor) /* Do all registers */
812 if (sequence.reg[i] != '.')
816 else /* just give likelihood 1 to offset_to_use. */
817 if ( (sequence.reg[i] != '.') && (j == (i+offset_to_use)%POSNUM))
821 if (!sequence.reg) scores[i][7]=0;
822 else if (sequence.reg[i] != '.')
835 /********************************************************/
837 /*** This averages score over coil, where a register shift does not start **/
838 /*** a new coil if start_at_reg_shift =0. ****/
839 /*** This defines whether in a coil by if likelihood is > bound. ***/
840 /** Value of 0 for avg_max means do average, 1 means do max. ***/
841 /*** For max score, find the maximum scoring residue in the CoilLike ***/
842 /*** and take the corresponding StructureScore. ***/
843 /*** To compute the max version, average over all residues in the coil **/
844 /*** that score the max coil likelihood. ***/
845 /*** BE VERY CAREFUL in making changes!!! i.e. need to count down on k **/
846 /*** since offset 7 needs to make use of previous offsets, so don't want **/
847 /*** them to have been changed to coil score yet. Also, we sometimes ***/
848 /*** pass in CoilLike as the StructureLike too, in which case have to be**/
849 /*** doubly careful. There is also ONE place where make a reference **/
850 /** to offset 7 when dealing with other registers, so had to save an **/
851 /** old version of offset 7 in case it changes (when StructureLike **/
852 /** is the same as CoilLike. **/
853 /*** ACTUALLY, all that should be moot now that we changed it so residue **/
854 /*** scores are always saved in tablenum +3, so pass those in and save **/
855 /*** coil scores in tablenum. ***/
856 void average_score_over_coils3(Sequence seq,
857 double StructureResLike[MAXSEQLEN][POSNUM+1],
858 double StructureCoilLike[MAXSEQLEN][POSNUM+1],
859 double CoilLike[MAXSEQLEN][POSNUM+1],
861 char offsets[MAXSEQLEN],
862 int start_at_reg_shift, double bound,
865 int i=0, j,k; /* k is the offset. */
866 int start_coil[POSNUM+1];
867 double total_coil_score[POSNUM+1];
868 double avg_coil_score[POSNUM+1];
869 double max_CoilLike[POSNUM+1];
870 double total_max_StructLike[POSNUM+1];
871 int number_max_residues[POSNUM+1];
872 double max_StructLike[POSNUM+1];
875 double total_coil_length[POSNUM +1]; /* A weighted length. */
876 double CoilLike_MaxOffset[MAXSEQLEN];
880 for (j=0; j<POSNUM+1; j++) {
882 total_coil_score[j]=0;
884 total_coil_length[j]=0;
886 total_max_StructLike[j]=0;
887 number_max_residues[j]=0;
890 while (i<=seq.seqlen) {
891 if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];
893 for (k=POSNUM; k>=0; k--) { /** Count down since k=7 needs to use **/
894 /** scores from other offset (i.e see **/
895 /** use of regist_to_use). **/
896 if (k!=7) regist = (i+k)%POSNUM;
899 /* Consider <= bound to be out of coil */
900 if ( ( CoilLike[i][regist] <= bound) || (i==seq.seqlen) ||
901 ( (i>0) && (offsets[i] != offsets[i-1]) && start_at_reg_shift &&
902 (CoilLike_MaxOffset[i-1]> bound) ) ) {
903 /** If not in a coil (value of -HUGE_VAL signals not in coil). */
904 if (start_coil[k] != -1) { /* Ended coil at i-1. */
905 avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
906 max_StructLike[k] = total_max_StructLike[k]/number_max_residues[k];
907 for (j=start_coil[k]; j<i; j++) {
910 StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
912 StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
915 StructureCoilLike[j][7]= avg_coil_score[k];
917 StructureCoilLike[j][7]= max_StructLike[k];
920 total_coil_score[k] =0; /***Resets for "not in coil" ***/
923 if (i != seq.seqlen) /* Consider <= bound to be out of coil */
924 if ( CoilLike[i][regist] <= bound) { /* Not in coil anymore. */
925 StructureCoilLike[i][regist]=0; /** Make coil score under bound */
928 /*** StructureCoilLike[i][regist]= StructureResLike[i][regist]; **/
929 /** If not in coil above bound make **/
930 /** the coil like the same as the **/
933 total_coil_length[k] =0;
936 else { /* Changing registers. */
938 if (CoilLike[i][regist] > max_CoilLike[k]) {
939 max_CoilLike[k] = CoilLike[i][regist];
940 total_max_StructLike[k] = StructureResLike[i][regist];
941 number_max_residues[k] =1;
943 else if (CoilLike[i][regist] == max_CoilLike[k]) {
944 total_max_StructLike[k] += StructureResLike[i][regist];
945 number_max_residues[k]++;
948 if (1) { /** Do weighted_avg **/
949 int reg_to_use=regist;
951 reg_to_use= (i+offsets[i])%POSNUM;
953 total_coil_length[k]= CoilLike[i][reg_to_use];
954 total_coil_score[k] = StructureResLike[i][regist]*
955 CoilLike[i][reg_to_use];
957 if (total_coil_length[k] ==0)
958 fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
963 } /** Ends if ending a coil. **/
965 else { /* If in a coil. */
966 if (CoilLike[i][regist] > max_CoilLike[k]) {
967 max_CoilLike[k] = CoilLike[i][regist];
968 total_max_StructLike[k] = StructureResLike[i][regist];
969 number_max_residues[k] =1;
971 else if (CoilLike[i][regist] == max_CoilLike[k]) {
972 total_max_StructLike[k] += StructureResLike[i][regist];
973 number_max_residues[k]++;
976 if (1) { /* Weighted average scores. */
977 int reg_to_use=regist;
979 reg_to_use= (i+offsets[i])%POSNUM;
981 total_coil_score[k] += StructureResLike[i][regist]*
982 CoilLike[i][reg_to_use];
983 total_coil_length[k] += CoilLike[i][reg_to_use];
985 if (start_coil[k] == -1)
993 /* Need the following hack to get the correct register in combo off method.*/
995 if ( (offset_to_use == 7) && (!start_at_reg_shift))
996 for (i=0; i<seq.seqlen; i++) {
997 for (j=0; j<POSNUM; j++)
998 StructureCoilLike[i][j]= -HUGE_VAL;
999 StructureCoilLike[i][(i+offsets[i]) %POSNUM] = StructureCoilLike[i][7];
1004 /** Need the following hack to get scores for max offset method. **/
1005 if (offset_to_use == -1)
1006 for (i=0; i<seq.seqlen; i++)
1007 StructureCoilLike[i][7]= StructureCoilLike[i][(i+offsets[i]) %POSNUM];