--- /dev/null
+
+/*****************************************************************************/
+/*** (genwin.c) ***/
+/*****************************************************************************/
+
+/*--------------------------------------------------------------(includes)---*/
+
+#include "genwin.h"
+
+/*---------------------------------------------------------------(defines)---*/
+
+#define STRSIZE 100
+
+/*----------------------------------------------------------------(protos)---*/
+
+struct Sequence *readentry();
+
+/*---------------------------------------------------------------(globals)---*/
+
+char *blastdbs[] =
+ {"bba", "bbn", "embl", "gbupdate", "genbank", "genpept", "gpupdate",
+ "nr", "nrdb", "nrdb.shuf", "pir", "pseq", "swissprot", "tfdaa"};
+
+int nblastdbs = 14;
+
+#ifndef MIN
+#define MIN(a,b) ((a) <= (b) ? (a) : (b))
+#endif
+
+char *blastdir = "/net/cruncher/usr/ncbi/db/fasta/";
+char *indexdir = "/net/cruncher/usr/ncbi/db/index/";
+
+int nabets;
+struct Alphabet **abets;
+int ntvecs;
+struct TransVector **tvecs;
+int nsvecs;
+struct ScoreVector **svecs;
+int nsmats;
+struct ScoreMatrix **smats;
+
+int aaindex[128];
+unsigned char aaflag[128];
+char aachar[20];
+
+struct strlist
+ {
+ char string[STRSIZE];
+ struct strlist *next;
+ } *str, *curstr;
+
+/*---------------------------------------------------------------(tmalloc)---*/
+
+#define TESTMAX 1000
+void *tmalloc();
+int record_ptrs[TESTMAX] = {0,0,0,0};
+int rptr = 0;
+
+/*------------------------------------------------------------(genwininit)---*/
+
+genwininit()
+{
+ char *cp, *cp0;
+ int i;
+ char c;
+
+ for (i = 0; i < sizeof(aaindex)/sizeof(aaindex[0]); ++i) {
+ aaindex[i] = 20;
+ aaflag[i] = TRUE;
+ }
+
+ for (cp = cp0 = "ACDEFGHIKLMNPQRSTVWY"; (c = *cp) != '\0'; ++cp) {
+ i = cp - cp0;
+ aaindex[c] = i;
+ aaindex[tolower(c)] = i;
+ aachar[i] = tolower(c);
+ aaflag[c] = FALSE;
+ aaflag[tolower(c)] = FALSE;
+ }
+ return;
+}
+
+/*-------------------------------------------------------------(opendbase)---*/
+
+extern struct Database *opendbase(name)
+ char *name;
+
+ {struct Database *dbase;
+
+ dbase = (struct Database *) malloc(sizeof(struct Database));
+
+ if (blastdb(name))
+ {
+ dbase->filename = (char *) malloc(strlen(blastdir)+strlen(name)+1);
+ dbase->indexname = (char *) malloc(strlen(indexdir)+strlen(name)+1);
+ strcpy(dbase->filename, blastdir);
+ strcat(dbase->filename, name);
+ strcpy(dbase->indexname, indexdir);
+ strcat(dbase->indexname, name);
+ }
+ else
+ {
+ dbase->filename = (char *) malloc(strlen(name)+1);
+ dbase->indexname = (char *) malloc(strlen(name)+1);
+ strcpy(dbase->filename, name);
+ strcpy(dbase->indexname, name);
+ }
+
+ if (strcmp(dbase->filename, "-")==0)
+ {
+ dbase->fp = stdin;
+ }
+ else if ((dbase->fp=fopen(dbase->filename, "r"))==NULL)
+ {
+ free(dbase->filename);
+ free(dbase->indexname);
+ free(dbase);
+ return((struct Database *) NULL);
+ }
+
+ dbase->filepos = 0L;
+
+ return(dbase);
+ }
+
+/*---------------------------------------------------------------(blastdb)---*/
+
+int blastdb(name)
+ char *name;
+
+ {int i;
+
+ for (i=0; i<nblastdbs; i++)
+ {
+ if (strcmp(name, blastdbs[i])==0) {return(TRUE);}
+ }
+
+ return(FALSE);
+ }
+
+/*------------------------------------------------------------(closedbase)---*/
+
+extern closedbase(dbase)
+ struct Database *dbase;
+
+ {
+ fclose(dbase->fp);
+ free(dbase->filename);
+ free(dbase->indexname);
+ free(dbase);
+
+ return;
+ }
+
+/*--------------------------------------------------------------(firstseq)---*/
+
+extern struct Sequence *firstseq(dbase)
+ struct Database *dbase;
+
+ {
+ if (dbase->filepos!=0L)
+ {
+ dbase->filepos = 0L;
+ if (fseek(dbase->fp, dbase->filepos, 0)!=0)
+ {fprintf(stderr, "Error positioning file %s for firstseq.\n",
+ dbase->filename);
+ exit(1);}
+ }
+
+ return(readentry(dbase));
+ }
+
+/*---------------------------------------------------------------(nextseq)---*/
+
+extern struct Sequence *nextseq(dbase)
+ struct Database *dbase;
+
+ {
+ return(readentry(dbase));
+ }
+
+/*--------------------------------------------------------------(closeseq)---*/
+
+extern closeseq(seq)
+ struct Sequence *seq;
+
+ {
+ if (seq==NULL) return;
+
+ if (seq->id!=NULL) free(seq->id);
+ if (seq->name!=NULL) free(seq->name);
+ if (seq->organism!=NULL) free(seq->organism);
+ if (seq->header!=NULL) free(seq->header);
+ if (seq->state!=NULL) free(seq->state);
+ if (seq->composition!=NULL) free(seq->composition);
+ free(seq->seq);
+ free(seq);
+ return;
+ }
+
+/*---------------------------------------------------------------(openwin)---*/
+
+extern struct Sequence *openwin(parent, start, length)
+ struct Sequence *parent;
+ int start, length;
+
+ {struct Sequence *win;
+ int i;
+
+ if (start<0 || length<0 || start+length>parent->length)
+ {
+ return((struct Sequence *) NULL);
+ }
+
+ win = (struct Sequence *) malloc(sizeof(struct Sequence));
+
+/*--- ---[set links, up and down]---*/
+
+ win->parent = parent;
+ if (parent->root==NULL)
+ {win->root = parent;}
+ else
+ {win->root = parent->root;}
+ win->children = (struct Sequence **) NULL;
+
+/* parent->children = ***foo*** ---[not yet implemented]---*/
+
+ win->id = (char *) NULL;
+ win->name = (char *) NULL;
+ win->organism = (char *) NULL;
+ win->header = (char *) NULL;
+
+/*--- ---[install the local copy of the sequence]---*/
+
+ win->start = start;
+ win->length = length;
+#if 0
+ win->seq = (char *) malloc(sizeof(char)*length + 1);
+ memcpy(win->seq, (parent->seq)+start, length);
+ win->seq[length] = '\0';
+#else
+ win->seq = parent->seq + start;
+#endif
+
+/*--- ---[setup window implementation parameters]---*/
+
+/*--- ---[set local flags]---*/
+
+ win->rubberwin = FALSE;
+ win->floatwin = FALSE;
+ win->punctuation = FALSE;
+
+/*--- ---[initially unconfiguerd window]---*/
+
+ win->entropy = -2.;
+ win->state = (int *) NULL;
+ win->composition = (int *) NULL;
+ win->classvec = (char *) NULL;
+ win->scorevec = (double *) NULL;
+
+ stateon(win);
+
+ return win;
+}
+
+/*---------------------------------------------------------------(nextwin)---*/
+
+extern struct Sequence *nextwin(win, shift)
+ struct Sequence *win;
+ int shift;
+
+ {
+ if ((win->start+shift)<0 ||
+ (win->start+win->length+shift)>win->parent->length)
+ {
+ return((struct Sequence *) NULL);
+ }
+ else
+ {
+ return(openwin(win->parent, win->start+shift, win->length));
+ }
+ }
+
+/*--------------------------------------------------------------(shiftwin1)---*/
+static void decrementsv(), incrementsv();
+
+extern int shiftwin1(win)
+ struct Sequence *win;
+{
+ register int j, length;
+ register int *comp;
+
+ length = win->length;
+ comp = win->composition;
+
+ if ((++win->start + length) > win->parent->length) {
+ --win->start;
+ return FALSE;
+ }
+
+ if (!aaflag[j = win->seq[0]])
+ decrementsv(win->state, comp[aaindex[j]]--);
+
+ j = win->seq[length];
+ ++win->seq;
+
+ if (!aaflag[j])
+ incrementsv(win->state, comp[aaindex[j]]++);
+
+ if (win->entropy > -2.)
+ win->entropy = entropy(win->state);
+
+ return TRUE;
+}
+
+/*--------------------------------------------------------------(closewin)---*/
+
+extern closewin(win)
+ struct Sequence *win;
+
+ {
+ if (win==NULL) return;
+
+ if (win->state!=NULL) free(win->state);
+ if (win->composition!=NULL) free(win->composition);
+ if (win->classvec!=NULL) free(win->classvec);
+ if (win->scorevec!=NULL) free(win->scorevec);
+
+ free(win);
+ return;
+ }
+
+/*----------------------------------------------------------------(compon)---*/
+
+extern compon(win)
+ struct Sequence *win;
+{
+ register int *comp;
+ register int aa;
+ register char *seq, *seqmax;
+
+ win->composition = comp = (int *) calloc(20*sizeof(*comp), 1);
+ seq = win->seq;
+ seqmax = seq + win->length;
+
+ while (seq < seqmax) {
+ aa = *seq++;
+ if (!aaflag[aa])
+ comp[aaindex[aa]]++;
+ }
+
+ return;
+}
+
+/*---------------------------------------------------------------(stateon)---*/
+
+static int state_cmp(s1, s2)
+ int *s1, *s2;
+{
+ return *s2 - *s1;
+}
+
+extern stateon(win)
+ struct Sequence *win;
+{
+ register int aa, nel, c;
+
+ if (win->composition == NULL)
+ compon(win);
+
+ win->state = (int *) malloc(21*sizeof(win->state[0]));
+
+ for (aa = nel = 0; aa < 20; ++aa) {
+ if ((c = win->composition[aa]) == 0)
+ continue;
+ win->state[nel++] = c;
+ }
+ for (aa = nel; aa < 21; ++aa)
+ win->state[aa] = 0;
+
+ qsort(win->state, nel, sizeof(win->state[0]), state_cmp);
+
+ return;
+}
+
+/*-----------------------------------------------------------------(enton)---*/
+
+extern enton(win)
+ struct Sequence *win;
+
+ {
+ if (win->state==NULL) {stateon(win);}
+
+ win->entropy = entropy(win->state);
+
+ return;
+ }
+
+/*---------------------------------------------------------------(entropy)---*/
+static int thewindow;
+static double *entray;
+
+#define LN2 0.69314718055994530941723212145818
+
+void
+entropy_init(window)
+ int window;
+{
+ int i;
+ double x, xw;
+
+ entray = (double *)malloc((window+1) * sizeof(*entray));
+ xw = window;
+ for (i = 1; i <= window; ++i) {
+ x = i / xw;
+ entray[i] = -x * log(x) / LN2;
+ }
+
+ thewindow = window;
+}
+
+extern double entropy(sv)
+ register int *sv;
+{
+ int *sv0 = sv;
+ register double ent;
+ register int i, total;
+ register int *svmax;
+ register double xtotrecip, xsv;
+
+ for (total = 0; (i = *sv) != 0; ++sv)
+ total += i;
+ svmax = sv;
+ ent = 0.0;
+ if (total == thewindow) {
+ for (sv = sv0; sv < svmax; ) {
+ ent += entray[*sv++];
+ }
+ return ent;
+ }
+ if (total == 0)
+ return 0.;
+
+ xtotrecip = 1./(double)total;
+ for (sv = sv0; sv < svmax; ) {
+ xsv = *sv++;
+ ent += xsv * log(xsv * xtotrecip);
+ }
+ return -ent * xtotrecip / LN2;
+}
+
+/*-----------------------------------------------------------(decrementsv)---*/
+
+static void
+decrementsv(sv, class)
+ register int *sv;
+ register int class;
+{
+ register int svi;
+
+ while ((svi = *sv++) != 0) {
+ if (svi == class && *sv < class) {
+ sv[-1] = svi - 1;
+ break;
+ }
+ }
+}
+
+/*-----------------------------------------------------------(incrementsv)---*/
+
+static void
+incrementsv(sv, class)
+ register int *sv;
+ int class;
+{
+ for (;;) {
+ if (*sv++ == class) {
+ sv[-1]++;
+ break;
+ }
+ }
+}
+
+/*-------------------------------------------------------------(readentry)---*/
+
+struct Sequence *readentry(dbase)
+ struct Database *dbase;
+
+ {struct Sequence *seq;
+ int c;
+
+ seq = (struct Sequence *) malloc(sizeof(struct Sequence));
+
+ seq->db = dbase;
+
+/*--- ---[backpointers null at the top]---*/
+
+ seq->parent = (struct Sequence *) NULL;
+ seq->root = (struct Sequence *) NULL;
+ seq->children = (struct Sequence **) NULL;
+
+/*--- ---[set flags]---*/
+
+ seq->rubberwin = FALSE;
+ seq->floatwin = FALSE;
+
+/*--- ---[read from file]---*/
+
+ if (!readhdr(seq))
+ {
+ return((struct Sequence *) NULL);
+ }
+ while (1) /*---[skip multiple headers]---*/
+ {
+ c = getc(dbase->fp);
+ if (c == EOF)
+ break;
+ if (c != '>') {
+ ungetc(c, dbase->fp);
+ break;
+ }
+ while ((c=getc(dbase->fp)) != EOF && c !='\n')
+ ;
+ if (c == EOF)
+ break;
+ }
+ readseq(seq);
+
+/*--- ---[set implementation parameters]---*/
+
+/*--- ---[initially unconfigured]---*/
+
+ seq->entropy = -2.;
+ seq->state = (int *) NULL;
+ seq->composition = (int *) NULL;
+ seq->classvec = (char *) NULL;
+ seq->scorevec = (double *) NULL;
+
+ return(seq);
+ }
+
+/*---------------------------------------------------------------(readhdr)---*/
+
+readhdr(seq)
+ struct Sequence *seq;
+
+ {FILE *fp;
+ char *bptr, *curpos;
+ int c, i, itotal;
+ int idend, namend, orgend;
+
+ fp = seq->db->fp;
+
+ if ((c=getc(fp)) == EOF)
+ {
+ free(seq);
+ return(FALSE);
+ }
+
+ while (c != EOF && isspace(c))
+ {
+ c = getc(fp);
+ }
+
+ if (c!='>')
+ {fprintf(stderr, "Error reading fasta format - '>' not found.\n");
+ exit(1);}
+ ungetc(c, fp);
+/* ---[read the header line]---*/
+ str = (struct strlist *) malloc (sizeof(struct strlist));
+ str->next = NULL;
+ curstr = str;
+
+ for (i=0,itotal=0,c=getc(fp); c != EOF; c=getc(fp))
+ {
+ if (c=='\n') break;
+
+ if (i==STRSIZE-1)
+ {curstr->string[i] = '\0';
+ curstr->next = (struct strlist *) malloc (sizeof(struct strlist));
+ curstr = curstr->next;
+ curstr->next = NULL;
+ i = 0;}
+
+ curstr->string[i] = c;
+ itotal++;
+ i++;
+ }
+
+ curstr->string[i] = '\0';
+ seq->header = (char *) malloc (itotal+2);
+ seq->header[0] = '\0';
+
+ for (curstr=str, curpos=seq->header; curstr!=NULL;)
+ {
+ if (curstr->next==NULL)
+ {memccpy(curpos, curstr->string, '\0', STRSIZE);}
+ else
+ {memccpy(curpos, curstr->string, '\0', STRSIZE-1);}
+
+ str = curstr;
+ curstr = curstr->next;
+ free (str);
+
+ if (curstr!=NULL) {curpos = curpos+STRSIZE-1;}
+ }
+
+ bptr = (seq->header)+1;
+ seq->name = (char *) NULL;
+ seq->organism = (char *) NULL;
+/* ---[parse out the id]---*/
+ idend = findchar(bptr, ' ');
+ if (idend==-1) {idend = findchar(bptr, '\n');}
+ if (idend==-1) {idend = findchar(bptr, '\0');}
+ if (idend==-1)
+ {fprintf(stderr, "Error parsing header line - id.\n");
+ fputs(seq->header, fp);
+ exit(1);}
+
+ seq->id = (char *) malloc((idend+1)*sizeof(char));
+ memcpy(seq->id, bptr, idend);
+ seq->id[idend] = '\0';
+
+ if (bptr[idend]=='\n' || bptr[idend]=='\0') {return(TRUE);}
+
+/* ---[parse out the protein name]---*/
+ bptr = bptr + idend + 1;
+ while (bptr[0]==' ') {bptr++;}
+
+ namend = findchar(bptr, '-');
+ if (namend==-1) {namend = findchar(bptr, '\n');}
+ if (namend==-1) {namend = findchar(bptr, '\0');}
+ if (namend==-1)
+ {fprintf(stderr, "Error parsing header line - name.\n");
+ fputs(seq->header, fp);
+ return(TRUE);}
+
+ seq->name = (char *) malloc((namend+1)*sizeof(char));
+ memcpy(seq->name, bptr, namend);
+ seq->name[namend] = '\0';
+
+ if (bptr[namend]=='\n' || bptr[namend]=='\0') {return(TRUE);}
+
+/* ---[parse out organism]---*/
+ bptr = bptr + namend + 1;
+ while (bptr[0]==' ') {bptr++;}
+
+ orgend = findchar(bptr, '|');
+ if (orgend==-1) {orgend = findchar(bptr, '#');}
+ if (orgend==-1) {orgend = findchar(bptr, '\n');}
+ if (orgend==-1) {orgend = findchar(bptr, '\0');}
+ if (orgend==-1)
+ {fprintf(stderr, "Error parsing header line - organism.\n");
+ fputs(seq->header, fp);
+ return(TRUE);}
+
+ seq->organism = (char *) malloc((orgend+1)*sizeof(char));
+ memcpy(seq->organism, bptr, orgend);
+ seq->organism[orgend] = '\0';
+
+/* ---[skip over multiple header lines]---*/
+ while (TRUE)
+ {
+ c = getc(fp);
+ if (c == EOF)
+ return(TRUE);
+ if (c=='>')
+ {
+ skipline(fp);
+ }
+ else
+ {
+ ungetc(c,fp);
+ break;
+ }
+ }
+
+ return(TRUE);
+ }
+
+/*--------------------------------------------------------------(skipline)---*/
+
+skipline(fp)
+ FILE *fp;
+
+ {int c;
+
+ while ((c=getc(fp))!='\n' && c!=EOF)
+ ;
+
+ return;
+ }
+
+/*--------------------------------------------------------------(findchar)---*/
+
+extern int findchar(str, chr)
+ char *str;
+ char chr;
+
+ {int i;
+
+ for (i=0; ; i++)
+ {
+ if (str[i]==chr)
+ {
+ return(i);
+ }
+ if (str[i]=='\0')
+ {
+ return(-1);
+ }
+ }
+ }
+
+/*---------------------------------------------------------------(readseq)---*/
+
+readseq(seq)
+ struct Sequence *seq;
+
+{FILE *fp;
+ int i, itotal;
+ int c;
+ char *curpos;
+
+ fp = seq->db->fp;
+
+ seq->punctuation = FALSE;
+
+ str = (struct strlist *) malloc (sizeof(struct strlist));
+ str->next = NULL;
+ curstr = str;
+
+ for (i = 0, itotal = 0, c = getc(fp); c != EOF; c = getc(fp)) {
+ if (!aaflag[c]) {
+Keep:
+ if (i < STRSIZE-1) {
+ curstr->string[i++] = c;
+ continue;
+ }
+ itotal += STRSIZE-1;
+ curstr->string[STRSIZE-1] = '\0';
+ curstr->next = (struct strlist *) malloc(sizeof(*curstr));
+ curstr = curstr->next;
+ curstr->next = NULL;
+ curstr->string[0] = c;
+ i = 1;
+ continue;
+ }
+
+ switch (c) {
+ case '>':
+ ungetc(c, fp);
+ goto EndLoop;
+ case '*': case '-':
+ seq->punctuation = TRUE;
+ goto Keep;
+ case 'b': case 'B':
+ case 'u': case 'U': /* selenocysteine */
+ case 'x': case 'X':
+ case 'z': case 'Z':
+ goto Keep;
+ default:
+ continue;
+ }
+ }
+EndLoop:
+ itotal += i;
+
+ curstr->string[i] = '\0';
+ seq->seq = (char *) malloc (itotal+2);
+ seq->seq[0] = '\0';
+
+ for (curstr = str, curpos = seq->seq; curstr != NULL;) {
+ if (curstr->next == NULL)
+ memccpy(curpos, curstr->string, '\0', STRSIZE);
+ else
+ memccpy(curpos, curstr->string, '\0', STRSIZE-1);
+
+ str = curstr;
+ curstr = curstr->next;
+ free(str);
+
+ if (curstr != NULL)
+ curpos = curpos+STRSIZE-1;
+ }
+
+ seq->length = strlen(seq->seq);
+
+ return;
+}
+/*-----------------------------------------------------------------(upper)---*/
+
+extern upper(string, len)
+ register char *string;
+ size_t len;
+{
+ register char *stringmax, c;
+
+ for (stringmax = string + len; string < stringmax; ++string)
+ if (islower(c = *string))
+ *string = toupper(c);
+}
+
+/*-----------------------------------------------------------------(lower)---*/
+
+extern lower(string, len)
+ char *string;
+ size_t len;
+{
+ register char *stringmax, c;
+
+ for (stringmax = string + len; string < stringmax; ++string)
+ if (isupper(c = *string))
+ *string = tolower(c);
+}
+
+/*-------------------------------------------------------------------(min)---*/
+
+int min(a, b)
+ int a, b;
+
+ {
+ if (a<b) {return(a);}
+ else {return(b);}
+ }
+
+/*-------------------------------------------------------------------(max)---*/
+
+int max(a, b)
+ int a, b;
+
+ {
+ if (a<b) {return(b);}
+ else {return(a);}
+ }
+
+/*---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------(tmalloc)---*/
+
+void *tmalloc(size)
+ size_t size;
+
+ {void *ptr;
+
+ ptr = (void *) malloc(size);
+
+ if (rptr>TESTMAX)
+ {
+ exit(2);
+ }
+
+ record_ptrs[rptr] = (int) ptr;
+ rptr++;
+
+ return(ptr);
+ }
+
+/*-----------------------------------------------------------------(tfree)---*/
+
+tfree(ptr)
+ void *ptr;
+
+ {int i;
+
+ for (i=0; i<rptr; i++)
+ {
+ if (record_ptrs[i]==(int)ptr)
+ {
+ record_ptrs[i] = 0;
+ break;
+ }
+ }
+
+ free(ptr);
+ }
+
+/*---------------------------------------------------------------------------*/
git://source.jalview.org / jpred.git / blobdiff