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