/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "encode.p" */ #include /* encode a book of sequences into strings of integers. by Gary Stormo. modified by: Dr. Thomas D. Schneider National Institutes of Health National Cancer Institute Center for Cancer Research Nanobiology Program Molecular Information Theory Group Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu (use only if first address fails) http://www.ccrnp.ncifcrf.gov/~toms/ modules needed: delmods, delman, matmods */ /* the end of the program */ /* begin module version */ #define version 1.42 /* of encode.p 2007 Jun 22 2007 Jun 22, 1.42: handle mutated names; load from delmod 2002 Nov 1, 1.41: fix output message 2000 Nov 20, 1.40: upgrade documentation and to gpc 2000 Nov 17, 1.38: upgrade to new delmod for reading name inst. 1999 Jun 11, 1.37: upgrade to theline in alignedlist stuff 1995 Jan 29: empty inst implies book alignment origin july 7, 1983 */ /* end module version */ /* begin module describe.encode */ /* name encode: encodes a book of sequences into strings of integers synopsis encode(inst: in, book: in, encseq: out, encodep: in, output: out) files inst: the instructions generating the book; for aligning the sequences If the inst file is empty, then the sequences are aligned by the zero coordinate of the book (this allows the use of the "default coordinate zero" option of Delila) or by the first base of the piece, as defined by the first parameter. book: the sequences to be encoded encseq: the encoded sequences encodep: parameter file for describing how the sequences are to be encoded. The first parameter, the first character on the first line, defines how to align the pieces. See the alist program for the detailed logic. There are three choices, as in alist: 'f' (for 'first') then the sequences are always aligned by their first base. 'i' then the sequences are aligned by the delila instructions. If the inst file is empty, alignment is forced to the 'b' mode. 'b' (for 'internal') then the alignment is on the internal zero of the book's sequence. This option is to be used when "default coordinate zero" is used in the Delila instructions. The remaining parameters are stored as a list of parameter records, of which there may be any number. Each parameter record has five lines of information which it must include (all i's and j's are integers): 1. i j specify the nucleotides, relative to the aligned base, over which this parameter record is to operate; these may be any integers, but i <= j is required; 2. i is the size of the windows to be encoded; within the window the number of each oligonucleotide of length 'coding' are determined and printed as part of the total sequence vector; 3. i is the shift to the next window to be encoded; 4. i : j1 j2 j3 ... is the 'coding'-level and arrangement; the 'coding'-level, i, is the number of nucleotides in the oligos we are counting, i.e., 1 means monos, 2 means dis, ...; if i > 1 then we can also skip bases between the ones we are encoding; if the i is followed next by a colon, there must be i-1 integers (j1..j(i-1)) which specify the number of bases to be skipped between the ones which are encoded; for example, if we have the sequence xyz and we are interested in the di-nucleotides we can get the xy by the parameter '2 : 0', or we could get the xz by parameter '2 : 1'; if there is no colon all the skips are assumed to be zero; 5. i is the shift to the next coding site within the window; this allows us to encode only some of the oligos within a window, such as only those that are in-frame; multiple parameter records can be concatenated in the encodep file and then each sequence in the book will be encoded according to each parameter record into a single vector of integers. output: for messages to the user description This program is used to encode a book of sequences into a string of integers. Each sequence in the book is encoded into a single string of integers (ended by an 'end of sequence' symbol) according to the user specified parameters, which are in the file 'encodep'. examples documentation @article{Schneider1984, author = "T. D. Schneider and G. D. Stormo and M. A. Yarus and L. Gold", title = "Delila system tools", journal = "Nucl. Acids Res.", volume = "12", pages = "129-140", year = "1984"} see also {Example parameter file:} encodep {delman.use.encode: } http://www.lecb.ncifcrf.gov/~toms/delman1.html#delman.use.encode.1 {delman.use.aligned.books:} http://www.lecb.ncifcrf.gov/~toms/delman1.html#delman.use.aligned.books {Before using encode, one should always check the sequences by looking at them as an aligned list with} alist.p {The output of this program is used by} rseq.p author Gary Stormo bugs none known technical notes */ /* end module describe.encode */ /* begin module book.const */ /* constants needed for book manipulations */ #define dnamax 1024 /* length of dna arrays */ #define namelength 100 /* maximum key name length */ #define linelength 200 /* maximum line readable in book */ /* end module book.const version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module vector.const */ #define maxvecpart 64 /* the number of elements in a 'part' of a vector */ #define vpagewidth 64 /* the width (in characters) of the output vector file from procedure writevector */ /* end module vector.const version = 'matmod 1.95 85 apr 18 tds/gds'; */ #define pagewidth 32 /* number of encoding elements written per line */ #define codingmax 6 /* the maximum coding-level allowed */ /* begin module interact.const */ /* begin module string.const */ #define maxstring 2000 /* the maximum string */ /* end module string.const version = 4.86; (@ of prgmod.p 2004 Sep 8 */ /* end module interact.const version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module filler.const */ #define fillermax 50 /* the size of the filler array for a string */ /* end module filler.const version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module interact.type */ /* begin module string.type */ /* pointer to a string */ typedef struct string { /* a string of characters */ Char letters[maxstring]; /* the letters in the string */ long length; /* the number of characters in the string */ long current; /* the letter we are working on */ Char *next; /* the next string in a series */ } string; /* end module string.type version = 4.86; (@ of prgmod.p 2004 Sep 8 */ /* end module interact.type version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module filler.type */ /* the following is an array used to fill a string. it is convenient to have it much shorter than the maxstring, so that it is easy to fill the string using procedure fillstring. the user must declare the value of constant fillermax. */ typedef Char filler[fillermax]; /* end module filler.type version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module trigger.type */ typedef struct trigger { /* an object to be searched for */ string seek; /* the characters looked for */ long state; /* how close to triggering we are */ boolean skip; /* trigger not found- skip the line */ /* the trigger was found */ boolean found; } trigger; /* end module trigger.type version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module book.type */ /* types needed for book manipulations */ typedef long chset[5]; /* types defined in book definition */ typedef Char alpha[namelength]; /* this is not alfa */ /* name is a left justified string with blanks following the characters */ typedef struct name { alpha letters; /* zero means an unspecified structure */ char length; } name; typedef struct line { /* a line of characters */ Char letters[linelength]; uchar length; struct line *next; } line; typedef enum { plus, minus, dircomplement, dirhomologous } direction; typedef enum { linear, circular } configuration; typedef enum { on, off } state; typedef struct header { /* header of key */ name keynam; /* key name of structure */ line *fulnam; /* full name of structure */ /* note key */ line *note; } header; /* begin module base.type */ /* define the four nucleotide bases */ typedef enum { a, c, g, t } base; /* end module base.type version = 7.72; {of delmod.p 2007 Jul 23} */ /* sequence types */ typedef short dnarange; /* p2c: encode.p, line 247: * Note: Field width for seq assumes enum base has 4 elements [105] */ typedef uchar seq[(dnamax + 3) / 4]; typedef struct dnastring { seq part; dnarange length; struct dnastring *next; } dnastring; typedef struct orgkey { /* organism key */ header hea; /* genetic map units */ line *mapunit; } orgkey; typedef struct chrkey { /* chromosome key */ header hea; double mapbeg; /* number of genetic map beginning */ /* number of genetic map ending */ double mapend; } chrkey; typedef struct piekey { /* piece key */ header hea; double mapbeg; /* genetic map beginning */ configuration coocon; /* configruation (circular/linear) */ direction coodir; /* direction (+/-) relative to genetic map */ long coobeg; /* beginning nucleotide */ long cooend; /* ending nucleotide */ configuration piecon; /* configruation (circular/linear) */ direction piedir; /* direction (+/-) relative to coordinates */ long piebeg; /* beginning nucleotide */ long pieend; /* ending nucleotide */ } piekey; typedef struct piece { piekey key; dnastring *dna; } piece; typedef struct reference { name pienam; /* name of piece referred to */ double mapbeg; /* genetic map beginning */ direction refdir; /* direction relative to coordinates */ long refbeg; /* beginning nucleotide */ long refend; /* ending nucleotide */ } reference; typedef struct genkey { /* gene key */ header hea; reference ref; } genkey; typedef struct trakey { /* transcript key */ header hea; reference ref; } trakey; typedef struct markey { /* marker key */ header hea; reference ref; state sta; line *phenotype; struct marker *next; } markey; typedef struct marker { markey key; dnastring *dna; } marker; /* end module book.type version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module encode.type.param */ /* the following types allow the user to specify parameters which will be used to encode the sequences in the book. */ /* spaces are allowed between the encoded bases, and the number of bases to be skipped between each encoded pair of bases are kept in this linked list of integers */ typedef struct spacelist { long skips; /* bases skipped to next coded base */ struct spacelist *next; /* points to next spacing number */ } spacelist; /* the encoding parameters for each region are stored in these records, and the records for each region are connected into a linked list of all the encoding parameters for the entire sequences */ typedef enum { start, stop } endpoints; /* end points of a coding region */ typedef struct parameter { /* these are the instructions for doing the coding */ long range[2]; /* the bases to be coded by these instructions, relative to the alignedbase */ long window; /* the number of bases included in a coding vector */ long wshift; /* movement of the window to the next site */ long coding; /* the 'level' at which the coding is being done, i.e., 1: mononucleotides; 2: dinucleotides; ... */ spacelist *spaces; /* the spacing between the encoded bases */ long cshift; /* shift to the next coding unit in the window */ /* values calculated at read-in time */ long wvlength; /* length of a window vector. this is coding raised to the 4th power */ long pvlength; /* parameter vector length. the vector made up of a series of window vectors. */ /* note: pvlength/wvlength is the number of windows in the parameter range */ struct parameter *next; /* set of instructions for coding another region */ } parameter; /* end module encode.type.param version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.type */ /* vectors are useful for many applications, including the encoding of sequences. this vector type is designed to be flexible enough to be used whenever one needs a vector. it is designed as a linked list so it can contain as many elements as are ever needed. the actual 'vector' is a record containing the length and a pointer to the first 'vectorpart'. that vectorpart is a record containing an array of the first 'maxvecpart' elements (maxvecpart is a constant, from module vector.const, which must be included) and a pointer to the next 'vectorpart'. the elements are of type real so that either integers or reals may be used. */ typedef struct vectorpart { double numbers[maxvecpart]; struct vectorpart *next; } vectorpart; typedef struct vector { long length; vectorpart *part; } vector; /* end module vector.type version = 'matmod 1.95 85 apr 18 tds/gds'; */ Static _TEXT inst; /* the instructions that generated the book, for alignment of the sequences */ Static _TEXT book; /* the book of sequences */ Static _TEXT encseq; /* the encoded sequences, which are arrays of integers */ Static _TEXT encodep; /* the file of parameters which specifies how to encode the sequences */ Static long regions; /* the number of independently coded regions */ Static long length_; /* the length of the sequence pointed to */ Static long alignedbase; /* base in each sequence to be aligned by */ Static piece *apiece; /* a pointer the the piece currently being encoded */ Static parameter *firstparam; /* the beginning of the list of parameters */ Static long fourpowers[codingmax + 1]; /* for storing powers of 4 */ Static vector seqvector; /* the encoded sequence vector */ Static boolean noinst; /* true if the inst file is empty */ Static long theline; /* location in book */ Static Char alignmenttype; /* 'f' means alignment by First internal coordinate base, 'b' means alignment by Book, 'i' means alignment by Instructions */ /* begin module book.var */ /* ************************************************************************ */ /* global variables needed for book manipulations */ /* free storage: */ Static line *freeline; /* unused lines */ Static dnastring *freedna; /* unused dnas */ Static boolean readnumber; /* whether to read a number from the notes, or to read in the notes */ Static long number; /* the number of the item just read */ Static boolean numbered; /* true when the item just read is numbered */ Static boolean skipunnum; Static jmp_buf _JL1; /* a control variable to allow skipping of un-numbered items in the book */ /* ************************************************************************ */ /* end module book.var version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module package.primitive */ /* ************************************************************************ */ /* begin module halt */ Static Void halt() { /* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. */ printf(" program halt.\n"); longjmp(_JL1, 1); } /* end module halt version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module copyaline */ Static Void copyaline(fin, fout) _TEXT *fin, *fout; { /* copy a line from file fin to file fout */ while (!P_eoln(fin->f)) { putc(P_peek(fin->f), fout->f); getc(fin->f); } fscanf(fin->f, "%*[^\n]"); getc(fin->f); putc('\n', fout->f); } /* copyaline */ /* end module copyaline version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module copylines */ Static long copylines(fin, fout, n) _TEXT *fin, *fout; long n; { /* copy n lines of file fin to file fout. the actual number of lines copied is returned. */ long index = 0; /* the current line number */ while (!BUFEOF(fin->f) && index < n) { copyaline(fin, fout); index++; } return index; } /* copylines */ /* end module copylines version = 7.72; {of delmod.p 2007 Jul 23} */ /* ************************************************************************ */ /* end module package.primitive version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module vector.functions */ /* *********************************************************************** */ /* begin module vector.io */ /* *********************************************************************** */ /* begin module vector.primitives */ /* these functions and procedures do manipulations on vectors. */ Static double vget(v, pos) vector v; long pos; { /* this returns from vector 'v' the value of the element at position 'pos' */ long i; if (pos > v.length || pos < 1) { printf( " error in call to function vget: position %ld is beyond the end of the vector\n", pos); halt(); } /* move to the correct 'vectorpart' */ for (i = 1; i <= (pos - 1) / maxvecpart; i++) v.part = v.part->next; /* get the proper array element from this part */ return (v.part->numbers[(pos - 1) & (maxvecpart - 1)]); } Static Void vput(v, pos, number) vector *v; long pos; double number; { /* this puts into vector 'v' the value of 'number' at position 'pos' */ long i; vectorpart *firstpart; /* of the vector */ if (pos > v->length || pos < 1) { printf( " error in call to function vput: position %ld is beyond the end of the vector\n", pos); halt(); } /* move to the correct 'vectorpart' */ firstpart = v->part; for (i = 1; i <= (pos - 1) / maxvecpart; i++) v->part = v->part->next; /* put the 'number' into the proper array element for this part */ v->part->numbers[(pos - 1) & (maxvecpart - 1)] = number; v->part = firstpart; } Static Void makevector(v, l) vector *v; long l; { /* create the linked list of vector-parts which are needed for a vector of length 'l' */ long numparts; /* number of parts needed for the vector */ long i; /* index to the parts of the vector */ vectorpart *firstpart; /* of the vector */ vectorpart *newpart; /* of the vector */ if (l < 1) { printf(" makevector: positive length required\n"); halt(); } v->length = l; v->part = (vectorpart *)Malloc(sizeof(vectorpart)); firstpart = v->part; numparts = (v->length - 1) / maxvecpart + 1; for (i = 1; i < numparts; i++) { newpart = (vectorpart *)Malloc(sizeof(vectorpart)); v->part->next = newpart; v->part = newpart; } v->part->next = NULL; v->part = firstpart; } /* end module vector.primitives version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.readvector */ Static Void readvector(thefile, v) _TEXT *thefile; vector *v; { /* read the elements of v into v from 'thefile'. v must already be set up (all the parts created and linked together) before calling. 'thefile' must contain, from the cursor point until the end of the vector, only numbers, either integers or reals; otherwise it will bomb. the vector ends with an end of line so that end of file can be tested for. */ long i, j; /* indecies */ long numparts; /* number of parts in the vector */ long lastpart; /* the number of elements in the last vector part */ vectorpart *firstpart; /* pointer to the first vector part */ firstpart = v->part; numparts = (v->length - 1) / maxvecpart + 1; lastpart = ((v->length - 1) & (maxvecpart - 1)) + 1; for (i = 1; i < numparts; i++) { for (j = 0; j < maxvecpart; j++) fscanf(thefile->f, "%lg", &v->part->numbers[j]); v->part = v->part->next; } for (j = 0; j < lastpart; j++) fscanf(thefile->f, "%lg", &v->part->numbers[j]); fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); v->part = firstpart; } /* end module vector.readvector version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.writevector */ Static Void writevector(thefile, v, y, z) _TEXT *thefile; vector v; long y, z; { /* writes the elements of 'v' to 'thefile'. the integers y and z are: y: the size of the field for printing the number; z: the size of the decimal part of the field, if z = 0 then integers will be printed. */ long pos = 0; /* posititon in the vector */ long i, j; /* indecies */ long numparts; /* number of parts in the vector */ long lastpart; /* the number of elements in the last vector part */ vectorpart *firstpart; /* pointer to the first vector part */ long x; /* the number of elements to write on a line */ x = (long)(vpagewidth / (y + 1.0)); firstpart = v.part; numparts = (v.length - 1) / maxvecpart + 1; lastpart = ((v.length - 1) & (maxvecpart - 1)) + 1; if (z == 0) { for (i = 1; i < numparts; i++) { for (j = 0; j < maxvecpart; j++) { fprintf(thefile->f, " %*ld", (int)y, (long)floor(v.part->numbers[j] + 0.5)); pos++; if (pos % x == 0) putc('\n', thefile->f); /* p2c: encode.p, line 591: * Note: Using % for possibly-negative arguments [317] */ } v.part = v.part->next; } for (j = 0; j < lastpart; j++) { fprintf(thefile->f, " %*ld", (int)y, (long)floor(v.part->numbers[j] + 0.5)); pos++; if (pos < v.length && pos % x == 0) putc('\n', thefile->f); /* p2c: encode.p, line 598: * Note: Using % for possibly-negative arguments [317] */ } } else { for (i = 1; i < numparts; i++) { for (j = 0; j < maxvecpart; j++) { fprintf(thefile->f, " %*.*f", (int)y, (int)z, v.part->numbers[j]); pos++; if (pos % x == 0) putc('\n', thefile->f); /* p2c: encode.p, line 606: * Note: Using % for possibly-negative arguments [317] */ } v.part = v.part->next; } for (j = 0; j < lastpart; j++) { fprintf(thefile->f, " %*.*f", (int)y, (int)z, v.part->numbers[j]); pos++; if (pos < v.length && pos % x == 0) putc('\n', thefile->f); /* p2c: encode.p, line 613: * Note: Using % for possibly-negative arguments [317] */ } } putc('\n', thefile->f); v.part = firstpart; } /* end module vector.writevector version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* *********************************************************************** */ /* end module vector.io version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.dotproduct */ Static double dotproduct(vectora, vectorb) vector vectora, vectorb; { /* this returns the dotproduct of 'vectora' and 'vectorb' */ long i; /* an index */ double j = 0.0; /* partial products */ if (vectora.length != vectorb.length) { printf(" function dotproduct: vector lengths must be equal\n"); halt(); } for (i = 1; i <= vectora.length; i++) j += vget(vectora, i) * vget(vectorb, i); return j; } /* end module vector.dotproduct version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.magnitude */ Static double magnitude(v) vector *v; { /* find the magnitude (length) of the vector v */ return sqrt(dotproduct(*v, *v)); } /* end module vector.magnitude version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module vector.normalize */ Static Void normalize(v) vector *v; { /* this replaces the vector v with the congruent vector of unit length, i.e., the resulting vector v has the property that v.v = 1 */ long i; /* index */ double length; /* of the unnormalized vector */ long FORLIM; length = magnitude(v); FORLIM = v->length; for (i = 1; i <= FORLIM; i++) vput(v, i, vget(*v, i) / length); } /* end module vector.normalize version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* *********************************************************************** */ /* end module vector.functions version = 'matmod 1.95 85 apr 18 tds/gds'; */ /* begin module encode.vectorsize */ Static long vectorsize(param) parameter *param; { /* determine the size of the vector generated by a string of parameter records, begin with 'param' */ long size = 0; /* of the vector, so far */ while (param != NULL) { size += param->pvlength; param = param->next; } return size; } /* end module encode.vectorsize version = 'matmod 1.65 85 apr 9 tds/gds'; */ /* begin module package.getpiece */ /* ************************************************************************ */ /* begin module package.brpiece */ /* ************************************************************************ */ /* begin module book.basis */ /* procedures needed for book manipulations */ /* get procedures should be used for all linked lists of records */ Static Void getline(l) line **l; { /* obtain a line from the free line list or by making a new one */ if (freeline != NULL) { *l = freeline; freeline = freeline->next; } else *l = (line *)Malloc(sizeof(line)); (*l)->length = 0; (*l)->next = NULL; } Static Void getdna(l) dnastring **l; { if (freedna != NULL) { *l = freedna; freedna = freedna->next; } else *l = (dnastring *)Malloc(sizeof(dnastring)); (*l)->length = 0; (*l)->next = NULL; } /* clear procedures should be called each time the records are no longer needed failure to do this may result in a stack overflow. */ Static Void clearline(l) line **l; { /* return a line to the free line list */ line *lptr; if (*l == NULL) return; lptr = *l; *l = (*l)->next; lptr->next = freeline; freeline = lptr; } Static Void writeline(afile, l, carriagereturn) _TEXT *afile; line *l; boolean carriagereturn; { /* write a line to a file, with carriage return if carriagereturn is true. */ long index; /* index to characters in l */ long FORLIM; FORLIM = l->length; for (index = 0; index < FORLIM; index++) putc(l->letters[index], afile->f); if (carriagereturn) putc('\n', afile->f); } Static Void showfreedna() { /* show the freedna list */ long counter = 0; /* count of freedna list */ dnastring *l; /* pointer into freedna list */ l = freedna; while (l != NULL) { counter++; printf("%ld", counter); printf(", length = %d\n", l->length); /* This is illegal according to gpc because one cannot write a pointer to a text file. It can be unearthed for debugging. write(output, ', pointer id: ',l:1); */ l = l->next; } } Static Void cleardna(l) dnastring **l; { /* clear the dna strutures to the free list */ dnastring *lptr; if (*l == NULL) return; lptr = *l; *l = (*l)->next; lptr->next = freedna; freedna = lptr; } Static Void clearheader(h) header *h; { /* clear the header h (remove lines to free storage) */ clearline(&h->fulnam); while (h->note != NULL) clearline(&h->note); } Static Void clearpiece(p) piece **p; { /* clear the dna of the piece */ while ((*p)->dna != NULL) cleardna(&(*p)->dna); clearheader(&(*p)->key.hea); } Static base chartobase(ch) Char ch; { /* convert a character into a base */ base Result; switch (ch) { case 'a': Result = a; break; case 'c': Result = c; break; case 'g': Result = g; break; case 't': Result = t; break; } return Result; } Static Char basetochar(ba) base ba; { /* convert a base into a character */ Char Result; switch (ba) { case a: Result = 'a'; break; case c: Result = 'c'; break; case g: Result = 'g'; break; case t: Result = 't'; break; } return Result; } Static base complement(ba) base ba; { /* take the complement of ba */ base Result; switch (ba) { case a: Result = t; break; case c: Result = g; break; case g: Result = c; break; case t: Result = a; break; } return Result; } Static Char chomplement(b) Char b; { /* create the character complement of base b. I must be getting hungry! */ return (basetochar(complement(chartobase(b)))); } Static long pietoint(p, pie) long p; piece *pie; { /* p is a coordinate on the piece. we want to transform p into a number from 1 to n: an internal coordinate system for easy manipulation of piece coordinates */ /* Note: the dirhomologous and dircomplement are treated as plus and minus directions, which MIGHT NOT BE RIGHT! */ long i; /* an intermediate value */ piekey *WITH; WITH = &pie->key; switch (WITH->piedir) { case dirhomologous: case plus: if (p >= WITH->piebeg) i = p - WITH->piebeg + 1; else i = p - WITH->coobeg + WITH->cooend - WITH->piebeg + 2; break; case dircomplement: case minus: if (p <= WITH->piebeg) i = WITH->piebeg - p + 1; else i = WITH->cooend - p + WITH->piebeg - WITH->coobeg + 2; break; } return i; } Static long inttopie(i, pie) long i; piece *pie; { /* i is in the range 1 to some maximum. it is an internal coordinate system for the program. we want to do a coordinate transformation to obtain a value in the range of the piece called pie: i=1 corresponds to piebeg and i=its maximum corresponds to pieend */ /* Note: the dirhomologous and dircomplement are treated as plus and minus directions, which MIGHT NOT BE RIGHT! */ long p; /* an intermediate value */ piekey *WITH; WITH = &pie->key; switch (WITH->piedir) { case dirhomologous: case plus: p = WITH->piebeg + i - 1; if (p > WITH->cooend) { if (WITH->coocon == circular) p += WITH->coobeg - WITH->cooend - 1; } break; case dircomplement: case minus: p = WITH->piebeg - i + 1; if (p < WITH->coobeg) { if (WITH->coocon == circular) p += WITH->cooend - WITH->coobeg + 1; } break; } return p; } Static long piecelength(pie) piece *pie; { /* return the length of the dna in pie */ return (pietoint(pie->key.pieend, pie)); } /* end module book.basis version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.getto */ Static Char getto(thefile, theline, ch) _TEXT *thefile; long *theline; long *ch; { /* search the file for a character in the first line which is a member of the set ch. Note: on 1999 March 10 the definition of this function was cleaned up. Instead of putting thefile on the line AFTER the charcter ch has been found, it puts thefile ON the line. Other routines like brdna and brpiece have to move to the next line themselves. This makes getto give the OBJECT. */ Char achar = ' '; /* a character in thefile */ boolean done = false; /* done finding achar */ while (!done) { if (BUFEOF(thefile->f)) { done = true; break; } achar = P_peek(thefile->f); if (P_inset(achar, ch)) { done = true; break; } fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } if (P_inset(achar, ch)) return achar; else { return ' '; /* The old method - while (not(achar in ch)) and (not eof(thefile)) do begin readln(thefile,achar); theline := succ(theline) end; if (achar in ch) then getto:=achar else getto:=' ' */ } } /* end module book.getto version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.skipstar */ Static Void skipstar(thefile) _TEXT *thefile; { /* skip start of line (or star = '*'). */ if (BUFEOF(thefile->f)) { printf(" procedure skipstar: end of book found\n"); halt(); return; } if (P_peek(thefile->f) != '*') { printf(" procedure skipstar: bad book\n"); printf(" \"*\" expected as first character on the line, but \"%c\" was found\n", P_peek(thefile->f)); halt(); } getc(thefile->f); /* skip the star */ if (P_peek(thefile->f) != ' ') { /* skip the blank */ printf(" procedure skipstar: bad book\n"); printf(" \"* \" expected on a line but \"*%c\" was found\n", P_peek(thefile->f)); halt(); } getc(thefile->f); } /* skipstar */ /* end module book.skipstar version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brreanum */ Static Void brreanum(thefile, theline, reanum) _TEXT *thefile; long *theline; double *reanum; { /* read a real number from the file */ skipstar(thefile); fscanf(thefile->f, "%lg%*[^\n]", reanum); getc(thefile->f); (*theline)++; } /* end module book.brreanum version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brnumber */ Static Void brnumber(thefile, theline, num) _TEXT *thefile; long *theline, *num; { /* read a number from the file */ skipstar(thefile); fscanf(thefile->f, "%ld%*[^\n]", num); getc(thefile->f); (*theline)++; } /* end module book.brnumber version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brname */ Static Void brname(thefile, theline, nam) _TEXT *thefile; long *theline; name *nam; { /* read a name from the file */ long i; /* an index to the name */ Char c_; /* a character read */ skipstar(thefile); nam->length = 0; do { nam->length++; c_ = getc(thefile->f); if (c_ == '\n') c_ = ' '; nam->letters[nam->length - 1] = c_; } while (!(P_eoln(thefile->f) || nam->length >= namelength || nam->letters[nam->length - 1] == ' ')); if (nam->letters[nam->length - 1] == ' ') nam->length--; if (nam->length < namelength) { for (i = nam->length; i < namelength; i++) nam->letters[i] = ' '; } fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } /* brname */ /* end module book.brname version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brline */ Static Void brline(thefile, theline, l) _TEXT *thefile; long *theline; line **l; { /* read a line from the file */ long i = 0; Char acharacter; skipstar(thefile); /* protection added 2005 Sep 15 */ while (!P_eoln(thefile->f) && i < linelength) { i++; acharacter = getc(thefile->f); if (acharacter == '\n') acharacter = ' '; (*l)->letters[i-1] = acharacter; } /* protection added 2005 Sep 15 */ if (!P_eoln(thefile->f)) { printf("***********************************************\n"); printf("* WARNING: brline: book line length exceeded\n"); printf("* linelength > %ld characters\n", (long)linelength); printf("* Only %ld characters read from book\n", (long)linelength); printf("***********************************************\n"); } (*l)->length = i; (*l)->next = NULL; fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } /* end module book.brline version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brdirect */ Static Void brdirect(thefile, theline, direct) _TEXT *thefile; long *theline; direction *direct; { /* read a direction */ Char ch; skipstar(thefile); fscanf(thefile->f, "%c%*[^\n]", &ch); getc(thefile->f); if (ch == '\n') ch = ' '; (*theline)++; if (ch == '+') *direct = plus; else *direct = minus; } /* end module book.brdirect version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brconfig */ Static Void brconfig(thefile, theline, config) _TEXT *thefile; long *theline; configuration *config; { /* read a configuration */ Char ch; skipstar(thefile); fscanf(thefile->f, "%c%*[^\n]", &ch); getc(thefile->f); if (ch == '\n') ch = ' '; (*theline)++; if (ch == 'l') *config = linear; else *config = circular; } /* end module book.brconfig version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brnotenumber */ Static Void brnotenumber(thefile, theline, note) _TEXT *thefile; long *theline; line **note; { /* book note reading to obtain the number of the object. the procedure returns the value of the number as a global. (this is not such a good practice, but we are stuck with it for now.) */ *note = NULL; numbered = false; number = 0; /* force number to zero if there is no number at all */ /* the next character is n or * depending on whether there are notes */ if (P_peek(thefile->f) != 'n') return; fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; if (P_peek(thefile->f) != 'n') { skipstar(thefile); if (!P_eoln(thefile->f)) { if (P_peek(thefile->f) == '#') { numbered = true; getc(thefile->f); /* move past the number symbol */ fscanf(thefile->f, "%ld", &number); } } do { fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } while (P_peek(thefile->f) != 'n'); fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; return; } fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } /* brnotenumber */ /* end module book.brnotenumber version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brnote */ Static Void brnote(thefile, theline, note) _TEXT *thefile; long *theline; line **note; { /* read note key */ line *newnote; /* the new note */ line *previousnote; /* the last line of the notes */ *note = NULL; if (P_peek(thefile->f) != 'n') /* enter note */ return; fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; if (P_peek(thefile->f) != 'n') { /* abort null note (n/n) */ getline(note); newnote = *note; while (P_peek(thefile->f) != 'n') { /* wait until end of note */ brline(thefile, theline, &newnote); previousnote = newnote; /* get next note */ getline(&newnote->next); newnote = newnote->next; } /* last note was not used, so: */ clearline(&newnote); previousnote->next = NULL; fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; return; } fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; } /* brnote */ /* end module book.brnote version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brheader */ Static Void brheader(thefile, theline, hea) _TEXT *thefile; long *theline; header *hea; { /* read the header of a key. */ fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); /* move past the object name - new definition 1999 Mar 13 */ (*theline)++; /*bbb*/ /* read key name */ brname(thefile, theline, &hea->keynam); /* read full name */ getline(&hea->fulnam); brline(thefile, theline, &hea->fulnam); /* read note key */ if (readnumber) brnotenumber(thefile, theline, &hea->note); else brnote(thefile, theline, &hea->note); } /* end module book.brheader version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.copyheader */ Static Void copyheader(fromhea, tohea) header fromhea, *tohea; { /* copy the header fromhea into tohea. Note that the linked objects are NOT copied, but merely pointed to. */ memcpy(tohea->keynam.letters, fromhea.keynam.letters, sizeof(alpha)); tohea->keynam.length = fromhea.keynam.length; tohea->note = fromhea.note; tohea->fulnam = fromhea.fulnam; } /* end module book.copyheader version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brpiekey */ Static Void brpiekey(thefile, theline, pie) _TEXT *thefile; long *theline; piekey *pie; { /* read piece key, track the line number */ brheader(thefile, theline, &pie->hea); brreanum(thefile, theline, &pie->mapbeg); brconfig(thefile, theline, &pie->coocon); brdirect(thefile, theline, &pie->coodir); brnumber(thefile, theline, &pie->coobeg); brnumber(thefile, theline, &pie->cooend); brconfig(thefile, theline, &pie->piecon); brdirect(thefile, theline, &pie->piedir); brnumber(thefile, theline, &pie->piebeg); brnumber(thefile, theline, &pie->pieend); } /* end module book.brpiekey version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brdna */ Static Void brdna(thefile, theline, dna) _TEXT *thefile; long *theline; dnastring **dna; { /* read in dna from thefile, track the line */ /* note: if the dna were circularized, by linking the last dnastring to the first, then the cleardna routine could not clear properly, and would loop forever... there is no reason to do that, since a simple mod function will allow one to access the circle. */ Char ch; dnastring *workdna; long SET[5]; long TEMP; getdna(dna); workdna = *dna; ch = getto(thefile, theline, P_addset(P_expset(SET, 0L), 'd')); fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); (*theline)++; ch = getc(thefile->f); /* skipstar */ if (ch == '\n') ch = ' '; while (ch == '*') { ch = getc(thefile->f); /* skip blank */ if (ch == '\n') ch = ' '; do { ch = getc(thefile->f); if (ch == '\n') ch = ' '; if (ch == 't' || ch == 'g' || ch == 'c' || ch == 'a') { if (workdna->length == dnamax) { getdna(&workdna->next); workdna = workdna->next; } workdna->length++; TEMP = workdna->length - 1; P_clrbits_B(workdna->part, TEMP, 1, 3); P_putbits_UB(workdna->part, TEMP, (int)chartobase(ch), 1, 3); } } while (!P_eoln(thefile->f)); fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); /* go to next line */ (*theline)++; ch = getc(thefile->f); /* ch is either '*' or 'd' */ if (ch == '\n') ch = ' '; } fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); /* read past the d */ (*theline)++; } /* end module book.brdna version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brpiece */ Static Void brpiece(thefile, theline, pie) _TEXT *thefile; long *theline; piece **pie; { /* read in a piece, change theline to reflect the lines traversed */ /* readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *) theline := succ(theline); (* BUG: was below! *) bbb*/ brpiekey(thefile, theline, &(*pie)->key); if (numbered || !skipunnum) brdna(thefile, theline, &(*pie)->dna); fscanf(thefile->f, "%*[^\n]"); getc(thefile->f); /* move past the word 'piece' - new definition 1999 Mar 13 */ (*theline)++; } /* end module book.brpiece version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.brinit */ Static Void brinit(book, theline) _TEXT *book; long *theline; { /* check that the book is ok to read, and set up the global variables for br routines */ /* halt if the book is bad (first word is 'halt') or the first character is not * */ if (*book->name != '\0') { if (book->f != NULL) book->f = freopen(book->name, "r", book->f); else book->f = fopen(book->name, "r"); } else rewind(book->f); if (book->f == NULL) _EscIO2(FileNotFound, book->name); RESETBUF(book->f, Char); if (!BUFEOF(book->f)) { /* check for the date line */ if (P_peek(book->f) != '*') { if (P_peek(book->f) != 'h') printf(" this is not the first line of a book:\n"); else printf(" bad book:\n"); putchar(' '); while (!(P_eoln(book->f) | BUFEOF(book->f))) { putchar(P_peek(book->f)); getc(book->f); } putchar('\n'); halt(); } } else { printf(" book is empty\n"); halt(); } /* initialize free storage */ freeline = NULL; freedna = NULL; readnumber = true; /* usually we read in numbers for items */ number = 0; /* arbitrary value */ numbered = false; /* the piece has no number (none yet read in) */ skipunnum = false; *theline = 1; } /* brinit */ /* end module book.brinit version = 7.72; {of delmod.p 2007 Jul 23} */ /* ************************************************************************ */ /* end module package.brpiece version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.getpiece */ Static Void getpiece(thefile, theline, pie) _TEXT *thefile; long *theline; piece **pie; { /* move to and read in the next piece in the book */ Char ch; long SET[5]; ch = getto(thefile, theline, P_addset(P_expset(SET, 0L), 'p')); /* get to the next p(iece) in the book */ if (ch != ' ') { brpiece(thefile, theline, pie); /* 1999 june 2: removed this: ch:=getto(thefile,theline,['p']); (* read to end of p *) */ /* bbb - now done in brpiece readln(thefile); (* read past piece *) theline := succ(theline); */ } else clearpiece(pie); } /* end module book.getpiece version = 7.72; {of delmod.p 2007 Jul 23} */ /* ************************************************************************ */ /* end module package.getpiece version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module package.trigger */ /* ************************************************************************ */ /* begin module interact.clearstring */ Static Void clearstring(ribbon) string *ribbon; { /* empty the string */ long index; /* to the ribbon */ for (index = 0; index < maxstring; index++) ribbon->letters[index] = ' '; ribbon->length = 0; ribbon->current = 0; } /* clearstring */ /* end module interact.clearstring version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module interact.writestring */ Static Void writestring(tofile, s) _TEXT *tofile; string *s; { /* write the string s to file tofile, no writeln */ long i; /* index to s */ long FORLIM; FORLIM = s->length; for (i = 0; i < FORLIM; i++) putc(s->letters[i], tofile->f); } /* writestring */ /* end module interact.writestring version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module filler.fillstring */ Static Void fillstring(s, a_) string *s; Char *a_; { /* this procedure makes it reasonably easy to fill the string s with characters. one calls the procedure as: */ /* 1 2 3 4 5 */ /* 12345678901234567890123456789012345678901234567890 */ /* fillstring(s, 'this-is-the-string '); the two comments make it easy to line the characters up. also, for this example, it was assumed that the length of filler as defined by the constant fillermax was 50. */ long length = fillermax; /* of the string without trailing blanks */ long index; /* of s */ clearstring(s); while (length > 1 && a_[length-1] == ' ') length--; if (length == 1 && a_[length-1] == ' ') { printf("fillstring: the string is empty\n"); halt(); } for (index = 0; index < length; index++) s->letters[index] = a_[index]; s->length = length; s->current = 1; } /* fillstring */ /* end module filler.fillstring version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module filler.filltrigger */ Static Void filltrigger(t_, a_) trigger *t_; Char *a_; { /* fill the trigger t */ fillstring(&t_->seek, a_); } /* fillstring */ /* end module filler.filltrigger version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module trigger.proc */ /* this module allows one to scan a series of characters, as from an array or a file, and to "trigger" or detect a simple string in the series. the advantage of the trigger is that several triggers can "observe" a stream of characters at once, each looking for a different thing. some other modules required: interact.const, interact.type */ Static Void resettrigger(t_) trigger *t_; { /* reset the trigger to ground state */ t_->state = 0; t_->skip = false; t_->found = false; } /* resettrigger */ Static Void testfortrigger(ch, t_) Char ch; trigger *t_; { /* look at the character ch. if it is part of the trigger (at the current trigger state), then the trigger state goes higher. if it is not part of the trigger then the trigger state is reset, skip is true and one should skip onward to find the trigger. if the trigger is found, found is true. 1996 Sep 12: Bug found! In the case of a trigger "ab", the program used to miss it for situations like "aab". This was because at the first a it would step up. Then it would see the second a and recognize that was not part of ab. It would fail to realize that it could be the start of a new one. The code now accounts for that possibility. */ t_->state++; /* writestring(list,seek); writeln(list,'testfortrigger seek.letters[',state:1,']:', seek.letters[state],' ch:',ch); */ if (t_->seek.letters[t_->state - 1] == ch) { t_->skip = false; if (t_->state == t_->seek.length) t_->found = true; else t_->found = false; return; } /* it failed. But wait! It could be the beginning of a NEW trigger string! */ if (t_->seek.letters[0] == ch) { t_->state = 1; t_->skip = false; t_->found = false; return; } t_->state = 0; t_->skip = true; t_->found = false; /* reset trigger */ } /* testfortrigger */ /* end module trigger.proc version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module skipblanks */ Static Void skipblanks(thefile) _TEXT *thefile; { /* skip over blanks until a non-blank, or end of line, is found */ while ((P_peek(thefile->f) == ' ') & (!P_eoln(thefile->f))) getc(thefile->f); } Static Void skipnonblanks(thefile) _TEXT *thefile; { /* skip over nonblanks until a blank, or end of line, is found */ while ((P_peek(thefile->f) != ' ') & (!P_eoln(thefile->f))) getc(thefile->f); } Static Void skipcolumn(thefile) _TEXT *thefile; { /* skip over a data column */ skipblanks(thefile); skipnonblanks(thefile); } /* end module skipblanks version = 4.18; (@ of prgmod.p 1996 September 12 */ /* ************************************************************************ */ /* end module package.trigger version = 4.18; (@ of prgmod.p 1996 September 12 */ /* begin module findblank */ Static Void findblank(afile) _TEXT *afile; { /* read a file to find the next blank character */ Char ch; do { ch = getc(afile->f); if (ch == '\n') ch = ' '; } while (ch != ' '); } /* end module findblank version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module findnonblank */ Static Void findnonblank(afile, ch) _TEXT *afile; Char *ch; { /* find the next non blank character in a file, return it in ch. */ *ch = ' '; while (!BUFEOF(afile->f) && *ch == ' ') { *ch = getc(afile->f); if (*ch == '\n') *ch = ' '; if (P_eoln(afile->f)) { fscanf(afile->f, "%*[^\n]"); getc(afile->f); } } } #define maximumrange 10000 /* if the alignment point is more than this distance from the piece ends, the program halts in an attempt to catch the alignment bug... 1991 Jan 11 It appears that the rewrite of the code has removed the bug, but the check will be kept. */ #define semicolon ';' /* end of delila instruction */ /* Local variables for align: */ struct LOC_align { _TEXT *inst; Char ch; /* a character in inst */ trigger endcomment; /* trigger to find '*-)' (ignore the dash!) */ trigger endcurly; /* trigger to find comments: '}' */ } ; /* a dot '.' has been found in the name - ignore the rest of the name - for comparisons with mutations. */ /* procedure rd(var f: text; var ch: char); (* read ch from f allowing inspection of the result *) begin read(f,ch); write(output,ch); write(list,ch); write(output,'<',ch,'>'); end; procedure rdln(var f: text); (* readln f allowing inspection of the result *) begin readln(f); writeln(output); writeln(list); end; */ Local Void skipcomment(f, LINK) _TEXT *f; struct LOC_align *LINK; { /* skip an entire comment */ boolean comment = true; /* true means we are inside a comment */ /* skip to end of comment */ resettrigger(&LINK->endcomment); while (comment) { if (BUFEOF(f->f)) { printf("A comment does not end!\n"); halt(); } if (P_eoln(f->f)) { fscanf(f->f, "%*[^\n]"); getc(f->f); continue; } /* rdln(f) */ LINK->ch = getc(f->f); if (LINK->ch == '\n') LINK->ch = ' '; testfortrigger(LINK->ch, &LINK->endcomment); if (LINK->endcomment.found) { comment = false; /*write(output,'<'); rd(f,ch); write(output,'>');*/ } } } Local Void skipcurly(f, LINK) _TEXT *f; struct LOC_align *LINK; { /* skip an entire comment made by {}*/ boolean comment = true; /* true means we are inside a comment */ /* skip to end of comment */ resettrigger(&LINK->endcurly); while (comment) { if (BUFEOF(f->f)) { printf("A comment does not end!\n"); halt(); } if (P_eoln(f->f)) { fscanf(f->f, "%*[^\n]"); getc(f->f); continue; } /* rdln(f) */ LINK->ch = getc(f->f); if (LINK->ch == '\n') LINK->ch = ' '; testfortrigger(LINK->ch, &LINK->endcurly); if (LINK->endcurly.found) { comment = false; /*write(output,'<'); rd(f,ch); write(output,'>');*/ } } } Local Void skipquote(quote, LINK) trigger quote; struct LOC_align *LINK; { /* skip an entire quote of either the ' or " persuasion */ Char kind; /* the kind of quote, ' or " */ kind = quote.seek.letters[0]; /*writeln(output,'skipquote ',kind);*/ do { findnonblank(LINK->inst, &LINK->ch); /* get to the quote */ } while (!((LINK->ch == kind) | BUFEOF(LINK->inst->f))); if (LINK->ch != kind) { printf("end of quote starting with %c not found\n", kind); halt(); } } /* end module findnonblank version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module align.align */ Static Void align(inst_, book, theline, pie, length, alignedbase) _TEXT *inst_, *book; long *theline; piece **pie; long *length, *alignedbase; { /* documentation on align is in module info.align and delman.use.aligned.books. 1996 Sep 12: The routine now uses the trigger functions found in prgmod. The bug in the oldalign routine (that it misses the end of comments that end in a series of asterisks) has been fixed. It now checks that the piece corresponds to the book. */ struct LOC_align V; long p; /* index to a piece name */ long p1; /* another index to a piece name */ boolean done = false; /* done finding an aligning get */ long thebase; /* the base read in */ boolean indefault = false; /* true when within a default statement. These can contain the word 'piece', which must be ignored. */ trigger gettrigger; /* trigger to find 'get' */ trigger defaulttrigger; /* trigger to find 'default' */ trigger nametrigger; /* trigger to find 'name' */ trigger piecetrigger; /* trigger to find 'piece' */ trigger settrigger; /* trigger to find 'set' */ trigger begincomment; /* trigger to find '(-*' (ignore the dash!) */ trigger begincurly; /* trigger to find comments: '{' */ trigger quote1trigger; /* trigger to find single quote ' */ trigger quote2trigger; /* trigger to find double quote " */ boolean dotteddone; name *WITH; V.inst = inst_; filltrigger(&defaulttrigger, "default "); filltrigger(&gettrigger, "get "); filltrigger(&nametrigger, "name "); filltrigger(&piecetrigger, "piece "); filltrigger(&settrigger, "set "); filltrigger(&begincomment, "(* "); filltrigger(&V.endcomment, "*) "); filltrigger(&begincurly, "{ "); filltrigger(&V.endcurly, "} "); filltrigger("e1trigger, "' "); filltrigger("e2trigger, "\" "); resettrigger(&defaulttrigger); resettrigger(&gettrigger); resettrigger(&nametrigger); resettrigger(&piecetrigger); resettrigger(&settrigger); resettrigger(&begincomment); resettrigger(&begincurly); resettrigger("e1trigger); resettrigger("e2trigger); if (BUFEOF(book->f)) /* if there is still more to the book ... */ return; getpiece(book, theline, pie); /* read in the piece */ if (BUFEOF(book->f)) /* if we found a piece ... */ return; *length = pietoint((*pie)->key.pieend, *pie); /* calculate piece length */ /* now find in inst the next occurance of 'get' */ while (!done) { if (BUFEOF(V.inst->f)) { /* no instructions? */ *alignedbase = 1; /* simply align by the first base */ done = true; break; } if (P_eoln(V.inst->f)) { fscanf(V.inst->f, "%*[^\n]"); getc(V.inst->f); continue; } /*then rdln(inst)*/ V.ch = getc(V.inst->f); if (V.ch == '\n') V.ch = ' '; testfortrigger(V.ch, &begincomment); testfortrigger(V.ch, &begincurly); if (begincomment.found || begincurly.found) { if (V.ch == '*') { skipcomment(V.inst, &V); resettrigger(&begincomment); } else { resettrigger(&begincurly); skipcurly(V.inst, &V); } continue; } testfortrigger(V.ch, &gettrigger); if (gettrigger.found) { findnonblank(V.inst, &V.ch); /* get to "from" */ findblank(V.inst); /* get past "from" */ fscanf(V.inst->f, "%ld", &thebase); /* read in the alignedbase */ /*writeln(output);writeln(output,'thebase = ',thebase:1);*/ *alignedbase = pietoint(thebase, *pie); /*writeln(output,'alignedbase=',alignedbase:1);*/ done = true; } testfortrigger(V.ch, "e1trigger); if (quote1trigger.found) skipquote(quote1trigger, &V); testfortrigger(V.ch, "e2trigger); if (quote2trigger.found) skipquote(quote2trigger, &V); testfortrigger(V.ch, &defaulttrigger); if (defaulttrigger.found) { indefault = true; resettrigger(&defaulttrigger); } if (V.ch == semicolon) indefault = false; testfortrigger(V.ch, &settrigger); if (settrigger.found) { indefault = true; resettrigger(&settrigger); } if (V.ch == semicolon) indefault = false; /* check that piece names are correct */ testfortrigger(V.ch, &piecetrigger); if (indefault) continue; if (!piecetrigger.found) continue; skipblanks(V.inst); /* get to name */ WITH = &(*pie)->key.hea.keynam; /* for p := 1 to length do begin */ /* 2007 Jun 22: replace loop with while so that we can drop out when dotted names are detected. */ p = 1; dotteddone = false; while (!dotteddone) { if (P_eoln(V.inst->f)) { dotteddone = true; break; } V.ch = getc(V.inst->f); if (V.ch == '\n') V.ch = ' '; /* ignore names after a dot */ /* if ch = '.' then writeln(output,'inst dotteddone'); */ if (V.ch == '.') dotteddone = true; if (WITH->letters[p-1] == '.') dotteddone = true; /* if ch = '.' then writeln(output,'book dotteddone'); writeln(output,'BUBBa ch = ',ch,' ',p:1); */ /*zzz*/ if (WITH->letters[p-1] != V.ch && !dotteddone && V.ch != ';') { printf("The piece name in the book: \n"); /* p2c: encode.p, line 1750: Note: * Format for packed-array-of-char will work only if width < length [321] */ printf("%.*s\n", WITH->length, WITH->letters); printf("does not match the inst file piece name:\n"); /* write the letters that matched: */ for (p1 = 0; p1 <= p - 2; p1++) putchar(WITH->letters[p1]); /* write the offending letter: */ putchar(V.ch); /* get the rest of the name and show it: */ done = P_eoln(V.inst->f); while (!done) { done = P_eoln(V.inst->f); if (done) break; V.ch = getc(V.inst->f); if (V.ch == '\n') V.ch = ' '; if (V.ch == ' ' || V.ch == ';') done = true; if (!done) putchar(V.ch); } putchar('\n'); /* mark the first letter that does not match: */ for (p1 = 1; p1 < p; p1++) putchar(' '); printf("^\n"); halt(); } p++; if (p > WITH->length) { dotteddone = true; /* we are not inside a comment */ } } } /*rd(inst,ch);*/ if (*alignedbase > -maximumrange && *alignedbase <= *length + maximumrange) return; printf(" In procedure align:\n"); printf(" read in base was %ld\n", thebase); printf(" in internal coordinates: %ld\n", *alignedbase); printf(" maximum range was %ld\n", (long)maximumrange); printf(" piece length was %ld\n", *length); WITH = &(*pie)->key.hea.keynam; /* p2c: encode.p, line 1797: Note: * Format for packed-array-of-char will work only if width < length [321] */ printf(" piece name: %.*s\n", WITH->length, WITH->letters); printf(" piece number: %ld\n", number); printf(" aligned base is too far away... see the code\n"); halt(); } #undef maximumrange #undef semicolon /* end module align.align version = 7.72; {of delmod.p 2007 Jul 23} */ /* begin module book.getbase */ Static base getbase(position, pie) long position; piece *pie; { /* Get a base from the position (internal coordinates) of the piece. Protection is made against positions outside the piece. In the case of circles it would be convenient to wrap around when requests are off the end. So the routine will do a modular wrap for positions outside the range 1 to the length. This is a new feature as of 2000 March 22. */ dnastring *workdna; /* pointer to the dna part of pie */ long p; /* current count of bases into the workdna */ long spot; /* the last base of the dna part */ long thelength; /* the length of the piece */ /* writeln(output,'NEW getbase: position=',position:1,'^^^^^^^^^^^^^^^^^^^^'); */ /* handle cases of position out of range by circular wrapping */ thelength = piecelength(pie); while (position < 1) position += thelength; while (position > thelength) position -= thelength; workdna = pie->dna; p = workdna->length; while (position > p) { /* writeln(output,' workdna^.length=',workdna^.length:1); */ workdna = workdna->next; if (workdna == NULL) { printf("error in function getbase!\n"); halt(); } p += workdna->length; } /* writeln(output,'p=',p:1); */ if (true) { spot = workdna->length - p + position; /* writeln(output,'spot=',spot:1); showdnasegment(output,workdna, spot); */ if (spot <= 0) { printf("error in getbase, spot (= %ld) must be positive\n", spot); halt(); } if (spot > workdna->length) { printf("error in getbase, spot (=%ld) must be less than length (=%d)\n", spot, workdna->length); halt(); } /* writeln(output,'base = ', workdna^.part[spot]); */ return ((base)P_getbits_UB(workdna->part, spot - 1, 1, 3)); } printf("error in getbase: request off end of piece\n"); halt(); } /* end module book.getbase version = 7.72; {of delmod.p 2007 Jul 23} */ Static Void initialize() { long i; /* an index */ printf(" encode, version %4.2f\n", version); if (*inst.name != '\0') { if (inst.f != NULL) inst.f = freopen(inst.name, "r", inst.f); else inst.f = fopen(inst.name, "r"); } else rewind(inst.f); if (inst.f == NULL) _EscIO2(FileNotFound, inst.name); RESETBUF(inst.f, Char); if (BUFEOF(inst.f)) noinst = true; else noinst = false; brinit(&book, &theline); if (*encseq.name != '\0') { if (encseq.f != NULL) encseq.f = freopen(encseq.name, "w", encseq.f); else encseq.f = fopen(encseq.name, "w"); } else { if (encseq.f != NULL) rewind(encseq.f); else encseq.f = tmpfile(); } if (encseq.f == NULL) _EscIO2(FileNotFound, encseq.name); SETUPBUF(encseq.f, Char); firstparam = (parameter *)Malloc(sizeof(parameter)); apiece = (piece *)Malloc(sizeof(piece)); fourpowers[0] = 1; /* 4**0 = 1 */ for (i = 1; i <= codingmax; i++) fourpowers[i] = fourpowers[i-1] * 4; } Static Void setparam() { /* get the parameter list from the file encodep and put the information into the linked parameter list */ parameter *newparam; /* the new parameter list from 'encodep' */ parameter *param; /* a parameter record being read in */ spacelist *firstspaces; /* in the spaces list */ spacelist *newspaces; /* for the spaces list */ long i; /* an index */ Char ch; /* to check for ':' */ long FORLIM; if (*encodep.name != '\0') { if (encodep.f != NULL) encodep.f = freopen(encodep.name, "r", encodep.f); else encodep.f = fopen(encodep.name, "r"); } else rewind(encodep.f); if (encodep.f == NULL) _EscIO2(FileNotFound, encodep.name); RESETBUF(encodep.f, Char); if (BUFEOF(encodep.f)) { printf(" error: empty parameter file (encodep)\n"); halt(); return; } /* read the alignment type */ fscanf(encodep.f, "%c%*[^\n]", &alignmenttype); getc(encodep.f); if (alignmenttype == '\n') alignmenttype = ' '; if (alignmenttype != 'b' && alignmenttype != 'i' && alignmenttype != 'f') { printf("alignment type inst must be f, b, or i\n"); halt(); } param = firstparam; regions = 0; do { regions++; fscanf(encodep.f, "%ld%ld%*[^\n]", param->range, ¶m->range[(long)stop]); getc(encodep.f); if (param->range[(long)start] > param->range[(long)stop]) { printf( " error in parameter file: end of range cannot be less than beginning of range\n"); halt(); } fscanf(encodep.f, "%ld%*[^\n]", ¶m->window); getc(encodep.f); if (param->window <= 0) { printf(" error in parameter file: window size must be positive\n"); halt(); } fscanf(encodep.f, "%ld%*[^\n]", ¶m->wshift); getc(encodep.f); if (param->wshift <= 0) { printf(" error in parameter file: window shift must be positive\n"); halt(); } fscanf(encodep.f, "%ld", ¶m->coding); if (param->coding <= 0) { printf(" error in parameter file: coding must be positive\n"); halt(); } if (param->coding > codingmax) { printf(" error in parameter file: coding can not be greater than %ld\n", (long)codingmax); halt(); } if (param->coding > codingmax) { printf(" error in coding file: requested coding level too large\n"); halt(); } param->spaces = (spacelist *)Malloc(sizeof(spacelist)); ch = ' '; while (!P_eoln(encodep.f) && ch == ' ') { ch = getc(encodep.f); if (ch == '\n') ch = ' '; } if (ch == ':') { firstspaces = (spacelist *)Malloc(sizeof(spacelist)); param->spaces = firstspaces; fscanf(encodep.f, "%ld", ¶m->spaces->skips); FORLIM = param->coding; for (i = 2; i < FORLIM; i++) { newspaces = (spacelist *)Malloc(sizeof(spacelist)); param->spaces->next = newspaces; param->spaces = newspaces; fscanf(encodep.f, "%ld", ¶m->spaces->skips); } param->spaces->next = NULL; param->spaces = firstspaces; } else param->spaces = NULL; fscanf(encodep.f, "%*[^\n]"); getc(encodep.f); fscanf(encodep.f, "%ld%*[^\n]", ¶m->cshift); getc(encodep.f); if (param->cshift <= 0) { printf(" error in parameter file: coding shift must be positive\n"); halt(); } /* wvlength is 4**coding */ param->wvlength = (long)floor(exp(param->coding * log(4.0)) + 0.5); /* pvlength is wvlength * (the number of windows in the range) */ param->pvlength = param->wvlength * (long)((double)(param->range[(long)stop] - param->range[(long)start]) / param->wshift + 1); if (BUFEOF(encodep.f)) param->next = NULL; else { newparam = (parameter *)Malloc(sizeof(parameter)); param->next = newparam; param = newparam; } } while (!BUFEOF(encodep.f)); } Static Void encheader() { /* write the header information to encseq */ parameter *param; /* points to a parameter record */ spacelist *aspace; /* in the spaces list */ long i; /* an index */ parameter *WITH; long FORLIM; fprintf(encseq.f, " encode %4.2f; encoding of sequences in\n", version); putc(' ', encseq.f); copyaline(&book, &encseq); fprintf(encseq.f, "\n %ld independently coded regions\n\n", regions); param = firstparam; while (param != NULL) { WITH = param; fprintf(encseq.f, " %ld to %ld", WITH->range[(long)start], WITH->range[(long)stop]); fprintf(encseq.f, " is encoded as:\n"); fprintf(encseq.f, "%6c%ld long windows\n", ' ', WITH->window); fprintf(encseq.f, "%6c%ld bases shift to new window\n", ' ', WITH->wshift); fprintf(encseq.f, "%6c%ld", ' ', WITH->coding); if (WITH->coding > 1) { fprintf(encseq.f, " :"); if (WITH->spaces == NULL) { FORLIM = WITH->coding; for (i = 1; i < FORLIM; i++) fprintf(encseq.f, " 0"); } else { aspace = WITH->spaces; FORLIM = WITH->coding; for (i = 1; i < FORLIM; i++) { fprintf(encseq.f, " %ld", aspace->skips); aspace = aspace->next; } } } fprintf(encseq.f, " are the coding units\n"); fprintf(encseq.f, "%6c%ld bases shift to new coding site\n", ' ', WITH->cshift); param = param->next; putc('\n', encseq.f); } fprintf(encseq.f, " %ld is the vector length\n\n", vectorsize(firstparam)); } Static Void codeit(apiece, alignedbase, v) piece *apiece; long alignedbase; vector *v; { /*aparam: paramptr;*/ /* this procedure takes one sequence (which has been aligned) and encodes it into a string of integers according to the parameters the user has specified. each window of the sequence causes a vector to be printed, the size of which is determined by the coding level specified. for example, if the coding level is 1 (monos are being counted) then the vector has four elements. if the coding level is 2 (dis are being counted) then the vector has 16 elements. the vectors for each window are concatenated to give the vector for each sequence, which is ended with the 'end of sequence' symbol. the procedure takes advantage of the fact that the type 'base' is an ordered set of the nucleotides, with a,c,g,t being assigned internally the values 0,1,2,3. the total number of elements of a vector is 4**coding-level, and these correspond to the number of oligos of each type from all a's to all t's. this makes it easy to increase the vector element for the oligo which exists at the coding site. for example, there are 64 tri-nucleotides and the vector which encodes them has elements 0 (for the oligo aaa) through 63 (for the oligo ttt). the number of the oligo cgt would be stored in element 27, as seen from 16 * 1 (16 for 4**(coding-level - 1) and 1 for the c) + 4 * 2 ( 4 for 4**(coding-level - 2) and 2 for the g) + 1 * 3 ( 1 for 4**(coding-level - 3) and 3 for the t) --------- 27 another example is finding the position of the oligo taac in the vector of all the 256 tetramers, from position 0 to 255. 64 * 3 (64 = 4**(coding-level - 1) and 3 for the t) + 16 * 0 (16 = 4**(coding-level - 2) and 0 for the a) + 4 * 0 ( 4 = 4**(coding-level - 3) and 0 for the a) + 1 * 1 ( 1 = 4**(coding-level - 4) and 1 for the c) --------- 193 ( is the vector element corresponding to the oligo taac ) */ parameter *param; /* list of parameters for encoding */ spacelist *aspace; /* in the spaces list */ long startsite; /* of a window on the sequence */ long firstpos; /* of a coding site within a window */ long pos; /* of a particular base to be encoded */ long sitesize; /* total number of bases in a coding site, both those that are encoded and those that are skipped */ long element = 1; /* of the encoding vector, from 1 to v.length */ long welement; /* of the window vector, from 0 to 4**coding - 1 */ long length; /* of the sequence */ long i; /* an index */ long FORLIM; length = piecelength(apiece); param = firstparam; FORLIM = v->length; /* clear the vector */ for (i = 1; i <= FORLIM; i++) vput(v, i, 0.0); /* start at the beginning fo the vector */ do { /* for each parameter record in the list */ /* determine the number of bases in a coding site */ aspace = param->spaces; sitesize = param->coding; while (aspace != NULL) { sitesize += aspace->skips; aspace = aspace->next; } /* set the window beginning at the range start */ startsite = alignedbase + param->range[(long)start]; do { /* move the window across the range */ /* set coding site beginning at window beginning */ firstpos = startsite; do { /* if the whole coding site is on the piece */ if (firstpos > 0) { if (firstpos + sitesize - 1 <= length || apiece->key.piecon == circular && firstpos <= length) { pos = firstpos; welement = 0; aspace = param->spaces; for (i = param->coding - 1; i >= 0; i--) { welement += fourpowers[i] * (int)getbase(pos, apiece); if (aspace != NULL) { pos += aspace->skips + 1; aspace = aspace->next; } else pos++; /* for circular sequences */ if (pos > length) pos -= length; } /* increment the vector at position 'element + welement' */ vput(v, element + welement, vget(*v, element + welement) + 1); } } firstpos += param->cshift; } while (firstpos < startsite + param->window); startsite += param->wshift; /* advance element beyond window */ element += param->wvlength; } while (startsite <= alignedbase + param->range[(long)stop]); param = param->next; } while (param != NULL); } /* codeit */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; encodep.f = NULL; strcpy(encodep.name, "encodep"); encseq.f = NULL; strcpy(encseq.name, "encseq"); book.f = NULL; strcpy(book.name, "book"); inst.f = NULL; strcpy(inst.name, "inst"); initialize(); setparam(); encheader(); seqvector.length = vectorsize(firstparam); makevector(&seqvector, seqvector.length); if (noinst) { if (alignmenttype == 'i') { alignmenttype = 'b'; printf("There are no instructions so alignment type is forced to b mode\n"); } } while (!BUFEOF(book.f)) { if (alignmenttype != 'i') { getpiece(&book, &theline, &apiece); /* read in the piece */ if (!BUFEOF(book.f)) { length_ = piecelength(apiece); switch (alignmenttype) { case 'i': printf("prgm error\n"); halt(); break; case 'b': alignedbase = pietoint(0L, apiece); break; case 'f': alignedbase = 1; break; } } } else align(&inst, &book, &theline, &apiece, &length_, &alignedbase); if (BUFEOF(book.f)) /*firstparam,*/ break; codeit(apiece, alignedbase, &seqvector); writevector(&encseq, seqvector, 3L, 0L); clearpiece(&apiece); } _L1: if (inst.f != NULL) fclose(inst.f); if (book.f != NULL) fclose(book.f); if (encseq.f != NULL) fclose(encseq.f); if (encodep.f != NULL) fclose(encodep.f); exit(EXIT_SUCCESS); } /* encode */ /* End. */