+++ /dev/null
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <ctype.h>
-#include <math.h>
-
-#include <ctype.h>
-#include "getquery.h"
-#include "gjutil.h"
-
-
-#define AA "GAVLIFPSTCMWYNQDEKRH"
-#define AAN 20
-#define MSL 40000
-#define ML 1000
-
-#define DLC 1
-#define DUC 100
-#define DWS 10
-#define DMin_Ene 0.3
-#define DJOIN 45
-#define DDEL 35
-
-
-#define MAX(a,b) ((a)>(b)?(a):(b))
-#define MIN(a,b) ((a)<(b)?(a):(b))
-
-
-typedef struct
-{
- char *name;
- int le;
- char *seq;
- double expscore;
- double *eprof;
- double *smp;
- double *en;
- int ngr;
- int **gr;
-} SEQ_STR;
-
-
-typedef struct
-{
- double **CC;
- double *distro;
- double min, max;
- double step;
- double cutoff;
- int nb;
-} P_STR;
-
-
-void read_mat(char *path, char *fn, double **MAT, double *matave);
-int getargs(char *line,char *args[],int max);
-void read_ref(char *path, char *fn, double **REF);
-void Get_Histo( P_STR *PARAMS, char *path, char *);
-double **DMatrix(int n_rows, int n_cols);
-void *my_malloc(size_t size);
-void IUPred(SEQ_STR *SEQ, P_STR *PARAMS);
-void getRegions(SEQ_STR *SEQ );
-void Get_Seq(char *fn, SEQ_STR *SEQ);
-void createOutputFiles(int type);
-void closeOutputFiles(int type);
-
-
-int LC, UC, WS;
-double Min_Ene;
-int JOIN, DEL;
-int Flag_EP;
-double EP;
-FILE *shortout;
-FILE *longout;
-FILE *globout;
-
-
-int main(int argc, char **argv)
-{
- P_STR *PARAMS;
- SEQ_STR *SEQ;
- int i,j;
- int type;
- char *path;
-
- if (argc<2) {
- printf(" Usage: %s seqfile type \n",argv[0]);
- printf(" where type stands for one of the options of \n");
- printf(" \"long\", \"short\", \"glob\". If no type is provided the program assumes all types\n");
- exit(1);
- }
-
- if ((path=getenv("IUPred_PATH"))==NULL) {
- fprintf(stderr,"IUPred_PATH environment variable is not set\n");
- path="./";
- }
-
- printf("# IUPred \n");
- printf("# Copyright (c) Zsuzsanna Dosztanyi, 2005\n");
- printf("#\n");
- printf("# Z. Dosztanyi, V. Csizmok, P. Tompa and I. Simon\n");
- printf("# J. Mol. Biol. (2005) 347, 827-839. \n");
- printf("#\n");
- printf("# Modified to work within JABAWS (http://www.compbio.dundee.ac.uk/jabaws) framework by \n");
- printf("# Peter Troshin (pvtroshin@dundee.ac.uk) and Geoff Barton gjbarton@dundee.ac.uk\n");
- printf("# June, 2011\n");
- printf("#\n");
- printf("#\n");
-
- if(argv[2]==NULL) {
- type=3;
- }
- else if ((strncmp(argv[2],"long",4))==0) {
- type=0;
- }
- else if ((strncmp(argv[2],"short",5))==0) {
- type=1;
- }
- else if ((strncmp(argv[2],"glob",4))==0) {
- type=2;
- }
- else {
- printf("Cannot recognise disorder type assuming long\n");
- type=0;
- }
-
- FILE *fasta;
-
- if ((fasta=fopen(argv[1],"r"))==NULL) {
- printf("Could not open %s\n",argv[1]),exit(1);
- }
-
-/* Creating output files depending on the type */
- createOutputFiles(type);
-
-/* Read input file sequence by sequence */
-
- SEQS *fastaseq;
-
- do {
-
- fastaseq = gseq_fasta(fasta);
- if(fastaseq==NULL) {
- break;
- }
- SEQ=malloc(sizeof(SEQ_STR));
-
- //printf("No: %s\n",fastaseq->id );
- //printf("L: %d\n",fastaseq->slen);
-
- SEQ->name=fastaseq->id;
- SEQ->seq=fastaseq->seq;
- SEQ->le=fastaseq->slen;
-
-#ifdef DEBUG
- printf("N: %s\n",SEQ->name);
- printf("S: %s\n",SEQ->seq);
- printf("L: %d\n",SEQ->le);
-#endif
-
-
- if (SEQ->le==0) {printf(" Sequence length 0\n");exit(1);}
-
-#ifdef DEBUG
- printf("%s %d\n%s\n",SEQ->name,SEQ->le,SEQ->seq);
-#endif
-
-
- PARAMS=malloc(sizeof(P_STR));
- PARAMS->CC= DMatrix(AAN,AAN);
-
- /* LONG - DISORDER */
- if (type==0 || type==3) {
- LC=1;
- UC=100;
- WS=10;
- Flag_EP=0;
-
- read_ref(path,"ss",PARAMS->CC);
- Get_Histo(PARAMS, path, "histo");
-
- IUPred(SEQ,PARAMS);
-
- fprintf(longout, "# %s\n",SEQ->name);
-
- for (i=0;i<SEQ->le;i++)
- fprintf(longout, "%5d %c %10.4f\n",i+1,SEQ->seq[i],SEQ->en[i]);
-
-}
-/* SHORT - DISORDER */
- if (type==1 || type==3) {
- LC=1;
- UC=25;
- WS=10;
- Flag_EP=1;
- EP=-1.26;
-
- read_ref(path,"ss_casp",PARAMS->CC);
- Get_Histo(PARAMS, path, "histo_casp");
-
- IUPred(SEQ,PARAMS);
-
- fprintf(shortout, "# %s\n",SEQ->name);
- for (i=0;i<SEQ->le;i++)
- fprintf(shortout, "%5d %c %10.4f\n",i+1,SEQ->seq[i],SEQ->en[i]);
-}
-
-/* GLOB - GLOBULAR DOMAINS */
- if (type==2 || type==3) {
- char *globseq;
-
- LC=1;
- UC=100;
- WS=15;
- Flag_EP=0;
-
- read_ref(path,"ss",PARAMS->CC);
- Get_Histo(PARAMS,path,"histo");
-
- IUPred(SEQ,PARAMS);
-
- Min_Ene=DMin_Ene;
- JOIN=DJOIN;
- DEL=DDEL;
-
- getRegions(SEQ);
- globseq=malloc((SEQ->le+1)*sizeof(char));
- for (i=0;i<SEQ->le;i++) globseq[i]=tolower(SEQ->seq[i]);
-
- fprintf(globout,"# %s\n",SEQ->name);
- fprintf(globout,"Number of globular domains: %5d \n",SEQ->ngr);
- for (i=0;i<SEQ->ngr;i++) {
- fprintf(globout," globular domain %5d. %d - %d \n",
- i+1,SEQ->gr[i][0]+1,SEQ->gr[i][1]+1);
- for (j=SEQ->gr[i][0];j<SEQ->gr[i][1]+1;j++) {
- globseq[j]=toupper(globseq[j]);
- }
- }
- fprintf(globout,">%s\n",SEQ->name);
- for (i=0;i<SEQ->le;i++) {
- if ((i>0)&&(i%60==0)) fprintf(globout,"\n");
- else if ((i>0)&&(i%10==0)) fprintf(globout," ");
- fprintf(globout,"%c",globseq[i]);
- }
- fprintf(globout,"\n");
- free(globseq);
- }
-
-#ifdef DEBUG
- for (i=0;i<SEQ->le;i++)
- printf("%5d %c %10.4f\n",i,SEQ->seq[i],SEQ->en[i]);
-#endif
-
-
- free(SEQ->name);
- free(SEQ->seq);
- free(SEQ);
- free(fastaseq);
-
- } while(fastaseq != NULL);
-
- fclose(fasta);
- closeOutputFiles(type);
-
- return 0;
-}
-
-
-void closeOutputFiles(int type) {
- if (type==0) {
- fclose(longout);
- }
- else if (type==1) {
- fclose(shortout);
- }
- else if (type==2) {
- fclose(globout);
- }
- else if (type==3) {
- fclose(longout);
- fclose(shortout);
- fclose(globout);
- }
-}
-
-void createOutputFiles(int type) {
- if (type==0) {
- longout = fopen("out.long", "w");
- }
- else if (type==1) {
- shortout = fopen("out.short", "w");
- }
- else if (type==2) {
- globout = fopen("out.glob", "w");
- }
- else if (type==3) {
- longout = fopen("out.long", "w");
- shortout = fopen("out.short", "w");
- globout = fopen("out.glob", "w");
- }
-}
-
-void IUPred(SEQ_STR *SEQ, P_STR *PARAMS)
-{
- int i,j, a1, a2, p;
- int naa;
- double n2;
- double min, max, step;
-
- naa=SEQ->le;
- min=PARAMS->min; max=PARAMS->max;step=PARAMS->step;
-
- SEQ->eprof=malloc(naa*sizeof(double));
- for (i=0;i<naa;i++) SEQ->eprof[i]=0;
- SEQ->en=malloc(naa*sizeof(double));
- for (i=0;i<naa;i++) SEQ->en[i]=0;
- SEQ->smp=malloc(naa*sizeof(double));
-
- for (i=0;i<naa;i++) SEQ->smp[i]=0;
-
- SEQ->expscore=0;
- for (i=0;i<naa;i++) {
- a1=strchr(AA,((toupper(SEQ->seq[i]))))-AA;
- if ((a1<0) || (a1>=AAN)) continue;
- n2=0;
- for (j=0;j<naa;j++) if (((abs(i-j))>LC)&&((abs(i-j))<UC)) {
- a2=strchr(AA,((toupper(SEQ->seq[j]))))-AA;
- if ((a2<0) || (a2>=AAN)) continue;
- SEQ->eprof[i]+=PARAMS->CC[a1][a2];
- n2++;
- }
- SEQ->expscore+=SEQ->eprof[i]/(naa*n2);
- SEQ->eprof[i]/=n2;
-
- }
-
- if (Flag_EP==0) {
- for (i=0;i<naa;i++) {
- n2=0;
- for (j=MAX(0,i-WS);j<=MIN(naa,i+WS+1);j++) {
- SEQ->smp[i]+=SEQ->eprof[j];
- n2++;
- }
- SEQ->smp[i]/=n2;
- }
- } else {
- for (i=0;i<naa;i++) {
- n2=0;
- for (j=i-WS;j<i+WS;j++) {
- if ((j<0)||(j>=naa)) SEQ->smp[i]+=EP;
- else SEQ->smp[i]+=SEQ->eprof[j];
-
- n2++;
- }
- SEQ->smp[i]/=n2;
- }
- }
-
- for (i=0;i<naa;i++) {
-
- if (SEQ->smp[i]<=min+2*step) SEQ->en[i]=1;
- if (SEQ->smp[i]>=max-2*step) SEQ->en[i]=0;
- if ((SEQ->smp[i]>min+2*step)&&(SEQ->smp[i]<max-2*step)) {
- p=(int)((SEQ->smp[i]-min)*(1.0/step));
- SEQ->en[i]=PARAMS->distro[p];
- }
-
-#ifdef DEBUG
- printf("%5d %10.4f %10.4f %10.4f\n",
- i,SEQ->eprof[i], SEQ->smp[i],SEQ->en[i]);
-
-#endif
- }
-
-}
-
-
-
-void getRegions(SEQ_STR *SEQ )
-{
- int naa;
- int i, k,kk;
- int **GR, **mGR;
- int in_GR;
- int nr, mnr;
- int beg_GR, end_GR;
- int beg,end;
-
-
- naa=SEQ->le;
-
-
- GR=NULL;
- nr=0;
- in_GR=0;
- beg_GR=end_GR=0;
-
-
- for (i=0;i<naa;i++) {
-
- if ((in_GR==1)&&(SEQ->smp[i]<=Min_Ene)) {
- GR=realloc(GR,(nr+1)*sizeof(int *));
- GR[nr]=malloc(2*sizeof(int));
- GR[nr][0]=beg_GR;
- GR[nr][1]=end_GR;
- in_GR=0;
- nr++;
- }
- else if (in_GR==1)
- end_GR++;
- if ((SEQ->smp[i]>Min_Ene)&&(in_GR==0)) {
- beg_GR=i;
- end_GR=i;
- in_GR=1;
- }
- }
- if (in_GR==1) {
- GR=realloc(GR,(nr+1)*sizeof(int *));
- GR[nr]=malloc(2*sizeof(int));
- GR[nr][0]=beg_GR;
- GR[nr][1]=end_GR;
- in_GR=0;
- nr++;
- }
-
-
- mnr=0;
- k=0;mGR=NULL;
- kk=k+1;
- if (nr>0) {
- beg=GR[0][0];end=GR[0][1];
- }
- while (k<nr) {
- if ((kk<nr)&&(GR[kk][0]-end)<JOIN) {
- beg=GR[k][0];
- end=GR[kk][1];
- kk++;
- }
- else if ((end-beg+1)<DEL) {
- k++;
- if (k<nr) {
- beg=GR[k][0];
- end=GR[k][1];
- }
- }
- else {
- mGR=realloc(mGR,(mnr+1)*sizeof(int*));
- mGR[mnr]=malloc(2*sizeof(int));
- mGR[mnr][0]=beg;
- mGR[mnr][1]=end;
- mnr++;
- k=kk;
- kk++;
- if (k<nr) {
- beg=GR[k][0];
- end=GR[k][1];
- }
- }
- }
-
-
- for (i=0;i<nr;i++) free(GR[i]);
- free(GR);
-
- SEQ->ngr=mnr;
- SEQ->gr=mGR;
-
-
-
-}
-
-
-void Get_Histo(P_STR *PARAMS, char *path, char *fn)
-{
- FILE *f;
- char ln[ML];
- int i,nb, set;
- double v, min, max, cutoff, c;
- char *fullfn;
- int sl;
-
- sl=strlen(path)+strlen(fn)+2;
- fullfn=malloc(sl*sizeof(char));
- sprintf(fullfn,"%s/%s",path,fn);
-
- if ((f=fopen(fullfn,"r"))==NULL) {
- printf("Could not open %s\n",fullfn);
- exit(1);
- }
-
-
-
- fscanf(f,"%*s %lf %lf %d\n",&min, &max, &nb);
-
- PARAMS->distro=malloc(nb*sizeof(double ));
- for (i=0;i<nb;i++) PARAMS->distro[i]=0;
-
-
- for (i=0,set=0;i<nb;i++) {
- fgets(ln,ML,f);
- if (feof(f)) break;
- if (ln[0]=='#') continue;
-
- sscanf(ln,"%*s %lf %*s %*s %lf\n", &c,&v);
- if ((set==0)&&(v<=0.5)) {set=1;cutoff=c;}
- PARAMS->distro[i]=v;
- }
-
- fclose(f);
- PARAMS->max=max;
- PARAMS->min=min;
- PARAMS->nb=nb;
- PARAMS->cutoff=cutoff;
-
-
- PARAMS->step=(max-min)/nb;
- PARAMS->cutoff-=PARAMS->step;
-
-}
-
-
-void read_ref(char *path, char *fn, double **REF)
-{
-
- FILE *f;
- int p1,p2;
- double v;
- char line[1000],s[20];
- int i,j;
- char *fullfn;
- int sl;
-
- sl=strlen(path)+strlen(fn)+2;
- fullfn=malloc(sl*sizeof(char));
- sprintf(fullfn,"%s/%s",path,fn);
-
- if ((f=fopen(fullfn,"r"))==NULL) {
- printf("Could not open %s\n",fullfn);
- exit(1);
- }
-
-
- while (!feof(f)) {
- fgets(line,1000,f);
-
- sscanf(line,"%d",&p1);
- sscanf(&line[8],"%d",&p2);
- sscanf(&line[17],"%s",s);
- v=atof(s);
- REF[p1][p2]=v;
-
- }
-
-
- if (REF[9][9]<0) {
- for (i=0;i<AAN;i++) for (j=0;j<AAN;j++)
- REF[i][j]*=-1;
- }
-
-
- fclose(f);
-
-}
-
-
-
-void read_mat(char *path, char *fn, double **MAT, double *matave)
-{
- FILE *f;
- char ln[ML];
- int numargs;
- char *args[AAN+2], AAL[AAN],a1,a2;
- int i,j,k ,p, q;
- double val;
- int sl;
- char *fullfn;
-
- sl=strlen(path)+strlen(fn)+2;
- fullfn=malloc(sl*sizeof(char));
- sprintf(fullfn,"%s/%s",path,fn);
-
- if ((f=fopen(fullfn,"r"))==NULL) {
- printf("Could not open %s\n",fullfn);
- exit(1);
- }
-
-
- fgets(ln,ML,f);
- sprintf(AAL,"%s",ln);
- i=0;
- while (!feof(f)) {
- fgets(ln,ML,f);
- k=0;j=0;
- if (feof(f)) break;
-
-
- if (ln[0] == '\n') continue;
- if (ln[0] == '#') continue;
- numargs = getargs(ln,args,AAN+2);
-
-
-
- for (j=0;j<numargs;j++) {
- val=atof(args[j]);
- a1=AAL[i];
- a2=AAL[j];
- p=strchr(AA,a1)-AA;
- q=strchr(AA,a2)-AA;
-
- MAT[p][q]=val;
-
- }
- i++;
-
- }
-
- *matave=0;
- for (i=0;i<AAN;i++) for (j=0;j<AAN;j++) *matave+=MAT[i][j];
- *matave/=(AAN*AAN);
-
-
- fclose(f);
-
-
-}
-
-
-int getargs(char *line,char *args[],int max)
-{
- char *inptr;
- int i;
-
- inptr=line;
- for (i=0;i<max;i++) {
- if ((args[i]=strtok(inptr," \t\n"))==NULL)
- break;
- inptr=NULL;
-
- }
-
- return(i);
-
-}
-
-void *my_malloc(size_t size)
-{
-
- void *new_mem;
-
- if (size == 0)
- return NULL;
- new_mem = malloc(size);
- if (new_mem == NULL) {
- fprintf(stderr, "can't allocate enough memory: %d bytes\n", size);
- }
- return new_mem;
-}
-
-
-
-double **DMatrix(int n_rows, int n_cols)
-{
- double **matrix;
- int i,j;
-
- matrix = (double **) my_malloc(n_rows*sizeof(double *));
- matrix[0] = (double *) my_malloc(n_rows*n_cols*sizeof(double));
-
- for (i = 1; i < n_rows; i++)
- matrix[i] = matrix[i-1] + n_cols;
- for (i=0;i<n_rows;i++) for (j=0;j<n_cols;j++) matrix[i][j]=0;
-
- return matrix;
-
-}
-
-
-/*
-void Get_Seq(char *fn, SEQ_STR *SEQ)
-{
- char line[ML];
- char c=0;
- int j;
- FILE *f;
-
- if ((fn==NULL)||(strlen(fn)==0)) {
- printf("No sequence filename\n"),exit(1);
- }
-
- if ((f=fopen(fn,"r"))==NULL) {
- printf("Could not open %s\n",fn),exit(1);
- }
-
-
- SEQ->seq=calloc(MSL,sizeof(char));
-
- fgets(line,ML,f);
- sscanf(&line[1],"%s",SEQ->name);
-
-
- j=0;
- while ((c!='>') && (!feof(f))) {
- c=fgetc(f);
- if (isalpha((int)(c))) {
- SEQ->seq[j]=c;
- j++;
- if ((j>0)&&(j%MSL==0)) {
- SEQ->seq=realloc(SEQ->seq,(j+MSL)*sizeof(char));
- }
- }
- }
- SEQ->le=j;
-
-#ifdef DEBUG
- printf("%s %5d\n%s\n",SEQ->name,SEQ->le,SEQ->seq);
-#endif
-
-}
-*/
-
git://source.jalview.org / jabaws.git / blobdiff