--- /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