sources/multicoil/output_cuts.c

   1 changed
   2 void finish_log(int mode, double bound, FILE *fin, FILE *flog,
   3                 FILE *flog_coil_conflicts,
   4                 FILE *fout_coils,
   5                 FILE *fout,
   6                 /* FILE *fout_raw[MAX_TABLE_NUMBER+1], */
   7                 double *m, double *b, double *m_single, double *b_single,
   8                 char lib[MAXFUNCTNUM],
   9                 char lib_differ[MAXFUNCTNUM],
  10                 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
  11                 int main_method, int main_table,
  12                 int main_preprocessor_method,
  13                 int *seqcnt, int avg_max, int which_priors, int which_priorp,
  14                 double prior_freq_single[MAX_TABLE_NUMBER],
  15                 double prior_freq_pair[MAX_TABLE_NUMBER],
  16                 int good_turing_fixed_disc, int structural_pos[POSNUM+1])
  17 to
  18 void finish_log(int mode, double bound, FILE *fin, FILE *flog,
  19                 FILE *fout_coils,
  20                 FILE *fout,
  21                 double *m, double *b, double *m_single, double *b_single,
  22                 char lib[MAXFUNCTNUM],
  23                 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
  24                 int main_method, int main_table,
  25                 int main_preprocessor_method,
  26                 int *seqcnt, int avg_max)
  27
  28 changed
  29     output_seq(lib,differ_lib, multi_lib,m,b,m_single, b_single,
  30                mode,bound,flog,flog_coil_conflicts,
  31                fout_coils,fout,  avg_max,
  32                main_method, main_preprocessor_method, main_table);
  33 to
  34     output_seq(lib, multi_lib,m,b,m_single, b_single,
  35                mode,bound,flog,
  36                fout_coils,fout,  avg_max,
  37                main_method, main_preprocessor_method, main_table);
  38
  39
  40
  41
  42
  43 changed
  44     NewSeq_nonX (mode, sequence, which_priors, which_priorp, prior_freq_single,
  45                  prior_freq_pair, good_turing_fixed_disc, structural_pos);
  46 to
  47     NewSeq_nonX (mode, sequence);
  48
  49
  50
  51
  52 cut
  53 /*** In TableDiffCoil method and PairDiffAvg method ***/
  54 /*** this computes the distribution of coil scores. ***/
  55 /*** If in pos mode it only counts coils that are   ***/
  56 /*** in the pos file.                               ***/
  57 /*** It is called from output_seq.  The results     ***/
  58 /*** are printed from finish_log to log file.       ***/
  59
  60 void count_positives(int *number_coils_at_most_half,
  61                      int *number_coils_above_half,
  62                      int *number_coils_above_3quarters,
  63                      int *number_coils_below_1quarter,
  64                      int *number_coils_missed_by_preproc,
  65                      Sequence seq, double scores[MAXSEQLEN][POSNUM+1],
  66                      double preproc_like[MAXSEQLEN][POSNUM+1],
  67                      int mode, double bound, int start_at_reg_shift,
  68                      FILE *flog,
  69                      double *minpositive)
  70      /* minpositive[0] is min differentiator score on a pos coil. */
  71      /* minpositive[1] is min differentiator score above .5 on pos coil. */
  72      /* minpositive[2] is min preproc score above bound on a pos coil.  */
  73 {
  74   int i,j;
  75   char offsets[MAXSEQLEN];
  76   int not_scored_current_coil =1;
  77   int in_pos_coil =0;
  78   int printed_to_log =0;
  79   double preproc_score_current_coil =0;
  80
  81   for (i=0; i<seq.seqlen; i++) {
  82     for (j=0; j<POSNUM; j++)
  83       if (preproc_like[i][j] == preproc_like[i][7]) {
  84         offsets[i] = (j-i) % POSNUM;
  85         if (offsets[i]<0) offsets[i] += POSNUM;
  86       }
  87   /* This is a new thing so that if an offset change is not real */
  88   /* we keep the old offset. **********/
  89     if ( (i>0) && (offsets[i-1] != -1) && (offsets[i] != offsets[i-1]) &&
  90                 (preproc_like[i][(i+offsets[i-1]) %POSNUM]
  91                           == preproc_like[i][7]) )
  92       offsets[i] = offsets[i-1];
  93   }
  94
  95
  96   if (mode & POS_MODE) {
  97     if (seq.reg)                  /* Only count coils if in posfile */
  98       for (i=0; i<seq.seqlen; i++) {
  99         if (seq.reg[i] != '.') {
 100           if (preproc_like[i][7] > preproc_score_current_coil)
 101             preproc_score_current_coil = preproc_like[i][7];  /* Score of */
 102                                                  /* coil is max over coil */
 103           in_pos_coil=1;
 104
 105           if (preproc_like[i][7]>bound)  {
 106             if ( not_scored_current_coil ||                /* A new coil. */
 107               ( (offsets[i] != offsets[i-1]) && start_at_reg_shift)) {
 108
 109               if (scores[i][7] < minpositive[0]) minpositive[0] = scores[i][7];
 110               if ((scores[i][7] < minpositive[1]) && (scores[i][7]>.5))
 111                 minpositive[1] = scores[i][7];
 112
 113               if (scores[i][7]>.5) {
 114                 (*number_coils_above_half)++;
 115                 if (scores[i][7] > .75)
 116                   (*number_coils_above_3quarters)++;
 117               }
 118               else  {
 119                 printed_to_log=1;
 120                 fprintf(flog,"ALERT: Low score at res %d in %s\n",i,seq.code);
 121                 (*number_coils_at_most_half)++;
 122                 if (scores[i][7] < .25)
 123                   (*number_coils_below_1quarter)++;
 124               }
 125               not_scored_current_coil =0;
 126             }
 127           }
 128           else { /* Coil in pos file, but not found by preproc. */
 129             if ( start_at_reg_shift && not_scored_current_coil
 130                 && (seq.reg[i]%POSNUM != ((seq.reg[i-1]+1) %POSNUM)) &&
 131                 (seq.reg[i-1] != '.') ) {
 132                                    /* Register shift in pos */
 133               (*number_coils_missed_by_preproc)++;
 134               printed_to_log=1;
 135               fprintf(flog,"ALERT:  Missed Coil at res %d in %s\n",i,seq.code);
 136               not_scored_current_coil =1;
 137             }
 138           }
 139
 140         }
 141
 142         else {   /* Not a posfile coil. */
 143           if ((in_pos_coil) && (preproc_score_current_coil < minpositive[2])
 144               && (preproc_score_current_coil>bound) )
 145             minpositive[2] = preproc_score_current_coil;
 146
 147           if ( (not_scored_current_coil) && (in_pos_coil)) {
 148                                           /* Never found current coil. */
 149               (*number_coils_missed_by_preproc)++;
 150               printed_to_log=1;
 151               fprintf(flog,"ALERT:  Missed Coil at res %d in %s\n",i,seq.code);
 152             }
 153           not_scored_current_coil =1;
 154           in_pos_coil =0;
 155         }
 156       }
 157
 158     if ((in_pos_coil) && (preproc_score_current_coil < minpositive[2])
 159         && (preproc_score_current_coil>bound) )  /* If last pos coilended at */
 160                                              /* end of seq. */
 161       minpositive[2] = preproc_score_current_coil;
 162
 163
 164     if (not_scored_current_coil && in_pos_coil) {
 165                       /* missed last coil. */
 166       (*number_coils_missed_by_preproc)++;
 167       printed_to_log=1;
 168       fprintf(flog,"ALERT:  Missed Coil at res %d in %s\n",i,seq.code);
 169     }
 170
 171     if (printed_to_log) fprintf(flog,"\n\n\n");
 172   }
 173
 174   else  /*  Not in pos mode, so score all coils found.  */
 175
 176     for (i=0; i<seq.seqlen; i++) {
 177       if (preproc_like[i][7]>bound)  {
 178         if ( not_scored_current_coil ||                /* A new coil. */
 179             ( (offsets[i] != offsets[i-1]) && start_at_reg_shift)) {
 180
 181           if (scores[i][7] < minpositive[0]) minpositive[0] = scores[i][7];
 182           if ((scores[i][7] < minpositive[1]) && (scores[i][7]>.5))
 183             minpositive[1] = scores[i][7];
 184
 185           if (scores[i][7]>.5) {
 186             (*number_coils_above_half)++;
 187             if (scores[i][7] > .75)
 188               (*number_coils_above_3quarters)++;
 189               }
 190           else  {
 191             (*number_coils_at_most_half)++;
 192             if (scores[i][7] < .25)
 193               (*number_coils_below_1quarter)++;
 194           }
 195           not_scored_current_coil =0;
 196         }
 197       }
 198       else  /* Nothing found by preproc. */
 199         not_scored_current_coil =1;
 200
 201
 202     }
 203
 204 }
 205
 206
 207 /*********************************ROUND 2*********************************/
 208
 209 /** Used to print out for deciding between 2 and 3 stranded. */
 210 void log_output3_prn(int mode, double score, double scs[MAXSEQLEN][POSNUM+1],
 211                     int seqlen, int cnt,
 212                     Sequence sequence, double bound, double upper_bound,
 213                     FILE *flog)
 214 {
 215   char *title, *code, *reg, *seq;
 216   int i, st;
 217   double max_pos, max_neg;
 218  double the_score;
 219   char regist, offset;
 220   char offsets[MAXSEQLEN];
 221
 222   title=sequence.title;
 223   code=sequence.code;
 224   reg=sequence.reg;
 225   seq=sequence.seq;
 226
 227    /*  Print out sequence number, code, title, and maxscore if > bound. */
 228   /* Note that score represents the max-likelihood by paircoil. */
 229   if (score > 0) {  /* As long as there was a coil found by paircoil.  */
 230     fprintf(flog, "%d \t%s \n\t%s \n\t%6.2lf \n", cnt, code,
 231              title, score);
 232
 233     if (sequence.reg && (mode & DEBUG_MODE))  /* Print out coil predictions like a pos file. */
 234       print_coil2(mode, seqlen, reg, scs, bound,
 235                   upper_bound, seq, flog); /* Gives register predicted by */
 236                                            /*  paircoil.  */
 237
 238
 239
 240    /* Now put scores for all the regions that score above bound.      */
 241    /* Note that if st<0 then we are not currently in a coiled region. */
 242
 243     fprintf(flog, " [");
 244     st = -1;
 245     max_pos = -HUGE_VAL; /* max==-HUGE_VAL and st==-1 equivalently mean that */
 246     max_neg = HUGE_VAL;  /* aren't currently in a coil.                      */
 247
 248     for (i = 0; i < seqlen; i++)  {
 249       the_score = -HUGE_VAL;
 250       for (offset=0; offset<POSNUM; offset++) {  /* Deals with scs[i][7]=0
 251                                                    for non-coiled region */
 252         regist= (i+ offset)%POSNUM;
 253         if (scs[i][regist]==scs[i][7]) {
 254           the_score = scs[i][7];
 255           offsets[i]=offset;
 256         }
 257       }
 258       /* This is a new thing so that if an offset change is not real */
 259       /* we keep the old offset. **********/
 260       if ( (i>0) && (offsets[i-1] != -1) && (offsets[i] != offsets[i-1]) &&
 261           (scs[i][(i+offsets[i-1]) %POSNUM] == scs[i][7]) )
 262         offsets[i] = offsets[i-1];
 263
 264
 265       if (the_score != -HUGE_VAL )  {   /* In a coil*/
 266         if (st < 0) st = i;  /* Start new coil at i after non-coiled region.*/
 267         else           /* else if i starts a new coil by changing register, */
 268                        /* then print previous coil and start new coil.      */
 269           if (offsets[i] != offsets[i-1]) {
 270             if (max_pos > -HUGE_VAL) { /* ended a coil at i-1 to print out. */
 271               if (max_pos != max_neg)
 272                 fprintf(flog,"%6.2lf,%6.2lf@%4d-%4d:%c;  ",
 273                         max_neg,max_pos,st,i-1,
 274                         numpos( (offsets[st] + st) %POSNUM));
 275               else
 276                 fprintf(flog,"%6.2lf@%4d-%4d:%c;  ",
 277                         max_pos,st,i-1, numpos( (offsets[st] + st) %POSNUM));
 278             }
 279             max_pos = scs[i][7];  /*Reset for new coil starting at i. */
 280             max_neg = scs[i][7];
 281             st = i;               /*  Start new coil at i.            */
 282           }
 283
 284         if (max_pos < scs[i][7]) max_pos = scs[i][7];
 285         if (max_neg > scs[i][7]) max_neg = scs[i][7];
 286          /* New max for current coil. */
 287       }
 288      else if (max_pos > -HUGE_VAL) { /* If just ended a coil by scoring low. */
 289        if (max_pos != max_neg)
 290          fprintf(flog,"%6.2lf@,%6.2lf%4d-%4d:%c;  ",max_neg,max_pos,st,i-1,
 291                  numpos((offsets[st]+st) %POSNUM) );
 292        else
 293          fprintf(flog,"%6.2lf%4d-%4d:%c;  ",max_pos,st,i-1,
 294                  numpos((offsets[st]+st) %POSNUM) );
 295        st = -1;
 296        max_pos = -HUGE_VAL;
 297        max_neg = HUGE_VAL;
 298      }
 299
 300     }  /* ends count on i. */
 301
 302     if (st >= 0)     /** If the end of the seq. ended a coiled region. */
 303       if (max_pos > -HUGE_VAL) {
 304         if (max_pos != max_neg)
 305           fprintf(flog,"%6.2lf,%6.2lf@%4d-%4d:%c;  ",max_neg,max_pos,st,i-1,
 306                   numpos((offsets[st]+st) %POSNUM));
 307         else
 308           fprintf(flog,"%6.2lf@%4d-%4d:%c;  ",max_pos,st,i-1,
 309                   numpos((offsets[st]+st) %POSNUM));
 310       }
 311     fprintf(flog, "End: %d]\n\n", seqlen-1);
 312
 313   }       /* Ends "if (score > bound)" */
 314 }
 315
 316
 317
 318
 319
 320
 321 /****************************************************************************/
 322
 323
 324 void tuple_output(Sequence sequence, int mode, FILE *fout,
 325                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM],
 326                   double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
 327                   char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
 328                   int number_tables, int number_multi_lib,
 329                   double bound)
 330 {
 331   int start, end = -1;
 332   int j;
 333   char *seq= sequence.seq;
 334   char *reg= sequence.reg;
 335   int seqlen= sequence.seqlen;
 336   int table;
 337   int posp = mode & POS_MODE;
 338   int above_bound_mode = mode & ABOVE_BOUND_MODE;
 339   int res_above_bound;
 340
 341   extern Sequence sequence;
 342 /*  extern FILE *debug_out; */
 343
 344   int dist;
 345
 346   if (posp && (!reg)) {
 347     /***  fprintf(stderr,"\nIgnored a sequence without posfile registers in txt_output."); **/
 348     return; }
 349
 350
 351 /****************************************************************************/
 352
 353
 354 void tuple_output(Sequence sequence, int mode, FILE *fout,
 355                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM],
 356                   double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
 357                   char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
 358                   int number_tables, int number_multi_lib,
 359                   double bound)
 360 {
 361   int start, end = -1;
 362   int j;
 363   char *seq= sequence.seq;
 364   char *reg= sequence.reg;
 365   int seqlen= sequence.seqlen;
 366   int table;
 367   int posp = mode & POS_MODE;
 368   int above_bound_mode = mode & ABOVE_BOUND_MODE;
 369   int res_above_bound;
 370
 371   extern Sequence sequence;
 372 /*  extern FILE *debug_out; */
 373
 374   int dist;
 375
 376   if (posp && (!reg)) {
 377     /***  fprintf(stderr,"\nIgnored a sequence without posfile registers in txt_output."); **/
 378     return; }
 379
 380
 381   while (nondot_advance(&start, &end, seq, reg, seqlen, posp)) {
 382     if (end-start < MINWINDOW) {
 383       continue; }                         /* If a coil is shorter than */
 384                                           /* MINWINDOW, it does not get*/
 385                                           /* output to the histogram.  */
 386
 387 /*** Note -HUGE_VAL means was a proline, should be -HUGE_VAL for all tables */
 388     for (j=start;j<end;j++) {
 389       res_above_bound = (!above_bound_mode) &&
 390                         (all_scs[0][j][offsets[0][j]] != -HUGE_VAL);
 391       if (!res_above_bound)
 392         for (table=0; table <number_tables; table++)
 393           for (dist=0; dist < number_multi_lib;dist++)
 394             res_above_bound |= (all_like[dist*number_tables+ table][j][7]
 395                                          >bound);
 396
 397
 398       if (res_above_bound && (!posp || reg[j] != '.')) {
 399 /***************
 400         if ( (all_scs[0][j][offsets[0][j]] < 0) &&
 401             (all_scs[1][j][offsets[1][j]] < 8) )
 402           fprintf(debug_out,"In sequence %s %s residue %d scored %lf, %lf\n",
 403                   sequence.code, sequence.title, j,
 404                   all_scs[0][j][offsets[0][j]],all_scs[1][j][offsets[1][j]]);
 405 ****************/
 406         fprintf(fout,"\n");
 407         for (dist=0; dist < number_multi_lib; dist++)
 408           for (table =0; table<number_tables; table++)
 409             fprintf(fout,"%lf ",all_scs[dist*number_tables + table][j]
 410                     [offsets[table+number_tables*dist][j]]);
 411       }
 412
 413     }
 414   }
 415 }
 416
 417
 418 /****************************************************************************/
 419
 420
 421
 422 /* Makes the view graphs (-g option) for paircoil. */
 423 /* If it reads the regions to be graphed from a file the format is:   */
 424 /*                                                                    */
 425 /*   >code_for_sequence1    start1 end1                               */
 426 /*   >code_for_sequence2    start2 end2                               */
 427 /*   ....                                                             */
 428 res_compare (FILE *ps, double *sc1, double *sc2, Sequence *sequence,
 429              FILE *regionfile)
 430 {
 431   char c;
 432   int  i;
 433   int reset;
 434   char codebuff[MAXLINE];
 435   int rstart, rend, seqlen;
 436   double max1=-HUGE_VAL, max2=-HUGE_VAL, rmax1=-HUGE_VAL, rmax2=-HUGE_VAL;
 437
 438   seqlen=sequence->seqlen;
 439
 440   if (regionfile) {
 441     while (1) {
 442       c=fgetc(regionfile);
 443       while (c!='>' && c!=EOF) c=fgetc(regionfile);
 444
 445       /* scan for the sequence code name, and a start and end value */
 446       if(c=='>')  {
 447         fscanf (regionfile, "%s %d %d", codebuff, &rstart, &rend);
 448         if ( !strcmp (codebuff, sequence->code))
 449           break;    /** We got the code for the sequence we are working on */
 450       }
 451       else {
 452         rstart=1;
 453         rend=seqlen;
 454         break; /** No more codes in the graph file to check, so just  */
 455         /*  break and graph the whole sequence.                */
 456       }
 457     }
 458
 459     --rstart;
 460   }
 461   else {                      /* Just graph the whole seq if no region file. */
 462     rstart = 0;
 463     rend = seqlen;
 464   }
 465
 466   for  (i=rstart; i < rend; ++i)
 467     if (rmax1<sc1[i])               /* sc1 is Stock score.  */
 468       rmax1=sc1[i];
 469   for  (i=rstart; i < rend; ++i)
 470     if (rmax2<sc2[i])              /* sc2 is PairCoil score. */
 471       rmax2=sc2[i];
 472   for  (i=0; i < seqlen; ++i)
 473     if (max1<sc1[i])
 474       max1=sc1[i];
 475   for  (i=0; i < seqlen; ++i)
 476     if (max2<sc2[i])
 477       max2=sc2[i];
 478   fputs ("%! \n", ps);
 479   fputs("90 rotate 50 -200 translate\n", ps);
 480   fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
 481   fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
 482   fputs("/Times-Roman findfont 18 scalefont setfont\n",ps);
 483   fprintf (ps, "50 150 moveto (%s) show 50 100 moveto(%s) show\n", sequence->code, sequence->title);
 484   fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps);
 485   fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps);
 486
 487   fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps);
 488   fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps);
 489   fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2);  fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2);
 490   fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend );
 491   --rend;
 492
 493     --rstart;
 494   }
 495   else {                      /* Just graph the whole seq if no region file. */
 496     rstart = 0;
 497     rend = seqlen;
 498   }
 499
 500   for  (i=rstart; i < rend; ++i)
 501     if (rmax1<sc1[i])               /* sc1 is Stock score.  */
 502       rmax1=sc1[i];
 503   for  (i=rstart; i < rend; ++i)
 504     if (rmax2<sc2[i])              /* sc2 is PairCoil score. */
 505       rmax2=sc2[i];
 506   for  (i=0; i < seqlen; ++i)
 507     if (max1<sc1[i])
 508       max1=sc1[i];
 509   for  (i=0; i < seqlen; ++i)
 510     if (max2<sc2[i])
 511       max2=sc2[i];
 512   fputs ("%! \n", ps);
 513   fputs("90 rotate 50 -200 translate\n", ps);
 514   fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
 515   fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
 516   fputs("/Times-Roman findfont 18 scalefont setfont\n",ps);
 517   fprintf (ps, "50 150 moveto (%s) show 50 100 moveto(%s) show\n", sequence->code, sequence->title);
 518   fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps);
 519   fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps);
 520
 521   fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps);
 522   fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps);
 523   fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2);  fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2);
 524   fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend );
 525   --rend;
 526
 527     --rstart;
 528   }
 529   else {                      /* Just graph the whole seq if no region file. */
 530     rstart = 0;
 531     rend = seqlen;
 532   }
 533
 534   for  (i=rstart; i < rend; ++i)
 535     if (rmax1<sc1[i])               /* sc1 is Stock score.  */
 536       rmax1=sc1[i];
 537   for  (i=rstart; i < rend; ++i)
 538     if (rmax2<sc2[i])              /* sc2 is PairCoil score. */
 539       rmax2=sc2[i];
 540   for  (i=0; i < seqlen; ++i)
 541     if (max1<sc1[i])
 542       max1=sc1[i];
 543   for  (i=0; i < seqlen; ++i)
 544     if (max2<sc2[i])
 545       max2=sc2[i];
 546   fputs ("%! \n", ps);
 547   fputs("90 rotate 50 -200 translate\n", ps);
 548   fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
 549   fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
 550   fputs("/Times-Roman findfont 18 scalefont setfont\n",ps);
 551   fprintf (ps, "50 150 moveto (%s) show 50 100 moveto(%s) show\n", sequence->code, sequence->title);
 552   fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps);
 553   fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps);
 554
 555   fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps);
 556   fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps);
 557   fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2);  fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2);
 558   fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend );
 559   --rend;
 560
 561   fprintf(ps, "/vx %lf def /vy 5 def\n",(11*72-100)/(double) sequence->seqlen);
 562   fputs("/adjust {exch vx mul exch vy mul} def\n",ps);
 563   fputs("/-Inf -150 def\n",ps);
 564   fputs("/m {adjust moveto vx 0 rlineto} def\n",ps);
 565   fputs("/l {adjust lineto vx 0 rlineto} def\n",ps);
 566   fputs("/min {1 index 1 index lt {pop} {exch pop} ifelse} def\n",ps);
 567   fputs("/bar {exch vx mul 0 moveto 0 exch vy mul rlineto stroke} def\n",ps);
 568   fputs("/foo {2 index 2 index 2 index\
 569       1 index 1 index gt\
 570       {0.3 0.7 1 setrgbcolor pop bar} {0.5 1 0 setrgbcolor exch pop bar} ifelse\n      min 1 0 0 setrgbcolor bar} def\n",ps);
 571
 572   /* Add axes to PS output */
 573   fputs("/Times-Roman findfont 10 scalefont setfont\n",ps);
 574   fputs("1 setlinewidth 0.7 setgray 0 15 adjust moveto 0 -50 adjust lineto stroke \n",ps);
 575   for(i=-50; i<=15; i+= 5) {
 576     fprintf(ps,"0 %d adjust moveto -5 0 rmoveto 10 0 rlineto stroke \n",i);
 577     fprintf(ps,"0 %d adjust moveto \n", i);
 578     fprintf(ps,"7 0 rmoveto (%d) show \n",i);
 579   }
 580   fprintf(ps,"1 setlinewidth 0 setgray %d 15 adjust moveto %d -50 adjust lineto stroke \n",sequence->seqlen, sequence->seqlen);
 581   for(i=-50; i<=15; i+= 5) {
 582     fprintf(ps,"%d %d adjust moveto -5 0 rmoveto 10 0 rlineto stroke \n",sequence->seqlen, i);
 583     fprintf(ps, "%d %d adjust moveto \n",sequence->seqlen, i);
 584     fprintf(ps,"-7 0 rmoveto (%5.2f) rshow \n",exp(i/(double)WINDOW));
 585   }
 586
 587   fputs("vx 2 div setlinewidth 0.5 setgray\n",ps);
 588   fputs("2 setlinewidth \n",ps);
 589   fputs("0 setgray [2] 0 setdash \n",ps);
 590   for(reset=1 ,i = 0; i <sequence->seqlen; ++i){
 591     int inf;
 592     inf=(sc1 [i] == -HUGE_VAL);
 593     if(!inf)
 594       fprintf (ps, "%d %lf %s\n" , i, sc1 [i], reset?"m":"l");
 595     reset=inf;
 596   }
 597   fputs("stroke \n",ps);
 598   fputs("0.7 setgray [] 0 setdash\n",ps);
 599   for(reset=1 ,i = 0; i <sequence->seqlen; ++i){
 600     int inf;
 601     inf=(sc2 [i] == -HUGE_VAL);
 602     if(!inf)
 603       fprintf (ps, "%d %lf %s\n" , i, sc2 [i], reset?"m":"l");
 604     reset=inf;
 605   }
 606   fputs("stroke \n",ps);
 607   fprintf(ps,"0 setlinewidth 0 setgray 0 0 adjust moveto %d 0 adjust lineto stroke \n",sequence->seqlen);
 608   fprintf(ps,"0 setlinewidth 0 setgray 0 10 adjust moveto %d 10 adjust lineto stroke \n",sequence->seqlen);
 609   fprintf(ps,"0 setlinewidth 0 setgray 0 -10 adjust moveto %d -10 adjust lineto stroke \n",sequence->seqlen);
 610   for (i=49;  i <sequence->seqlen; i+=50)
 611     fprintf(ps, "%d 0 adjust moveto vx 2 div  -5 rmoveto  0 10  rlineto stroke\n", i);
 612   fprintf(ps, "%d -15 adjust moveto vx 2 div 0 rmoveto 0 30 vy mul rlineto stroke\n", rstart);
 613   fprintf(ps, "%d -15 adjust moveto vx 2 div 0 rmoveto 0 30 vy mul rlineto stroke\n", rend);
 614
 615   fputs("showpage\n\n", ps);
 616 }
 617
 618
 619
 620
 621
 622 /**************************************************************************/
 623 /************The following function is intended to give a user friendly ***/
 624 /************output of where the coil predictors and coil distinguishers **/
 625 /************come up with different predictions using different tables.  **/
 626 /************If all_scores=preproc_scores, then we are just concerned    **/
 627 /************with the coil predictions.  Otherwise we want the 2-3       **/
 628 /************stranded predictions too. **/
 629 /***  Note that for 2-3 stranded distinguishing, this is only called if  **/
 630 /***  doing coil scores, not residue scores.  **/
 631 /*** Note right now in the printing we only deal with the max scoring **/
 632 /*** offset.  Maybe make an interactive thing to deal with others?    **/
 633 /*** NOTE that we use MAX_TABLE_NUMBER +1 for most things, since we use **/
 634 /*** the number_tables +1 entry for pos file info if in pos mode.     **/
 635
 636 /** If avg_or_max is 0 then average residue scores over coil.  ***/
 637 /** Otherwise do max residue score for the coil score.   ***/
 638
 639 void coil_conflicts(int mode,
 640                     double all_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1],
 641              double all_preproc_like[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1],
 642                     Sequence sequence, int seqnum, double bound,
 643                     double bound_preprocessor, int number_tables,
 644                     FILE *flog_coil_conflicts, int avg_or_max)
 645
 646 {
 647   double all_preproc_coil_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1];
 648    /* This will give each coil found the avg_max_score over that coil. */
 649  int i,j, tablenum, offset, row;
 650   char offsets[MAX_TABLE_NUMBER +1 ][MAXSEQLEN];
 651
 652
 653   int coil_starts[MAX_TABLE_NUMBER +1][MAXNUMCOILS][POSNUM +1];
 654   int coil_ends[MAX_TABLE_NUMBER +1][MAXNUMCOILS][POSNUM +1];
 655   /* coil_...[table][i][offset] holds the start and end of the ith coil */
 656   /* in offset. The end of the coil list is marked by a -1. */
 657
 658   int current_coil_number[MAX_TABLE_NUMBER+1][POSNUM+1];
 659   int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1];
 660   int current_coil;
 661
 662   int label_for_coil[MAX_TABLE_NUMBER+1][MAXNUMCOILS];
 663                /* Gives each coil in each table a number, so overlapping */
 664                /* coils have same range of numbers. */
 665
 666   int max_coil_number;
 667
 668   int coil_above_bound =0;
 669   int label_in_row=0;
 670   int printed_blank_line;
 671   int number_tables_here;
 672
 673
 674   if ( (mode & POS_MODE) && sequence.reg)
 675     number_tables_here= number_tables+1;  /* Locally add 1 to number_tables */
 676                                           /* since number_tables entry will */
 677                                           /* contain pos file info.         */
 678   else number_tables_here = number_tables;
 679
 680   for (tablenum=0; tablenum<number_tables_here; tablenum++)
 681     for (offset =0; offset<=POSNUM; offset++) {
 682       current_coil_number[tablenum][offset] = 0;
 683       number_coils[tablenum][offset]=0;
 684     }
 685
 686   get_coil_score(sequence, bound_preprocessor, all_preproc_like,
 687                  all_preproc_coil_scores, offsets, number_tables,
 688                  coil_starts, coil_ends, number_coils, mode, avg_or_max,
 689                  all_scores,bound);
 690
 691
 692   max_coil_number = get_coil_order(coil_starts, coil_ends, label_for_coil,
 693                                    number_tables_here, number_coils);
 694
 695
 696
 697   for (tablenum=0; tablenum<number_tables_here; tablenum++)
 698     if (number_coils[tablenum][7] >0) coil_above_bound =1;
 699
 700
 701   if (coil_above_bound)
 702     fprintf(flog_coil_conflicts, "\n\n%d \t%s \n\t%s", seqnum, sequence.code,
 703                               sequence.title);
 704
 705
 706
 707 /*  if (all_scores == all_preproc_like)*/   /* Just concerned with finding coils */
 708
 709
 710
 711
 712 /*** Put 3  coils per row.  ***/
 713   for (row =0; 3*row <= max_coil_number; row++) {
 714     printed_blank_line = 0;
 715     for (tablenum=0; tablenum<number_tables_here; tablenum++) {
 716       current_coil = current_coil_number[tablenum][7];
 717
 718       if (current_coil < number_coils[tablenum][7])
 719         for (i=0; i<3; i++)
 720           if (label_for_coil[tablenum][current_coil]  == 3*row + i)
 721             label_in_row =1;
 722
 723       if ( (label_in_row) && (current_coil < number_coils[tablenum][7]) )  {
 724                                                   /*Label if there is a coil */
 725
 726         if ((row>0) && (!printed_blank_line)) {
 727           fprintf(flog_coil_conflicts,"\n");  /* To print pretty. */
 728           printed_blank_line=1;
 729         }
 730
 731         if (tablenum != number_tables)  /* Regular table. */
 732           fprintf(flog_coil_conflicts,"\nTab%d", tablenum);
 733         else     /* Posfile coils.  */
 734           fprintf(flog_coil_conflicts,"\nPos ");
 735
 736         /** Print score@start-end:offset,reg **/
 737         for (i=0; (i<3) && (current_coil < number_coils[tablenum][7]) ; i++)
 738                                /* 3 coils per row...coil label is 3*row+i */
 739           if ( label_for_coil[tablenum][current_coil] == 3*row + i) {
 740
 741             if (tablenum != number_tables) {  /* Just a regular table. */
 742               if (all_scores == all_preproc_like)
 743                 fprintf(flog_coil_conflicts,"      %5.2lf",
 744                      all_preproc_coil_scores[tablenum]
 745                                 [coil_starts[tablenum][current_coil][7]] [7]);
 746               else                  /* Both preproc. and coil_diff score */
 747                 fprintf(flog_coil_conflicts," %5.2lf,%4.2lf",
 748                      all_preproc_coil_scores[tablenum]
 749                                 [coil_starts[tablenum][current_coil][7]] [7],
 750                      all_scores[tablenum]
 751                         [coil_starts[tablenum][current_coil][7]] [7]);
 752
 753               /* Now print coil location, etc. */
 754               fprintf(flog_coil_conflicts,"@%4d-%4d:%d,%c",
 755                      coil_starts[tablenum][current_coil][7],
 756                      coil_ends[tablenum][current_coil][7],
 757                      offsets[tablenum][coil_starts[tablenum][current_coil][7]],
 758                      'a' + (( coil_starts[tablenum][current_coil][7] +
 759                      offsets[tablenum]
 760                           [coil_starts[tablenum][current_coil][7]]) %POSNUM) );
 761             }
 762             else    /*  Posfile coil so don't print score.  */
 763               fprintf(flog_coil_conflicts,"            %4d-%4d:%d,%c",
 764                      coil_starts[tablenum][current_coil][7],
 765                      coil_ends[tablenum][current_coil][7],
 766                      offsets[tablenum][coil_starts[tablenum][current_coil][7]],
 767                      'a' + (( coil_starts[tablenum][current_coil][7] +
 768                      offsets[tablenum]
 769                           [coil_starts[tablenum][current_coil][7]]) %POSNUM) );
 770
 771
 772             current_coil_number[tablenum][7]++;
 773             current_coil++;
 774
 775           }
 776           else   /** The current_coil does not have current_label.  **/
 777                  /** Made each coil take (19) 24 spaces + 1 at start.     **/
 778             fprintf(flog_coil_conflicts, "                         ");
 779       }
 780       label_in_row =0;
 781     }
 782   }
 783 }
 784 /**************************************************************************/
 785
 786
 787 /*** This will compute the max scoring coils and their offsets... AND **/
 788 /*** compute the scores for all offsets.....  *****/
 789 /*** Note that coils start and end at reg shifts only for max offset. */
 790 /*** offset -1 indicates non-coiled region.  */
 791 /*** Also computes the coil_starts[table][i][offset] for the i-th coil **/
 792 /*** found   in offset by the table. Also computes the coil ends.      **/
 793 /** If in pos mode than store info in slot for number_tables   **/
 794
 795 /** If avg_or_max is 0 then average residue scores over coil.  ***/
 796 /** Otherwise do max residue score for the coil score.   ***/
 797
 798 get_coil_score(Sequence sequence, double bound_preprocessor,
 799                 double all_preproc_like[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1],
 800          double all_preproc_coil_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1],
 801                 char offsets[MAX_TABLE_NUMBER +1][MAXSEQLEN],
 802                 int number_tables,
 803                 int coil_starts[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1],
 804                 int coil_ends[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1],
 805                 int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1],
 806                 int mode, int avg_or_max,
 807                 double all_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1],
 808                 double bound)
 809 {
 810   int i, j, tablenum, offset;
 811   double max_coil_score[POSNUM+1];
 812   double total_coil_score[POSNUM+1];
 813   int start_coil[POSNUM+1];
 814   int current_coil[POSNUM +1];
 815   int first_error=1;
 816
 817   int started_coil = 0;
 818   int pos_coil_number = 0;
 819
 820   extern int offset_to_use;
 821   int offset_to_use_here;
 822
 823   if (offset_to_use ==-1) offset_to_use_here = 7;
 824   else offset_to_use_here = offset_to_use;
 825 /*********************************************************************/
 826
 827   if ((mode & POS_MODE) && sequence.reg) { /** Compute pos coil info in **/
 828                             /** index number_tables +1.              **/
 829     for (i=0; i<=sequence.seqlen; i++)  {
 830       if ( (i<sequence.seqlen) && (sequence.reg[i] != '.')) {
 831         offsets[number_tables][i] = (sequence.reg[i]-i) % POSNUM;
 832         if (offsets[number_tables][i]<0)
 833            offsets[number_tables][i] += POSNUM;
 834
 835         if (!started_coil) {  /** Starting new coil **/
 836           coil_starts[number_tables][pos_coil_number][7] = i;
 837           started_coil =1;
 838         }
 839         else
 840           if ( offsets[number_tables][i] != offsets[number_tables][i-1])  {
 841                                     /** reg shift. **/
 842             coil_ends[number_tables][pos_coil_number][7] = i-1;
 843             pos_coil_number++;
 844             coil_starts[number_tables][pos_coil_number][7] = i;
 845           }
 846       }
 847       else {
 848         if (i<sequence.seqlen)
 849           offsets[number_tables][i] = -1;   /* Not in coil. **/
 850         if (started_coil) { /** Ending a coil. **/
 851           coil_ends[number_tables][pos_coil_number][7] = i-1;
 852           pos_coil_number++;
 853           started_coil =0;
 854         }
 855       }
 856     }
 857     number_coils[number_tables][7] = pos_coil_number;
 858   }   /*  Ends computation for pos file coils.  **/
 859
 860 /***********************************************************************/
 861 /* Note... offset = -1 indicates non-coiled region. **/
 862 /** Note that a "non-coiled" region is non-coiled either because it does **/
 863 /** not meet the preproc_bound, OR it does not meet the score (2-3stranded) **/
 864 /** bound if we are in a trimer_dimer distinguisher mode. **/
 865
 866   for (tablenum =0; tablenum < number_tables; tablenum++) {
 867
 868     for (offset =0; offset<=POSNUM; offset++) {  /*  Initialization */
 869       max_coil_score[offset] =0;
 870       total_coil_score[offset]=0;
 871       start_coil[offset] = -1;
 872       current_coil[offset] = 0;
 873     }
 874
 875     for(i=0; i<sequence.seqlen; i++) {
 876
 877       /** Get offset. **/
 878       if ( (all_preproc_like[tablenum][i][7] <= bound_preprocessor ) ||
 879            ( (all_scores != all_preproc_like) &&
 880                (all_scores[tablenum][i][7] <= bound) )     )/*  Not coiled. */
 881         offsets[tablenum][i] = -1;
 882       else
 883         for (j=0; j<POSNUM; j++)
 884           if (all_preproc_like[tablenum][i][j] ==
 885                     all_preproc_like[tablenum][i][7]) {
 886             offsets[tablenum][i] = (j-i) % POSNUM;
 887             if (offsets[tablenum][i]<0) offsets[tablenum][i] += POSNUM;
 888
 889           /* This is a new thing so that if an offset change is not real */
 890           /* we keep the old offset. **********/
 891             if ( (i>0) && (offsets[tablenum][i-1] != -1) &&
 892                 (offsets[tablenum][i] != offsets[tablenum][i-1]) &&
 893                 (all_preproc_like[tablenum][i]
 894                             [(i+offsets[tablenum][i-1]) %POSNUM]
 895                           == all_preproc_like[tablenum][i][7]) )
 896               offsets[tablenum][i] = offsets[tablenum][i-1];
 897           }
 898       /******************/
 899       /*** Find coils for given offset (7 is max offset).  ***/
 900       for (offset=0; offset<=POSNUM; offset++) {
 901
 902
 903         if ( (all_preproc_like[tablenum][i][offset] > bound_preprocessor) &&
 904                ( (all_scores == all_preproc_like)  ||
 905                  (all_scores[tablenum][i][offset] > bound) )   ) {
 906        /* In coil that meets both bounds, or only has to meet preproc bound. */
 907         /*** For  offset 7,  max register change so score prev coil.  ***/
 908           if (  (start_coil[7] != -1) && (offset == 7) &&
 909               (i>0) && (offsets[tablenum][i] != offsets[tablenum][i-1])) {
 910             for (j = start_coil[7]; j<i; j++)
 911               if (avg_or_max == 0)   /* Do average score. */
 912                 all_preproc_coil_scores[tablenum][j][7] =
 913                        total_coil_score[7]/(i-start_coil[7]);
 914               else                  /* Do max score over coil. */
 915                 all_preproc_coil_scores[tablenum][j][7] = max_coil_score[7];
 916
 917
 918 /****** Make note of the start and end of the coil just found. ****/
 919             if (current_coil[7] > MAXNUMCOILS) {
 920               if (first_error) {
 921                 fprintf(stderr,
 922                         "\n\nERROR: Found more than max number of coils.  Found %d coils in seq %s.",current_coil[7],sequence.title);
 923                 first_error =0;
 924               }
 925             }
 926             else {
 927               coil_starts[tablenum][current_coil[7]][7] = start_coil[7];
 928               coil_ends[tablenum][current_coil[7]][7] = i-1;
 929               number_coils[tablenum][7]++;
 930               current_coil[7]++;
 931             }
 932             max_coil_score[7]=0;
 933             total_coil_score[7]=0;
 934             start_coil[7]=-1;
 935           }
 936
 937           if (start_coil[offset] == -1)  /* Starting a coil. */
 938             start_coil[offset] =  i;
 939
 940           if (all_preproc_like[tablenum][i][offset] >max_coil_score[offset])
 941             max_coil_score[offset] = all_preproc_like[tablenum][i][offset];
 942
 943           total_coil_score[offset] += all_preproc_like[tablenum][i][offset];
 944         }
 945
 946         else     {   /** Not in a coil.  */
 947           all_preproc_coil_scores[tablenum][i][offset]=0;
 948           if (start_coil[offset] != -1)  {  /* ended a coil, so score it. */
 949             for (j = start_coil[offset]; j<i; j++)
 950               if (avg_or_max == 0)   /* Do average score. */
 951                 all_preproc_coil_scores[tablenum][j][7] =
 952                                 total_coil_score[7]/(i-start_coil[7]);
 953               else                  /* Do max score over coil. */
 954                 all_preproc_coil_scores[tablenum][j][offset] =
 955                                              max_coil_score[offset];
 956
 957 /****** Make note of the start and end of the coil just found. ****/
 958             if (current_coil[offset] > MAXNUMCOILS) {
 959               if (first_error && (offset == offset_to_use_here)) {
 960                 fprintf(stderr,
 961                         "\n\nERROR: Found more than max number of coils in offset %d.  Found %d coils in %s using table %d",current_coil[offset],sequence.title,tablenum, offset_to_use);
 962                 first_error=0;
 963               }
 964             }
 965             else {
 966               coil_starts[tablenum][current_coil[offset]][offset] =
 967                           start_coil[offset];
 968               coil_ends[tablenum][current_coil[offset]][offset] = i-1;
 969               number_coils[tablenum][offset]++;
 970               current_coil[offset]++;
 971             }
 972             max_coil_score[offset]=0;
 973             total_coil_score[offset]=0;
 974             start_coil[offset]=-1;
 975           }
 976         }
 977
 978       } /* Count on offset. */
 979     }   /* Count on i (position in sequence) */
 980
 981
 982
 983  /************** Now deal with coils ended by end of sequence.   *****/
 984     for (offset =0; offset<=POSNUM; offset++)
 985       if (start_coil[offset] != -1)  {  /* ended a coil, with end of seq. */
 986         for (j = start_coil[offset]; j<i; j++)
 987           if (avg_or_max == 0)   /* Do average score. */
 988             all_preproc_coil_scores[tablenum][j][7] =
 989                                 total_coil_score[7]/(i-start_coil[7]);
 990           else                  /* Do max score over coil. */
 991             all_preproc_coil_scores[tablenum][j][offset] =
 992                                            max_coil_score[offset];
 993
 994 /****** Make note of the start and end of the coil just found. ****/
 995         if (current_coil[7] > MAXNUMCOILS) {
 996           if (first_error) {
 997             fprintf(stderr,"\n\nERROR: Found more than max number of coils.");
 998             first_error =0;
 999           }
1000         }
1001         else {
1002           coil_starts[tablenum][current_coil[offset]][offset] =
1003                                             start_coil[offset];
1004           coil_ends[tablenum][current_coil[offset]][offset] = i-1;
1005           number_coils[tablenum][offset]++;
1006           current_coil[offset]++;
1007         }
1008       }   /* Ends if started a coil in offset.  Also ends count on offset. */
1009
1010  /***********************************************************************/
1011     for (offset= 0; offset<=POSNUM; offset++) {
1012                              /* Indicate end of coil list with -1 */
1013             coil_starts[tablenum][current_coil[offset]][offset] = -1;
1014             coil_ends[tablenum][current_coil[offset]][offset] = -1;
1015             current_coil[offset]++;
1016           }
1017
1018   }  /* Count on table. */
1019 }
1020
1021
1022
1023
1024
1025
1026
1027 int get_coil_order(int coil_starts[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1],
1028                     int coil_ends[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1],
1029                     int number_for_coil[MAX_TABLE_NUMBER+1][MAXNUMCOILS],
1030                     int number_tables,
1031                     int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1])
1032 {
1033   int current_coil_number=0;
1034   int current_coil[MAX_TABLE_NUMBER];
1035   int tablenum;
1036   int min_start_current=HUGE_VAL;
1037   int max_end_current= -HUGE_VAL;
1038   int last_end = -HUGE_VAL;
1039   int coil_table;
1040   int found_coil =1;
1041
1042
1043   for (tablenum=0; tablenum < number_tables; tablenum++)
1044     current_coil[tablenum]=0;
1045
1046 /*** Find the minimum starting unlabeled coil. For coils with same start **/
1047 /**  Take the max length one.  For a coil completely covered by the prev.**/
1048 /** coil, give it the same number.  */
1049
1050   while (found_coil)  {
1051     found_coil =0;  /*  Need to verify are still coils to continue. */
1052     min_start_current=HUGE_VAL;
1053     max_end_current= -HUGE_VAL;
1054
1055     for (tablenum=0; tablenum< number_tables; tablenum++) {
1056
1057       if (current_coil[tablenum] < number_coils[tablenum][7]) {
1058
1059         found_coil=1;
1060
1061         if (coil_starts[tablenum][current_coil[tablenum]][7] ==
1062                                      min_start_current)
1063           if (coil_ends[tablenum][current_coil[tablenum]][7] >
1064                                      max_end_current) {
1065             coil_table= tablenum;
1066             max_end_current= coil_ends[tablenum][current_coil[tablenum]][7];
1067           }
1068
1069         if (coil_starts[tablenum][current_coil[tablenum]][7]
1070                                                     < min_start_current) {
1071           min_start_current=coil_starts[tablenum][current_coil[tablenum]][7];
1072           coil_table= tablenum;
1073           max_end_current= coil_ends[tablenum][current_coil[tablenum]][7];
1074         }
1075       }
1076
1077     }     /* Ends count on tablenum */
1078
1079     if (found_coil) {
1080       if (max_end_current <= last_end) {
1081         number_for_coil[coil_table][current_coil[coil_table]] =
1082                                 current_coil_number -1;
1083
1084 /*
1085         fprintf(stderr,"\nFound coil in table %d. Coil number is %d",
1086                   coil_table, current_coil_number-1);
1087         fprintf(stderr,"\n    start= %d, end=%d", min_start_current,
1088                   max_end_current);
1089 */
1090       }
1091
1092       else {
1093         number_for_coil[coil_table][current_coil[coil_table]]=
1094                                       current_coil_number;
1095         last_end = max_end_current;
1096 /*
1097         fprintf(stderr,"\nFound coil2 in table %d. Coil number is %d",
1098                   coil_table, current_coil_number);
1099         fprintf(stderr,"\n    start= %d, end=%d", min_start_current,
1100                   max_end_current);
1101 */
1102         current_coil_number++;
1103       }
1104
1105       current_coil[coil_table]++;
1106     }   /* If found coil*/
1107   }       /* Ends while loop.  */
1108
1109
1110   return(current_coil_number -1);
1111 }
1112
1113
1114
1115  Changed:
1116 void tuple_output2(Sequence sequence, int mode, FILE *fout,
1117                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1118                    int number_tables, int number_multi_lib[MAX_TABLE_NUMBER],
1119                    double bound,
1120                    int number_score_dim)
1121 to
1122 void tuple_output2(Sequence sequence, int mode, FILE *fout,
1123                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1124                    int number_tables, double bound,int num_score_dim)
1125
1126 and cut
1127
1128   int num_score_dim;
1129   if (number_score_dim) num_score_dim = number_score_dim;
1130   else {
1131     num_score_dim =0;
1132     for (table=0; table<number_tables; table++)
1133       num_score_dim += number_multi_lib[table];
1134   }
1135
1136
1137
1138 Changed
1139
1140 void tuple_output_for_max(Sequence sequence, int mode, FILE *fout,
1141                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1142                    int number_tables, int number_multi_lib[MAX_TABLE_NUMBER],
1143                           double bound,
1144                           int number_score_dim)
1145 to
1146 void tuple_output_for_max(Sequence sequence, int mode, FILE *fout,
1147                 double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1148                    int number_tables, double bound,int num_score_dim)
1149
1150 and cut
1151   int num_score_dim;
1152
1153
1154   if (number_score_dim) num_score_dim = number_score_dim;
1155   else { /*  Old order for tuple, gaussian param. */
1156     num_score_dim =0;
1157     for (table=0; table<number_tables; table++)
1158       num_score_dim += number_multi_lib[table];
1159   }
1160
1161
1162
1163 In tuple_output2()
1164 right after
1165       if ( res_above_bound && (!posp || sequence.reg[j] != '.')  ) {
1166         fprintf(fout,"\n");
1167 cut
1168         if (!number_score_dim) {  /** Old style order.  **/
1169           for (dist=0; dist<7; dist++)
1170             for (table =0; table<number_tables; table++)
1171               if (dist < number_multi_lib[table])
1172                 fprintf(fout,"%lf ",all_scs[dist*number_tables + table][j][7]);
1173         }
1174         else
1175
1176 In tuple_output_max()
1177 right after
1178           for (reg=0; reg<7; reg++) {
1179             fprintf(fout,"\n");
1180 cut
1181             if (!number_score_dim) { /** Old order for tuple. **/
1182               for (dist=0; dist<7; dist++)
1183                 for (table =0; table<number_tables; table++)
1184                   if (dist<number_multi_lib[table])
1185                     fprintf(fout,"%lf ",
1186                           all_scs[dist*number_tables + table][j][reg]);
1187             }
1188             else
1189
1190 and after
1191         else if (sequence.reg[j] != '.') { /* Output correct reg scores. */
1192           fprintf(fout,"\n");
1193 cut
1194           if (!number_score_dim) {   /*** Old order of tuples, scores. **/
1195             for (dist=0; dist<7; dist++)
1196               for (table =0; table<number_tables; table++)
1197                 if (dist < number_multi_lib[table])
1198                   fprintf(fout,"%lf ",
1199                       all_scs[dist*number_tables + table][j][sequence.reg[j]]);
1200           }
1201           else