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