/******************************************************************************/
/*
   FoldIndex 
   Prilusky, J., et al. (2005) Bioinformatics, 21, 3435-3438.
                                                                              */
/******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define MAXSEQ  10000

/******************************************************************************/
/******************************************************************************/
static char amino[21] = "ARNDCQEGHLIKMFPSTWYV";
double hydropathy[20] = {0.7000,0.0000,0.1111,0.1111,0.7778,
                         0.1111,0.1111,0.4556,0.1444,0.9222,
                         1.0000,0.0667,0.7111,0.8111,0.3222,
                         0.4111,0.4222,0.4000,0.3556,0.9667};
int charge[20] = {+0,+1,+0,-1,+0,+0,-1,+0,+0,+0,
                  +0,+1,+0,+0,+0,+0,+0,+0,+0,+0};

/******************************************************************************/
/******************************************************************************/
FILE *openfile(filename)
char *filename;
{
	FILE *fp;

        if ((fp = fopen(filename,"r")) == NULL) {
             fprintf(stderr,"Cannot open %s\n",filename);
             exit (-1);
        }
	return(fp);
}

/******************************************************************************/
/******************************************************************************/
int readseq(fp,seq)
FILE *fp;
char seq[];
{
	int c;
	int i;

	i = 0;
	while((c = getc(fp)) != EOF) {
	     if (c == ' ') continue;
	     if (c == '\n') continue;

	     seq[i++] = c;
	     if (i  >  MAXSEQ) {
	          fprintf(stderr,"Increase MAXSEQ\n");
		  exit (-1);
	     }
	}
	return(i);
}

/******************************************************************************/
/******************************************************************************/
double foldindex(len, seq)
int len;
char seq[MAXSEQ];
{
	int 	    i, j, k;
	double     htot, postot, negtot, ctot;
	double     y;

	htot = postot = negtot = 0.0;
	for (i = 0; i < len; i++){
	     for (j = 0; j < 20; j++)
	          if (seq[i] == amino[j]){
	               htot += hydropathy[j];
	               if (charge[j] > 0)
	                    postot++;
	               else if (charge[j] < 0)
	                    negtot++;
	               break;
	          }
	}

	htot /= (double)len;
	postot /= (double)len;
	negtot /= (double)len;
	     
	ctot = fabs(postot - negtot);
	y = 2.785*htot - ctot - 1.151;

	return (y);
}

/******************************************************************************/
/******************************************************************************/
int main(argc, argv)
int argc;
char *argv[];
{
	int        len;
	char	   seq[MAXSEQ];
	double     fi;
	FILE      *fp;
	FILE      *openfile();
	int        readseq();
	double     foldindex();	

	if (argv[1] == NULL){
	     fprintf(stderr,"%s   (seq)\n", argv[0]);
	     exit (-1);
	}

	fp = openfile(argv[1]);
	len = readseq(fp,seq);
	fclose(fp);

	fi = foldindex(len,seq);

	printf("%s\n",seq);
	printf("foldindex = %f ",fi);

	if (fi < 0)
	     printf(" < 0 : disordered\n");
	else
	     printf(" > 0 : ordered\n");
	
	return 0;
}
/* EOF */