sources/multicoil/sc2seq_interface_always_likeline.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 #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"
  22
  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];
  25
  26
  27 double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];
  28
  29 double (*pprobs)[POSNUM];
  30 double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];
  31
  32
  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];
  37
  38
  39 #define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
  40
  41
  42
  43 void *Malloc (n)
  44      int n;
  45 {
  46   void *ptr;
  47   ptr = (void*)malloc(n);
  48   if (ptr) return ptr;
  49   fprintf(stderr,"Memory shortage!!\n");
  50   exit(-1);
  51 }
  52
  53
  54
  55
  56 void switch_tables(int table)
  57 {
  58   extern int WINDOW, window_length[MAX_TABLE_NUMBER];
  59
  60   pprobs= many_pprobs[table];
  61   pprobp= many_pprobp[table];
  62   WINDOW = window_length[table];
  63 }
  64
  65
  66
  67 /******** Used to subtract off a sequence from the table.  ***/
  68
  69 void recalc_prob(Sequence sequence, int input_table, int mode)
  70 {
  71   int ProlineFreeWindow = mode & PROLINE_FREE_MODE;
  72
  73   /* Remove the sequence from the positive table before computing on it. */
  74   add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
  75              ptotalp[input_table],
  76              sequence.reg, sequence.seq, sequence.seqlen, -1);
  77
  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],
  85              ptotalp[input_table],
  86              sequence.reg, sequence.seq, sequence.seqlen, 1);
  87
  88   #ifdef AVG_WEIRD_PROBS
  89     calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
  90   #else
  91     calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
  92   #endif
  93
  94 }
  95
  96
  97
  98
  99
 100
 101
 102 void PairCoilData (FILE *fgin, FILE **fpin,
 103                    int ProlineFree, int table_to_remove_from,
 104                    char *input_filename, int mode)
 105 {
 106   long gfreqs[AANUM]; long gtotals;
 107   long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM];
 108
 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;
 114
 115
 116   int i=0;
 117
 118   if (fgin) {
 119     readgenfreq(fgin, &gtotals, gfreqs, gtotalp, gfreqp, stderr);
 120     calcgprobs(gfreqs, gtotals, gprobs);
 121     calcgprobp(gfreqp, gtotalp, gprobp);
 122   }
 123   else setgprob_to_uniform(gprobs, gprobp);
 124
 125   #ifdef AVG_WEIRD_PROBS
 126     calc_weird_gprob_avg(gprobp, gprobs);
 127   #else
 128     calc_weird_gprob_sum(gprobp, gprobs);
 129   #endif
 130
 131
 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);
 135
 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);
 141
 142       estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
 143                                 pfreqp[i], ptotalp[i], many_pprobp,
 144                                 ProlineFree);
 145     }
 146     else {
 147       fprintf(stderr,"\nCould not make input table number %d.",i);
 148       exit(-1);
 149     }
 150
 151
 152 #ifdef AVG_WEIRD_PROBS        /* make table.     */
 153     calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);
 154 #else
 155     calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);
 156 #endif
 157
 158   }
 159
 160   switch_tables(0);
 161
 162 }
 163
 164
 165 NEWCOILSinit ()
 166 {
 167   initialize_relative();
 168 }
 169
 170 #define NEGCCRATIO 30
 171 #define NEGMEAN 0.77
 172 #define NEGSD 0.20
 173 #define CCMEAN 1.63
 174 #define CCSD 0.24
 175 #define SQRT2PI 2.5066282746310002
 176
 177
 178
 179 /** To compute the paircoil likelihood line. */
 180
 181
 182 void likelihood_from_txt(char *txt_filename, double *m, double *b1)
 183 {
 184   FILE *txt_file;
 185   int mode2=0;
 186   double percent_coil=PERCENT_COILED;
 187   double score1;
 188   int filenum=1;
 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];
 195   long *counts[2];
 196   double  *means[2];
 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;
 201   int place;
 202
 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]);
 207
 208
 209   mode2 |= RATIO_TRAP_METHOD;
 210   mode2 |= PERCENT_COIL_MODE;
 211
 212
 213
 214    for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
 215      counts[0][i] = 0;
 216      means[0][i]=0;
 217     }
 218
 219
 220   txt_file = sopen (txt_filename, "r");
 221   while (fread(&score1,sizeof(double),1,txt_file)) {
 222     if (score1<MINSCORE || score1>-(MINSCORE)) {
 223 /*
 224   fprintf(stderr,"Strange score: %lf\n",score1);
 225   exit(-1);
 226 */
 227       number_strange++;
 228     }
 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);
 233       /* exit(-1); */
 234       counts[0][-2*MINSCORE]++;
 235       means[0][-2*MINSCORE]++;
 236       number_high++;
 237       maxplaces[0]=-2*MINSCORE;
 238     }
 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;
 245         /* exit(-1); */
 246     }
 247
 248     else {
 249       counts[0][place]++;
 250       means[0][place] += score1;
 251       if (place>maxplaces[0]) maxplaces[0]=place;
 252       if (place<minplaces[0]) minplaces[0]=place;
 253       ++n[0];
 254     }
 255   }
 256
 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);
 259    sclose(txt_file);
 260
 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*/
 265
 266
 267    /*** We need a mean***/
 268
 269    mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
 270                          maxplaces[0],n[0]);
 271
 272
 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])) {
 277       ++n_halfgauss[0];
 278       sum2[0] += (score1-mean[0])*(score1-mean[0]);
 279     }
 280   sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
 281   down[0]= (double) 2*n_halfgauss[0]/n[0];
 282
 283   fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
 284   sclose(txt_file);
 285
 286   minplace=minplaces[0];  maxplace=maxplaces[0];
 287
 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,
 294              a,b,.20,.80,.03);
 295   sclose(outfile);
 296   sclose(outpsfile);
 297   {
 298     char clean_command[200];
 299
 300     strcpy(clean_command, "rm -f ");
 301
 302 /*** The following removes the pir.txt output file. */
 303     system(strcat(clean_command,txt_filename));
 304
 305     }                            /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
 306                                  /** WRITE TO SAME FILESNAMES.             **/
 307   {
 308     char clean_command[200];
 309
 310     strcpy(clean_command, "rm -f ");
 311
 312     system(strcat(clean_command,PIR_LIKES));
 313
 314     }
 315
 316   *m=a[0]; *b1=b[0];
 317 }
 318
 319
 320
 321
 322
 323 /*  Look in file .likelihood for likelihood line.   If not there, compute */
 324 /*  it from pir.                                               */
 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)
 326 {
 327  char *tmp;
 328  int libnumber=0;
 329  int  i, cnt=0;
 330  int current_lib;
 331  char buf[80];
 332
 333  FILE *pir;
 334  FILE *linefile;
 335  FILE *txt_file1, *txt_file2;
 336  double score, scs[MAXSEQLEN][POSNUM];
 337  Sequence sequence;
 338  char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
 339  char title[500], code[500];
 340  int seqlen;
 341  int likefile;
 342  int found_like_lib=0;  int found_single=0;
 343
 344  sequence.seq=seq;
 345  sequence.reg=reg;
 346  sequence.title=title;
 347  sequence.code=code;
 348
 349
 350  for (i=0; i<functnum; i++)
 351    libnumber |= (int) pow(2,lib[i]);
 352  fprintf(stderr,"\nlibnumber= %d",libnumber);
 353
 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++) {
 357    current_lib = -1;
 358    found_like_lib =0;
 359    found_single =0;
 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);
 364      else {
 365        strcpy(likelihoods[likefile],getenv("HOME"));
 366        strcat(likelihoods[likefile],"/.likelihood");
 367      }
 368    }
 369
 370
 371
 372    if (linefile=sopen(likelihoods[likefile],"r")) {
 373      while ( (fgets(buf,80,linefile)) && !(found_single && found_like_lib))  {
 374        if (buf[0] == '0') {
 375          sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
 376          found_single=1;
 377        }
 378        else if (!found_like_lib) {
 379          sscanf(buf, "%d %lf %lf",&current_lib, &m[likefile], &b[likefile]);
 380          if (current_lib == libnumber)  found_like_lib=1;;
 381        }
 382      }
 383      sclose(linefile);
 384
 385      if (!(found_like_lib && found_single)) {
 386        if (!found_single)
 387          fprintf(stderr,"\nSingles likelihood function %d unknown.",likefile);
 388
 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);
 392
 393
 394        if (pir_name[0]==',') {
 395          fprintf(stderr,"\nUnable to continue.  Please input the location of the pir in the .paircoil file.");
 396          exit(-1);
 397        }
 398        else {
 399          fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
 400
 401          pir = sopen(pir_name, "r");
 402
 403          if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
 404          if (!found_single)   txt_file2 = sopen(SINGLE_PIR_TXT,"w");
 405
 406          for (cnt=1;getseq2(pir,&sequence); cnt++)
 407            if (sequence.seqlen >= MINWINDOW) {
 408              seqlen= sequence.seqlen;
 409              if (!(cnt %10)) {
 410                putc('.',stderr);
 411                fflush(stderr);
 412              }
 413
 414              WINDOW=PAIRWINDOW;
 415              if (!found_like_lib) {
 416                switch_tables(likefile);
 417                scseqadj(lib,seq,seqlen,scs,offsets, &score,
 418                         mode,0);
 419                txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
 420              }
 421              if (!found_single) {
 422                switch_tables(likefile);
 423
 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);
 428              }
 429
 430            }  /* Ends if sequence > MINWINDOW, and hence the for */
 431               /* loop that reads from pir. */
 432          sclose(pir);
 433
 434
 435          if (!found_like_lib) {    /* Now compute the lib likelihood line. */
 436            sclose(txt_file1);
 437
 438            likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
 439
 440            linefile=sopen(likelihoods[likefile],"a");
 441            fprintf(linefile,"%d %lf %lf\n",libnumber,m[likefile], b[likefile]);
 442
 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");
 448
 449            sclose(linefile);
 450          }
 451
 452          if (!found_single) {  /*  Now compute the singles method like line */
 453            sclose(txt_file2);
 454
 455            likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
 456                                &b_single[likefile]);
 457
 458            linefile=sopen(likelihoods[likefile],"a");
 459            fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile],
 460                    b_single[likefile]);
 461
 462            fprintf(linefile,"\t (lib =");
 463                 /* Now print human readable lib number. */
 464            fprintf(linefile," Singles Method");
 465            fprintf(linefile," )\n");
 466
 467            sclose(linefile);
 468          }
 469        }      /* Ends else from if (pir_name == ',') */
 470
 471      }       /*  Ends section for file not having lib or singles like line */
 472
 473    }         /*  Ends section where search likelihoods[likefile] file. */
 474    else   {
 475      fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
 476      exit(-1);
 477    }
 478  }          /*  Ends count through tables=linefiles.  */
 479 }
 480
 481
 482
 483 double stocksc2likelihood (double sc)
 484 {
 485   double gg, gcc;
 486
 487   sc = exp(sc/28);
 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));
 493 }
 494
 495
 496
 497 /*  Likelihood line is mx +b */
 498 double sc2likelihood (double sc, double m, double b)
 499 {
 500   double likelihood;
 501
 502   likelihood = m*sc + b;
 503
 504   if (likelihood>1) return(1);
 505   else if (likelihood<0) return(0);
 506   else return(likelihood);
 507 }
 508
 509
 510 static int current_seq_number=0;
 511
 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)
 513 {
 514   int i,j;
 515   double sc2scores[MAXSEQLEN][POSNUM];
 516
 517   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 518
 519   double  frac_coiled_this_table;
 520   extern double init_class_prob[MAX_CLASS_NUMBER];
 521
 522   frac_coiled_this_table = init_class_prob[table];
 523
 524
 525   switch_tables(table);
 526
 527   /* score the sequence using sc2seq */
 528   scseqadj(lib,
 529             sequence.seq, sequence.seqlen,
 530             sc2scores,
 531             offsets,
 532             maxscore, mode, 0);
 533
 534
 535   if (!raw_score)
 536     if (mode & USE_LIKE_LINE)
 537       *maxscore= sc2likelihood(*maxscore,m,b);
 538     else {
 539       *maxscore= frac_coiled_this_table * exp(*maxscore);
 540     }
 541
 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]];
 547       else
 548         if (mode & USE_LIKE_LINE)
 549           scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
 550         else   {
 551           scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
 552           if (scores[i][7] >1) scores[i][7]=1;
 553         }
 554
 555     else    /*  Offset to use is 0-6. */
 556       if (raw_score)
 557         scores[i][7] = sc2scores[i][offset_to_use];
 558       else
 559         if (mode & USE_LIKE_LINE)
 560           scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
 561         else {
 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;
 565         }
 566
 567     /* do other registers */
 568     for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 569                            i+j%POSNUM  **/
 570       if (raw_score)
 571         scores[i][(i+j)%POSNUM] = sc2scores[i][j];
 572       else
 573         if (mode & USE_LIKE_LINE)
 574           scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
 575         else  {
 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;
 579         }
 580     }
 581   }
 582
 583 }
 584
 585
 586
 587
 588
 589 /***************************************************************************/
 590
 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,
 595                         int raw_score)
 596 {
 597   int i,j;
 598   double singlescores[MAXSEQLEN][POSNUM];
 599   char offsets[MAXSEQLEN];
 600
 601   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 602
 603   double  frac_coiled_this_table;
 604
 605
 606   switch_tables(table);
 607
 608   /* score the thing */
 609   scorestock(
 610             sequence.seq, sequence.seqlen,
 611             singlescores,
 612             offsets,
 613             maxscore,ProlineFreeWindow, 0);  /** Score using our table. */
 614
 615
 616   if (!raw_score)
 617     if (mode & USE_LIKE_LINE)
 618       *maxscore= sc2likelihood(*maxscore,m,b);
 619     else {
 620       /*    fprintf(stderr,"\nmaxscore = %lf", *maxscore);
 621        */
 622       *maxscore= frac_coiled_this_table * exp(*maxscore);
 623       /*    fprintf(stderr,"  maxlike  = %lf",  *maxscore);
 624        */
 625     }
 626
 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. */
 630       if (!raw_score)
 631         if (mode & USE_LIKE_LINE)
 632           scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b);
 633         else {
 634           scores[i][7] = frac_coiled_this_table*
 635                       exp(singlescores[i][offsets[i]]);
 636           if (scores[i][7] >1) scores[i][7]=1;
 637         }
 638       else scores[i][7] = singlescores[i][offsets[i]];
 639
 640     else                  /*  Use offset_to_use. */
 641       if (!raw_score)
 642         if (mode & USE_LIKE_LINE)
 643           scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
 644         else {
 645           scores[i][7] = frac_coiled_this_table*
 646                      exp(singlescores[i][offsets[i]]);
 647           if (scores[i][7] >1) scores[i][7]=1;
 648         }
 649       else scores[i][7] = singlescores[i][offset_to_use];
 650
 651     /* do the other registers */
 652     for (j=0; j<7; ++j)
 653       if (!raw_score)
 654         if (mode & USE_LIKE_LINE)
 655           scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b);
 656         else {
 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;
 660         }
 661       else scores[i][(i+j)%POSNUM] = singlescores[i][j];  /*Raw_score*/
 662
 663   }
 664
 665
 666 }
 667 /***************************************************************************/
 668
 669
 670
 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)
 672 {
 673   int i,j;
 674   double sc2scores[MAXSEQLEN][POSNUM];
 675
 676   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 677
 678
 679   switch_tables(table);
 680
 681   /* score the sequence using sc2seq */
 682     scseqadj(lib,
 683             sequence.seq, sequence.seqlen,
 684             sc2scores,
 685             offsets,
 686             maxscore, mode, 0);
 687
 688 /*** Do a SingleCoil type score. ****/
 689 /***
 690     scorestock(sequence.seq, sequence.seqlen,
 691             sc2scores,offsets,
 692             maxscore,ProlineFreeWindow, 0);  **/ /** Score using our table. */
 693
 694
 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];
 701
 702     /* do other registers */
 703     for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 704                            i+j%POSNUM  **/
 705       scores[i][(i+j)%POSNUM] = sc2scores[i][j];
 706     }
 707   }
 708
 709 }
 710
 711
 712
 713
 714
 715 /***************************************************************************/
 716
 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],
 720                    int offset_to_use)
 721 {
 722   int i,j;
 723   double stockscores[MAXSEQLEN][POSNUM];
 724   char offsets[MAXSEQLEN];
 725
 726
 727   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 728
 729   /* score the thing */
 730   scorestock(
 731             sequence.seq, sequence.seqlen,
 732             stockscores,
 733             offsets,
 734             maxscore,ProlineFreeWindow, 1);  /** Score using STOCK table. */
 735
 736
 737   *maxscore= stocksc2likelihood(*maxscore);
 738
 739
 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 */
 747     for (j=0; j<7; ++j)
 748       scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]);
 749   }
 750
 751
 752 }
 753
 754
 755 /***********************************************************/
 756
 757
 758
 759
 760
 761
 762 void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1],
 763                   int offset_to_use, int preprocessor)
 764 {
 765   int i,j;
 766
 767   if (!sequence.reg)
 768     fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");
 769
 770
 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. */
 774
 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])
 778           scores[i][j] = 1;
 779         else scores[i][j]=0;
 780
 781       else if ( (offset_to_use == -1) &&  preprocessor) /* Do all registers */
 782         if (sequence.reg[i] != '.')
 783           scores[i][j] = 1;
 784         else scores[i][j]=0;
 785
 786       else  /* just give likelihood 1 to offset_to_use. */
 787         if ( (sequence.reg[i] != '.') && (j ==  (i+offset_to_use)%POSNUM))
 788           scores[i][j] = 1;
 789         else scores[i][j]=0;
 790
 791     if (!sequence.reg) scores[i][7]=0;
 792     else if (sequence.reg[i] != '.')
 793       scores[i][7]=1;
 794     else scores[i][7]=0;
 795   }
 796
 797 }
 798
 799
 800
 801
 802
 803
 804
 805 /********************************************************/
 806
 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],
 830                                int offset_to_use,
 831                                char offsets[MAXSEQLEN],
 832                                int start_at_reg_shift, double bound,
 833                                int avg_max)
 834 {
 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];
 843   int regist;
 844
 845   double total_coil_length[POSNUM +1];  /* A weighted length. */
 846   double CoilLike_MaxOffset[MAXSEQLEN];
 847
 848
 849
 850   for (j=0; j<POSNUM+1; j++) {
 851     start_coil[j] = -1;
 852     total_coil_score[j]=0;
 853     avg_coil_score[j]=0;
 854     total_coil_length[j]=0;
 855     max_CoilLike[j]=0;
 856     total_max_StructLike[j]=0;
 857     number_max_residues[j]=0;
 858   }
 859
 860   while (i<=seq.seqlen) {
 861     if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];
 862
 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;
 867       else regist =7;
 868
 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++) {
 878             if (k!=7)
 879               if (avg_max == 0)
 880                 StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
 881               else
 882                 StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
 883             else
 884               if (avg_max == 0)
 885                 StructureCoilLike[j][7]= avg_coil_score[k];
 886               else
 887                 StructureCoilLike[j][7]= max_StructLike[k];
 888           }
 889         }
 890         total_coil_score[k] =0;      /***Resets for "not in coil"  ***/
 891         max_CoilLike[k]=0;
 892
 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 */
 896                                              /** be 0**/
 897
 898 /***        StructureCoilLike[i][regist]= StructureResLike[i][regist];  **/
 899                                       /** If not in coil above bound make   **/
 900                                       /** the coil like the same as the     **/
 901                                       /** regist like. **/
 902             start_coil[k] = -1;
 903             total_coil_length[k] =0;
 904             max_CoilLike[k]=0;
 905           }
 906           else {  /* Changing registers. */
 907             start_coil[k]= i;
 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;
 912             }
 913             else if (CoilLike[i][regist] == max_CoilLike[k]) {
 914               total_max_StructLike[k] += StructureResLike[i][regist];
 915               number_max_residues[k]++;
 916             }
 917
 918             if (1) {   /** Do weighted_avg **/
 919               int reg_to_use=regist;
 920               if (regist ==7)
 921                 reg_to_use= (i+offsets[i])%POSNUM;
 922
 923               total_coil_length[k]= CoilLike[i][reg_to_use];
 924               total_coil_score[k] = StructureResLike[i][regist]*
 925                 CoilLike[i][reg_to_use];
 926 /*********
 927               if (total_coil_length[k] ==0)
 928                 fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
 929                         i, 'a'+reg_to_use);
 930 *************/
 931             }
 932           }
 933       }   /** Ends if ending a coil.  **/
 934
 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;
 940         }
 941         else  if (CoilLike[i][regist] == max_CoilLike[k]) {
 942           total_max_StructLike[k] += StructureResLike[i][regist];
 943           number_max_residues[k]++;
 944         }
 945
 946         if (1) {  /* Weighted average scores.  */
 947           int reg_to_use=regist;
 948           if (regist ==7)
 949             reg_to_use= (i+offsets[i])%POSNUM;
 950
 951           total_coil_score[k] += StructureResLike[i][regist]*
 952                                            CoilLike[i][reg_to_use];
 953           total_coil_length[k] += CoilLike[i][reg_to_use];
 954         }
 955         if (start_coil[k] == -1)
 956           start_coil[k] =i;
 957       }
 958     }
 959
 960     i++;
 961   }
 962
 963  /*  Need the following hack to get the correct register in combo off method.*/
 964 /************
 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];
 970     }
 971   else
 972 *************/
 973
 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];
 978
 979 }
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007