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