sources/ncoils/coils_old.c

   1
   2 /* Program COILS version 2.1 */
   3   /* written by J.Lupas, 22 JUN 1993 */
   4   /* edited on 15 JAN 1994 */
   5   /* revised and corrected 4 APR 1994 */
   6   /* incorporates the option to use either the old MTK chart or
   7      the new MTIDK chart AND the choice of weighting positions
   8      a & d more (2.5 times) */
   9   /* 4 output options:
  10         - probabilities for window sizes 14, 21 & 28 in columns,
  11         - probabilities for window sizes 14, 21 & 28 in rows,
  12         - scores only for a user-input window size, or
  13         - the probabilities above a user-input cutoff.
  14      Requests input and output files from user. */
  15   /* transferred to c++ by Larry Harvie, harvie@owl.WPI.EDU */
  16   /* adapted to C by Reinhard Schneider, SCHNEIDER@EMBL-Heidelberg.DE */
  17 #include <math.h>
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21
  22
  23   typedef FILE *text;
  24   typedef double chartstore[21][8];
  25   typedef char residuestore[10001];
  26   typedef double calcstore[4][10001];
  27   typedef int heptadstore[4][10001];
  28
  29   int window;
  30   char inputfile[80],outputfile[80];
  31   text protein,pickup;
  32   chartstore chartvalue;   /* probability values */
  33   residuestore residue;    /* stores amino acid sequence */
  34   heptadstore heptnum;     /* stores heptad position number */
  35   calcstore calcnumb;      /* stores best scores */
  36   int res_total;       /* total no. of residues in protein */
  37   char junkline[9000];     /* reads opening protein info. */
  38   char nameprot[9000];
  39   int nomore;      /* in case there's junk after last protein */
  40   short windows;
  41   double winheight;
  42   char option,chartchar,weightchar;
  43   double ad_weight;
  44   int hept_weight;
  45   calcstore tempcalc;
  46   heptadstore temphept;
  47   char fr14[10001],fr21[10001],fr28[10001];
  48   int pr14[10001],pr21[10001],pr28[10001];
  49
  50 short mod(short x,short y)
  51 {
  52    double x1=(double) x,y1=(double) y;
  53
  54         return (short) floor(fmod(x,y));
  55 }
  56 void  Init_newchart()
  57   /* store the probability values for each of the 20
  58      amino acids to occur in each of the heptad positions */
  59 {
  60 chartvalue[1][1]=2.998; chartvalue[1][2]=0.269; chartvalue[1][3]=0.367;
  61 chartvalue[1][4]=3.852; chartvalue[1][5]=0.510; chartvalue[1][6]=0.514;
  62 chartvalue[1][7]=0.562;
  63 chartvalue[2][1]=2.408; chartvalue[2][2]=0.261; chartvalue[2][3]=0.345;
  64 chartvalue[2][4]=0.931; chartvalue[2][5]=0.402; chartvalue[2][6]=0.440;
  65 chartvalue[2][7]=0.289;
  66 chartvalue[3][1]=1.525; chartvalue[3][2]=0.479; chartvalue[3][3]=0.350;
  67 chartvalue[3][4]=0.887; chartvalue[3][5]=0.286; chartvalue[3][6]=0.350;
  68 chartvalue[3][7]=0.362;
  69 chartvalue[4][1]=2.161; chartvalue[4][2]=0.605; chartvalue[4][3]=0.442;
  70 chartvalue[4][4]=1.441; chartvalue[4][5]=0.607; chartvalue[4][6]=0.457;
  71 chartvalue[4][7]=0.570;
  72 chartvalue[5][1]=0.490; chartvalue[5][2]=0.075; chartvalue[5][3]=0.391;
  73 chartvalue[5][4]=0.639; chartvalue[5][5]=0.125; chartvalue[5][6]=0.081;
  74 chartvalue[5][7]=0.038;
  75 chartvalue[6][1]=1.319; chartvalue[6][2]=0.064; chartvalue[6][3]=0.081;
  76 chartvalue[6][4]=1.526; chartvalue[6][5]=0.204; chartvalue[6][6]=0.118;
  77 chartvalue[6][7]=0.096;
  78 chartvalue[7][1]=0.084; chartvalue[7][2]=0.215; chartvalue[7][3]=0.432;
  79 chartvalue[7][4]=0.111; chartvalue[7][5]=0.153; chartvalue[7][6]=0.367;
  80 chartvalue[7][7]=0.125;
  81 chartvalue[8][1]=1.283; chartvalue[8][2]=1.364; chartvalue[8][3]=1.077;
  82 chartvalue[8][4]=2.219; chartvalue[8][5]=0.490; chartvalue[8][6]=1.265;
  83 chartvalue[8][7]=0.903;
  84 chartvalue[9][1]=1.233; chartvalue[9][2]=2.194; chartvalue[9][3]=1.817;
  85 chartvalue[9][4]=0.611; chartvalue[9][5]=2.095; chartvalue[9][6]=1.686;
  86 chartvalue[9][7]=2.027;
  87 chartvalue[10][1]=1.014; chartvalue[10][2]=1.476;
  88 chartvalue[10][3]=1.771;
  89 chartvalue[10][4]=0.114; chartvalue[10][5]=1.667;
  90 chartvalue[10][6]=2.006;
  91 chartvalue[10][7]=1.844;
  92 chartvalue[11][1]=0.590; chartvalue[11][2]=0.646;
  93 chartvalue[11][3]=0.584;
  94 chartvalue[11][4]=0.842; chartvalue[11][5]=0.307;
  95 chartvalue[11][6]=0.611;
  96 chartvalue[11][7]=0.396;
  97 chartvalue[12][1]=0.281; chartvalue[12][2]=3.351;
  98 chartvalue[12][3]=2.998;
  99 chartvalue[12][4]=0.789; chartvalue[12][5]=4.868;
 100 chartvalue[12][6]=2.735;
 101 chartvalue[12][7]=3.812;
 102 chartvalue[13][1]=0.068; chartvalue[13][2]=2.103;
 103 chartvalue[13][3]=1.646;
 104 chartvalue[13][4]=0.182; chartvalue[13][5]=0.664;
 105 chartvalue[13][6]=1.581;
 106 chartvalue[13][7]=1.401;
 107 chartvalue[14][1]=0.311; chartvalue[14][2]=2.290;
 108 chartvalue[14][3]=2.330;
 109 chartvalue[14][4]=0.811; chartvalue[14][5]=2.596;
 110 chartvalue[14][6]=2.155;
 111 chartvalue[14][7]=2.585;
 112 chartvalue[15][1]=1.231; chartvalue[15][2]=1.683;
 113 chartvalue[15][3]=2.157;
 114 chartvalue[15][4]=0.197; chartvalue[15][5]=1.653;
 115 chartvalue[15][6]=2.430;
 116 chartvalue[15][7]=2.065;
 117 chartvalue[16][1]=0.332; chartvalue[16][2]=0.753;
 118 chartvalue[16][3]=0.930;
 119 chartvalue[16][4]=0.424; chartvalue[16][5]=0.734;
 120 chartvalue[16][6]=0.801;
 121 chartvalue[16][7]=0.518;
 122 chartvalue[17][1]=0.197; chartvalue[17][2]=0.543;
 123 chartvalue[17][3]=0.647;
 124 chartvalue[17][4]=0.680; chartvalue[17][5]=0.905;
 125 chartvalue[17][6]=0.643;
 126 chartvalue[17][7]=0.808;
 127 chartvalue[18][1]=0.918; chartvalue[18][2]=0.002;
 128 chartvalue[18][3]=0.385;
 129 chartvalue[18][4]=0.440; chartvalue[18][5]=0.138;
 130 chartvalue[18][6]=0.432;
 131 chartvalue[18][7]=0.079;
 132 chartvalue[19][1]=0.066; chartvalue[19][2]=0.064;
 133 chartvalue[19][3]=0.065;
 134 chartvalue[19][4]=0.747; chartvalue[19][5]=0.006;
 135 chartvalue[19][6]=0.115;
 136 chartvalue[19][7]=0.014;
 137 chartvalue[20][1]=0.004; chartvalue[20][2]=0.108;
 138 chartvalue[20][3]=0.018;
 139 chartvalue[20][4]=0.006; chartvalue[20][5]=0.010;
 140 chartvalue[20][6]=0.004;
 141 chartvalue[20][7]=0.007;
 142 }
 143 void  Init_oldchart()
 144 {
 145   chartvalue[1][1]=3.167; chartvalue[1][2]=0.297;
 146 chartvalue[1][3]=0.398;
 147   chartvalue[1][4]=3.902; chartvalue[1][5]=0.585;
 148 chartvalue[1][6]=0.501;
 149   chartvalue[1][7]=0.483;
 150   chartvalue[2][1]=2.597; chartvalue[2][2]=0.098;
 151 chartvalue[2][3]=0.345;
 152   chartvalue[2][4]=0.894; chartvalue[2][5]=0.514;
 153 chartvalue[2][6]=0.471;
 154   chartvalue[2][7]=0.431;
 155   chartvalue[3][1]=1.665; chartvalue[3][2]=0.403;
 156 chartvalue[3][3]=0.386;
 157   chartvalue[3][4]=0.949; chartvalue[3][5]=0.211;
 158 chartvalue[3][6]=0.342;
 159   chartvalue[3][7]=0.360;
 160   chartvalue[4][1]=2.240; chartvalue[4][2]=0.37;
 161 chartvalue[4][3]=0.480;
 162   chartvalue[4][4]=1.409; chartvalue[4][5]=0.541;
 163 chartvalue[4][6]=0.772;
 164   chartvalue[4][7]=0.663;
 165   chartvalue[5][1]=0.531; chartvalue[5][2]=0.076;
 166 chartvalue[5][3]=0.403;
 167   chartvalue[5][4]=0.662; chartvalue[5][5]=0.189;
 168 chartvalue[5][6]=0.106;
 169   chartvalue[5][7]=0.013;
 170   chartvalue[6][1]=1.417; chartvalue[6][2]=0.090;
 171 chartvalue[6][3]=0.122;
 172   chartvalue[6][4]=1.659; chartvalue[6][5]=0.19;  chartvalue[6][6]=0.13;
 173   chartvalue[6][7]=0.1550;
 174   chartvalue[7][1]=0.045; chartvalue[7][2]=0.275;
 175 chartvalue[7][3]=0.578;
 176   chartvalue[7][4]=0.216; chartvalue[7][5]=0.211;
 177 chartvalue[7][6]=0.426;
 178   chartvalue[7][7]=0.156;
 179   chartvalue[8][1]=1.297; chartvalue[8][2]=1.551;
 180 chartvalue[8][3]=1.084;
 181   chartvalue[8][4]=2.612; chartvalue[8][5]=0.377;
 182 chartvalue[8][6]=1.248;
 183   chartvalue[8][7]=0.877;
 184   chartvalue[9][1]=1.375; chartvalue[9][2]=2.639;
 185 chartvalue[9][3]=1.763;
 186   chartvalue[9][4]=0.191; chartvalue[9][5]=1.815;
 187 chartvalue[9][6]=1.961;
 188   chartvalue[9][7]=2.795;
 189   chartvalue[10][1]=0.659; chartvalue[10][2]=1.163;
 190 chartvalue[10][3]=1.210;
 191   chartvalue[10][4]=0.031; chartvalue[10][5]=1.358;
 192 chartvalue[10][6]=1.937;
 193   chartvalue[10][7]=1.798;
 194   chartvalue[11][1]=0.347; chartvalue[11][2]=0.275;
 195 chartvalue[11][3]=0.679;
 196   chartvalue[11][4]=0.395; chartvalue[11][5]=0.294;
 197 chartvalue[11][6]=0.579;
 198   chartvalue[11][7]=0.213;
 199   chartvalue[12][1]=0.262; chartvalue[12][2]=3.496;
 200 chartvalue[12][3]=3.108;
 201   chartvalue[12][4]=0.998; chartvalue[12][5]=5.685;
 202 chartvalue[12][6]=2.494;
 203   chartvalue[12][7]=3.048;
 204   chartvalue[13][1]=0.03;  chartvalue[13][2]=2.352;
 205 chartvalue[13][3]=2.268;
 206   chartvalue[13][4]=0.237; chartvalue[13][5]=0.663;
 207 chartvalue[13][6]=1.62;
 208   chartvalue[13][7]=1.448;
 209   chartvalue[14][1]=0.179; chartvalue[14][2]=2.114;
 210 chartvalue[14][3]=1.778;
 211   chartvalue[14][4]=0.631; chartvalue[14][5]=2.55;
 212 chartvalue[14][6]=1.578;
 213   chartvalue[14][7]=2.526;
 214   chartvalue[15][1]=0.835; chartvalue[15][2]=1.475;
 215 chartvalue[15][3]=1.534;
 216   chartvalue[15][4]=0.039; chartvalue[15][5]=1.722;
 217 chartvalue[15][6]=2.456;
 218   chartvalue[15][7]=2.280;
 219   chartvalue[16][1]=0.382; chartvalue[16][2]=0.583;
 220 chartvalue[16][3]=1.052;
 221   chartvalue[16][4]=0.419; chartvalue[16][5]=0.525;
 222 chartvalue[16][6]=0.916;
 223   chartvalue[16][7]=0.628;
 224   chartvalue[17][1]=0.169; chartvalue[17][2]=0.702;
 225 chartvalue[17][3]=0.955;
 226   chartvalue[17][4]=0.654; chartvalue[17][5]=0.791;
 227 chartvalue[17][6]=0.843;
 228   chartvalue[17][7]=0.647;
 229   chartvalue[18][1]=0.824; chartvalue[18][2]=0.022;
 230 chartvalue[18][3]=0.308;
 231   chartvalue[18][4]=0.152; chartvalue[18][5]=0.180;
 232 chartvalue[18][6]=0.156;
 233   chartvalue[18][7]=0.044;
 234   chartvalue[19][1]=0.24;  chartvalue[19][2]=0;
 235 chartvalue[19][3]=0.00;
 236   chartvalue[19][4]=0.456; chartvalue[19][5]=0.019;
 237 chartvalue[19][6]=0.00;
 238   chartvalue[19][7]=0.00;
 239   chartvalue[20][1]=0.00;  chartvalue[20][2]=0.008; chartvalue[20][3]=0;
 240   chartvalue[20][4]=0.013; chartvalue[20][5]=0.0;   chartvalue[20][6]=0;
 241   chartvalue[20][7]=0;
 242 }
 243
 244
 245 void  Get_info()
 246 {
 247   char junk;
 248
 249   printf("COILS version 2.1\n");
 250   printf("ENTER INPUT FILE:  ");
 251   fflush(stdout);
 252   scanf("%s",inputfile);
 253   printf("-->%s\n",inputfile);
 254   protein=fopen(inputfile,"r");  /* modify for reading */
 255   if (protein==NULL) {
 256      perror("problem");
 257      printf("'%s' can not be found.  Try again after verifying the file name.\n",inputfile);
 258      exit(0);
 259   } /* endif */
 260
 261   printf("ENTER OUTPUT FILE:  ");
 262   scanf("%s",outputfile);
 263   pickup=fopen(outputfile,"w");
 264   if (pickup==NULL) {
 265      printf("The system could not open the output file\n");
 266      printf(".  Check the free space on disk\n");
 267      exit(0);
 268   } /* endif */
 269   printf("Two scoring matrices are available:\n");
 270   printf("   1 - MTK\n");
 271   printf("   2 - MTIDK\n");
 272   printf("Which matix? (1|2) <enter>? ");
 273   chartchar='\0';
 274   while ((chartchar!='1') && (chartchar!='2') && (chartchar!=(char) 13))
 275      scanf("%c",&chartchar);
 276   printf("Do you want a weight of 2.5 for positions a & d? (Y|N) <enter>? ");
 277   weightchar='\0';
 278   while ((weightchar!='n') && (weightchar!='y') && (weightchar!='Y') &&
 279 (weightchar!='N')  && (weightchar!=(char) 13))
 280      scanf("%c",&weightchar);
 281
 282   fprintf(pickup,"COILS version 2.1\n");
 283   if (chartchar=='2')
 284   {
 285     fprintf(pickup,"using MTIDK matrix.\n");
 286     Init_newchart();
 287   } else {
 288     fprintf(pickup,"using MTK matrix.\n");
 289     Init_oldchart();
 290   }
 291   if ((weightchar=='y') || (weightchar=='Y'))
 292   {
 293     fprintf(pickup,"weights: a,d=2.5 and b,c,e,f,g=1.0\n");
 294     hept_weight=10;
 295     ad_weight=2.5;
 296   } else {
 297     fprintf(pickup,"no weights\n");
 298     hept_weight=7;
 299     ad_weight=1.0;
 300   }
 301   fprintf(pickup,"Input file is %s\n",inputfile);
 302 }
 303 void  Select_option()
 304 {
 305     printf("OUTPUT OPTIONS:\n");
 306     printf("   p - probabilities in columns\n");
 307     printf("   a - probabilities in rows, abbreviated to the first digit\n");
 308     printf("   b - scores (size of the scanning window defined by the user)\n");
 309     printf("   c - only probabilities above a user-defined cutoff\n");
 310     printf("ENTER OPTION: ");
 311   while ((option!='a') && (option!='b') && (option!='c') && (option!='A')
 312 && (option!='B') && (option!='C') && (option!='p') && (option!='P')) {
 313      scanf("%c",&option);
 314   } /* endwhile */
 315     printf("\n");
 316   if ((option=='b') || (option=='B'))
 317   {
 318     printf("ENTER WINDOW SIZE:  ");
 319     scanf("%d",&window);
 320  } /* end */
 321   else window=14;
 322 }
 323
 324
 325 int checkforend()
 326 {
 327   int junklen;
 328   char *last4;
 329   int checkforend=0;
 330
 331   junkline[strlen(junkline)-1]='\0';
 332   junklen=strlen(junkline);
 333   if (junklen > 0)
 334   {
 335     if (junkline[0] == '>')
 336        checkforend=1;
 337     strcpy(nameprot,junkline);
 338   }
 339   if (junklen > 3)
 340   {
 341     last4=(char *) &(junkline[junklen-4]);
 342     if (strcmp(last4,"  ..")==0)
 343     {
 344       checkforend=1;
 345       strcpy(nameprot,junkline);
 346     }
 347   }
 348   return checkforend;
 349 }
 350 void  Skip_opening_info()
 351   /* Skip opening lines of protein to get to amino acid sequence */
 352 {
 353   int endingfound;
 354
 355    endingfound=checkforend();
 356    while (!endingfound)
 357     {
 358       if (!feof(protein))
 359       {
 360         fgets(junkline,9999,protein);
 361         endingfound=checkforend();
 362      } else {
 363         strcpy(nameprot," ");
 364         fclose(protein);
 365         protein=fopen(inputfile,"r");
 366         endingfound=1;
 367      } /* endif */
 368    } /* endwhile */
 369 }
 370 void  Store_residues()
 371   /* stores only the amino acids of the sequence, skipping numbers,
 372 blanks,
 373      and carriage returns and ending with '/' || '*' */
 374 {
 375   char tempstore;
 376   int validaa;
 377
 378   res_total=0;
 379   tempstore=fgetc(protein);
 380   while ((tempstore != '/') && (tempstore != '*') && (!feof(protein)))
 381   {
 382      switch (tempstore & 95) {
 383       case 'L':
 384       case 'I':
 385       case 'V':
 386       case 'M':
 387       case 'F':
 388       case 'Y':
 389       case 'G':
 390       case 'A':
 391       case 'K':
 392       case 'P':
 393       case 'R':
 394       case 'H':
 395       case 'E':
 396       case 'D':
 397       case 'Q':
 398       case 'N':
 399       case 'S':
 400       case 'T':
 401       case 'C':
 402       case 'W':
 403          validaa=1;
 404         break;
 405      default:
 406         validaa=0;
 407        break;
 408      } /* endswitch */
 409
 410     if (validaa)
 411     {
 412       res_total=res_total+1;
 413       residue[res_total]=tempstore;
 414     } /* endif */
 415     tempstore=fgetc(protein);
 416   } /* endwhile */
 417   if (tempstore == '/') fgets(junkline,9999,protein); /* read other '/'
 418 at end */
 419 }
 420
 421
 422 int Getx(char res)
 423 {
 424    int X;
 425
 426   switch (res & (char) 95) {
 427   case 'L':
 428      X=1;     break;
 429   case 'I':
 430      X=2;     break;
 431   case 'V':
 432      X=3;     break;
 433   case 'M':
 434      X=4;     break;
 435   case 'F':
 436      X=5;     break;
 437   case 'Y':
 438      X=6;     break;
 439   case 'G':
 440      X=7;     break;
 441   case 'A':
 442      X=8;     break;
 443   case 'K':
 444      X=9;     break;
 445   case 'R':
 446      X=10;     break;
 447   case 'H':
 448      X=11;     break;
 449   case 'E':
 450      X=12;     break;
 451   case 'D':
 452      X=13;     break;
 453   case 'Q':
 454      X=14;     break;
 455   case 'N':
 456      X=15;     break;
 457   case 'S':
 458      X=16;     break;
 459   case 'T':
 460      X=17;     break;
 461   case 'C':
 462      X=18;     break;
 463   case 'W':
 464      X=19;     break;
 465   case 'P':
 466      X=20;     break;
 467   default:
 468      X=-1;    break;
 469   } /* endswitch */
 470   return X;
 471 }
 472 void Calculate(int startwindow,int endindex)
 473   /* calculates best scores for each window frame */
 474 {
 475   int X,x,y,extras,heptad_pos,window_pos,res_pos;
 476   double root_inverse,scores,misc;
 477   int hept,index,window;
 478   double weight;
 479   int root;
 480
 481   printf("Calculating...\n");
 482   index=1;
 483   window=startwindow;
 484   while (index<=endindex)
 485   {
 486    printf(".");
 487    fflush(stdout);
 488    res_pos=0;
 489    root=(window / 7) * hept_weight;
 490    do {
 491     res_pos=res_pos+1;
 492     tempcalc[index][res_pos]=0;
 493     for (heptad_pos=0; heptad_pos<=6; heptad_pos++)  /* go through each
 494 residue in each heptad pos */
 495     {
 496       scores=1.0;
 497       hept=1;
 498       for (window_pos=0;window_pos<=window-1;window_pos++)  /* get
 499 values
 500 at all 21 positions */
 501       {
 502         switch (residue[window_pos + res_pos] & 95){
 503           case 'L':
 504              x=1;     break;
 505           case 'I':
 506              x=2;     break;
 507           case 'V':
 508              x=3;     break;
 509           case 'M':
 510              x=4;     break;
 511           case 'F':
 512              x=5;     break;
 513           case 'Y':
 514              x=6;     break;
 515           case 'G':
 516              x=7;     break;
 517           case 'A':
 518              x=8;     break;
 519           case 'K':
 520              x=9;     break;
 521           case 'R':
 522              x=10;     break;
 523           case 'H':
 524              x=11;     break;
 525           case 'E':
 526              x=12;     break;
 527           case 'D':
 528              x=13;     break;
 529           case 'Q':
 530              x=14;     break;
 531           case 'N':
 532              x=15;     break;
 533           case 'S':
 534              x=16;     break;
 535           case 'T':
 536              x=17;     break;
 537           case 'C':
 538              x=18;     break;
 539           case 'W':
 540              x=19;     break;
 541           case 'P':
 542              x=20;     break;
 543           default:
 544              x=-1;    break;
 545        } /* endswitch */
 546         y=(int) fmod((double)(window_pos + res_pos + heptad_pos),7.0);
 547         if (y==0) y=7;
 548         if (window_pos==0) hept=y;
 549         if (y==0) y=7;
 550         if ((y==1) || (y==4)) weight=ad_weight;
 551         else weight=1.0;
 552         root_inverse=1.0/(double) root;
 553         misc=pow(chartvalue[x][y],weight);
 554         scores=scores*(pow(misc,root_inverse));
 555       } /* end of window_pos loop */
 556       if (scores>tempcalc[index][res_pos])
 557       {
 558         tempcalc[index][res_pos]=scores;
 559         temphept[index][res_pos]=(int) fmod((double) (hept-1),7.0);
 560       }
 561     }  /* end of heptad_pos loop*/
 562    } while (res_pos+window != res_total+1);
 563    for (extras=1; extras<=window-1;extras++)
 564    {
 565     tempcalc[index][res_pos+extras]=tempcalc[index][res_pos];
 566     temphept[index][res_pos+extras]=(int) fmod((double)
 567 (temphept[index][res_pos]+extras),7.0);
 568    }
 569    window=window+7;
 570    index=index+1;
 571   }  /* for window sizes 14, 21 & 28 */
 572   /*maximize loop*/
 573   index=1;
 574   window=startwindow;
 575   while (index<=endindex)
 576   {
 577    res_pos=0;
 578    do {
 579     res_pos=res_pos+1;
 580     calcnumb[index][res_pos]=tempcalc[index][res_pos];
 581     heptnum[index][res_pos]=temphept[index][res_pos];
 582     window_pos=0;
 583     do {
 584       window_pos=window_pos+1;
 585       if (res_pos-window_pos<1) window_pos=window-1;
 586       else
 587        if (tempcalc[index][res_pos-window_pos]>calcnumb[index][res_pos])
 588        {
 589          calcnumb[index][res_pos]=tempcalc[index][res_pos-window_pos];
 590          heptnum[index][res_pos]=(int)
 591 fmod((temphept[index][res_pos-window_pos]+window_pos),7.0);
 592        } /* endif */
 593     } while (window_pos!=window-1); /* enddo */
 594    } while (res_pos!=res_total); /* enddo */
 595    index=index+1;
 596    window=window+7;
 597   }
 598   printf("\n");
 599 }
 600 double Calcprob(double x,double meancc,double stddevcc,double
 601 meangl,double stddevgl,double ratio_gl_cc)
 602 {
 603   double prob1,prob2,prob3,prob4;
 604
 605   prob1=(0.5) * pow(((x-meancc) / stddevcc),2.0);
 606   prob2=(0.5) * pow(((x-meangl) / stddevgl),2.0);
 607   prob3=stddevgl * exp(-prob1);
 608   prob4=ratio_gl_cc * stddevcc * exp(-prob2);
 609   return (prob3) / (prob3+prob4);
 610 }
 611 char frame(int heptnum)
 612 {
 613    switch (heptnum) {
 614    case 0:
 615       return 'a';      break;
 616    case 1:
 617       return 'b';      break;
 618    case 2:
 619       return 'c';      break;
 620    case 3:
 621       return 'd';      break;
 622    case 4:
 623       return 'e';      break;
 624    case 5:
 625       return 'f';      break;
 626    case 6:
 627       return 'g';     break;
 628    } /* endswitch */
 629    return 'x';
 630 }
 631 int i_trunc(double n)
 632 {
 633         return (int) floor(n);
 634 }
 635 void  Column_probs(double peakmin)
 636 {
 637   int res_pos;
 638   char fr14,fr21,fr28;
 639   double prob14,prob21,prob28,old14,old21,old28,xx=1,lastxx=1;
 640   int prevres,color=8,c14=9,c21=15,c28=12;
 641    char t1[]="14",t2[]="21",t3[]="28";
 642
 643   fprintf(pickup,"%s\n",nameprot);
 644   fprintf(pickup,"  Residue       Window=14             Window=21        Window=28\n");
 645   fprintf(pickup,"             Score  Probability     Score Probability  Score  Probability\n");
 646   prevres=1;
 647   for (res_pos=1;res_pos<=res_total;res_pos++)
 648   {
 649     fr14=frame(heptnum[1][res_pos]);
 650     fr21=frame(heptnum[2][res_pos]);
 651     fr28=frame(heptnum[3][res_pos]);
 652     if (chartchar=='2')
 653     {
 654      if ((weightchar=='y') || (weightchar=='Y'))
 655      {
 656       prob28=Calcprob(calcnumb[3][res_pos],1.74,0.20,0.86,0.18,30);
 657       prob21=Calcprob(calcnumb[2][res_pos],1.79,0.24,0.92,0.22,25);
 658       prob14=Calcprob(calcnumb[1][res_pos],1.89,0.30,1.04,0.27,20);
 659      } else {
 660       prob28=Calcprob(calcnumb[3][res_pos],1.69,0.18,0.80,0.18,30);
 661       prob21=Calcprob(calcnumb[2][res_pos],1.74,0.23,0.86,0.21,25);
 662       prob14=Calcprob(calcnumb[1][res_pos],1.82,0.28,0.95,0.26,20);
 663     } /* end */
 664    } /* end */
 665     else
 666     {
 667      if ((weightchar=='y') || (weightchar=='Y'))
 668      {
 669       if (calcnumb[3][res_pos] > 0.0)
 670         prob28=Calcprob(calcnumb[3][res_pos],1.70,0.24,0.79,0.23,30);
 671       else {
 672         prob28=0.0;
 673         fr28='x';
 674       } /* endif */
 675       if (calcnumb[2][res_pos] > 0.0)
 676         prob21=Calcprob(calcnumb[2][res_pos],1.76,0.28,0.86,0.26,25);
 677       else { prob21=0.0; fr21='x'; }
 678       if (calcnumb[1][res_pos] > 0.0)
 679         prob14=Calcprob(calcnumb[1][res_pos],1.88,0.34,1.00,0.33,20);
 680       else { prob14=0.0; fr14='x'; }
 681     } /* end */
 682      else
 683      {
 684       if (calcnumb[3][res_pos] > 0.0)
 685
 686 prob28=Calcprob(calcnumb[3][res_pos],1.628,0.243,0.770,0.202,30);
 687       else { prob28=0.0; fr28='x'; }
 688       if (calcnumb[2][res_pos] > 0.0)
 689
 690 prob21=Calcprob(calcnumb[2][res_pos],1.683,0.285,0.828,0.236,25);
 691       else { prob21=0.0; fr21='x'; }
 692       if (calcnumb[1][res_pos] > 0.0)
 693
 694 prob14=Calcprob(calcnumb[1][res_pos],1.782,0.328,0.936,0.289,20);
 695       else { prob14=0.0; fr14='x'; }
 696     } /* end */
 697     }
 698     if ((prob14>=peakmin) || (prob21>=peakmin) || (prob28>=peakmin))
 699     {
 700       if (res_pos-1!=prevres) fprintf(pickup,"...\n");
 701       prevres=res_pos;
 702       fprintf(pickup,"%5d %c",res_pos,residue[res_pos]);
 703       fprintf(pickup,"     %c %5.3lf    %5.3lf",fr14,calcnumb[1][res_pos],prob14);
 704       fprintf(pickup,"       %c %5.3lf    %5.3lf",fr21,calcnumb[2][res_pos],prob21);
 705       fprintf(pickup,"       %c %5.3lf    %5.3lf",fr28,calcnumb[3][res_pos],prob28);
 706       fprintf(pickup,"\n");
 707     }
 708   }
 709 }
 710 void  Row_probs()
 711 {
 712   int res_pos,maxline,startpos;
 713
 714   fprintf(pickup,"%s\n",nameprot);
 715   for (res_pos=1;res_pos<=res_total;res_pos++)
 716   {
 717     fr14[res_pos]=frame(heptnum[1][res_pos]);
 718     fr21[res_pos]=frame(heptnum[2][res_pos]);
 719     fr28[res_pos]=frame(heptnum[3][res_pos]);
 720     if (chartchar=='2')
 721     {
 722      if ((weightchar=='y') || (weightchar=='Y'))
 723      {
 724
 725 pr14[res_pos]=i_trunc(Calcprob(calcnumb[1][res_pos],1.89,0.30,1.04,0.27,20)*10.0);
 726
 727 pr21[res_pos]=i_trunc(Calcprob(calcnumb[2][res_pos],1.79,0.24,0.92,0.22,25)*10.0);
 728
 729 pr28[res_pos]=i_trunc(Calcprob(calcnumb[3][res_pos],1.74,0.20,0.86,0.18,30)*10.0);
 730     } /* end */
 731      else
 732      {
 733
 734 pr14[res_pos]=i_trunc(Calcprob(calcnumb[1][res_pos],1.82,0.28,0.95,0.26,20)*10.0);
 735
 736 pr21[res_pos]=i_trunc(Calcprob(calcnumb[2][res_pos],1.74,0.23,0.86,0.21,25)*10.0);
 737
 738 pr28[res_pos]=i_trunc(Calcprob(calcnumb[3][res_pos],1.69,0.18,0.80,0.18,30)*10.0);
 739     } /* end */
 740    } /* end */
 741     else
 742     {
 743      if ((weightchar=='y') || (weightchar=='Y'))
 744      {
 745       if (calcnumb[3][res_pos] > 0)
 746
 747 pr28[res_pos]=i_trunc(Calcprob(calcnumb[3][res_pos],1.70,0.24,0.79,0.23,30)*10.0);
 748       else { pr28[res_pos]=0; fr28[res_pos]='x'; }
 749       if (calcnumb[2][res_pos] > 0)
 750
 751 pr21[res_pos]=i_trunc(Calcprob(calcnumb[2][res_pos],1.76,0.28,0.86,0.26,25)*10.0);
 752       else { pr21[res_pos]=0; fr21[res_pos]='x'; }
 753       if (calcnumb[1][res_pos] > 0)
 754
 755 pr14[res_pos]=i_trunc(Calcprob(calcnumb[1][res_pos],1.88,0.34,1.00,0.33,20)*10.0);
 756       else { pr14[res_pos]=0; fr14[res_pos]='x'; }
 757     } /* end */
 758      else
 759      {
 760       if (calcnumb[3][res_pos] > 0)
 761
 762 pr28[res_pos]=i_trunc(Calcprob(calcnumb[3][res_pos],1.628,0.243,0.770,0.202,30)*10.0);
 763       else { pr28[res_pos]=0; fr28[res_pos]='x'; }
 764       if (calcnumb[2][res_pos] > 0)
 765
 766 pr21[res_pos]=i_trunc(Calcprob(calcnumb[2][res_pos],1.683,0.285,0.828,0.236,25)*10.0);
 767       else { pr21[res_pos]=0; fr21[res_pos]='x'; }
 768       if (calcnumb[1][res_pos] > 0)
 769
 770 pr14[res_pos]=i_trunc(Calcprob(calcnumb[1][res_pos],1.782,0.328,0.936,0.289,20)*10.0);
 771       else { pr14[res_pos]=0; fr14[res_pos]='x'; }
 772     } /* end */
 773     }
 774     if (pr14[res_pos]>9) pr14[res_pos]=9;
 775     if (pr21[res_pos]>9) pr21[res_pos]=9;
 776     if (pr28[res_pos]>9) pr28[res_pos]=9;
 777   }
 778   maxline=60;
 779   startpos=1;
 780   while (startpos+maxline-1 <= res_total)
 781   {
 782     fprintf(pickup,"%d\n",startpos);
 783     fprintf(pickup,"    .    |    .    |    .    |    .    |    .   |    .    |\n");
 784     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 785       fprintf(pickup,"%c",residue[res_pos]);
 786     fprintf(pickup,"\n\n");
 787     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 788       fprintf(pickup,"%d",pr14[res_pos]);
 789     fprintf(pickup,"\n");
 790     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 791       fprintf(pickup,"%d",pr21[res_pos]);
 792     fprintf(pickup,"\n");
 793     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 794       fprintf(pickup,"%d",pr28[res_pos]);
 795     fprintf(pickup,"\n");
 796     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 797       fprintf(pickup,"%c",fr14[res_pos]);
 798     fprintf(pickup,"\n");
 799     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 800       fprintf(pickup,"%c",fr21[res_pos]);
 801     fprintf(pickup,"\n");
 802     for (res_pos=startpos;res_pos<=startpos+maxline-1;res_pos++)
 803       fprintf(pickup,"%c",fr28[res_pos]);
 804     fprintf(pickup,"\n");
 805     fprintf(pickup,"\n");
 806     startpos=startpos+60;
 807   }
 808   /* finish up! */
 809   if (startpos<=res_total)
 810   {
 811     fprintf(pickup,"%d\n",startpos);
 812     fprintf(pickup,"    .    |    .    |    .    |    .    |    .   |    .    |\n");
 813     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 814       fprintf(pickup,"%c",residue[res_pos]);
 815     fprintf(pickup,"\n");
 816     fprintf(pickup,"\n");
 817     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 818       fprintf(pickup,"%d",pr14[res_pos]);
 819     fprintf(pickup,"\n");
 820     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 821       fprintf(pickup,"%d",pr21[res_pos]);
 822     fprintf(pickup,"\n");
 823     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 824       fprintf(pickup,"%d",pr28[res_pos]);
 825     fprintf(pickup,"\n");
 826     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 827       fprintf(pickup,"%c",fr14[res_pos]);
 828     fprintf(pickup,"\n");
 829     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 830       fprintf(pickup,"%c",fr21[res_pos]);
 831     fprintf(pickup,"\n");
 832     for (res_pos=startpos;res_pos<=res_total;res_pos++)
 833       fprintf(pickup,"%c",fr28[res_pos]);
 834     fprintf(pickup,"\n");
 835     fprintf(pickup,"\n");
 836   }
 837 }
 838
 839
 840 void  Scores_only()
 841 {
 842   int res_pos;
 843   char fr;
 844
 845   fprintf(pickup,"%s\n",nameprot);
 846   fprintf(pickup,"       Residue  Frame  Score\n");
 847   for (res_pos=1;res_pos<=res_total;res_pos++)
 848   {
 849     if (calcnumb[1][res_pos] > 0.0)
 850     {
 851       fr=frame(heptnum[1][res_pos]);
 852       fprintf(pickup,"    %6d %c      %c    %2.6lf\n",res_pos,residue[res_pos],fr,
 853           calcnumb[1][res_pos]);
 854    } else
 855       fprintf(pickup,"    %6d %c      x    %2.6lf\n",res_pos,residue[res_pos],
 856           calcnumb[1][res_pos]);
 857   }
 858 }
 859
 860 void  Calc_print()
 861 {
 862   double peakmin=0.0;
 863   int num_wind;
 864
 865   if (res_total >= window)
 866   {
 867     num_wind=3;
 868     if ((option=='a') || (option=='A'))
 869     {
 870       Calculate(window,num_wind);
 871       Row_probs();
 872    } else
 873     if ((option=='b') || (option=='B'))
 874     {
 875       num_wind=1;
 876       Calculate(window,num_wind);
 877       Scores_only();
 878    } /* end */
 879     else
 880     {
 881       if ((option=='c') || (option=='C'))
 882       {
 883         printf("ENTER MINIMUM PROBABILITY (as .##):  ");
 884         scanf("%lf",&peakmin);
 885       }
 886       Calculate(window,num_wind);
 887       Column_probs(peakmin);
 888     }
 889  } else
 890     fprintf(pickup,"Protein too short. Length=%d\n",res_total);
 891 }
 892
 893 main(int argc, char *argv[], char *envp[])
 894 {    /* Main Program */
 895   Get_info();
 896   Select_option();
 897   fgets(junkline,9999,protein);
 898
 899   do {
 900     Skip_opening_info();
 901     if (!nomore)
 902     {
 903       Store_residues();
 904       Calc_print();
 905       if (!feof(protein)) fgets(junkline,9999,protein);
 906     }
 907   } while (!feof(protein)); /* enddo */
 908   fclose(protein);
 909   fclose(pickup);
 910 }
 911