sources/multicoil/sc2seq_interface.c

   1 /* Interface to sc2seq created by David Wilson and modified by Ethan Wolf */
   2 /* Other interfaces by Ethan Wolf 1995. */
   3
   4 #include <stdio.h>
   5 #include <stdlib.h>
   6 #include <math.h>
   7 #include <string.h>
   8 #include "scio.h"
   9 #include "switches.h"
  10 #include "likelihood.h"
  11 #include "sc.h"
  12 #include "scconst.h"
  13 #include "options.h"
  14 #include "compute_like.h"
  15 #include "interface.h"
  16 #include "stats.h"
  17
  18 #include "debug_printout_flag.h"
  19
  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"
  24
  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];
  27
  28
  29 double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];
  30
  31 double (*pprobs)[POSNUM];
  32 double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];
  33
  34
  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];
  39
  40 double getdlog();
  41
  42 #define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
  43
  44
  45
  46 void *Malloc (n)
  47      int n;
  48 {
  49   void *ptr;
  50   ptr = (void*)malloc(n);
  51   if (ptr) return ptr;
  52   fprintf(stderr,"Memory shortage!!\n");
  53   exit(-1);
  54 }
  55
  56
  57
  58
  59 void switch_tables(int table)
  60 {
  61   extern int WINDOW, window_length[MAX_TABLE_NUMBER];
  62
  63   pprobs= many_pprobs[table];
  64   pprobp= many_pprobp[table];
  65   WINDOW = window_length[table];
  66 }
  67
  68
  69
  70 /******** Used to subtract off a sequence from the table.  ***/
  71
  72 void recalc_prob(Sequence sequence, int input_table, int mode, long pfreqp2[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
  73 {
  74   int ProlineFreeWindow = mode & PROLINE_FREE_MODE;
  75
  76   /* Remove the sequence from the positive table before computing on it. */
  77   add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
  78              ptotalp[input_table],
  79              sequence.reg, sequence.seq, sequence.seqlen, -1);
  80
  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],
  88              ptotalp[input_table],
  89              sequence.reg, sequence.seq, sequence.seqlen, 1);
  90
  91   #ifdef AVG_WEIRD_PROBS
  92     calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
  93   #else
  94     calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
  95   #endif
  96
  97 }
  98
  99
 100
 101 void normalize_gen(double gprobs[NUM_RES_TYPE],
 102                    double gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM])
 103 {
 104   int res1, res2, dist;
 105
 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]))) );
 115       }
 116   }
 117
 118 }
 119
 120
 121
 122
 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])
 127 {
 128   long gfreqs[AANUM]; long gtotals;
 129   long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM];
 130
 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;
 136
 137
 138   int i=0;
 139
 140   if (fgin) {
 141     readgenfreq(fgin, &gtotals, gfreqs, gtotalp, gfreqp, stderr);
 142     calcgprobs(gfreqs, gtotals, gprobs);
 143     calcgprobp(gfreqp, gtotalp, gprobp);
 144 #ifdef NoGenProline
 145     normalize_gen(gprobs, gprobp);
 146 #endif
 147   }
 148   else { setgprob_to_uniform(gprobs, gprobp); }
 149
 150   #ifdef AVG_WEIRD_PROBS
 151     calc_weird_gprob_avg(gprobp, gprobs);
 152   #else
 153     calc_weird_gprob_sum(gprobp, gprobs);
 154   #endif
 155
 156
 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);
 160
 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);
 166
 167       estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
 168                                 pfreqp[i], ptotalp[i], many_pprobp,
 169                                 ProlineFree);
 170     }
 171     else {
 172       fprintf(stderr,"\nCould not make input table number %d.",i);
 173       exit(-1);
 174     }
 175
 176
 177 #ifdef AVG_WEIRD_PROBS        /* make table.     */
 178     calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);
 179 #else
 180     calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);
 181 #endif
 182
 183   }
 184
 185   switch_tables(0);
 186
 187 }
 188
 189
 190 NEWCOILSinit ()
 191 {
 192   initialize_relative();
 193 }
 194
 195 #define NEGCCRATIO 30
 196 #define NEGMEAN 0.77
 197 #define NEGSD 0.20
 198 #define CCMEAN 1.63
 199 #define CCSD 0.24
 200 #define SQRT2PI 2.5066282746310002
 201
 202
 203
 204 /** To compute the paircoil likelihood line. */
 205
 206
 207 void likelihood_from_txt(char *txt_filename, double *m, double *b1)
 208 {
 209   FILE *txt_file;
 210   int mode2=0;
 211   double percent_coil=PERCENT_COILED;
 212   double score1;
 213   int filenum=1;
 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];
 220   long *counts[2];
 221   double  *means[2];
 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;
 226   int place;
 227
 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]);
 232
 233
 234   mode2 |= RATIO_TRAP_METHOD;
 235   mode2 |= PERCENT_COIL_MODE;
 236
 237
 238
 239    for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
 240      counts[0][i] = 0;
 241      means[0][i]=0;
 242     }
 243
 244
 245   txt_file = sopen (txt_filename, "r");
 246   while (fread(&score1,sizeof(double),1,txt_file)) {
 247     if (score1<MINSCORE || score1>-(MINSCORE)) {
 248 /*
 249   fprintf(stderr,"Strange score: %lf\n",score1);
 250   exit(-1);
 251 */
 252       number_strange++;
 253     }
 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);
 258       /* exit(-1); */
 259       counts[0][-2*MINSCORE]++;
 260       means[0][-2*MINSCORE]++;
 261       number_high++;
 262       maxplaces[0]=-2*MINSCORE;
 263     }
 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;
 270         /* exit(-1); */
 271     }
 272
 273     else {
 274       counts[0][place]++;
 275       means[0][place] += score1;
 276       if (place>maxplaces[0]) maxplaces[0]=place;
 277       if (place<minplaces[0]) minplaces[0]=place;
 278       ++n[0];
 279     }
 280   }
 281
 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);
 284    sclose(txt_file);
 285
 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*/
 290
 291
 292    /*** We need a mean***/
 293
 294    mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
 295                          maxplaces[0],n[0]);
 296
 297
 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])) {
 302       ++n_halfgauss[0];
 303       sum2[0] += (score1-mean[0])*(score1-mean[0]);
 304     }
 305   sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
 306   down[0]= (double) 2*n_halfgauss[0]/n[0];
 307
 308   fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
 309   sclose(txt_file);
 310
 311   minplace=minplaces[0];  maxplace=maxplaces[0];
 312
 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,
 319              a,b,.20,.80,.03);
 320   sclose(outfile);
 321   sclose(outpsfile);
 322   {
 323     char clean_command[200];
 324
 325     strcpy(clean_command, "rm -f ");
 326
 327 /*** The following removes the pir.txt output file. */
 328     system(strcat(clean_command,txt_filename));
 329
 330     }                            /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
 331                                  /** WRITE TO SAME FILESNAMES.             **/
 332   {
 333     char clean_command[200];
 334
 335     strcpy(clean_command, "rm -f ");
 336
 337     system(strcat(clean_command,PIR_LIKES));
 338
 339     }
 340
 341   *m=a[0]; *b1=b[0];
 342 }
 343
 344
 345
 346
 347
 348 /*  Look in file .likelihood for likelihood line.   If not there, compute */
 349 /*  it from pir.                                               */
 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)
 351 {
 352  char *tmp;
 353  int libnumber=0;
 354  int  i, cnt=0;
 355  int current_lib;
 356  char buf[80];
 357
 358  FILE *pir;
 359  FILE *linefile;
 360  FILE *txt_file1, *txt_file2;
 361  double score, scs[MAXSEQLEN][POSNUM];
 362  Sequence sequence;
 363  char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
 364  char title[500], code[500];
 365  int seqlen;
 366  int likefile;
 367  int found_like_lib=0;  int found_single=0;
 368
 369  sequence.seq=seq;
 370  sequence.reg=reg;
 371  sequence.title=title;
 372  sequence.code=code;
 373
 374
 375  for (i=0; i<functnum; i++)
 376    libnumber |= (int) pow(2,lib[i]);
 377 #ifdef DEBUG_PRINTOUT
 378  fprintf(stderr,"\nlibnumber= %d",libnumber);
 379 #endif
 380
 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++) {
 384    current_lib = -1;
 385    found_like_lib =0;
 386    found_single =0;
 387 #ifdef DEBUG_PRINTOUT
 388    fprintf(stderr,"\nLikelihood %d=%s",likefile,likelihoods[likefile]);
 389 #endif
 390    if (likelihoods[likefile][0] == ',') { /* Don't bother loading likefile */
 391 #ifdef DEBUG_PRINTOUT
 392      fprintf(stderr,
 393              "\n Not doing likelihood conversion for paircoil with table %d\n",
 394              likefile);
 395 #endif
 396    }
 397    else {
 398      if (linefile=sopen(likelihoods[likefile],"r")) {
 399        while ((fgets(buf,80,linefile)) && !(found_single && found_like_lib))  {
 400          if (buf[0] == '0') {
 401            sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
 402            found_single=1;
 403          }
 404          else if (!found_like_lib) {
 405            sscanf(buf, "%d %lf %lf",&current_lib, &m[likefile], &b[likefile]);
 406            if (current_lib == libnumber)  found_like_lib=1;;
 407          }
 408        }
 409        sclose(linefile);
 410
 411        if (!(found_like_lib && found_single)) {
 412          if (!found_single)
 413            fprintf(stderr,"\nSingles likelihood function %d unknown.",
 414                    likefile);
 415
 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);
 419
 420
 421          if (pir_name[0]==',') {
 422            fprintf(stderr,"\nUnable to continue.  Please input the location of the pir in the .paircoil file.");
 423            exit(-1);
 424          }
 425          else {
 426            fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
 427
 428            pir = sopen(pir_name, "r");
 429
 430            if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
 431            if (!found_single)   txt_file2 = sopen(SINGLE_PIR_TXT,"w");
 432
 433            for (cnt=1;getseq2(pir,&sequence); cnt++)
 434              if (sequence.seqlen >= MINWINDOW) {
 435                seqlen= sequence.seqlen;
 436                if (!(cnt %10)) {
 437                  putc('.',stderr);
 438                  fflush(stderr);
 439                }
 440
 441                WINDOW=PAIRWINDOW;
 442                if (!found_like_lib) {
 443                  switch_tables(likefile);
 444                  scseqadj(lib,seq,seqlen,scs,offsets, &score,
 445                           mode,0);
 446                  txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
 447                }
 448                if (!found_single) {
 449                  switch_tables(likefile);
 450
 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);
 455                }
 456
 457              }  /* Ends if sequence > MINWINDOW, and hence the for */
 458                  /* loop that reads from pir. */
 459            sclose(pir);
 460
 461
 462            if (!found_like_lib) {    /* Now compute the lib likelihood line. */
 463              sclose(txt_file1);
 464
 465              likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
 466
 467              linefile=sopen(likelihoods[likefile],"a");
 468              fprintf(linefile,"%d %lf %lf\n",
 469                      libnumber,m[likefile], b[likefile]);
 470
 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");
 476
 477              sclose(linefile);
 478            }
 479
 480            if (!found_single) { /*  Now compute the singles method like line */
 481              sclose(txt_file2);
 482
 483              likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
 484                                &b_single[likefile]);
 485
 486              linefile=sopen(likelihoods[likefile],"a");
 487              fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile],
 488                      b_single[likefile]);
 489
 490              fprintf(linefile,"\t (lib =");
 491              /* Now print human readable lib number. */
 492              fprintf(linefile," Singles Method");
 493              fprintf(linefile," )\n");
 494
 495              sclose(linefile);
 496            }
 497          }      /* Ends else from if (pir_name == ',') */
 498
 499        }       /*  Ends section for file not having lib or singles like line */
 500
 501      }         /*  Ends section where search likelihoods[likefile] file. */
 502      else   {
 503        fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
 504        exit(-1);
 505      }
 506    }  /** Ends else the likelihood file exists to load a likelihood line... **/
 507
 508  }          /*  Ends count through tables=linefiles.  */
 509 }
 510
 511
 512
 513 double stocksc2likelihood (double sc)
 514 {
 515   double gg, gcc;
 516
 517   sc = exp(sc/28);
 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));
 523 }
 524
 525
 526
 527 /*  Likelihood line is mx +b */
 528 double sc2likelihood (double sc, double m, double b)
 529 {
 530   double likelihood;
 531
 532   likelihood = m*sc + b;
 533
 534   if (likelihood>1) return(1);
 535   else if (likelihood<0) return(0);
 536   else return(likelihood);
 537 }
 538
 539
 540 static int current_seq_number=0;
 541
 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)
 543 {
 544   int i,j;
 545   double sc2scores[MAXSEQLEN][POSNUM];
 546
 547   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 548
 549   double  frac_coiled_this_table;
 550   extern double init_class_prob[MAX_CLASS_NUMBER];
 551
 552   frac_coiled_this_table = init_class_prob[table];
 553
 554
 555   switch_tables(table);
 556
 557   /* score the sequence using sc2seq */
 558   scseqadj(lib,
 559             sequence.seq, sequence.seqlen,
 560             sc2scores,
 561             offsets,
 562             maxscore, mode, 0);
 563
 564
 565   if (!raw_score)
 566     if (mode & USE_LIKE_LINE)
 567       *maxscore= sc2likelihood(*maxscore,m,b);
 568     else {
 569       *maxscore= frac_coiled_this_table * exp(*maxscore);
 570     }
 571
 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]];
 577       else
 578         if (mode & USE_LIKE_LINE)
 579           scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
 580         else   {
 581           scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
 582           if (scores[i][7] >1) scores[i][7]=1;
 583         }
 584
 585     else    /*  Offset to use is 0-6. */
 586       if (raw_score)
 587         scores[i][7] = sc2scores[i][offset_to_use];
 588       else
 589         if (mode & USE_LIKE_LINE)
 590           scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
 591         else {
 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;
 595         }
 596
 597     /* do other registers */
 598     for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 599                            i+j%POSNUM  **/
 600       if (raw_score)
 601         scores[i][(i+j)%POSNUM] = sc2scores[i][j];
 602       else
 603         if (mode & USE_LIKE_LINE)
 604           scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
 605         else  {
 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;
 609         }
 610     }
 611   }
 612
 613 }
 614
 615
 616
 617
 618
 619 /***************************************************************************/
 620
 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,
 625                         int raw_score)
 626 {
 627   int i,j;
 628   double singlescores[MAXSEQLEN][POSNUM];
 629   char offsets[MAXSEQLEN];
 630
 631   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 632
 633   double  frac_coiled_this_table;
 634
 635
 636   switch_tables(table);
 637
 638   /* score the thing */
 639   scorestock(
 640             sequence.seq, sequence.seqlen,
 641             singlescores,
 642             offsets,
 643             maxscore,ProlineFreeWindow, 0);  /** Score using our table. */
 644
 645
 646   if (!raw_score)
 647     if (mode & USE_LIKE_LINE)
 648       *maxscore= sc2likelihood(*maxscore,m,b);
 649     else {
 650       /*    fprintf(stderr,"\nmaxscore = %lf", *maxscore);
 651        */
 652       *maxscore= frac_coiled_this_table * exp(*maxscore);
 653       /*    fprintf(stderr,"  maxlike  = %lf",  *maxscore);
 654        */
 655     }
 656
 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. */
 660       if (!raw_score)
 661         if (mode & USE_LIKE_LINE)
 662           scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b);
 663         else {
 664           scores[i][7] = frac_coiled_this_table*
 665                       exp(singlescores[i][offsets[i]]);
 666           if (scores[i][7] >1) scores[i][7]=1;
 667         }
 668       else scores[i][7] = singlescores[i][offsets[i]];
 669
 670     else                  /*  Use offset_to_use. */
 671       if (!raw_score)
 672         if (mode & USE_LIKE_LINE)
 673           scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
 674         else {
 675           scores[i][7] = frac_coiled_this_table*
 676                      exp(singlescores[i][offsets[i]]);
 677           if (scores[i][7] >1) scores[i][7]=1;
 678         }
 679       else scores[i][7] = singlescores[i][offset_to_use];
 680
 681     /* do the other registers */
 682     for (j=0; j<7; ++j)
 683       if (!raw_score)
 684         if (mode & USE_LIKE_LINE)
 685           scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b);
 686         else {
 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;
 690         }
 691       else scores[i][(i+j)%POSNUM] = singlescores[i][j];  /*Raw_score*/
 692
 693   }
 694
 695
 696 }
 697 /***************************************************************************/
 698
 699
 700
 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)
 702 {
 703   int i,j;
 704   double sc2scores[MAXSEQLEN][POSNUM];
 705
 706   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 707
 708
 709   switch_tables(table);
 710
 711   /* score the sequence using sc2seq */
 712     scseqadj(lib,
 713             sequence.seq, sequence.seqlen,
 714             sc2scores,
 715             offsets,
 716             maxscore, mode, 0);
 717
 718 /*** Do a SingleCoil type score. ****/
 719 /***
 720     scorestock(sequence.seq, sequence.seqlen,
 721             sc2scores,offsets,
 722             maxscore,ProlineFreeWindow, 0);  **/ /** Score using our table. */
 723
 724
 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];
 731
 732     /* do other registers */
 733     for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 734                            i+j%POSNUM  **/
 735       scores[i][(i+j)%POSNUM] = sc2scores[i][j];
 736     }
 737   }
 738
 739 }
 740
 741
 742
 743
 744
 745 /***************************************************************************/
 746
 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],
 750                    int offset_to_use)
 751 {
 752   int i,j;
 753   double stockscores[MAXSEQLEN][POSNUM];
 754   char offsets[MAXSEQLEN];
 755
 756
 757   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 758
 759   /* score the thing */
 760   scorestock(
 761             sequence.seq, sequence.seqlen,
 762             stockscores,
 763             offsets,
 764             maxscore,ProlineFreeWindow, 1);  /** Score using STOCK table. */
 765
 766
 767   *maxscore= stocksc2likelihood(*maxscore);
 768
 769
 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 */
 777     for (j=0; j<7; ++j)
 778       scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]);
 779   }
 780
 781
 782 }
 783
 784
 785 /***********************************************************/
 786
 787
 788
 789
 790
 791
 792 void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1],
 793                   int offset_to_use, int preprocessor)
 794 {
 795   int i,j;
 796
 797   if (!sequence.reg)
 798     fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");
 799
 800
 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. */
 804
 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])
 808           scores[i][j] = 1;
 809         else scores[i][j]=0;
 810
 811       else if ( (offset_to_use == -1) &&  preprocessor) /* Do all registers */
 812         if (sequence.reg[i] != '.')
 813           scores[i][j] = 1;
 814         else scores[i][j]=0;
 815
 816       else  /* just give likelihood 1 to offset_to_use. */
 817         if ( (sequence.reg[i] != '.') && (j ==  (i+offset_to_use)%POSNUM))
 818           scores[i][j] = 1;
 819         else scores[i][j]=0;
 820
 821     if (!sequence.reg) scores[i][7]=0;
 822     else if (sequence.reg[i] != '.')
 823       scores[i][7]=1;
 824     else scores[i][7]=0;
 825   }
 826
 827 }
 828
 829
 830
 831
 832
 833
 834
 835 /********************************************************/
 836
 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],
 860                                int offset_to_use,
 861                                char offsets[MAXSEQLEN],
 862                                int start_at_reg_shift, double bound,
 863                                int avg_max)
 864 {
 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];
 873   int regist;
 874
 875   double total_coil_length[POSNUM +1];  /* A weighted length. */
 876   double CoilLike_MaxOffset[MAXSEQLEN];
 877
 878
 879
 880   for (j=0; j<POSNUM+1; j++) {
 881     start_coil[j] = -1;
 882     total_coil_score[j]=0;
 883     avg_coil_score[j]=0;
 884     total_coil_length[j]=0;
 885     max_CoilLike[j]=0;
 886     total_max_StructLike[j]=0;
 887     number_max_residues[j]=0;
 888   }
 889
 890   while (i<=seq.seqlen) {
 891     if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];
 892
 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;
 897       else regist =7;
 898
 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++) {
 908             if (k!=7)
 909               if (avg_max == 0)
 910                 StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
 911               else
 912                 StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
 913             else
 914               if (avg_max == 0)
 915                 StructureCoilLike[j][7]= avg_coil_score[k];
 916               else
 917                 StructureCoilLike[j][7]= max_StructLike[k];
 918           }
 919         }
 920         total_coil_score[k] =0;      /***Resets for "not in coil"  ***/
 921         max_CoilLike[k]=0;
 922
 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 */
 926                                              /** be 0**/
 927
 928 /***        StructureCoilLike[i][regist]= StructureResLike[i][regist];  **/
 929                                       /** If not in coil above bound make   **/
 930                                       /** the coil like the same as the     **/
 931                                       /** regist like. **/
 932             start_coil[k] = -1;
 933             total_coil_length[k] =0;
 934             max_CoilLike[k]=0;
 935           }
 936           else {  /* Changing registers. */
 937             start_coil[k]= i;
 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;
 942             }
 943             else if (CoilLike[i][regist] == max_CoilLike[k]) {
 944               total_max_StructLike[k] += StructureResLike[i][regist];
 945               number_max_residues[k]++;
 946             }
 947
 948             if (1) {   /** Do weighted_avg **/
 949               int reg_to_use=regist;
 950               if (regist ==7)
 951                 reg_to_use= (i+offsets[i])%POSNUM;
 952
 953               total_coil_length[k]= CoilLike[i][reg_to_use];
 954               total_coil_score[k] = StructureResLike[i][regist]*
 955                 CoilLike[i][reg_to_use];
 956 /*********
 957               if (total_coil_length[k] ==0)
 958                 fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
 959                         i, 'a'+reg_to_use);
 960 *************/
 961             }
 962           }
 963       }   /** Ends if ending a coil.  **/
 964
 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;
 970         }
 971         else  if (CoilLike[i][regist] == max_CoilLike[k]) {
 972           total_max_StructLike[k] += StructureResLike[i][regist];
 973           number_max_residues[k]++;
 974         }
 975
 976         if (1) {  /* Weighted average scores.  */
 977           int reg_to_use=regist;
 978           if (regist ==7)
 979             reg_to_use= (i+offsets[i])%POSNUM;
 980
 981           total_coil_score[k] += StructureResLike[i][regist]*
 982                                            CoilLike[i][reg_to_use];
 983           total_coil_length[k] += CoilLike[i][reg_to_use];
 984         }
 985         if (start_coil[k] == -1)
 986           start_coil[k] =i;
 987       }
 988     }
 989
 990     i++;
 991   }
 992
 993  /*  Need the following hack to get the correct register in combo off method.*/
 994 /************
 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];
1000     }
1001   else
1002 *************/
1003
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];
1008
1009 }