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 #define PIR_LIKES "/tmp/pir.likes"
19 #define PS_OF_LIKELIHOOD "~/likelihood.ps"
20 #define PAIR_PIR_TXT "/tmp/pir1.txt"
21 #define SINGLE_PIR_TXT "/tmp/pir2.txt"
23 double many_pprobs[MAX_TABLE_NUMBER][NUM_RES_TYPE][POSNUM],
24 many_pprobp[MAX_TABLE_NUMBER][NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM];
27 double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];
29 double (*pprobs)[POSNUM];
30 double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];
33 long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
34 long ptotals[MAX_TABLE_NUMBER][POSNUM];
35 long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
36 long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
39 #define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
47 ptr = (void*)malloc(n);
49 fprintf(stderr,"Memory shortage!!\n");
56 void switch_tables(int table)
58 extern int WINDOW, window_length[MAX_TABLE_NUMBER];
60 pprobs= many_pprobs[table];
61 pprobp= many_pprobp[table];
62 WINDOW = window_length[table];
67 /******** Used to subtract off a sequence from the table. ***/
69 void recalc_prob(Sequence sequence, int input_table, int mode)
71 int ProlineFreeWindow = mode & PROLINE_FREE_MODE;
73 /* Remove the sequence from the positive table before computing on it. */
74 add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
76 sequence.reg, sequence.seq, sequence.seqlen, -1);
78 if (!(mode & WEIGHTED_PROBS))
79 estimate_database_probs2(input_table, pfreqs[input_table],
80 ptotals[input_table], many_pprobs,
81 pfreqp[input_table], ptotalp[input_table],
82 many_pprobp, ProlineFreeWindow);
83 /** Add seq back into the frequence counts. **/
84 add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
86 sequence.reg, sequence.seq, sequence.seqlen, 1);
88 #ifdef AVG_WEIRD_PROBS
89 calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
91 calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
102 void PairCoilData (FILE *fgin, FILE **fpin,
103 int ProlineFree, int table_to_remove_from,
104 char *input_filename, int mode)
106 long gfreqs[AANUM]; long gtotals;
107 long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM];
109 extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
110 extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
111 extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
112 extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
113 extern double SCALE0;
119 readgenfreq(fgin, >otals, gfreqs, gtotalp, gfreqp, stderr);
120 calcgprobs(gfreqs, gtotals, gprobs);
121 calcgprobp(gfreqp, gtotalp, gprobp);
123 else setgprob_to_uniform(gprobs, gprobp);
125 #ifdef AVG_WEIRD_PROBS
126 calc_weird_gprob_avg(gprobp, gprobs);
128 calc_weird_gprob_sum(gprobp, gprobs);
132 for (i=0; i<NUMBER_TABLES; i++) { /* Go through all input tables. */
133 if (fpin[i] != NULL) { /* Calculate probabilities for that table. */
134 readposfreq(fpin[i],ptotals[i],pfreqs[i], ptotalp[i], pfreqp[i], stderr);
136 if (i == table_to_remove_from)
137 remove_all_input_seq_from_table(pfreqs[table_to_remove_from],
138 ptotals[table_to_remove_from],
139 pfreqp[table_to_remove_from],
140 ptotalp[table_to_remove_from], input_filename, -1);
142 estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
143 pfreqp[i], ptotalp[i], many_pprobp,
147 fprintf(stderr,"\nCould not make input table number %d.",i);
152 #ifdef AVG_WEIRD_PROBS /* make table. */
153 calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);
155 calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);
167 initialize_relative();
170 #define NEGCCRATIO 30
175 #define SQRT2PI 2.5066282746310002
179 /** To compute the paircoil likelihood line. */
182 void likelihood_from_txt(char *txt_filename, double *m, double *b1)
186 double percent_coil=PERCENT_COILED;
189 long n[2]={0,0}, n_halfgauss[2]={0,0};
190 double sum[2]={0,0}, sum2[2]={0,0};
191 double a[2],b[2]; /* ax+b is the likelihood line approximation */
192 double mean[2], sd[2], down[2];
193 long _counts[2][-4*MINSCORE +1];
194 double _means[2][-4*MINSCORE+1];
197 int i, maxplaces[2]= {2*MINSCORE, 2*MINSCORE}, minplaces[2]={-2*MINSCORE,-2*MINSCORE};
198 int minplace, maxplace;
199 int number_strange=0, number_high=0;
200 FILE *devnull, *outfile, *outpsfile;
203 counts[0] = &(_counts[0][-2*MINSCORE]);
204 counts[1] = &(_counts[1][-2*MINSCORE]);
205 means[0] = &(_means[0][-2*MINSCORE]);
206 means[1] = &(_means[1][-2*MINSCORE]);
209 mode2 |= RATIO_TRAP_METHOD;
210 mode2 |= PERCENT_COIL_MODE;
214 for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
220 txt_file = sopen (txt_filename, "r");
221 while (fread(&score1,sizeof(double),1,txt_file)) {
222 if (score1<MINSCORE || score1>-(MINSCORE)) {
224 fprintf(stderr,"Strange score: %lf\n",score1);
229 place = ( (int)((score1-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
230 if (place > -2*MINSCORE) {
231 fprintf(stderr, "Out of bounds: %d.\n", place);
232 fprintf(stderr, "Score: %lf\n",score1);
234 counts[0][-2*MINSCORE]++;
235 means[0][-2*MINSCORE]++;
237 maxplaces[0]=-2*MINSCORE;
239 else if (place < 2*MINSCORE) {
240 fprintf(stderr, "Out of bounds: %d.\n", place);
241 fprintf(stderr, "Score: %lf\n",score1);
242 counts[0][2*MINSCORE]++;
243 means[0][2*MINSCORE]++;
244 minplaces[0]=2*MINSCORE;
250 means[0][place] += score1;
251 if (place>maxplaces[0]) maxplaces[0]=place;
252 if (place<minplaces[0]) minplaces[0]=place;
257 fprintf(stderr,"\n%d out of bounds scores > %d were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE));
258 fprintf(stderr,"\n%d strange scores < %d were scored as %d\n",number_strange-number_high,MINSCORE, MINSCORE);
261 /** Compute the means of the scores in each box. **/
262 for (i= minplaces[0]; i<=maxplaces[0]; ++i)
263 if (counts[0][i] != 0)
264 means[0][i] /= counts[0][i]; /** Each bucket gets average score*/
267 /*** We need a mean***/
269 mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
273 /*** Calculate gaussian fit to scores below the mean. **/
274 txt_file=sopen(txt_filename, "r");
275 while (fread(&score1,sizeof(double),1,txt_file))
276 if ( (score1 > MINSCORE) && (score1 < mean[0])) {
278 sum2[0] += (score1-mean[0])*(score1-mean[0]);
280 sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
281 down[0]= (double) 2*n_halfgauss[0]/n[0];
283 fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
286 minplace=minplaces[0]; maxplace=maxplaces[0];
288 devnull=sopen("/dev/null","r+");
289 outfile=sopen(PIR_LIKES,"w");
290 outpsfile=sopen(PS_OF_LIKELIHOOD,"w");
291 /* outfile=sopen("~/pir.out","w"); */
292 likelihood(n,mode2, filenum, 0.0 , minplace, maxplace, means,
293 counts, mean, sd, down, outfile,outpsfile,devnull,devnull,
298 char clean_command[200];
300 strcpy(clean_command, "rm -f ");
302 /*** The following removes the pir.txt output file. */
303 system(strcat(clean_command,txt_filename));
305 } /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
306 /** WRITE TO SAME FILESNAMES. **/
308 char clean_command[200];
310 strcpy(clean_command, "rm -f ");
312 system(strcat(clean_command,PIR_LIKES));
323 /* Look in file .likelihood for likelihood line. If not there, compute */
325 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)
335 FILE *txt_file1, *txt_file2;
336 double score, scs[MAXSEQLEN][POSNUM];
338 char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
339 char title[500], code[500];
342 int found_like_lib=0; int found_single=0;
346 sequence.title=title;
350 for (i=0; i<functnum; i++)
351 libnumber |= (int) pow(2,lib[i]);
352 fprintf(stderr,"\nlibnumber= %d",libnumber);
354 /* What directory to look for paircoil defauts in. Set environment */
355 /* variable LIKELINE_DIR to be a path to this file. */
356 for (likefile=0; likefile <number_tables; likefile++) {
360 fprintf(stderr,"\nLikelihood %d=%s",likefile,likelihoods[likefile]);
361 if (likelihoods[likefile][0] == ',') {
362 if (tmp=getenv("LIKELINE_DIR"))
363 strcpy(likelihoods[likefile],tmp);
365 strcpy(likelihoods[likefile],getenv("HOME"));
366 strcat(likelihoods[likefile],"/.likelihood");
372 if (linefile=sopen(likelihoods[likefile],"r")) {
373 while ( (fgets(buf,80,linefile)) && !(found_single && found_like_lib)) {
375 sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
378 else if (!found_like_lib) {
379 sscanf(buf, "%d %lf %lf",¤t_lib, &m[likefile], &b[likefile]);
380 if (current_lib == libnumber) found_like_lib=1;;
385 if (!(found_like_lib && found_single)) {
387 fprintf(stderr,"\nSingles likelihood function %d unknown.",likefile);
389 if (current_lib != libnumber ) /* Wasn't in file so should create it
390 by running on pir. Add this later.*/
391 fprintf(stderr, "\nLikelihood function %d unknown...", likefile);
394 if (pir_name[0]==',') {
395 fprintf(stderr,"\nUnable to continue. Please input the location of the pir in the .paircoil file.");
399 fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
401 pir = sopen(pir_name, "r");
403 if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
404 if (!found_single) txt_file2 = sopen(SINGLE_PIR_TXT,"w");
406 for (cnt=1;getseq2(pir,&sequence); cnt++)
407 if (sequence.seqlen >= MINWINDOW) {
408 seqlen= sequence.seqlen;
415 if (!found_like_lib) {
416 switch_tables(likefile);
417 scseqadj(lib,seq,seqlen,scs,offsets, &score,
419 txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
422 switch_tables(likefile);
424 scorestock(seq,seqlen,scs,offsets, &score,
425 mode & PROLINE_FREE_MODE, 0);
426 /* Score stock = singles method using our table. */
427 txt_output(sequence,0,txt_file2, scs, NULL,offsets,0);
430 } /* Ends if sequence > MINWINDOW, and hence the for */
431 /* loop that reads from pir. */
435 if (!found_like_lib) { /* Now compute the lib likelihood line. */
438 likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
440 linefile=sopen(likelihoods[likefile],"a");
441 fprintf(linefile,"%d %lf %lf\n",libnumber,m[likefile], b[likefile]);
443 fprintf(linefile,"\t (lib =");
444 /* Now print human readable lib number. */
445 for (i=0; i<functnum; i++)
446 fprintf(linefile," %d",lib[i]);
447 fprintf(linefile," )\n");
452 if (!found_single) { /* Now compute the singles method like line */
455 likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
456 &b_single[likefile]);
458 linefile=sopen(likelihoods[likefile],"a");
459 fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile],
462 fprintf(linefile,"\t (lib =");
463 /* Now print human readable lib number. */
464 fprintf(linefile," Singles Method");
465 fprintf(linefile," )\n");
469 } /* Ends else from if (pir_name == ',') */
471 } /* Ends section for file not having lib or singles like line */
473 } /* Ends section where search likelihoods[likefile] file. */
475 fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
478 } /* Ends count through tables=linefiles. */
483 double stocksc2likelihood (double sc)
488 gg = (sc - NEGMEAN)/NEGSD;
489 gcc = (sc - CCMEAN)/CCSD;
490 gg = exp((-gg * gg) / 2) / NEGSD / SQRT2PI;
491 gcc = exp((-gcc * gcc) / 2) / CCSD / SQRT2PI;
492 return(gcc / (NEGCCRATIO * gg + gcc));
497 /* Likelihood line is mx +b */
498 double sc2likelihood (double sc, double m, double b)
502 likelihood = m*sc + b;
504 if (likelihood>1) return(1);
505 else if (likelihood<0) return(0);
506 else return(likelihood);
510 static int current_seq_number=0;
512 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)
515 double sc2scores[MAXSEQLEN][POSNUM];
517 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
519 double frac_coiled_this_table;
520 extern double init_class_prob[MAX_CLASS_NUMBER];
522 frac_coiled_this_table = init_class_prob[table];
525 switch_tables(table);
527 /* score the sequence using sc2seq */
529 sequence.seq, sequence.seqlen,
536 if (mode & USE_LIKE_LINE)
537 *maxscore= sc2likelihood(*maxscore,m,b);
539 *maxscore= frac_coiled_this_table * exp(*maxscore);
542 for (i=0; i<sequence.seqlen; ++i) {
543 /* Do the most probabable register if offset_to_use is -1 or 7. */
544 if ((offset_to_use == -1) || (offset_to_use ==7))
545 if (raw_score) /*** Return rawscore. **/
546 scores[i][7] = sc2scores[i][offsets[i]];
548 if (mode & USE_LIKE_LINE)
549 scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
551 scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
552 if (scores[i][7] >1) scores[i][7]=1;
555 else /* Offset to use is 0-6. */
557 scores[i][7] = sc2scores[i][offset_to_use];
559 if (mode & USE_LIKE_LINE)
560 scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
562 scores[i][7] = frac_coiled_this_table *
563 exp(sc2scores[i][offset_to_use]);
564 if (scores[i][7] >1) scores[i][7]=1;
567 /* do other registers */
568 for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
571 scores[i][(i+j)%POSNUM] = sc2scores[i][j];
573 if (mode & USE_LIKE_LINE)
574 scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
576 scores[i][(i+j)%POSNUM] = frac_coiled_this_table *
577 exp(sc2scores[i][j]);
578 if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
589 /***************************************************************************/
591 /** Score using the singles method on our tables. **/
592 double *SingleCoilScore(int mode, Sequence sequence, double m, double b,
593 double *maxscore, int table,
594 double scores[MAXSEQLEN][POSNUM+1], int offset_to_use,
598 double singlescores[MAXSEQLEN][POSNUM];
599 char offsets[MAXSEQLEN];
601 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
603 double frac_coiled_this_table;
606 switch_tables(table);
608 /* score the thing */
610 sequence.seq, sequence.seqlen,
613 maxscore,ProlineFreeWindow, 0); /** Score using our table. */
617 if (mode & USE_LIKE_LINE)
618 *maxscore= sc2likelihood(*maxscore,m,b);
620 /* fprintf(stderr,"\nmaxscore = %lf", *maxscore);
622 *maxscore= frac_coiled_this_table * exp(*maxscore);
623 /* fprintf(stderr," maxlike = %lf", *maxscore);
627 for (i=0; i<sequence.seqlen; ++i) {
628 /* Do the most likely register. */
629 if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
631 if (mode & USE_LIKE_LINE)
632 scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b);
634 scores[i][7] = frac_coiled_this_table*
635 exp(singlescores[i][offsets[i]]);
636 if (scores[i][7] >1) scores[i][7]=1;
638 else scores[i][7] = singlescores[i][offsets[i]];
640 else /* Use offset_to_use. */
642 if (mode & USE_LIKE_LINE)
643 scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
645 scores[i][7] = frac_coiled_this_table*
646 exp(singlescores[i][offsets[i]]);
647 if (scores[i][7] >1) scores[i][7]=1;
649 else scores[i][7] = singlescores[i][offset_to_use];
651 /* do the other registers */
654 if (mode & USE_LIKE_LINE)
655 scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b);
657 scores[i][(i+j)%POSNUM] = frac_coiled_this_table*
658 exp(singlescores[i][j]);
659 if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
661 else scores[i][(i+j)%POSNUM] = singlescores[i][j]; /*Raw_score*/
667 /***************************************************************************/
671 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)
674 double sc2scores[MAXSEQLEN][POSNUM];
676 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
679 switch_tables(table);
681 /* score the sequence using sc2seq */
683 sequence.seq, sequence.seqlen,
688 /*** Do a SingleCoil type score. ****/
690 scorestock(sequence.seq, sequence.seqlen,
692 maxscore,ProlineFreeWindow, 0); **/ /** Score using our table. */
695 for (i=0; i<sequence.seqlen; ++i) {
696 /* Do the most probabable register if offset_to_use is -1 or 7. */
697 if ((offset_to_use == -1) || (offset_to_use ==7))
698 scores[i][7] = sc2scores[i][offsets[i]];
699 else /* Offset to use is 0-6. */
700 scores[i][7] = sc2scores[i][offset_to_use];
702 /* do other registers */
703 for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
705 scores[i][(i+j)%POSNUM] = sc2scores[i][j];
715 /***************************************************************************/
717 /*** Score using the singles method on stock table and likelihood. **/
718 void NEWCOILSScore(int mode, Sequence sequence, double *maxscore,
719 double scores[MAXSEQLEN][POSNUM+1],
723 double stockscores[MAXSEQLEN][POSNUM];
724 char offsets[MAXSEQLEN];
727 int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
729 /* score the thing */
731 sequence.seq, sequence.seqlen,
734 maxscore,ProlineFreeWindow, 1); /** Score using STOCK table. */
737 *maxscore= stocksc2likelihood(*maxscore);
740 for (i=0; i<sequence.seqlen; ++i) {
741 /* Do the most likely register. */
742 if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
743 scores[i][7] = stocksc2likelihood(stockscores[i][offsets[i]]);
744 else /* Use offset_to_use. */
745 scores[i][7] = stocksc2likelihood(stockscores[i][offset_to_use]);
746 /* do the other registers */
748 scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]);
755 /***********************************************************/
762 void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1],
763 int offset_to_use, int preprocessor)
768 fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");
771 for (i=0; i<sequence.seqlen; i++) {
772 for (j=0; j<POSNUM; j++)
773 if (!sequence.reg) scores[i][j] =0; /* No coils in input file. */
775 else if ( (offset_to_use ==7) || /* max offset is pos file reg */
776 ((offset_to_use == -1) && !preprocessor) )
777 if (j== sequence.reg[i])
781 else if ( (offset_to_use == -1) && preprocessor) /* Do all registers */
782 if (sequence.reg[i] != '.')
786 else /* just give likelihood 1 to offset_to_use. */
787 if ( (sequence.reg[i] != '.') && (j == (i+offset_to_use)%POSNUM))
791 if (!sequence.reg) scores[i][7]=0;
792 else if (sequence.reg[i] != '.')
805 /********************************************************/
807 /*** This averages score over coil, where a register shift does not start **/
808 /*** a new coil if start_at_reg_shift =0. ****/
809 /*** This defines whether in a coil by if likelihood is > bound. ***/
810 /** Value of 0 for avg_max means do average, 1 means do max. ***/
811 /*** For max score, find the maximum scoring residue in the CoilLike ***/
812 /*** and take the corresponding StructureScore. ***/
813 /*** To compute the max version, average over all residues in the coil **/
814 /*** that score the max coil likelihood. ***/
815 /*** BE VERY CAREFUL in making changes!!! i.e. need to count down on k **/
816 /*** since offset 7 needs to make use of previous offsets, so don't want **/
817 /*** them to have been changed to coil score yet. Also, we sometimes ***/
818 /*** pass in CoilLike as the StructureLike too, in which case have to be**/
819 /*** doubly careful. There is also ONE place where make a reference **/
820 /** to offset 7 when dealing with other registers, so had to save an **/
821 /** old version of offset 7 in case it changes (when StructureLike **/
822 /** is the same as CoilLike. **/
823 /*** ACTUALLY, all that should be moot now that we changed it so residue **/
824 /*** scores are always saved in tablenum +3, so pass those in and save **/
825 /*** coil scores in tablenum. ***/
826 void average_score_over_coils3(Sequence seq,
827 double StructureResLike[MAXSEQLEN][POSNUM+1],
828 double StructureCoilLike[MAXSEQLEN][POSNUM+1],
829 double CoilLike[MAXSEQLEN][POSNUM+1],
831 char offsets[MAXSEQLEN],
832 int start_at_reg_shift, double bound,
835 int i=0, j,k; /* k is the offset. */
836 int start_coil[POSNUM+1];
837 double total_coil_score[POSNUM+1];
838 double avg_coil_score[POSNUM+1];
839 double max_CoilLike[POSNUM+1];
840 double total_max_StructLike[POSNUM+1];
841 int number_max_residues[POSNUM+1];
842 double max_StructLike[POSNUM+1];
845 double total_coil_length[POSNUM +1]; /* A weighted length. */
846 double CoilLike_MaxOffset[MAXSEQLEN];
850 for (j=0; j<POSNUM+1; j++) {
852 total_coil_score[j]=0;
854 total_coil_length[j]=0;
856 total_max_StructLike[j]=0;
857 number_max_residues[j]=0;
860 while (i<=seq.seqlen) {
861 if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];
863 for (k=POSNUM; k>=0; k--) { /** Count down since k=7 needs to use **/
864 /** scores from other offset (i.e see **/
865 /** use of regist_to_use). **/
866 if (k!=7) regist = (i+k)%POSNUM;
869 /* Consider <= bound to be out of coil */
870 if ( ( CoilLike[i][regist] <= bound) || (i==seq.seqlen) ||
871 ( (i>0) && (offsets[i] != offsets[i-1]) && start_at_reg_shift &&
872 (CoilLike_MaxOffset[i-1]> bound) ) ) {
873 /** If not in a coil (value of -HUGE_VAL signals not in coil). */
874 if (start_coil[k] != -1) { /* Ended coil at i-1. */
875 avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
876 max_StructLike[k] = total_max_StructLike[k]/number_max_residues[k];
877 for (j=start_coil[k]; j<i; j++) {
880 StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
882 StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
885 StructureCoilLike[j][7]= avg_coil_score[k];
887 StructureCoilLike[j][7]= max_StructLike[k];
890 total_coil_score[k] =0; /***Resets for "not in coil" ***/
893 if (i != seq.seqlen) /* Consider <= bound to be out of coil */
894 if ( CoilLike[i][regist] <= bound) { /* Not in coil anymore. */
895 StructureCoilLike[i][regist]=0; /** Make coil score under bound */
898 /*** StructureCoilLike[i][regist]= StructureResLike[i][regist]; **/
899 /** If not in coil above bound make **/
900 /** the coil like the same as the **/
903 total_coil_length[k] =0;
906 else { /* Changing registers. */
908 if (CoilLike[i][regist] > max_CoilLike[k]) {
909 max_CoilLike[k] = CoilLike[i][regist];
910 total_max_StructLike[k] = StructureResLike[i][regist];
911 number_max_residues[k] =1;
913 else if (CoilLike[i][regist] == max_CoilLike[k]) {
914 total_max_StructLike[k] += StructureResLike[i][regist];
915 number_max_residues[k]++;
918 if (1) { /** Do weighted_avg **/
919 int reg_to_use=regist;
921 reg_to_use= (i+offsets[i])%POSNUM;
923 total_coil_length[k]= CoilLike[i][reg_to_use];
924 total_coil_score[k] = StructureResLike[i][regist]*
925 CoilLike[i][reg_to_use];
927 if (total_coil_length[k] ==0)
928 fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
933 } /** Ends if ending a coil. **/
935 else { /* If in a coil. */
936 if (CoilLike[i][regist] > max_CoilLike[k]) {
937 max_CoilLike[k] = CoilLike[i][regist];
938 total_max_StructLike[k] = StructureResLike[i][regist];
939 number_max_residues[k] =1;
941 else if (CoilLike[i][regist] == max_CoilLike[k]) {
942 total_max_StructLike[k] += StructureResLike[i][regist];
943 number_max_residues[k]++;
946 if (1) { /* Weighted average scores. */
947 int reg_to_use=regist;
949 reg_to_use= (i+offsets[i])%POSNUM;
951 total_coil_score[k] += StructureResLike[i][regist]*
952 CoilLike[i][reg_to_use];
953 total_coil_length[k] += CoilLike[i][reg_to_use];
955 if (start_coil[k] == -1)
963 /* Need the following hack to get the correct register in combo off method.*/
965 if ( (offset_to_use == 7) && (!start_at_reg_shift))
966 for (i=0; i<seq.seqlen; i++) {
967 for (j=0; j<POSNUM; j++)
968 StructureCoilLike[i][j]= -HUGE_VAL;
969 StructureCoilLike[i][(i+offsets[i]) %POSNUM] = StructureCoilLike[i][7];
974 /** Need the following hack to get scores for max offset method. **/
975 if (offset_to_use == -1)
976 for (i=0; i<seq.seqlen; i++)
977 StructureCoilLike[i][7]= StructureCoilLike[i][(i+offsets[i]) %POSNUM];