binaries/src/iupred/iupred.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <ctype.h>
   5 #include <math.h>
   6
   7
   8
   9 #define AA "GAVLIFPSTCMWYNQDEKRH"
  10 #define AAN 20
  11 #define MSL 40000
  12 #define ML 1000
  13
  14 #define DLC 1
  15 #define DUC 100
  16 #define DWS 10
  17 #define DMin_Ene 0.3
  18 #define DJOIN 45
  19 #define DDEL 35
  20
  21
  22 #define MAX(a,b) ((a)>(b)?(a):(b))
  23 #define MIN(a,b) ((a)<(b)?(a):(b))
  24
  25
  26 typedef struct
  27 {
  28   char name[1000];
  29   int le;
  30   char *seq;
  31   double expscore;
  32   double *eprof;
  33   double *smp;
  34   double *en;
  35   int ngr;
  36   int **gr;
  37 } SEQ_STR;
  38
  39
  40 typedef struct
  41 {
  42   double **CC;
  43   double *distro;
  44   double min, max;
  45   double step;
  46   double cutoff;
  47   int nb;
  48 } P_STR;
  49
  50
  51 void            read_mat(char *path, char *fn, double **MAT, double *matave);
  52 int             getargs(char *line,char *args[],int max);
  53 void            read_ref(char *path, char *fn, double **REF);
  54 void            Get_Histo( P_STR *P, char *path, char *);
  55 double          **DMatrix(int n_rows, int n_cols);
  56 void            *my_malloc(size_t size);
  57 void            IUPred(SEQ_STR *SEQ, P_STR *P);
  58 void            getRegions(SEQ_STR *SEQ );
  59 void            Get_Seq(char *fn, SEQ_STR *SEQ);
  60
  61
  62 int LC, UC, WS;
  63 double Min_Ene;
  64 int JOIN, DEL;
  65 int Flag_EP;
  66 double EP;
  67
  68
  69 int main(int argc, char **argv)
  70 {
  71   P_STR *P;
  72   SEQ_STR *SEQ;
  73   int i,j;
  74   int type;
  75   char *path;
  76
  77
  78
  79   if (argc!=3) {
  80     printf(" Usage: %s seqfile type \n",argv[0]);
  81     printf("      where type stands for one of the options of \n");
  82     printf("      \"long\", \"short\" or \"glob\"\n");
  83     exit(1);
  84   }
  85   if ((path=getenv("IUPred_PATH"))==NULL) {
  86     fprintf(stderr,"IUPred_PATH environment variable is not set\n");
  87     path="./";
  88   }
  89
  90   printf("# IUPred \n");
  91   printf("# Copyright (c) Zsuzsanna Dosztanyi, 2005\n");
  92   printf("#\n");
  93   printf("# Z. Dosztanyi, V. Csizmok, P. Tompa and I. Simon\n");
  94   printf("# J. Mol. Biol. (2005) 347, 827-839. \n");
  95   printf("#\n");
  96   printf("#\n");
  97
  98
  99
 100
 101   if ((strncmp(argv[2],"long",4))==0) {
 102     type=0;
 103   }
 104   else if ((strncmp(argv[2],"short",5))==0) {
 105     type=1;
 106   }
 107   else if ((strncmp(argv[2],"glob",4))==0) {
 108     type=2;
 109   }
 110   else {
 111     printf("Wrong argument\n");exit(1);
 112   }
 113
 114
 115   SEQ=malloc(sizeof(SEQ_STR));
 116   Get_Seq(argv[1],SEQ);
 117   if (SEQ->le==0) {printf(" Sequence length 0\n");exit(1);}
 118
 119 #ifdef DEBUG
 120   printf("%s %d\n%s\n",SEQ->name,SEQ->le,SEQ->seq);
 121 #endif
 122
 123
 124   P=malloc(sizeof(P_STR));
 125   P->CC= DMatrix(AAN,AAN);
 126
 127   if (type==0) {
 128
 129     LC=1;
 130     UC=100;
 131     WS=10;
 132     Flag_EP=0;
 133
 134     read_ref(path,"ss",P->CC);
 135     Get_Histo(P, path, "histo");
 136
 137
 138     IUPred(SEQ,P);
 139
 140     printf("# Prediction output \n");
 141     printf("# %s\n",SEQ->name);
 142     for (i=0;i<SEQ->le;i++)
 143       printf("%5d %c %10.4f\n",i+1,SEQ->seq[i],SEQ->en[i]);
 144   }
 145   if (type==1) {
 146     LC=1;
 147     UC=25;
 148     WS=10;
 149     Flag_EP=1;
 150     EP=-1.26;
 151
 152     read_ref(path,"ss_casp",P->CC);
 153     Get_Histo(P, path, "histo_casp");
 154
 155     IUPred(SEQ,P);
 156
 157     printf("# Prediction output \n");
 158     printf("# %s\n",SEQ->name);
 159     for (i=0;i<SEQ->le;i++)
 160       printf("%5d %c %10.4f\n",i+1,SEQ->seq[i],SEQ->en[i]);
 161   }
 162   if (type==2) {
 163     char *globseq;
 164
 165     LC=1;
 166     UC=100;
 167     WS=15;
 168     Flag_EP=0;
 169
 170
 171
 172     read_ref(path,"ss",P->CC);
 173     Get_Histo(P,path,"histo");
 174
 175     IUPred(SEQ,P);
 176
 177     Min_Ene=DMin_Ene;
 178     JOIN=DJOIN;
 179     DEL=DDEL;
 180
 181     getRegions(SEQ);
 182     globseq=malloc((SEQ->le+1)*sizeof(char));
 183     for (i=0;i<SEQ->le;i++) globseq[i]=tolower(SEQ->seq[i]);
 184
 185     printf("# Prediction output \n");
 186     printf("# %s\n",SEQ->name);
 187     printf("Number of globular domains: %5d \n",SEQ->ngr);
 188     for (i=0;i<SEQ->ngr;i++) {
 189       printf("          globular domain   %5d.    %d - %d \n",
 190              i+1,SEQ->gr[i][0]+1,SEQ->gr[i][1]+1);
 191       for (j=SEQ->gr[i][0];j<SEQ->gr[i][1]+1;j++) {
 192         globseq[j]=toupper(globseq[j]);
 193       }
 194     }
 195     printf(">%s\n",SEQ->name);
 196     for (i=0;i<SEQ->le;i++) {
 197       if ((i>0)&&(i%60==0)) printf("\n");
 198       else if ((i>0)&&(i%10==0)) printf(" ");
 199       printf("%c",globseq[i]);
 200     }
 201     printf("\n");
 202     free(globseq);
 203
 204
 205 #ifdef DEBUG
 206     for (i=0;i<SEQ->le;i++)
 207       printf("%5d %c %10.4f\n",i,SEQ->seq[i],SEQ->en[i]);
 208 #endif
 209   }
 210
 211   free(SEQ->seq);
 212   free(SEQ->eprof);free(SEQ->en);free(SEQ->smp);
 213   free(SEQ);
 214
 215
 216   return 0;
 217 }
 218
 219
 220 void IUPred(SEQ_STR *SEQ, P_STR *P)
 221 {
 222   int i,j, a1, a2,   p;
 223   int naa;
 224   double n2;
 225   double min, max, step;
 226
 227   naa=SEQ->le;
 228   min=P->min; max=P->max;step=P->step;
 229
 230   SEQ->eprof=malloc(naa*sizeof(double));
 231   for (i=0;i<naa;i++) SEQ->eprof[i]=0;
 232   SEQ->en=malloc(naa*sizeof(double));
 233   for (i=0;i<naa;i++) SEQ->en[i]=0;
 234   SEQ->smp=malloc(naa*sizeof(double));
 235   for (i=0;i<naa;i++) SEQ->smp[i]=0;
 236
 237
 238   SEQ->expscore=0;
 239   for (i=0;i<naa;i++) {
 240     a1=strchr(AA,((toupper(SEQ->seq[i]))))-AA;
 241     if ((a1<0) || (a1>=AAN)) continue;
 242     n2=0;
 243     for (j=0;j<naa;j++) if (((abs(i-j))>LC)&&((abs(i-j))<UC)) {
 244       a2=strchr(AA,((toupper(SEQ->seq[j]))))-AA;
 245       if ((a2<0) || (a2>=AAN)) continue;
 246      SEQ->eprof[i]+=P->CC[a1][a2];
 247      n2++;
 248    }
 249    SEQ->expscore+=SEQ->eprof[i]/(naa*n2);
 250    SEQ->eprof[i]/=n2;
 251
 252   }
 253
 254   if (Flag_EP==0) {
 255     for (i=0;i<naa;i++) {
 256       n2=0;
 257       for (j=MAX(0,i-WS);j<=MIN(naa,i+WS+1);j++) {
 258         SEQ->smp[i]+=SEQ->eprof[j];
 259         n2++;
 260       }
 261       SEQ->smp[i]/=n2;
 262     }
 263   } else {
 264     for (i=0;i<naa;i++) {
 265       n2=0;
 266       for (j=i-WS;j<i+WS;j++) {
 267         if ((j<0)||(j>=naa)) SEQ->smp[i]+=EP;
 268         else SEQ->smp[i]+=SEQ->eprof[j];
 269
 270         n2++;
 271       }
 272       SEQ->smp[i]/=n2;
 273     }
 274   }
 275
 276   for (i=0;i<naa;i++) {
 277
 278     if (SEQ->smp[i]<=min+2*step) SEQ->en[i]=1;
 279     if (SEQ->smp[i]>=max-2*step) SEQ->en[i]=0;
 280     if ((SEQ->smp[i]>min+2*step)&&(SEQ->smp[i]<max-2*step)) {
 281       p=(int)((SEQ->smp[i]-min)*(1.0/step));
 282       SEQ->en[i]=P->distro[p];
 283     }
 284
 285 #ifdef DEBUG
 286     printf("%5d %10.4f %10.4f %10.4f\n",
 287            i,SEQ->eprof[i], SEQ->smp[i],SEQ->en[i]);
 288
 289 #endif
 290
 291   }
 292
 293 }
 294
 295
 296
 297 void getRegions(SEQ_STR *SEQ )
 298 {
 299   int naa;
 300   int i, k,kk;
 301   int **GR, **mGR;
 302   int in_GR;
 303   int nr, mnr;
 304   int beg_GR, end_GR;
 305   int beg,end;
 306
 307
 308   naa=SEQ->le;
 309
 310
 311   GR=NULL;
 312   nr=0;
 313   in_GR=0;
 314   beg_GR=end_GR=0;
 315
 316
 317   for (i=0;i<naa;i++) {
 318
 319     if ((in_GR==1)&&(SEQ->smp[i]<=Min_Ene)) {
 320       GR=realloc(GR,(nr+1)*sizeof(int *));
 321       GR[nr]=malloc(2*sizeof(int));
 322       GR[nr][0]=beg_GR;
 323       GR[nr][1]=end_GR;
 324       in_GR=0;
 325       nr++;
 326     }
 327     else if (in_GR==1)
 328       end_GR++;
 329     if ((SEQ->smp[i]>Min_Ene)&&(in_GR==0)) {
 330       beg_GR=i;
 331       end_GR=i;
 332       in_GR=1;
 333     }
 334   }
 335   if (in_GR==1) {
 336     GR=realloc(GR,(nr+1)*sizeof(int *));
 337     GR[nr]=malloc(2*sizeof(int));
 338     GR[nr][0]=beg_GR;
 339     GR[nr][1]=end_GR;
 340     in_GR=0;
 341     nr++;
 342   }
 343
 344
 345   mnr=0;
 346   k=0;mGR=NULL;
 347   kk=k+1;
 348   if (nr>0) {
 349     beg=GR[0][0];end=GR[0][1];
 350   }
 351   while (k<nr) {
 352     if ((kk<nr)&&(GR[kk][0]-end)<JOIN) {
 353       beg=GR[k][0];
 354       end=GR[kk][1];
 355       kk++;
 356     }
 357     else if ((end-beg+1)<DEL) {
 358       k++;
 359       if (k<nr) {
 360         beg=GR[k][0];
 361         end=GR[k][1];
 362       }
 363     }
 364     else {
 365       mGR=realloc(mGR,(mnr+1)*sizeof(int*));
 366       mGR[mnr]=malloc(2*sizeof(int));
 367       mGR[mnr][0]=beg;
 368       mGR[mnr][1]=end;
 369       mnr++;
 370       k=kk;
 371       kk++;
 372       if (k<nr) {
 373         beg=GR[k][0];
 374         end=GR[k][1];
 375       }
 376     }
 377   }
 378
 379
 380   for (i=0;i<nr;i++) free(GR[i]);
 381   free(GR);
 382
 383   SEQ->ngr=mnr;
 384   SEQ->gr=mGR;
 385
 386
 387
 388 }
 389
 390
 391 void Get_Histo(P_STR *P, char *path, char *fn)
 392 {
 393   FILE *f;
 394   char ln[ML];
 395   int i,nb, set;
 396   double v, min, max, cutoff, c;
 397   char *fullfn;
 398   int sl;
 399
 400   sl=strlen(path)+strlen(fn)+2;
 401   fullfn=malloc(sl*sizeof(char));
 402   sprintf(fullfn,"%s/%s",path,fn);
 403
 404   if ((f=fopen(fullfn,"r"))==NULL) {
 405     printf("Could not open %s\n",fullfn);
 406     exit(1);
 407   }
 408
 409
 410
 411   fscanf(f,"%*s %lf %lf %d\n",&min, &max, &nb);
 412
 413   P->distro=malloc(nb*sizeof(double ));
 414   for (i=0;i<nb;i++) P->distro[i]=0;
 415
 416
 417   for (i=0,set=0;i<nb;i++) {
 418     fgets(ln,ML,f);
 419     if (feof(f)) break;
 420     if (ln[0]=='#') continue;
 421
 422     sscanf(ln,"%*s %lf %*s %*s   %lf\n", &c,&v);
 423     if ((set==0)&&(v<=0.5)) {set=1;cutoff=c;}
 424     P->distro[i]=v;
 425   }
 426
 427   fclose(f);
 428   P->max=max;
 429   P->min=min;
 430   P->nb=nb;
 431   P->cutoff=cutoff;
 432
 433
 434
 435   P->step=(max-min)/nb;
 436   P->cutoff-=P->step;
 437
 438 }
 439
 440
 441 void read_ref(char *path, char *fn, double **REF)
 442 {
 443
 444   FILE *f;
 445   int p1,p2;
 446   double v;
 447   char line[1000],s[20];
 448   int i,j;
 449   char *fullfn;
 450   int sl;
 451
 452   sl=strlen(path)+strlen(fn)+2;
 453   fullfn=malloc(sl*sizeof(char));
 454   sprintf(fullfn,"%s/%s",path,fn);
 455
 456   if ((f=fopen(fullfn,"r"))==NULL) {
 457     printf("Could not open %s\n",fullfn);
 458     exit(1);
 459   }
 460
 461
 462   while (!feof(f)) {
 463     fgets(line,1000,f);
 464
 465     sscanf(line,"%d",&p1);
 466     sscanf(&line[8],"%d",&p2);
 467     sscanf(&line[17],"%s",s);
 468     v=atof(s);
 469     REF[p1][p2]=v;
 470
 471   }
 472
 473
 474   if (REF[9][9]<0) {
 475     for (i=0;i<AAN;i++) for (j=0;j<AAN;j++)
 476       REF[i][j]*=-1;
 477   }
 478
 479
 480   fclose(f);
 481
 482 }
 483
 484
 485
 486 void read_mat(char *path, char *fn, double **MAT, double *matave)
 487 {
 488   FILE *f;
 489   char ln[ML];
 490   int numargs;
 491   char *args[AAN+2], AAL[AAN],a1,a2;
 492   int i,j,k ,p, q;
 493   double val;
 494   int sl;
 495   char *fullfn;
 496
 497   sl=strlen(path)+strlen(fn)+2;
 498   fullfn=malloc(sl*sizeof(char));
 499   sprintf(fullfn,"%s/%s",path,fn);
 500
 501   if ((f=fopen(fullfn,"r"))==NULL) {
 502     printf("Could not open %s\n",fullfn);
 503     exit(1);
 504   }
 505
 506
 507   fgets(ln,ML,f);
 508   sprintf(AAL,"%s",ln);
 509   i=0;
 510   while (!feof(f))   {
 511     fgets(ln,ML,f);
 512     k=0;j=0;
 513     if (feof(f)) break;
 514
 515
 516     if (ln[0] == '\n') continue;
 517     if (ln[0] == '#') continue;
 518     numargs = getargs(ln,args,AAN+2);
 519
 520
 521
 522     for (j=0;j<numargs;j++)  {
 523       val=atof(args[j]);
 524       a1=AAL[i];
 525       a2=AAL[j];
 526       p=strchr(AA,a1)-AA;
 527       q=strchr(AA,a2)-AA;
 528
 529       MAT[p][q]=val;
 530
 531     }
 532     i++;
 533
 534   }
 535
 536   *matave=0;
 537   for (i=0;i<AAN;i++) for (j=0;j<AAN;j++) *matave+=MAT[i][j];
 538   *matave/=(AAN*AAN);
 539
 540
 541   fclose(f);
 542
 543
 544 }
 545
 546
 547 int getargs(char *line,char *args[],int max)
 548 {
 549   char *inptr;
 550   int i;
 551
 552   inptr=line;
 553   for (i=0;i<max;i++)  {
 554     if ((args[i]=strtok(inptr," \t\n"))==NULL)
 555       break;
 556     inptr=NULL;
 557
 558   }
 559
 560   return(i);
 561
 562 }
 563
 564 void *my_malloc(size_t size)
 565 {
 566
 567   void    *new_mem;
 568
 569   if (size == 0)
 570     return NULL;
 571   new_mem = malloc(size);
 572   if (new_mem == NULL) {
 573     fprintf(stderr, "can't allocate enough memory: %d bytes\n", size);
 574   }
 575   return new_mem;
 576 }
 577
 578
 579
 580 double **DMatrix(int n_rows, int n_cols)
 581 {
 582   double **matrix;
 583   int   i,j;
 584
 585   matrix = (double **) my_malloc(n_rows*sizeof(double *));
 586   matrix[0] = (double *) my_malloc(n_rows*n_cols*sizeof(double));
 587
 588   for (i = 1; i < n_rows; i++)
 589     matrix[i] = matrix[i-1] + n_cols;
 590   for (i=0;i<n_rows;i++) for (j=0;j<n_cols;j++) matrix[i][j]=0;
 591
 592   return matrix;
 593
 594 }
 595
 596
 597 void Get_Seq(char *fn, SEQ_STR *SEQ)
 598 {
 599   char line[ML];
 600   char c=0;
 601   int j;
 602   FILE *f;
 603
 604   if ((fn==NULL)||(strlen(fn)==0)) {
 605     printf("No sequence filename\n"),exit(1);
 606   }
 607
 608   if ((f=fopen(fn,"r"))==NULL) {
 609     printf("Could not open %s\n",fn),exit(1);
 610   }
 611
 612
 613   SEQ->seq=calloc(MSL,sizeof(char));
 614
 615   fgets(line,ML,f);
 616   sscanf(&line[1],"%s",SEQ->name);
 617
 618
 619   j=0;
 620   while ((c!='>') && (!feof(f))) {
 621     c=fgetc(f);
 622     if (isalpha((int)(c))) {
 623       SEQ->seq[j]=c;
 624       j++;
 625       if ((j>0)&&(j%MSL==0)) {
 626         SEQ->seq=realloc(SEQ->seq,(j+MSL)*sizeof(char));
 627       }
 628     }
 629   }
 630   SEQ->le=j;
 631
 632 #ifdef DEBUG
 633   printf("%s %5d\n%s\n",SEQ->name,SEQ->le,SEQ->seq);
 634 #endif
 635
 636
 637 }
 638
 639