program delila (inst, book, listing, marksdelila, lib1, cat1, lib2, cat2, lib3, cat3, debug, output); (* lll ccc *) (* delila: the librarian for sequence manipulation by Thomas Schneider Gary Stormo Paul Morrissett useful suggestions by Jeff haemer module libraries needed: delman, delmods. Dr. Thomas D. Schneider National Cancer Institute Laboratory of Experimental and Computational Biology Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu (use only if first address fails) http://www.lecb.ncifcrf.gov/~toms/ *) (* label ******************************************************************) label 1; (* end of the program used in procedure halt only *) (* constants **************************************************************) const (* begin module version *) version = 4.89; (* of delila.p 2007 Dec 06 2007 Dec 06: 4.89; always reduce coordinates if the piece is circular 2005 Sep 6: 4.88; demote error 220 to a warning: deleting a whole sequence merely means not to put it in the book! Implement by not writing to the book ... 2005 Jan 11: 4.87; cleanup 2005 Jan 11: 4.86; long names cause infnite loop in irlongname this was missing: longname^.next := nil; 2004 Jan 21: 4.85; finish up marks. 2004 Jan 21: 4.84; cleanup 2004 Jan 21: 4.83; gpctime.p upgrade. compiles but gives segmentation fault! needed to set pk^.dna := nil; libpie^.dna := nil; 2004 Jan 21: 4.82; marks. arrow -> marks.arrow 2003 Aug 18: 4.81; halt message includes the name 'delila' 2003 Apr 7: 4.80; bug fix - request for change before start of piece now ok. 2002 Apr 17: 4.79; tab and other odd ASCII characters disallowed: error 33. 2001 Oct 5: 4.78; mutations off piece are now warnings, not errors (215, 216) 2001 Mar 28: 4.77; pass 2 errors written independently of pass 1 2001 Mar 28: 4.75; bug: two titles bombs! 2001 Mar 16: 4.74; clean up 2001 Mar 16: 4.73; gene functions; specpiece must use libpie not pk! 2001 Mar 16: 4.72; upgrade bug documentation (gene no longer functions) 2000 Nov 15: 4.70; catch unclosed comments 2000 Oct 26: 4.69; set maximum book size 2000 Oct 18: 4.67; stmts w/o with still trigger marksdelila stepping 2000 Oct 18: 4.66; upgrade to gpc 2000 Oct 17: 4.63; 'direction + with a4021131c[cr]; caused parse bomb. 2000 Oct 17: 4.62; blank after 'direction +' caused parse bomb. 2000 Aug 17: 4.60; fixed withused bug; marksdelila only written in pass 2 2000 July 12: 4.59; withused control: marksdelila not touched unless needed. 2000 June 21: 4.58; upgrade the See Also section. 2000 May 26: 4.55; allow inserts to be complemented. zzz111 2000 May 25: 4.54; Fix insertion bug in complements (proc. changesequence)! 2000 May 25: 4.53; Remove duplication of error reports to output. 2000 Mar 29: 4.52; Report error numbers to output, saves checking listing. 2000 January 3: 4.50; {} comments for the html links 2000 January 3: 4.49; html link for delila instructions 1999 August 17: 4.48; changes not noted 1999 July 17: 4.47: Add message You_need_a_marks.arrow_definition to marksdelila file so user will see this if they forget. 1999 July 13: 4.46: putbase, getbase etc now allow any length in each dnasegment - much more robust. 1999 July 9: 4.20: bug fix in getdnasegment dnamax could not be small (b and bDNAptr need to be stepped if they exceed dnamax) 1999 May 27: 4.07: user does not propagate coordinates 1999 May 24: 4.03: mutations cannot overlap - makes all marks work! not implemented - need other changes first 1999 May 23: 4.00: a26g works as t26c on the complement 1999 May 14: 3.76: and always are in increasing order, following the new definition in Libdef. 1999 May 12: 3.66: new comment type, {} 1999 May 11: 3.56: ability of insertion/deletions to wrap around 1999 May 10: 3.50: upgraded error reporting for incorrect change commands. 1999 May 6: 3.21: added errors 215, 216; moved mutation routines below geteoinst so error output is cleaner. 1999 May 5: 3.17: fix delila position misreading 1999 April 27: 3.10: fix memory leak in circledna/invert 1999 April 26: 3.06: fix bug in: get from 6 to 1 with i4,3ggttgg; 1999 April 15: 3.03: fix bug in: get from 1 to 6 with i3,4ggttgg; 1999 April 14: 3.00: allow insert off ends of piece (fix bug) 1999 April 14: 2.99: fixed equalname: need to set i initially 1999 April 13: 2.98: fixed memory leak 1999 March 21: 2.90: marks are sorted so most displays will work. 1999 March 21: 2.85: Synonymes: default = set 1999 March 20: 2.84: Delila can now create mutation marks for lister! 1999 March 19: 2.70: Delila can now create mutations! 1999 March 17: 2.61: For completeness, Delila can now extract a single base from the complementary strand, as in "get from 1 to 1 direction -;" 1999 March 14: 2.44: dopiece functions entirely on internal variables! 1999 March 13: 2.40: conversion to standard delmod book reading routines. 1999 Mar 9: 2.34 The catal file gives a line number that each object starts on. The bookreading routine getto used to get to the line after the start of an object, and the routines like brdna, brpiece, getocp all took this into account. To make it more clear, getto in delmod and here now gets to the start of the object. The routines that then read the object may readln past the start if they want. 1999 Mar 9: 2.31 Convert to standard book reading routings: Procedure dopiece needs to read in the dna, make mutations, clip out the relevant part and then spit it out to the book. It needs to keep track of library lines if the catal is to be of any use. But the standard book reading routine brpiece does not track the lines read. Therefore it was necessary to alter delmod.p so that the standard book reading routines keep track of the line number. This was done. 1999 Mar 6: Delila absored the mutation mechanism from dbmutate. Not functional yet because it needs to read pieces in standard way. 1998 June 24: default coordinate 0 was being set in pass 1 but not reset to normal at the start of pass 2. 1998 January 26: namelength set to 100 to allow long names 1998 January 4: dnamax set to 10 million for faster grabs 1997 January 10: For convenience, the default is: only pieces are numbered. 1996 September 2: Two reductions off end of piece is now flagged. 1996 August 12: introduction of 'same'. 1995 December 7: objects can now be named in the long name 1995 Nov 13: default coordinate zero now allows default coordinate (number) last changes: 1989 November 14 origin: 1980 or so *) (* end module version *) (* begin module describe.delila *) (* name delila: the librarian for sequence manipulation synopsis delila(inst: in, book: out, listing: out, marksdelila: out, lib1: in, cat1: in, lib2: in, cat2: in, lib3: in, cat3: in, output: out, debug: out) files inst: instructions written in the language delila that tell the program delila what sequences to pull out of the library. book: the set of sequences pulled out of the library. listing: the instructions are listed along with errors found or actions taken. marksdelila: Colored marks for the lister program that indicate the locations of base changes, insertions and deletions. lib1: the first library from which to obtain sequences cat1: the first catalogue, corresponding to lib1 lib2: the second library cat2: the second catalogue, corresponding to lib2 lib3: the third library cat3: the third catalogue, corresponding to lib3 debug: traces through the actions taken, for debugging delila (only produced if variable debugging is true.) output: messages to the user description Delila is a data base manager for nucleic acid sequences. It takes a set of instructions, written in the language delila (DEoxyribonucleic acid LIbrary LAnguage) and a large set of sequences called a library. The output is a listing of the actions taken (or errors) corresponding to the instructions, and a "book" containing the sequences desired. examples see the documentation documentation libdef (defines delila), delman.intro, delman.use, delman.construction philgen.ps see also {Definition of the database system:} libdef {Introduction to Delila Instructions:} http://www.lecb.ncifcrf.gov/~toms/delilainstructions.html {Related programs:} alist.p, catal.p, loocat.p, dbbk.p, lister.p dbmutate.p {is deprecated: use the mutation method, 'with'} {(} http://www.m-w.com/cgi-bin/dictionary?deprecated {)} author Thomas D. Schneider, Gary D. Stormo and Paul Morrisett useful suggestions by Jeff Haemer bugs There used to be many known bugs in delila, related to extracting linear fragments of circular sequences. Delila was rebuilt in the spring of 1999 and is much more robust now. The following features are not available in this program: recognition classes and enzymes, markers, automatic printing to the book of structures that intersect a piece, get all (for org, chr, rec and enz), get every and if. The gene mechanism was revived on 2001 Mar 16, eventually it may be used to implement features. Delila programs use a packed array of bases. This means that a smart Pascal compiler can store DNA sequences in two bits per base. When delila was originally designed, I thought that everybody would complete their sequences and therefore there would never be unknown bases in a database. Silly me. GenBank has plenty of such cases. The dbbk program avoids this problem by converting such bases to 'a'. Fortunately these can now be displayed with lister. Someday Delila may be upgraded to handle this case, but it might be at the cost of reducing the maximum sequence that can be handled. *) (* end module describe.delila *) (* files used: inst = instructions book = the book that is printed listing = instruction listing lib1,lib2,...numlibfil = the files of the library cat1,cat2,...numcatfil = the files of the catalogue debug = listing for debugging the code output = for fatal error messages to the terminal procedure name conventions: cr catalogue read ir instruction read lr library read bw book write flow of information in the librarian [file name] (procedure) [library] [catalogue] [instructions] : : : v v v : : : : (catalogue : : procedures) : : : : :....... ......: : : : : (lr"s) (ir"s) : : v v : : : ................: : : : : (delila) : : v v : : .......: :................ : : : (ir"s) (bw"s) (writeerror"s) : (writevalue"s) : : v v : : [book] [listing] further documentation for this program is in: 'organism and recognition class library definition: a dna sequence data base' 1980 june 9 note.. the following features are not yet available in this program: recognition class and enzymes markers automatic printing to the book of structures that intersect a piece get all (for org, chr, rec and enz) get every if lll = places that must be changed when one changes the number of library files: numlibfil ccc = places that must be changed when one changes the number of catalogue files: numcatfil *) numlibfil = 3; (* number of library files lll *) numcatfil = 3; (* number of catalogue files ccc *) namelength = 100; (* maximum key name length *) linelength = 200; (* maximum line readable in library *) dnamax = 1024; (* maximum bases of dna in a part of a dna string *) (* (60 bits/word)*(1base/2bits)*100 words=3000 bases *) (* 1998 Jan 4: set to 10 million bytes - why not?? *) (* 1999 Jul 7: 10 thousand bytes - ran out of memory! *) dnalinemax = 60; (* bases of dna per line written in the book *) sitemax=20; (* bases of a site in a part of a site string *) (* (60 bits/word)*(1 sitebase/3bits) * (1 word) = 20 sitebases *) maxstep = 10; (* the maximum number of steps of the traversal chart *) instwidth = 6; (* width of left edge numbers on the listing *) decbase = 2; (* number of decimal places for bases in marksdelila *) widbase = 6; (* width of places for bases in marksdelila *) decbits = 2; (* number of decimal places for bits in marksdelila *) widbits = 6; (* width of places for bits in marksdelila *) maxbook = maxint; (* maximum book size *) (* begin module datetime.const *) datetimearraylength = 19; (* length of dataarray for dates, It is just long enough to include the 4 digit year - solving the year 2000 problem: 1980/06/09 18:49:11 123456789 123456789 1 2 *) (* end module datetime.const version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module changeset.const *) changesetmax = 20; (* maximum number of changes allowed *) insertmax = 1000; (* maximum insertion length allowed (bp) *) (* end module changeset.const version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* types *****************************************************************) type (* the nodes in the library tree *) node = (libnode,orgnode,chrnode,marnode,mardnanode,tranode,gennode, pienode,piednanode,recnode,enznode,sitenode); (* step is the sequence of characters that must be printed in the book to go from one node to another *) step = packed array[1..maxstep] of char; (* begin module book.type ****************) (* types needed for book manipulations *) chset = set of 'a'..'z'; (* types defined in book definition *) alpha = packed array[1..namelength] of char; (* this is not alfa *) (* name is a left justified string with blanks following the characters *) name = record letters: alpha; length: 0..namelength (* zero means an unspecified structure *) end; lineptr = ^line; line = record (* a line of characters *) letters: packed array [1..linelength] of char; length: 0..linelength; next: lineptr end; direction = (plus, minus, dircomplement, dirhomologous); configuration = (linear, circular); state = (on, off); header = record (* header of key *) keynam: name; (* key name of structure *) fulnam: lineptr; (* full name of structure *) note: lineptr (* note key *) end; (* begin module base.type *) (* define the four nucleotide bases *) base = (a,c,g,t); (* end module base.type version = 7.52; {of delmod.p 2000 Jul 30} *) (* sequence types *) dnaptr = ^dnastring; dnarange = 0..dnamax; seq = packed array[1..dnamax] of base; dnastring = record part: seq; length: dnarange; next: dnaptr end; orgkey = record (* organism key *) hea: header; mapunit: lineptr (* genetic map units *) end; chrkey = record (* chromosome key *) hea: header; mapbeg: real; (* number of genetic map beginning *) mapend: real (* number of genetic map ending *) end; pieceptr = ^piece; piekey = record (* piece key *) hea: header; mapbeg: real; (* genetic map beginning *) coocon: configuration; (* configruation (circular/linear) *) coodir: direction; (* direction (+/-) relative to genetic map *) coobeg: integer; (* beginning nucleotide *) cooend: integer; (* ending nucleotide *) piecon: configuration; (* configruation (circular/linear) *) piedir: direction; (* direction (+/-) relative to coordinates *) piebeg: integer; (* beginning nucleotide *) pieend: integer; (* ending nucleotide *) end; piece = record key: piekey; dna: dnaptr end; reference = record pienam : name; (* name of piece referred to *) mapbeg : real; (* genetic map beginning *) refdir : direction; (* direction relative to coordinates *) refbeg : integer; (* beginning nucleotide *) refend : integer; (* ending nucleotide *) end; genkey = record (* gene key *) hea : header; ref : reference; end; trakey = record (* transcript key *) hea : header; ref : reference; end; markerptr = ^marker; markey = record (* marker key *) hea : header; ref : reference; sta : state; phenotype : lineptr; next : markerptr; end; marker = record key : markey; dna : dnaptr; end; (* end module book.type version = 7.52; {of delmod.p 2000 Jul 30} *) (* types defined in library definition that are not in the standard book reading routines *) reckey = record (* dna recognition class key *) hea: header end; enzymeptr = ^enzyme; enzkey = record (* enzyme key *) hea: header end; sitebase = (sa,sc,sg,st,pu,py,sn, modification,cleaveage,unknown,alternative); siteptr = ^sitestring; sitestring = record (* a single recognition site *) part: packed array[1..sitemax] of sitebase; length: 0..sitemax; next: siteptr (* pointer to another recognition site *) end; (* c) storage of the recognition of a dna sequence *) enzyme = record key: enzkey; sites: siteptr end; (* default types *) defaultkey = record (* to print the key or not to print the key *) note: state; (* note *) mar: state; (* marker *) gen: state; (* gene *) tra: state (* transcript *) end; defaultsite = record (* to act on the site or not to act on the site *) expand: state; (* expand sites *) modify: state; (* modify sites *) cleave: state; (* cleave sites *) end; (* how to act when a request is out of range: *) rangeaction = (rreduce,rcontinue,rhalt); numberedstructure = (orgnum, chrnum, marnum, tranum, gennum, pienum, recnum, enznum); defaultnumber = record sta: state; (* whether or not to number in the book *) str: array[numberedstructure] of state; (* what can be numbered *) item: integer; (* the next item"s number *) end; (* coordinate type: *) coordinatetype = (coornormal,coorzero); default = record key: defaultkey; sit: defaultsite; defout: rangeaction; (* odd name to get around sun4 bug... *) num: defaultnumber; coo: coordinatetype; doubling: state; (* whether to make two pieces when mutating *) arrowlength: real; (* length of mark arrow in lines *) end; (* catalogue types *) item = record (* an item in the catalogue *) letter: char; (* type of structure *) nam: name; (* the structure"s key name *) line: integer (* location of the structure in the library *) end; catfile = file of item; (* begin module datetime.type *) (* array for dates *) datetimearray = packed array[1..datetimearraylength] of char; (* end module datetime.type version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module delila.changeset.type *) changedata = record changetype: char; (* the type of change: c(hange), i(nsertion), d(eletion) *) baseold: char; (* the old base given in a change instruction *) basenew: char; (* the new base given in a change instruction *) basecoo1: integer; (* the first coordinate, external coordinates *) basecoo2: integer; (* the second coordinate, external coordinates *) internal1: integer; (* the first coordinate, internal coordinates *) internal2: integer; (* the second coordinate, internal coordinates *) insertasdeletion: boolean; (* insertion is acting as a deletion *) {dddDDDddd} inserts: integer; (* number of bases to insert *) insert: array[1..insertmax] of char; (* bases to insert *) {zzz111} insertcomplement: boolean; (* insert the complement of the typed sequence *) end; changeset = record (* the complete set of changes for an entry *) data: array[1..changesetmax] of changedata; number: integer; (* number of changes *) end; (* from changeset.type version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* end module delila.changeset.type *) (* variables **************************************************************) var (* file variables *) lib1,lib2,lib3: text; (* lll *) cat1,cat2,cat3: catfile; (* ccc *) inst,listing,book,marksdelila,debug: text; libline: array[1..numlibfil] of integer; (* location in the libfiles *) firstlibrary: integer; (* the first library used. this is the one from which the parent date of the book originates in bwbookheader *) versioninbook: boolean; (* this variable allows delila to write its version into the book after the first organism name *) (* time variables *) datetime: datetimearray; (* date of book withdrawal *) (* catalogue search varibles *) catnumber : 1..numcatfil; (* number of catfile being searched *) catline: array[1..numcatfil] of integer; (* line in the catfile *) currento,currentc,currentr,currentl : integer; (* the line in the current catalogue of the previously found organism, recognition site, chromosome or other key name *) itemfound: boolean; (* whether the requested item was found or not *) (* default variables *) def: default; (* all of the default variables *) indnum: numberedstructure; (* an index for def.num.str *) (* traverse tree variables *) (* the traversal chart defines how one can move around the library tree. some "steps" are illegal (illegaltraversal). the others are the first characters of words to be printed in the book between nodes of the tree *) traversalchart : array[node,node] of step; illegaltraversal: step; pastlibrary: node; (* the currently specified node from the library *) pastbook: node; (* the currently specified node for the book *) pastcheck: node; (* node for instruction checking *) nodechar: array[node] of char; (* the characters for each node *) nodeletter: char; (* the current node character *) debugging: boolean; (* for debugging purposes *) withused: boolean; (* set true if the 'with' command was used. This determines whether we need to write the marksdelila file *) booksize: integer; (* current size of book in bases to compare to maxbook *) (* begin module book.var *) (* ************************************************************************ *) (* global variables needed for book manipulations *) (* free storage: *) freeline: lineptr; (* unused lines *) freedna: dnaptr; (* unused dnas *) readnumber: boolean; (* whether to read a number from the notes, or to read in the notes *) number: integer; (* the number of the item just read *) numbered: boolean; (* true when the item just read is numbered *) skipunnum: boolean; (* a control variable to allow skipping of un-numbered items in the book *) (* ************************************************************************ *) (* end module book.var version = 7.52; {of delmod.p 2000 Jul 30} *) (* free storage for unimplemented functions *) freemarker: markerptr; (* unused markers *) freesite: siteptr; (* unused sites *) (* character reading variables *) ch: char; (* spare character *) firstch: char; (* first character of each library line *) blank: char; (* second character of each library attribute *) skipping: boolean; (* whether to skip the first two characters of a library line in routine skipstar *) (* memory of DNA and piece key information in the library *) libpie: pieceptr; (* the entire piece as recorded in the library *) libpieit: item; (* used to record the location of the dna of a piece in the library *) libpiefi: integer; (* file for libpieit *) (* who points where lineptrs point to a linked list of lines markers point to a linked list of markers, with one (or few) pieces of dna per marker pieces point to a linked list of dna strings dnaptrs point to a linked list of dnastrings enzymes point to a linked list of sites siteptrs point to a linked list of sites *) (* variable for creating zeroed coordinate system *) zerobase: integer; zeroshift: integer; (* amount to shift away from zero coordinate *) zeroBS: integer; (* sum of zerobase and zeroshift, for efficiency *) (* variables to handle the NAME to be given to the next long name *) longname: lineptr; (* the longname for the next object *) longnameexists: boolean; (* true if a longname was read *) (* variables for mutating sequences *) mutations: changeset; mutposition1: integer; (* error: mutation position 1 *) mutposition2: integer; (* error: mutation position 2 *) mutnotchar: char; (* error: what the mutation character is not *) mutischar: char; (* error: what the mutation character or command is *) mutcd1,mutcd2: changedata; (* changes to a sequence *) (* variable to remember last piece to avoid rereading it if we can. The mechanism will reread if the org/chr is respecified or we star pass 2. *) lastpiecename: name; (* variables for making marksdelila file *) insertupperbits: real; (* upperbits for insertion symbol *) insertlowerbits: real; (* lowerbits for insertion symbol *) deleteupperbits: real; (* upperbits for deletion symbol *) deletelowerbits: real; (* lowerbits for deletion symbol *) changeupperbits: real; (* upperbits for change symbol *) changelowerbits: real; (* lowerbits for change symbol *) (* begin module crash *) procedure crash; (* Crash the program by trying to open a nonexistant file. This allows tracing by the dbx program. To use: insert call into the halt program or whereever a traceable stop is desired. *) var bogus: text; (* boghous internal file *) begin writeln(output,' program crash.'); reset(bogus); end; (* end module crash version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module halt *) procedure 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. *) begin writeln(output,' program halt.'); goto 1 end; (* end module halt version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* 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 *) procedure getline(var l: lineptr); (* obtain a line from the free line list or by making a new one *) begin if freeline<>nil then begin l:=freeline; freeline:=freeline^.next end else new(l); l^.length:=0; l^.next:=nil end; procedure getdna(var l: dnaptr); begin if freedna<>nil then begin l:=freedna; freedna:=freedna^.next end else new(l); l^.length:=0; l^.next:=nil end; (* clear procedures should be called each time the records are no longer needed failure to do this may result in a stack overflow. *) procedure clearline(var l: lineptr); (* return a line to the free line list *) var lptr: lineptr; begin if l<>nil then begin lptr:=l; l:=l^.next; lptr^.next:=freeline; freeline:=lptr end end; procedure writeline(var afile: text; l: lineptr; carriagereturn: boolean); (* write a line to a file, with carriage return if carriagereturn is true. *) var index: integer; (* index to characters in l *) begin with l^ do begin for index := 1 to length do write(afile, letters[index]); end; if carriagereturn then writeln(afile); end; procedure showfreedna; (* show the freedna list *) var counter: integer; (* count of freedna list *) l: dnaptr; (* pointer into freedna list *) begin l := freedna; counter := 0; while l <> nil do begin counter := succ(counter); write(output,counter:1); write(output, ', length = ',l^.length:1); { 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); } writeln(output); l := l^.next end; end; procedure cleardna(var l: dnaptr); var lptr: dnaptr; begin if l<>nil then begin lptr:=l; l:=l^.next; lptr^.next:=freedna; freedna:=lptr end end; procedure clearheader(var h: header); (* clear the header h (remove lines to free storage) *) begin with h do begin clearline(fulnam); while note<>nil do clearline(note) end end; procedure clearpiece(var p: pieceptr); (* clear the dna of the piece *) begin while p^.dna<>nil do cleardna(p^.dna); clearheader(p^.key.hea) end; function chartobase(ch:char):base; (* convert a character into a base *) begin case ch of 'a': chartobase:=a; 'c': chartobase:=c; 'g': chartobase:=g; 't': chartobase:=t end end; function basetochar(ba:base):char; (* convert a base into a character *) begin case ba of a: basetochar:='a'; c: basetochar:='c'; g: basetochar:='g'; t: basetochar:='t'; end end; function complement(ba:base):base; (* take the complement of ba *) begin case ba of a: complement:=t; c: complement:=g; g: complement:=c; t: complement:=a; end end; function chomplement(b: char): char; (* create the character complement of base b. I must be getting hungry! *) begin chomplement := basetochar(complement(chartobase(b))); end; function pietoint(p: integer; pie: pieceptr): integer; (* 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 *) var i: integer; (* an intermediate value *) begin with pie^.key do begin case piedir of plus: if p>=piebeg then i:=p-piebeg+1 else i:=(p-coobeg)+(cooend-piebeg)+2; minus: if p<=piebeg then i:=piebeg-p+1 else i:=(cooend-p)+(piebeg-coobeg)+2 end; pietoint:=i end end; function inttopie(i: integer; pie: pieceptr):integer; (* 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 *) var p: integer; (* an intermediate value *) begin with pie^.key do begin case piedir of plus: begin p:=piebeg+(i-1); if p>cooend then if coocon=circular then p:=p-(cooend-coobeg+1) end; minus: begin p:=piebeg-(i-1); if p '*' then begin writeln(output,' procedure skipstar: bad book'); writeln(output,' "*" expected as first character on the line, but "', thefile^,'" was found'); halt end; get(thefile); (* skip the star *) if thefile^ <> ' ' then begin writeln(output,' procedure skipstar: bad book'); writeln(output,' "* " expected on a line but "*', thefile^,'" was found'); halt end; get(thefile) (* skip the blank *) end end; (* skipstar *) (* end module book.skipstar version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brreanum *) procedure brreanum(var thefile: text; var theline: integer; var reanum: real); (* read a real number from the file *) begin skipstar(thefile); readln(thefile,reanum); theline := succ(theline) end; (* end module book.brreanum version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brnumber *) procedure brnumber(var thefile: text; var theline: integer; var num: integer); (* read a number from the file *) begin skipstar(thefile); readln(thefile,num); theline := succ(theline) end; (* end module book.brnumber version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brname *) procedure brname(var thefile: text; var theline: integer; var nam: name); (* read a name from the file *) var i: integer; (* an index to the name *) c: char; (* a character read *) begin (* brname *) skipstar(thefile); with nam do begin length:=0; repeat length:=succ(length); read(thefile,c); letters[length] := c until (eoln(thefile)) or (length>=namelength) or (letters[length]=' '); if letters[length]=' ' then length:=length-1; if length 'n' then begin skipstar(thefile); if not eoln(thefile) then begin if thefile^ = '#' then begin numbered := true; get(thefile); (* move past the number symbol *) read(thefile,number); end end; repeat readln(thefile); theline := succ(theline) until thefile^ = 'n'; readln(thefile); theline := succ(theline) end else begin readln(thefile); theline := succ(theline) end end end; (* brnotenumber *) (* end module book.brnotenumber version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brnote *) procedure brnote(var thefile: text; var theline: integer; var note: lineptr); (* read note key *) var newnote: lineptr; (* the new note *) previousnote: lineptr; (* the last line of the notes *) begin (* brnote *) note:=nil; if thefile^ = 'n' then begin (* enter note *) readln(thefile); theline := succ(theline); if thefile^ <> 'n' then begin (* abort null note (n/n) *) getline(note); newnote:=note; while thefile^ <> 'n' do begin (* wait until end of note *) brline(thefile,theline,newnote); previousnote:=newnote; (* get next note *) getline(newnote^.next); newnote:=newnote^.next; end; (* last note was not used, so: *) clearline(newnote); previousnote^.next:=nil; readln(thefile); theline := succ(theline); end else begin readln(thefile); theline := succ(theline); end; end end; (* brnote *) (* end module book.brnote version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brheader *) procedure brheader(var thefile: text; var theline: integer; var hea: header); (* read the header of a key. *) begin with hea do begin readln(thefile); (* move past the object name - new definition 1999 Mar 13 *) theline := succ(theline); {bbb} (* read key name *) brname(thefile,theline,keynam); (* read full name *) getline(fulnam); brline(thefile,theline,fulnam); (* read note key *) if readnumber then brnotenumber(thefile,theline,note) else brnote(thefile,theline,note) end end; (* end module book.brheader version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.copyheader *) procedure copyheader(fromhea: header; var tohea: header); (* copy the header fromhea into tohea. Note that the linked objects are NOT copied, but merely pointed to. *) begin tohea.keynam.letters := fromhea.keynam.letters; tohea.keynam.length := fromhea.keynam.length; tohea.note := fromhea.note; tohea.fulnam := fromhea.fulnam; end; (* end module book.copyheader version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brpiekey *) procedure brpiekey(var thefile: text; var theline: integer; var pie: piekey); (* read piece key, track the line number *) begin with pie do begin brheader(thefile,theline,hea); brreanum(thefile,theline,mapbeg); brconfig(thefile,theline,coocon); brdirect(thefile,theline,coodir); brnumber(thefile,theline,coobeg); brnumber(thefile,theline,cooend); brconfig(thefile,theline,piecon); brdirect(thefile,theline,piedir); brnumber(thefile,theline,piebeg); brnumber(thefile,theline,pieend); end end; (* end module book.brpiekey version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brdna *) procedure brdna(var thefile: text; var theline: integer; var dna: dnaptr); (* 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. *) var ch: char; workdna: dnaptr; begin getdna(dna); workdna:=dna; ch:=getto(thefile,theline,['d']); readln(thefile); theline := succ(theline); read(thefile,ch); (* skipstar *) while (ch = '*') do begin read(thefile,ch); (* skip blank *) repeat read(thefile,ch); if ch in ['a','c','g','t'] then begin if workdna^.length=dnamax then begin getdna(workdna^.next); workdna:=workdna^.next end; workdna^.length:=succ(workdna^.length); workdna^.part[workdna^.length]:=chartobase(ch) end until eoln(thefile); readln(thefile); (* go to next line *) theline := succ(theline); read(thefile,ch); (* ch is either '*' or 'd' *) end; readln(thefile); (* read past the d *) theline := succ(theline); end; (* end module book.brdna version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brpiece *) procedure brpiece(var thefile: text; var theline: integer; var pie: pieceptr); (* read in a piece, change theline to reflect the lines traversed *) begin { 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 or (not skipunnum) then brdna(thefile,theline,pie^.dna); readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *) theline := succ(theline); end; (* end module book.brpiece version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brinit *) procedure brinit(var book: text; var theline: integer); (* check that the book is ok to read, and set up the global variables for br routines *) begin (* brinit *) (* halt if the book is bad (first word is 'halt') or the first character is not * *) reset(book); if not eof(book) then begin (* check for the date line *) if book^ <> '*' then begin if book^ <> 'h' then writeln(output, ' this is not the first line of a book:') else writeln(output, ' bad book:'); write(output, ' '); while not (eoln(book) or eof(book)) do begin write(output, book^); get(book) end; writeln(output); halt end end else begin writeln(output, ' book is empty'); halt end; (* initialize free storage *) freeline:=nil; freedna:=nil; 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; end; (* brinit *) (* end module book.brinit version = 7.52; {of delmod.p 2000 Jul 30} *) (* ************************************************************************ *) (* end module package.brpiece version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.getpiece *) procedure getpiece(var thefile: text; var theline: integer; var pie: pieceptr); (* move to and read in the next piece in the book *) var ch: char; begin ch:=getto(thefile,theline,['p']); (* get to the next p(iece) in the book *) if ch<>' ' then begin 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); } end else clearpiece(pie); end; (* end module book.getpiece version = 7.52; {of delmod.p 2000 Jul 30} *) (* ************************************************************************ *) (* end module package.getpiece version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brorgkey *) procedure brorgkey(var thefile: text; var theline: integer; var org: orgkey); (* read organism key *) begin with org do begin {bbb} brheader(thefile,theline,hea); getline(mapunit); brline(thefile,theline,mapunit); end end; (* end module book.brorgkey version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.brchrkey *) procedure brchrkey(var thefile: text; var theline: integer; var chr: chrkey); (* read chromosome key *) begin with chr do begin {bbb} brheader(thefile,theline,hea); brreanum(thefile,theline,mapbeg); brreanum(thefile,theline,mapend); end end; (* end module book.brchrkey version = 7.52; {of delmod.p 2000 Jul 30} *) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* begin module package.bwrite *) (****************************************************************************) (* this is a package of procedures for writing books, by gary stormo, aug 17, 1982 *) (* begin module book.bwbasics *) procedure bwstartline(var book: text); (* start a line of output to the book *) begin write(book,'* ') end; procedure bwline(var book: text; l: lineptr); (* write a line to the book *) var i: integer; begin if l<>nil then begin bwstartline(book); if l^.length<>0 then for i:=1 to l^.length do write(book,l^.letters[i]); writeln(book); end end; procedure bwtext(var book: text; lines: lineptr); (* write a set of lines to the book *) var l: lineptr; begin l := lines; while l<>nil do begin bwline(book,l); l := l^.next; end; end; procedure bwnote(var book: text; note: lineptr); (* writes the notes pointed to by 'note' to 'book' *) begin if note<>nil then begin writeln(book,'note'); bwtext(book,note); writeln(book,'note'); end end; procedure bwnumber(var book: text; num: integer); (* write a number to the book *) begin bwstartline(book); writeln(book,num:1) (* pascal will expand the field as far as needed *) end; procedure bwreanum(var book: text; reanum: real); (* write a real number to the book *) begin bwstartline(book); writeln(book,reanum:1:2) (* pascal will expand the field *) end; procedure bwstate(var book: text; sta: state); (* write a state to the book *) begin bwstartline(book); case sta of on: writeln(book,'on'); off: writeln(book,'off') end end; procedure bwname(var book: text; nam: name); (* write a name to the book *) var i: integer; begin with nam do begin bwstartline(book); for i:=1 to length do write(book,letters[i]); writeln(book) end end; procedure bwdirect(var book: text; direct: direction); (* write a direction to the book *) begin bwstartline(book); case direct of plus: writeln(book,'+'); minus: writeln(book,'-') end end; procedure bwconfig(var book: text; config: configuration); (* write a configuration to the book *) begin bwstartline(book); case config of linear: writeln(book,'linear'); circular: writeln(book,'circular') end end; procedure bwheader(var book: text; hea: header); (* write a key header to the book *) begin with hea do begin bwname(book,keynam); bwline(book,fulnam); bwnote(book,note); end end; procedure bworgkey(var book: text; org: orgkey); (* write the organism key *) begin with org do begin bwheader(book,hea); bwline (book,mapunit) end end; procedure bwchrkey(var book: text; chr: chrkey); (* write the chromosome key *) begin with chr do begin bwheader(book,hea); bwreanum(book,mapbeg); bwreanum(book,mapend) end end; (* end module book.bwbasics version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bworg *) procedure bworg(var thefile: text; org: orgkey; var chropen, orgopen: boolean); (* this writes the organism key 'org' to 'thefile', and returns 'orgopen' as true. if there is already a chromosome or organism open they are closed. *) begin if chropen then begin writeln(thefile,'chromosome'); chropen := false; end; if orgopen then writeln(thefile,'organism'); writeln(thefile,'organism'); bworgkey(book,org); orgopen := true; end; (* end module book.bworg version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwchr *) procedure bwchr(var thefile: text; chr: chrkey; var chropen: boolean); (* write the chromosome key 'chr' to 'thefile' and return chropen as true. if a chromosome is already open it is closed. *) begin if chropen then writeln(thefile,'chromosome'); writeln(thefile,'chromosome'); bwchrkey(book,chr); chropen := true; end; (* end module book.bwchr version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwdna *) procedure bwdna(var thefile: text; d: dnaptr); (* write the dna pointed to by 'd' to 'thefile' *) var i: integer; (* index to the sequence *) l: integer; (* index to the number of bases on the line *) newline: boolean; (* true when a new line should be started *) begin writeln(thefile,'dna'); newline := true; while (d <> nil) do begin for i := 1 to d^.length do begin if newline then begin bwstartline(thefile); l := 0; newline := false end; write(thefile,basetochar(d^.part[i])); l := succ(l); if (l mod 60 = 0) or (* end of line *) ((i = d^.length) and (d^.next = nil)) (* last base *) then begin writeln(thefile); newline := true end; end; d := d^.next; end; if not newline then writeln(thefile); (* last base *) writeln(thefile,'dna'); end; (* end module book.bwdna version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwpie *) procedure bwpie(var thefile: text; pie: pieceptr); (* writes the information pointed to by 'pie' to 'thefile' *) begin writeln(thefile,'piece'); with pie^ do begin bwheader(thefile,key.hea); bwstartline(thefile); writeln(thefile,key.mapbeg:1:2); bwstartline(thefile); if (key.coocon = circular) then writeln(thefile,'circular') else writeln(thefile,'linear'); bwstartline(thefile); if (key.coodir = plus) then writeln(thefile,'+') else writeln(thefile,'-'); bwstartline(thefile); writeln(thefile,key.coobeg:1); bwstartline(thefile); writeln(thefile,key.cooend:1); bwstartline(thefile); if (key.piecon = circular) then writeln(thefile,'circular') else writeln(thefile,'linear'); bwstartline(thefile); if (key.piedir = plus) then writeln(thefile,'+') else writeln(thefile,'-'); bwstartline(thefile); writeln(thefile,key.piebeg:1); bwstartline(thefile); writeln(thefile,key.pieend:1); end; bwdna(thefile,pie^.dna); writeln(thefile,'piece'); end; (* end module book.bwpie version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwref *) procedure bwref(var book: text; ref: reference); (* write a key reference to the book *) begin with ref do begin bwname (book,pienam); bwreanum(book,mapbeg); bwdirect(book,refdir); bwnumber(book,refbeg); bwnumber(book,refend) end end; (* end module book.bwref version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwgen *) procedure bwgen(var thefile: text; gene: genkey); (* this proecdure writes to 'thefile' the information in 'gene', properly formatted; *) begin writeln(thefile,'gene'); with gene do begin bwheader(thefile,hea); bwref(thefile,ref); end; writeln(thefile,'gene'); end; (* end module book.bwgen version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwtra *) procedure bwtra(var thefile: text; trans: trakey); (* this proecdure writes to 'thefile' the information in 'trans', properly formatted; *) begin writeln(thefile,'transcript'); with trans do begin bwheader(thefile,hea); bwref(thefile,ref); end; writeln(thefile,'transcript'); end; (* end module book.bwtra version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.bwmar *) procedure bwmar(var thefile: text; mark: marker); (* this proecdure writes to 'thefile' the information in 'mark', properly formatted; *) var i : integer; begin writeln(thefile,'marker'); with mark.key do begin bwheader(thefile,hea); bwref(thefile,ref); bwstate(thefile,sta); bwstartline(thefile); for i := 1 to phenotype^.length do write(thefile,phenotype^.letters[i]); writeln(thefile); end; bwdna(thefile,mark.dna); writeln(thefile,'marker'); end; (* end module book.bwmar version = 7.52; {of delmod.p 2000 Jul 30} *) (****************************************************************************) (* end module package.bwrite version = 7.52; {of delmod.p 2000 Jul 30} *) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* begin module nwpietoint *) function nwpietoint(p: integer; pie: pieceptr): integer; (* no wrap version of pietoint *) (* 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. If the coordinate is off the end, bring it to just before the end. *) var i: integer; (* an intermediate value *) begin with pie^.key do begin case piedir of plus: i := p - piebeg + 1; minus: i := piebeg - p + 1; end; if i > piecelength(pie) then i := piecelength(pie) + 1; if i < 0 then i := 0; { writeln(output,'nwpietoint ---------------'); if piedir = minus then writeln(output,'nwpietoint: pidir is minus') else writeln(output,'nwpietoint: pidir is plus'); writeln(output,'nwpietoint: p: ',p:1); writeln(output,'nwpietoint: pieend: ',pieend:1); writeln(output,'nwpietoint: i: ',i:1); } nwpietoint:=i end end; (* end module nwpietoint version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* note: 3 bars indicate breaks between major sections of code *) (**************************************************************************) (* catalogue procedures ***************************************************) (**************************************************************************) procedure tvrslibrary(future: node); (* check that the library is traversed properly, set nodeletter *) var thisstep: step; begin thisstep:=traversalchart[pastlibrary,future]; if thisstep=illegaltraversal then begin writeln(output,' program error: illegal library traversal ', ord(pastlibrary),' to ',ord(future)); halt end; pastlibrary:=future; nodeletter:=nodechar[future] end; (**************************************************************************) procedure grabitem(var cat: catfile; var it: item); (* grab an item from the cat *) begin (* grabitem *) it := cat^; if not eof(cat) then get(cat) end; (* grabitem *) procedure crcatheader(var cat: catfile; catnum: integer; var alp: datetimearray); (* read the catalogue header (date of creation) from the catalogue file, cat *) var it: item; (* a catalogue item for unloading *) i: integer; (* index to date *) begin (* crcatheader *) (* the header of the catalogue resides in the first item *) reset(cat); if not eof(cat) then begin grabitem(cat,it); for i:=1 to datetimearraylength do alp[i]:=it.nam.letters[i]; catline[catnum]:=2 end (* else: objections to this case are made in checklibcat *) end; (* crcatheader *) procedure nextitem (var newitem: item); (* this procedure finds the next item in the files, returns that item and the new current line and makes itemfound = true. if itemfound is false after calling this routine it means that the last item was the end of last file. by gary stormo aug 28,1979 modified by tom schneider 79 sep 7 *) begin (* nextitem *) itemfound:= false; case catnumber of (* ccc *) 1: if not eof(cat1) then begin grabitem(cat1,newitem); itemfound:= true end; 2: if not eof(cat2) then begin grabitem(cat2,newitem); itemfound:= true end; 3: if not eof(cat3) then begin grabitem(cat3,newitem); itemfound:= true end end; if itemfound then catline[catnumber]:=succ(catline[catnumber]) else begin (* go to the next file *) case catnumber of (* ccc *) 1: begin catnumber:= 2; reset(cat2); if not eof(cat2) then begin grabitem(cat2,newitem); (* skip date *) grabitem(cat2,newitem); itemfound:= true end; end; 2: begin catnumber:= 3; reset(cat3); if not eof(cat3) then begin grabitem(cat3,newitem); (* skip date *) grabitem(cat3,newitem); itemfound:= true end; end; 3: end; catline[catnumber]:=2; (* 2 because first item was read *) end; end; (* nextitem *) procedure finditem (newnode: node; n: name; var newitem: item; var fi: integer); (* procedure to search the catalogue for the requested item, returns the file number (fi) and line number of the request in the library. by gary stormo, aug 28, 1979 modified by tom schneider 79 sep 5 *) (* this version assumes that no organism or recognition has been split in two between two files *) var numsearched: integer; (* the number of catalogues searched in this particular finditem call *) procedure search(var cat: catfile); (* search one of the catfiles *) var stuck: boolean; (* when the keyname is not found between the beginning of the search and the nodeletter of a higher structure - to reset and research *) procedure readcat; (* read a catalogue item *) begin grabitem(cat,newitem); catline[catnumber]:=succ(catline[catnumber]) (*;if debugging then writeln(debug,'readcat, catnumber=',catnumber:3,' catline=' ,catline[catnumber]:4);*) end; procedure find; (* find an item. if found set itemfound true *) begin (* find *) (*if debugging then writeln(debug,'find');*) readcat; if newitem.letter = nodeletter then begin if newitem.nam.letters = n.letters then begin if debugging then writeln(debug,'found ',newitem.letter,' ',newitem.nam.letters); itemfound:= true; (* the file is the same as the catalogue *) fi := catnumber; case nodeletter of 'o': currento:= catline[catnumber]; 'c': currentc:= catline[catnumber]; 'r': currentr:= catline[catnumber]; 't','m','g','p','e': currentl:= catline[catnumber] end end end end; (* find *) procedure searchbetween(start,stop: integer); (* search the cat between start and stop *) begin (* reset the cat *) reset(cat); catline[catnumber]:=1; readcat; (* skip header of cat *) (* get to start *) while catline[catnumber]= place then digit else acharacter := '0' end; pluckdigit:=acharacter end; (* pluckdigit *) (* end module pluckdigit *) begin (* according to: http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_109.html#SEC109 *) GetTimeStamp(t); (* Predefined time stamp: http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_389.html#SEC389 TimeStamp = {@@packed} record DateValid, TimeValid : Boolean; Year : Integer; Month : 1 .. 12; Day : 1 .. 31; DayOfWeek : 0 .. 6; { 0 means Sunday } Hour : 0 .. 23; Minute : 0 .. 59; Second : 0 .. 61; { to allow for leap seconds } MicroSecond : 0 .. 999999 end; *) with t do begin if TimeValid then begin { writeln(output,'valid time'); writeln(output,'year =',year:4); writeln(output,'month =',month:2); writeln(output,'day =',day:2); writeln(output,'hour =',hour:2); writeln(output,'minute =',minute:2); writeln(output,'second =',second:2); } adatetime := 'year/mm/dd hh:mm:ss'; adatetime[ 1] := pluckdigit(year,3); adatetime[ 2] := pluckdigit(year,2); adatetime[ 3] := pluckdigit(year,1); adatetime[ 4] := pluckdigit(year,0); adatetime[ 6] := pluckdigit(month,1); adatetime[ 7] := pluckdigit(month,0); adatetime[ 9] := pluckdigit(day,1); adatetime[10] := pluckdigit(day,0); adatetime[12] := pluckdigit(hour,1); adatetime[13] := pluckdigit(hour,0); adatetime[15] := pluckdigit(minute,1); adatetime[16] := pluckdigit(minute,0); adatetime[18] := pluckdigit(second,1); adatetime[19] := pluckdigit(second,0); end else begin writeln(output,'getdatetime: invalid time!'); halt; end end; { Sun compiler method: date(adatetime,'%Y/%m/%d %H:%M:%S'); } end; (* end module getdatetime version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module readdatetime *) procedure readdatetime (var thefile: text; var adatetime: datetimearray); (* read the date and time from the file *) var index: integer; (* to the udatetime *) (* the following is an unpacked date time array, to avoid reading into a packed array. reading into a packed array is not transportable *) udatetime: array[1..datetimearraylength] of char; begin for index:=1 to datetimearraylength do read(thefile,udatetime[index]); pack(udatetime, 1, adatetime); if (adatetime[3]='/') and (adatetime[12]=':') then begin writeln(output,' old datetime (only 2 year digits) read: ', adatetime:datetimearraylength); end; end; (* end module readdatetime version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module writedatetime *) procedure writedatetime(var thefile: text; adatetime: datetimearray); (* expand the date and time out and print in the file *) var index: integer; (* index of datetime *) begin for index:=1 to datetimearraylength do write(thefile,adatetime[index]) end; (* end module writedatetime version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module timeseed *) procedure addtoseed(var seed, power: real; c: char); (* add the digit represented by c to the seed at the power position *) var n: integer; (* the character represented by c *) begin (* addtoseed *) power := power/10; { writeln(output,'addtoseed, c = ',c); writeln(output,'addtoseed, ord(c) = ',ord(c)); } case c of ' ': begin writeln(output,'timeseed: error in datetime'); halt; end; '0': n := 0; '1': n := 1; '2': n := 2; '3': n := 3; '4': n := 4; '5': n := 5; '6': n := 6; '7': n := 7; '8': n := 8; '9': n := 9 end; (*writeln(output,'timeseed number is [',n:1,']'); (@ debug *) seed := seed + power*n end; (* addtoseed *) procedure makeseed(adatetime: datetimearray; var seed: real); (* convert adatetime to a real number in seed, reversed order *) var power: real; (* a digit of the seed such as 0.01 *) begin seed := 0.0; power := 1.0; addtoseed(seed, power, adatetime[19]); addtoseed(seed, power, adatetime[18]); (* : *) addtoseed(seed, power, adatetime[16]); addtoseed(seed, power, adatetime[15]); (* : *) addtoseed(seed, power, adatetime[13]); addtoseed(seed, power, adatetime[12]); (* *) addtoseed(seed, power, adatetime[10]); addtoseed(seed, power, adatetime[9]); (* / *) addtoseed(seed, power, adatetime[7]); addtoseed(seed, power, adatetime[6]); (* / *) addtoseed(seed, power, adatetime[4]); addtoseed(seed, power, adatetime[3]); addtoseed(seed, power, adatetime[2]); addtoseed(seed, power, adatetime[1]); end; procedure orderseed(adatetime: datetimearray; var seed: real); (* convert adatetime to a real number in seed, normal order *) var power: real; (* a digit of the seed such as 0.01 *) begin seed := 0.0; power := 1.0; addtoseed(seed, power, adatetime[1]); addtoseed(seed, power, adatetime[2]); addtoseed(seed, power, adatetime[3]); addtoseed(seed, power, adatetime[4]); (* / *) addtoseed(seed, power, adatetime[6]); addtoseed(seed, power, adatetime[7]); (* / *) addtoseed(seed, power, adatetime[9]); addtoseed(seed, power, adatetime[10]); (* *) addtoseed(seed, power, adatetime[12]); addtoseed(seed, power, adatetime[13]); (* : *) addtoseed(seed, power, adatetime[15]); addtoseed(seed, power, adatetime[16]); (* : *) addtoseed(seed, power, adatetime[18]); addtoseed(seed, power, adatetime[19]); end; procedure timeseed(var seed: real); (* read the computer date and time. reverse the order of the digits and put a decimal point in front. this gives a fraction between zero and one that varies quite quickly, and is always unique (if the computer has sufficient accuracy). it is to be used as a seed to a random number generator. *) var adatetime: datetimearray; (* a date and time *) begin (* timeseed *) getdatetime(adatetime); { writeln(output,'timeseed: adatetime: ',adatetime); } makeseed(adatetime, seed); end; (* timeseed *) (* end module timeseed version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* begin module limitdate *) (* end module limitdate *) (* end module package.datetime version = 1.08; (@ of gpctime.p 2004 Jan 21 *) (* ************************************************************************ *) (* ************************************************************************ *) (* ************************************************************************ *) (* begin module book.other *) (* other useful dna, line and name manipulating functions *) (* begin module book.name *) procedure clearname(n: name); (* clear the name n *) var i: integer; (* index to piece name *) begin n.length := 0; for i := 1 to namelength do n.letters[i] := ' '; end; procedure writename(var f: text; n: name); (* write the name n to file f *) var i: integer; (* index to piece name *) begin for i := 1 to n.length do write(f,n.letters[i]); end; procedure copyname(a: name; var b: name); (* copy name a to name b *) var i: integer; (* index to piece name *) begin for i := 1 to a.length do b.letters[i] := a.letters[i]; b.length := a.length; end; function equalname(a,b: name): boolean; (* is name a equal to name b? *) var i: integer; (* index to piece name *) same: boolean; (* temporary variable to hold the answer *) begin same := true; (* what optimism!! *) if b.length = a.length then begin i := 1; while same and (i <= a.length) do begin same := (b.letters[i] = a.letters[i]); i := succ(i); end end else same := false; equalname := same; end; (* end module book.name version = 7.52; {of delmod.p 2000 Jul 30} *) procedure emptydna(var l: dnaptr); (* empty all of the dna in l onto the freedna list *) begin while l<>nil do cleardna(l); end; procedure emptyline(var l: lineptr); (* empty all of the line in l onto the freeline list *) var count: integer; (* count of lines removed *) begin while l<>nil do clearline(l); end; procedure copyline(fromline: lineptr; var toline: lineptr); (* copy a set of lines *) var f,t : lineptr; (* internal pointers *) begin (* destroy previous to line *) emptyline(toline); if fromline <> nil then begin getline(toline); f:=fromline; t:= toline; while f<>nil do begin t^.letters:=f^.letters; t^.length:=f^.length; f:=f^.next; if f<>nil then begin getline(t^.next); t:=t^.next end else t^.next:=nil end end end; procedure copydna(var fromdna,todna: dnaptr); (* copy a set of dna *) var aDNAptr, memdna: dnaptr; begin (* destroy previous to dna *) while todna<>nil do cleardna(todna); aDNAptr:=fromdna; if aDNAptr<>nil then begin getdna(memdna); todna:=memdna; while aDNAptr<>nil do begin memdna^.part:=aDNAptr^.part; memdna^.length:=aDNAptr^.length; aDNAptr:=aDNAptr^.next; if aDNAptr<>nil then begin getdna(memdna^.next); memdna:=memdna^.next end else memdna^.next:=nil end end; end; procedure copypiece(Apiece: pieceptr; var Bpiece: pieceptr); (* copy the pice a to piece b *) begin Bpiece^.key.hea.keynam := Apiece^.key.hea.keynam; copyline(Apiece^.key.hea.fulnam, Bpiece^.key.hea.fulnam); copyline(Apiece^.key.hea.note, Bpiece^.key.hea.note); Bpiece^.key.mapbeg := Apiece^.key.mapbeg; Bpiece^.key.coocon := Apiece^.key.coocon; Bpiece^.key.coodir := Apiece^.key.coodir; Bpiece^.key.coobeg := Apiece^.key.coobeg; Bpiece^.key.cooend := Apiece^.key.cooend; Bpiece^.key.piecon := Apiece^.key.piecon; Bpiece^.key.piedir := Apiece^.key.piedir; Bpiece^.key.piebeg := Apiece^.key.piebeg; Bpiece^.key.pieend := Apiece^.key.pieend; copydna(Apiece^.dna,Bpiece^.dna); end; function between(a,b,c: integer): boolean; (* is b between a and c? *) (* this is an inclusive between *) begin between:=((a<=b) and (b<=c)) or ((c<=b) and (b<=a)) end; function within(pie: pieceptr; p: integer): boolean; (* is p (external coordinates) within the piece pie? *) (* note 1: if coocon is linear then piecon must be linear. note 2: does the piece not go over the coordinate boundaries? note 3: if coocon is circular and piecon is circular, then one has the entire piece, so we can ask if p is within the coordinate system. *) begin with pie^.key do begin case coocon of linear: within:=between(piebeg,p,pieend); (* note 1 *) circular: case piecon of linear: case piedir of plus: if pieend>=piebeg (* note 2 *) then within:=between(piebeg,p,pieend) else within:=between(piebeg,p,cooend) or between(coobeg,p,pieend); minus: if pieend<=piebeg (* note 2 *) then within:=between(piebeg,p,pieend) else within:=between(piebeg,p,coobeg) or between(cooend,p,pieend) end; circular: within:=between(coobeg,p,cooend) (* note 3 *) end end end end; function withininternal( pie: pieceptr; p: integer): boolean; (* Is the internal position p inside the piece pie? *) begin withininternal := (p >= 1) and (p <= piecelength(pie)) end; (* end module book.other version = 7.52; {of delmod.p 2000 Jul 30} *) procedure showdnasegment(var f: text; d: dnaptr; spot: integer); (* show the dna segment to file f (for debugging), and mark the given spot by surrouding it with parenthesis *) var j: integer; (* index to the segment *) begin for j := 1 to d^.length do begin if j = spot then write(f,'('); if d^.part[j] in [a,c,g,t] then write(f, basetochar(d^.part[j])) else write(f, ord(d^.part[j]):1); if j = spot then write(f,')'); end; write(f,' '); end; procedure showsegments(var f: text; d: dnaptr); (* show all dna segments to file f (for debugging) *) var i: dnaptr; (* index to the segments *) n: integer; (* count of the segments *) begin i := d; n := 0; while i <> nil do begin showdnasegment(f,i,i^.length); n := succ(n); writeln(f, ' segment ',n:1); i := i^.next end; writeln(f); end; (* begin module book.getbase *) function getbase(position: integer; pie: pieceptr):base; (* 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. *) var workdna: dnaptr; (* pointer to the dna part of pie *) p: integer; (* current count of bases into the workdna *) spot: integer; (* the last base of the dna part *) thelength: integer; (* the length of the piece *) begin { writeln(output,'NEW getbase: position=',position:1,'^^^^^^^^^^^^^^^^^^^^'); } (* handle cases of position out of range by circular wrapping *) thelength := piecelength(pie); while position < 1 do position := position + thelength; while position > thelength do position := position - thelength; workdna:=pie^.dna; p:=workdna^.length; while position > p do begin { writeln(output,' workdna^.length=',workdna^.length:1); } workdna := workdna^.next; if workdna = nil then begin writeln(output,'error in function getbase!'); halt end; p := p + workdna^.length; end; { writeln(output,'p=',p:1); } if workdna = nil then begin writeln(output,'error in getbase: request off end of piece'); halt end else begin spot := workdna^.length - (p-position); { writeln(output,'spot=',spot:1); showdnasegment(output,workdna, spot); } if (spot <= 0) then begin writeln(output,'error in getbase, spot (= ',spot:1, ') must be positive'); halt end; if (spot > workdna^.length) then begin writeln(output,'error in getbase, spot (=',spot:1, ') must be less than length (=',workdna^.length:1,')'); halt end; { writeln(output,'base = ', workdna^.part[spot]); } getbase:=workdna^.part[spot] end end; (* end module book.getbase version = 7.52; {of delmod.p 2000 Jul 30} *) procedure reduceposition(pie: pieceptr; var p: integer); (* reduce the position p into the piece pie *) var size: integer; (* size of a circle *) begin with pie^.key do begin if piecon = circular then begin (* circular correction of p *) size:=cooend-coobeg+1; while p < coobeg do p:=p+size; while p > cooend do p:=p-size end else case piedir of (* linear correction of p brings p to the closest edge of pie *) plus: if pnil then begin l:=freemarker; freemarker:=freemarker^.key.next end else new(l); l^.key.next:=nil end; procedure clearmarker(var l: markerptr); var lptr: markerptr; begin if l<>nil then begin lptr:=l; l:=l^.key.next; lptr^.key.next:=freemarker; freemarker:=lptr end end; procedure getsite(var l: siteptr); begin if freesite<>nil then begin l:=freesite; freesite:=freesite^.next end else new(l); l^.length:=0; l^.next:=nil end; procedure clearsite(var l: siteptr); var lptr: siteptr; begin if l<>nil then begin lptr:=l; l:=l^.next; lptr^.next:=freesite; freesite:=lptr end end; (* ************************************************************************ *) (* ************************************************************************ *) (* ************************************************************************ *) function copylibname(var tofile: text; var lib: text; libnumber: integer) : boolean; (* copy the library name to a file. true returned means that one line was copied *) var ch: char; (* reading character *) begin reset(lib); if not eof(lib) then begin write(tofile,libnumber:1,' '); while not eoln(lib) do begin read(lib,ch); write(tofile,ch) end; readln(lib); writeln(tofile); libline[libnumber]:=2; copylibname:=true end else copylibname:=false end; procedure checklib(var thefile: text; chexpected: char; fi, li: integer); (* check that the character found in the library was the one expected fi is the file number and li is the line number *) (* the global firstch is the character found by reading in the library *) (* no longer valid as of 1999 April 29: this routine must not be removed since it affects thefile. *) begin { readln(thefile,firstch); This readln is now handled by brheader! 1999 April 29 } firstch := thefile^; if chexpected <> firstch then begin writeln(output,' delila: library does not match the catalogue'); writeln(output,' library ',fi:1, ' line ',li:1, ' character expected: ',chexpected, ' character found: ', firstch); writeln(output,' (the order of the libraries and catalogues is', ' probably wrong,'); writeln(output,' or the file is not a library)'); halt end; { zzzbbb no longer valid as of 1999 April 29: libline[fi]:=succ(libline[fi]) } end; (* read library variables *************************************************) (* these procedure read attributes from the library. they are all prefixed by l to indicate this. *) {zzz many of these are now out of date (and not called) and should be removed. However, they are tied to the marker gene and transcript stuff, which is still hanging around. Those should all be replaced by 'features', or maybe gene and transcript by features and markers by changes. } procedure libskipstar(var thefile: text); (* skip start of line (or star = '*'). the first character read is the global firstch. *) var c: char; (* a character in the file *) begin if skipping then begin read(thefile,firstch); if firstch<>'*' then begin writeln(output,' delila: the file lib',catnumber:1,' is bad.'); writeln(output,' either line ',libline[catnumber]:1, ' is missing an asterisk (*) on the attribute'); writeln(output,' or the file is not a delila data base.'); writeln(output,'The library line is:'); write(output,firstch); while not eoln(thefile) do begin read(thefile,c); write(output,c); end; writeln(output); halt end end else skipping:=true; (* turn skipping back on for next time *) read(thefile,blank); (* skip the blank *) end; procedure lrreanum(var thefile: text; var reanum: real); (* read a real number from the file *) begin libskipstar(thefile); readln(thefile,reanum); end; procedure lrnumber(var thefile: text; var num: integer); (* read a number from the file *) begin libskipstar(thefile); readln(thefile,num) end; procedure lrname(var thefile: text; var nam: name); (* read a name from the file *) type (* unpacked alpha for reading from files: *) ualpha = array[1..namelength] of char; var i: integer; (* index to name *) unamlets: ualpha; (* unpacked name *) begin libskipstar(thefile); with nam do begin length:=0; while (not eoln(thefile)) and (length < namelength) do begin length:=succ(length); read(thefile,unamlets[length]) end; pack(unamlets,1,letters); if length'o' then begin writeln(output,'delila: in procedure lrstate: state does not start with o'); halt end; readln(thefile,ch); if ch='n' then sta:=on else sta:=off end; (* read partial keys of library *******************************************) procedure lrlibheader(var lib: text; libnum: integer; var alp: datetimearray); (* read the library header (date of creation) from the library file, lib *) begin (* the header of the library resides in the first line *) reset(lib); if not eof(lib) then begin libskipstar(lib); readdatetime(lib,alp); readln(lib); libline[libnum]:=2 end end; procedure lrheader(var thefile: text; fi: integer; var hea: header); (* read the header of a key. *) var lines: integer; procedure lrnote(var note: lineptr); (* read note key *) var newnote: lineptr; (* the new note *) previousnote: lineptr; (* the last line of the notes *) begin read(thefile,firstch); if firstch='n' then begin (* enter note *) readln(thefile); lines:=succ(lines); read(thefile,firstch); if firstch<>'n' then begin (* abort null note (n/n) *) getline(note); newnote:=note; while firstch<>'n' do begin (* wait until end of note *) skipping:=false; lrline(thefile,newnote); lines:=succ(lines); read(thefile,firstch); previousnote:=newnote; (* get next note *) getline(newnote^.next); newnote:=newnote^.next; end; (* last note was not used, so: *) clearline(newnote); previousnote^.next:=nil; readln(thefile); lines:=succ(lines); end else begin readln(thefile); lines:=succ(lines) end end else skipping:=false end; begin clearheader(hea); with hea do begin (* read key name *) lrname(thefile,keynam); (* read full name *) getline(fulnam); lrline(thefile,fulnam); lines:=2; (* read note key *) lrnote(note); libline[fi]:=libline[fi]+lines end end; procedure lrreference(var thefile: text; fi: integer; var ref: reference); begin with ref do begin lrname (thefile,pienam); lrreanum(thefile,mapbeg); lrdirect(thefile,refdir); lrnumber(thefile,refbeg); lrnumber(thefile,refend); libline[fi]:=libline[fi]+5 end end; (* library file location routine ***************************************) procedure libgetto(var thefile: text; fi, li: integer); (* move thefile to the location requested: file fi and line li *) var currentline: integer; dummyname: datetimearray; begin (* libgetto *) { writeln(output,'libgetto START'); } {zzzbbb} (* get to the line *) (*;if debugging then writeln(debug,'getfileprelibline[',fi:3,']=',libline[fi]:4);*) currentline:= libline [fi]; if currentline > li then begin (* reset and skip library header *) lrlibheader(thefile,fi,dummyname); currentline:=libline[fi] end; while currentline < li do begin { writeln(output,'libgetto: thefile^ = ',thefile^); } {zzzbbb} readln(thefile); currentline:= succ(currentline) end; libline[fi]:=currentline (*if debugging then ;writeln(debug,'getfilepostlibline[',fi:3,']=',libline[fi]:4);*) end; (* libgetto *) (* top level librarian functions ******************************************) (* begin module delila.marksautomate *) procedure marksautomate(var markspots: text); (* define the components necessary for markspots *) begin writeln(markspots,'* marksdelila: define markings for the lister program'); writeln(markspots,'* written by Delila version ',version:4:2); writeln(markspots,'* The standard marks.arrow must be used prior to this file.'); writeln(markspots); writeln(markspots,'u'); writeln(markspots,'% This message will appear if you forgot your'); writeln(markspots,'% marks.arrow definition:'); writeln(markspots,'You_need_a_marks.arrow_definition'); writeln(markspots,''); writeln(markspots,'/setmarkspotarrow{'); writeln(markspots,'/bodycolor {black} def'); writeln(markspots,'/strokecolor {black} def'); writeln(markspots,'/BodyThick 0.30 fs def'); writeln(markspots,'/HeadWidth 0.90 fs def'); writeln(markspots,'/HeadLength 1.50 fs def'); writeln(markspots,'} def'); writeln(markspots,'setmarkspotarrow'); writeln(markspots); writeln(markspots,'/change {% tailx taily headx heady shift change'); writeln(markspots,'% the head of an arrow'); writeln(markspots,'pop'); writeln(markspots,'setmarkspotarrow'); writeln(markspots,'fixedarrow'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/changeworra{% tailx taily headx heady shift changeworra'); writeln(markspots,'% the tail of an arrow is a ''worra'' (spelling backards)'); writeln(markspots,'pop'); writeln(markspots,'setmarkspotarrow'); writeln(markspots,'fixedworra'); writeln(markspots,'} def'); (* no longer needed: *) writeln(markspots); writeln(markspots,'/insertion{% tailx taily headx heady shift insertion'); writeln(markspots,'% an insertion is a green rectangle'); writeln(markspots,'pop'); writeln(markspots,'/bodycolor {lightgreen} def'); writeln(markspots,'boundrectangle'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/deletion {% tailx taily headx heady shift deletion'); writeln(markspots,'% a deletion is a red rectangle'); writeln(markspots,'pop'); writeln(markspots,'/bodycolor {lightred} def'); writeln(markspots,'boundrectangle'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/fullinsertion{% tailx taily headx heady shift insertion'); writeln(markspots,'% an insertion is a green rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightgreen} {black} fullbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/leftinsertion{% tailx taily headx heady shift insertion'); writeln(markspots,'% an insertion is a green rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightgreen} {black} leftbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/midinsertion{% tailx taily headx heady shift insertion'); writeln(markspots,'% an insertion is a green rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightgreen} {black} midbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/rightinsertion{% tailx taily headx heady shift insertion'); writeln(markspots,'% an insertion is a green rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightgreen} {black} rightbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/fulldeletion {% tailx taily headx heady shift deletion'); writeln(markspots,'% a deletion is a red rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightred} {black} fullbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/leftdeletion {% tailx taily headx heady shift deletion'); writeln(markspots,'% a deletion is a red rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightred} {black} leftbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/middeletion {% tailx taily headx heady shift deletion'); writeln(markspots,'% a deletion is a red rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightred} {black} midbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'/rightdeletion {% tailx taily headx heady shift deletion'); writeln(markspots,'% a deletion is a red rectangle'); writeln(markspots,'pop'); writeln(markspots,'{lightred} {black} rightbox'); writeln(markspots,'} def'); writeln(markspots); writeln(markspots,'!'); writeln(markspots); end; (* end module delila.marksautomate version = 1.89; (@ of dbmutate.p 1999 April 26 *) procedure initlibrarian; (* reset library files and set up variables *) var i: integer; { apparently not used! indnum: numberedstructure; (* an index for def.num.str *) } nodefrom, nodeto: node; libcount: integer; (* the number of libraries in use *) procedure checklibcat(var lib: text; libnum: integer; var cat: catfile; catnum: integer); (* check that a library corresponds to its catalogue *) var libheader, catheader: datetimearray; (* the dates of the two *) begin (* checklibcat *) if not eof(lib) then begin if firstlibrary = 0 then firstlibrary := libnum; catnumber := catnum; (* keep the global up to date *) lrlibheader(lib,libnum,libheader); crcatheader(cat,catnum,catheader); if eof(cat) then begin writeln(output, ' delila: missing catalogue ', catnum:1, ' although library ', libnum:1, ' is not empty'); halt end; if libheader <> catheader then begin writeln(output,' delila: library ', libnum:1, ' and catalogue ',catnum:1, ' dates are not identical'); writeln(output,' ',libheader,' <> ',catheader); halt end; libcount := succ(libcount); end else begin (* the library does not exist: does the cat? *) if not eof(cat) then begin writeln(output,' delila: library ',libnum:1, ' is empty, yet catalogue ',catnum:1, ' is not empty'); halt end end end; (* checklibcat *) begin (* initlibrarian *) (* announce ourselves *) writeln(output,'delila ',version:4:2); writeln(output,'sequence segment size (dnamax) is ',dnamax:1); (* get date and time *) getdatetime(datetime); (* initialize library files *) for i:=1 to numlibfil do libline[i]:=1; (* lll *) reset(lib1); reset(lib2); reset(lib3); (* initialize catalogue files *) for i := 1 to numcatfil do catline[i] := 1; catnumber:=1; (* first cat to search *) (* ccc *) reset(cat1); reset(cat2); reset(cat3); reset(inst); rewrite(listing); rewrite(book); if debugging then rewrite(debug); (* note: this variable must be set before the checklib calls *) (* skip both characters before an attribute in the library *) skipping:=true; (* are there instructions? *) if eof(inst) then begin writeln(output,' delila: no instructions'); halt end; (* check that all libraries correspond to their catalogues *) libcount := 0; (* how many libraries are we talking to? *) firstlibrary := 0; (* we have not found the first library yet *) checklibcat(lib1,1,cat1,1); checklibcat(lib2,2,cat2,2); checklibcat(lib3,3,cat3,3); if libcount = 0 then begin writeln(output,'delila: no libraries'); halt end; write(output,libcount:1,' librar'); (* and now for the tail... *) if libcount = 1 then writeln(output,'y') else writeln(output,'ies'); (* initialize free storage *) freeline:=nil; freemarker:=nil; freedna:=nil; freesite:=nil; (* clear memory of library piece *) new(libpie); libpie^.key.hea.fulnam:=nil; libpie^.key.hea.note:=nil; libpie^.dna := nil; getdna(libpie^.dna); (* set up traversal chart *) (* define illegal traversals: *) illegaltraversal:='illegal '; (* fill traversal chart with illegal traversals: *) for nodefrom:=libnode to sitenode do for nodeto:=libnode to sitenode do traversalchart[nodefrom,nodeto]:=illegaltraversal; (* refill the chart with all the legal traversals *) (* lib *) traversalchart[libnode,orgnode]:='o. '; traversalchart[libnode,recnode]:='r. '; (* org *) traversalchart[orgnode,libnode]:='o. '; traversalchart[orgnode,orgnode]:='oo. '; traversalchart[orgnode,chrnode]:='c. '; traversalchart[orgnode,recnode]:='or. '; (* chr *) traversalchart[chrnode,libnode]:='co. '; traversalchart[chrnode,orgnode]:='coo. '; traversalchart[chrnode,chrnode]:='cc. '; traversalchart[chrnode,marnode]:='m. '; traversalchart[chrnode,tranode]:='t. '; traversalchart[chrnode,gennode]:='g. '; traversalchart[chrnode,pienode]:='p. '; traversalchart[chrnode,recnode]:='cor. '; (* mar *) traversalchart[marnode,mardnanode]:='d. '; (* mardna *) traversalchart[mardnanode,libnode]:='dmco. '; traversalchart[mardnanode,orgnode]:='dmcoo. '; traversalchart[mardnanode,chrnode]:='dmcc. '; traversalchart[mardnanode,marnode]:='dmm. '; traversalchart[mardnanode,tranode]:='dmt. '; traversalchart[mardnanode,gennode]:='dmg. '; traversalchart[mardnanode,pienode]:='dmp. '; traversalchart[mardnanode,recnode]:='dmcor. '; (* tra *) traversalchart[tranode,libnode]:='tco. '; traversalchart[tranode,orgnode]:='tcoo. '; traversalchart[tranode,chrnode]:='tcc. '; traversalchart[tranode,marnode]:='tm. '; traversalchart[tranode,tranode]:='tt. '; traversalchart[tranode,gennode]:='tg. '; traversalchart[tranode,pienode]:='tp. '; traversalchart[tranode,recnode]:='tcor. '; (* gen *) traversalchart[gennode,libnode]:='gco. '; traversalchart[gennode,orgnode]:='gcoo. '; traversalchart[gennode,chrnode]:='gcc. '; traversalchart[gennode,marnode]:='gm. '; traversalchart[gennode,tranode]:='gt. '; traversalchart[gennode,gennode]:='gg. '; traversalchart[gennode,pienode]:='gp. '; traversalchart[gennode,recnode]:='gcor. '; (* pie *) traversalchart[pienode,piednanode]:='d. '; (* piedna *) traversalchart[piednanode,libnode]:='dpco. '; traversalchart[piednanode,orgnode]:='dpcoo. '; traversalchart[piednanode,chrnode]:='dpcc. '; traversalchart[piednanode,marnode]:='dpm. '; traversalchart[piednanode,tranode]:='dpt. '; traversalchart[piednanode,gennode]:='dpg. '; traversalchart[piednanode,pienode]:='dpp. '; traversalchart[piednanode,recnode]:='dpcor. '; (* rec *) traversalchart[recnode,orgnode]:='r. '; traversalchart[recnode,orgnode]:='ro. '; traversalchart[recnode,recnode]:='rr. '; traversalchart[recnode,enznode]:='e. '; (* enz *) traversalchart[enznode,sitenode]:='s. '; (* site *) traversalchart[sitenode,libnode]:='ser. '; traversalchart[sitenode,orgnode]:='sero. '; traversalchart[sitenode,recnode]:='serr. '; traversalchart[sitenode,enznode]:='see. '; traversalchart[sitenode,sitenode]:='ss. '; {experimental: allow piece directly since supposedly accession numbers are unique (ha!) zzz traversalchart[libnode,pienode]:='p. '; } (* start library at library (top level) node *) pastlibrary:=libnode; pastbook:= libnode; pastcheck:= libnode; (* set up the characters that correspond to each node *) nodechar[libnode]:= 'l'; nodechar[orgnode]:= 'o'; nodechar[chrnode]:= 'c'; nodechar[marnode]:= 'm'; nodechar[mardnanode]:='d'; nodechar[tranode]:= 't'; nodechar[gennode]:= 'g'; nodechar[pienode]:= 'p'; nodechar[piednanode]:='d'; nodechar[recnode]:= 'r'; nodechar[enznode]:= 'e'; nodechar[sitenode]:= 's'; nodeletter:=nodechar[libnode]; (* for zeroing the coordinate system *) zerobase := 0; zeroshift := 0; (* make sure that there is no previous name: *) clearname(lastpiecename); withused := false; (* control rewrite of marksdelila here *) { old code - delete later - rewrite(marksdelila); marksautomate(marksdelila); } booksize := 0; (* for comparing to maxbook *) end; (* initlibrarian *) (* library reading procedures *********************************************) procedure lrorgkey(nam: name; var org: orgkey); (* read organism key *) var it: item; fi: integer; procedure rorgkey(var thefile: text); begin (* rorgkey *) with org,it do begin libgetto(thefile,fi,line); {zzzbbb} { writeln(output,'lrorgkey'); writeln(output,'thefile^ = ',thefile^); writeln(output,'line = ',line:1); writeln(output,'fi = ',fi:1); } checklib(thefile,nodeletter,fi, line); { writeln(output,'thefile^ = ',thefile^); writeln(output,'about to brorgkey:'); } brorgkey(thefile,libline[fi],org); { writeln(output,'lrorgkey DONE'); } end end; (* rchrkey *) begin (* lrogkey *) (* a new chromosome means new pieces: *) clearname(lastpiecename); finditem(orgnode,nam,it,fi); if itemfound then case fi of (* lll *) 1: rorgkey(lib1); 2: rorgkey(lib2); 3: rorgkey(lib3); end end; (* lrorgkey *) (* upgraded: *) procedure lrchrkey(nam: name; var chr: chrkey); (* read chromosome key *) var it: item; fi: integer; procedure rchrkey(var thefile: text); begin (* rorgkey *) with chr,it do begin { writeln(output,'in rchrkey, about to libgetto'); } libgetto(thefile,fi,line); {zzzbbb} { writeln(output,'rCHRkey fi = ',fi:1); writeln(output,'rCHRkey line = ',line:1); writeln(output,'rCHRkey thefile^ = ',thefile^); writeln(output,'rCHRkey about to checklib'); } checklib(thefile,nodeletter,fi,line); { writeln(output,'in rchrkey, about to brchrkey'); } brchrkey(thefile,libline[fi],chr); { writeln(output,'DONE rCHRkey fi = ',fi:1); writeln(output,'DONE rCHRkey line = ',line:1); writeln(output,'DONE rCHRkey libline[fi] = ',libline[fi]:1); writeln(output,'DONE rCHRkey thefile^ = ',thefile^); } end end; (* rchrkey *) begin (* lrchrkey *) (* a new chromosome means new pieces: *) clearname(lastpiecename); finditem(chrnode,nam,it,fi); if itemfound then case fi of (* lll *) 1: rchrkey(lib1); 2: rchrkey(lib2); 3: rchrkey(lib3); end end; (* lrchrkey *) { not implemented procedure lrmarker(nam: name; var markers: markerptr); (* read a marker, string it into the markers linked list *) var it: item; fi: integer; procedure rmarker(var thefile: text); procedure lrmrky(nam: name; var mar: markey); begin (* lrmrky *) with mar,it do begin checklib(thefile,nodeletter,fi,line); clearline(phenotype); getline(phenotype); lrheader(thefile,fi,hea); lrreference(thefile,fi,ref); lrstate(thefile,sta); lrline (thefile,phenotype); libline[fi]:=libline[fi]+2 end end; (* lrmrky *) procedure lrmardna(var dna: dnaptr); (* read marker dna *) begin (* lrmardna *) tvrslibrary(mardnanode); dna:=nil; (* simple read of all the dna *) (* to write this code, refer to lrpiedna, but simplify it since the whole dna is read, and there is no complementing *) (* dummy *) end; (* lrmardna *) begin (* rmarker *) libgetto(thefile,fi, it.line); (* the linked list of markers should be scanned to find the last marker. the new marker should then be inserted end *) (* alternatively, the new marker should be placed in the marker list in its proper position with respect to the coordinate system *) with markers^ do begin lrmrky(nam,key); lrmardna(dna) end end; (* rmarker *) begin (* lrmarker *) finditem(marnode,nam,it,fi); if itemfound then case fi of (* lll *) 1: rmarker(lib1); 2: rmarker(lib2); 3: rmarker(lib3); end end; (* lrmarker *) } procedure lrtrakey(nam: name; var tra: trakey); (* read transcript key *) var it: item; fi: integer; procedure rtrakey(var thefile: text); begin (* lrtrakey *) with tra,it do begin libgetto(thefile,fi,line); checklib(thefile,nodeletter,fi,line); lrheader(thefile,fi,hea); lrreference(thefile,fi,ref); end end; (* rtrakey *) begin (* lrtrakey *) finditem(tranode,nam,it,fi); if itemfound then case fi of (* lll *) 1: rtrakey(lib1); 2: rtrakey(lib2); 3: rtrakey(lib3); end end; (* lrtrakey *) procedure lrgenkey(nam: name; var gen: genkey); (* read gene key *) var it: item; filenumber: integer; procedure rgenkey(var thefile: text); begin (* rgenkey *) with gen,it do begin libgetto(thefile,filenumber,line); checklib(thefile,nodeletter,filenumber,line); (* to match the new getto rules, readln the file here: *) readln(thefile); libline[filenumber] := libline[filenumber] + 1; lrheader(thefile,filenumber,hea); lrreference(thefile,filenumber,ref); end end; (* rgenkey *) begin (* lrgenkey *) finditem(gennode,nam,it,filenumber); if itemfound then case filenumber of (* lll *) 1: rgenkey(lib1); 2: rgenkey(lib2); 3: rgenkey(lib3); end end; (* lrgenkey *) procedure lrpiece(nam: name; var libpie: pieceptr); (* read a piece *) var it: item; filenumber: integer; (* the file to read *) procedure rpiece(var thefile: text); begin (* rpiece *) libgetto(thefile, filenumber, it.line); (* use standard book reading routines: *) { writeln(output,'lrpiece: line is ',it.line:1); } tvrslibrary(piednanode); (* record that we traversed the tree *) write(output,'Request for: '); writename(output,nam); if equalname(nam, lastpiecename) and (lastpiecename.length > 0) then begin writeln(output, ' - same name as last time, keeping the same piece'); end else begin writeln(output,' - reading in the piece'); brpiece(thefile,libline[filenumber],libpie); copyname(libpie^.key.hea.keynam, lastpiecename); end end; (* rpiece *) begin (* lrpiece *) finditem(pienode,nam,it,filenumber); if itemfound then case filenumber of (* lll *) 1: rpiece(lib1); 2: rpiece(lib2); 3: rpiece(lib3); end end; (* lrpiece *) {not implemented procedure lrreckey(nam: name; var rec: reckey); (* read rececognition key *) var it: item; fi: integer; procedure rreckey(var thefile: text); begin (* rreckey *) with rec,it do begin libgetto(thefile,fi,line); checklib(thefile,nodeletter,fi,line); lrheader(thefile,fi,hea); end end; (* rreckey *) begin (* lrreckey *) finditem(recnode,nam,it,fi); if itemfound then case fi of (* lll *) 1: rreckey(lib1); 2: rreckey(lib2); 3: rreckey(lib3); end end; (* lrreckey *) procedure lrenzyme(nam: name; var enz: enzymeptr); (* read enzyme *) var it: item; fi: integer; procedure renzyme(var thefile: text); procedure lrenzkey(nam:name; var enz: enzkey); begin (* lrenzkey *) with enz,it do begin checklib(thefile,nodeletter,fi,line); lrheader(thefile,fi,hea); end end; (* lrenzkey *) procedure lrsites(var sit: siteptr); (* read dna sites of this enzyme *) begin (* lrsites *) tvrslibrary(sitenode); sit:=nil; (* here we read in all the sites *) (* dummy *) end; (* lrsites *) begin (* renzyme *) with enz^,it do begin libgetto(thefile,fi,line); lrenzkey(nam,key); lrsites(sites) end end; (* renzyme *) begin (* lrenzyme *) finditem(enznode,nam,it,fi); if itemfound then case fi of (* lll *) 1: renzyme(lib1); 2: renzyme(lib2); 3: renzyme(lib3); end end; (* lrenzyme *) } (**************************************************************************) (* book writing procedures ************************************************) (**************************************************************************) (* these are called bw for book writes *) procedure tvrsbook(future: node); (* procedure to print in the book the appropriate tree traversal characters: traverse to the new (future) node *) var thisstep: step; steps: 1..maxstep; begin thisstep:=traversalchart[pastbook,future]; if thisstep=illegaltraversal then begin writeln(output,' delila: program error: illegal book traversal ', ord(pastbook), ' to ', ord(future)); halt end; steps:=1; while thisstep[steps]<>'.' do begin case thisstep[steps] of 'o': begin write(book,'organism'); if not versioninbook then begin write(book,' - book produced by delila ',version:4:2); versioninbook := true; end; writeln(book); end; 'c': writeln(book,'chromosome'); 'm': writeln(book,'marker'); 'd': writeln(book,'dna'); 't': writeln(book,'transcript'); 'g': writeln(book,'gene'); 'p': writeln(book,'piece'); 'r': writeln(book,'recognition-class'); 'e': writeln(book,'enzyme'); 's': writeln(book,'site') end; steps:=succ(steps) end; pastbook:=future end; procedure bwbookheader(var title: lineptr); (* write out the book title *) var libheader: datetimearray; i: integer; begin bwstartline(book); (* date of withdrawal *) writedatetime(book,datetime); (* date of library creation *) write(book,', '); (* finish date of withdrawal *) case firstlibrary of (* lll *) 1: lrlibheader(lib1,1,libheader); 2: lrlibheader(lib2,2,libheader); 3: lrlibheader(lib3,3,libheader); end; writedatetime(book,libheader); (* book title *) if title<>nil then begin if title^.length<>0 then begin write(book,', '); (* finish date of library creation *) for i:=1 to title^.length do write(book,title^.letters[i]); end; clearline(title); end; (* finish the title line *) writeln(book) end; { not implemented procedure bwmarkers(var mar: markerptr); (* write a set of marks (or only one) to the book *) procedure bwmrky(var mar: markey); (* write a marker key *) begin tvrsbook(marnode); with mar do begin bwheader(hea); bwreference(ref); bwstate(sta); bwline (phenotype); end end; procedure bwmardna(dna: dnaptr); (* write marker dna *) begin (* it may be possible to simply move several of the dna writing routines from bwpiedna (grabbase, writebase) into the low level library routines and use those *) (* dummy *) end; begin (* bwmarks *) while mar<>nil do begin if def.key.mar=on then begin with mar^ do begin bwmrky(key); bwmardna(dna) end end; clearmarker(mar) end end; } {not implemented procedure bwtrakey(var tra: trakey); (* write a transcript key *) begin if (def.key.tra = on) then begin tvrsbook(tranode); with tra do begin bwheader(book,hea); bwref(book,ref) end end end; procedure bwgenkey(var gen: genkey); (* write a gene key *) begin if (def.key.gen=on) then begin tvrsbook(gennode); with gen do begin bwheader(book,hea); bwref(book,ref) end end end; } {not implemented procedure bwreckey(var rec: reckey); (* write a recognition key to the book *) begin tvrsbook(recnode); with rec do begin bwheader(book,hea); end end; procedure bwenzyme(var enz: enzymeptr); (* write an enzyme to the book *) procedure bwenzkey(var enz: enzkey); begin tvrsbook(enznode); with enz do begin bwheader(hea); end end; procedure bwsites(var sit: siteptr); (* write a set of sites out to the book *) begin (* dummy *) end; procedure expander(var sit: siteptr); (* recursive copy of each site: when end of copy is reached on any one branch, print *) begin (* dummy *) end; begin (* bwenzyme *) with enz^ do begin bwenzkey(key); if def.sit.expand = on then expander(sites) else bwsites(sites) end end; procedure bwall(n: node; cut: integer); (* write all of the object *) (* first call must be specified, for this recursive routine *) begin (* dummy *) end; procedure bwevery(n: node); (* write out every n within the currently specified structure *) begin (* dummy *) end; } (**************************************************************************) (**************************************************************************) (**************************************************************************) {OOO} (* begin module delila.writechangeset.describechangeset *) procedure wchangedata(var f: text; c: changedata); (* write the changedata to f in shorthand notation *) var i: integer; (* index to insertion *) begin with c do case changetype of 'c': write(f,baseold,basecoo1:1,basenew); 'i': begin write(f,'i',basecoo1:1,',',basecoo2:1); if inserts > 0 then for i := 1 to inserts do write(f,insert[i]) (* Don't do this here because the changeset as written is like a name: else write(f,'NOTHING'); *) end; 'd': write(f,'d',basecoo1:1,',',basecoo2:1) end; end; procedure writechangeset(var f: text; changes: changeset); (* write the changeset to file f in shorthand notation *) var n: integer; (* index to changes *) begin with changes do for n := 1 to number do begin if n > 1 then write(f,'.'); wchangedata(f, data[n]); { with data[n] do case changetype of 'c': write(f,baseold,basecoo1:1,basenew); 'i': begin write(f,'i',basecoo1:1,',',basecoo2:1); if inserts > 0 then for i := 1 to inserts do write(f,insert[i]) (* Don't do this here because the changeset as written is like a name: else write(f,'NOTHING'); *) end; 'd': write(f,'d',basecoo1:1,',',basecoo2:1) end; } end; end; procedure checkchangeset(var f: text; changes: changeset); (* write the changeset to file f in shorthand notation for INTERNAL coordinates *) var i: integer; (* index to insertion *) n: integer; (* index to changes *) begin with changes do for n := 1 to number do begin if n > 1 then write(f,'.'); with data[n] do case changetype of 'c': write(f,baseold,internal1:1,basenew); 'i': begin write(f,'i',internal1:1,',',internal2:1); if inserts > 0 then for i := 1 to inserts do write(f,insert[i]) (* Don't do this here because the changeset as written is like a name: else write(f,'NOTHING'); *) end; 'd': write(f,'d',internal1:1,',',internal2:1) end; end; end; procedure describechangeset(var f: text; changes: changeset); (* describe in English the changeset to file f *) var i: integer; (* index to insertion *) n: integer; (* index to changes *) begin if changes.number = 0 then write(f,'no changes') else begin with changes do for n := 1 to number do begin if n > 1 then write(f,', '); with data[n] do case changetype of 'c': write(f,'at ',basecoo1:1,' ',baseold,'->',basenew); 'i': begin write(f,'insert '); if inserts > 0 then for i := 1 to inserts do write(f,insert[i]) else write(f,'NOTHING'); write(f,' between ',basecoo1:1,' and ',basecoo2:1); end; 'd': write(f,'delete ',basecoo1:1,' to ',basecoo2:1) end; end; end; end; (* end module delila.writechangeset.describechangeset *) (* NOTE: SEE ALSO: linechangeset for output to the book *) {OOO} {OOO} {OOO} (**************************************************************************) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* instruction interpretation procedures **********************************) (**************************************************************************) function tvrschecks(futurecheck: node):boolean; (* this routine is to be used to check that the instructions are in proper order: 'if tvrschecks(this-try) then continue-analysis-of-instructions;' in three cases we must advance (force) the node since no reads are done at this time. (library reads do a force for the other traversal routines) *) begin tvrschecks:= traversalchart[pastcheck,futurecheck]<>illegaltraversal; if futurecheck=marnode then pastcheck:=mardnanode else if futurecheck=pienode then pastcheck:=piednanode else if futurecheck=enznode then pastcheck:=sitenode else pastcheck:=futurecheck end; (**************************************************************************) (**************************************************************************) procedure librarian; (* by paul morrissett readinstruction modified, and supporting routines written by tom schnieder *) (* there are two passes made through the delila instructions. pass 1 reads the code using ir routines (which all rely on procedure ichread to read the instructions). nothing happens when delila gets to a 'peak' in the code (after passing all the if statements), since this represents a syntactically correct instruction. (in other words, the fact that one got to a certain spot in the code means that so far the instruction is correct.) ichread is responsible for writing the lines into the listing during this pass. instructions are written as they are scanned. all writes are done inside the read procedures. syntax errors are pointed to. the title is read, if it is in the instructions. pass 2 runs only if pass 1 had no errors. the instructions are reread, and variables are collected. calls to the library read, and the book write routines are then made. numerical values and numbers of items are pointed to. note on irs and the parser: these routines assume that the previous actions left the inst file in good order. this means that they have read past all they needed, and they have read their parse stop symbol. they need not have done a findword however, so this is often required as the first action of an ir. (extra findword calls will not hurt.) *) const pagelength = 57; (* non-header lines on a page of listing. *) widinst = namelength; (* maximum width in characters of an instruction word *) minword = 2; (* the minimum length delila word *) max1errors = 33; (* maximum error number in pass 1 *) {zzz111} {EEE} max2errors = 223; (* maximum error number in pass 2 *) maxerrors = 10; (* listed per line in listing *) maxpositions = 15; (* the maximum number of position numbers per line *) maxnumbers = 15; (* the maximum number of numberings per line *) numberlength = 2; (* the length of the array numberword *) type delilaword = ( (* delila instruction types *) alldelila, arrdelila, (* 1999 Mar 20 *) begdelila,chrdelila,cledelila,comdelila,condelila,coodelila, cutdelila,defdelila,dirdelila, doudelila, (* 1999 Mar 20 *) elsdelila,enddelila,enzdelila, evedelila,expdelila,frodelila,gendelila,getdelila, haldelila,homdelila, ifdelila,keydelila,mardelila,moddelila,namdelila, notdelila,numdelila,offdelila, ondelila,orgdelila,outdelila, piedelila, recdelila, reddelila, samdelila, (* 1996 Aug 12 *) setdelila, (* 1999 Mar 20 *) sitdelila,sizdelila,thedelila,titdelila, todelila,tradelila,witdelila, nordelila,zerdelila, eodelila); (* instruction word spelling: *) instword = packed array[1..widinst] of char; elnrecord = record (* errors on a line record *) pos: array[1..maxerrors] of integer; err: array[1..maxerrors] of integer end; plnrecord = record (* dna position numbers on a line record *) pos: array[1..maxpositions] of integer; psn: array[1..maxpositions] of integer; rea: array[1..maxpositions] of real; (* real numbers on the line *) end; nlnrecord = record (* numbers on a line record *) pos: array[1..maxnumbers] of integer; num: array[1..maxnumbers] of integer end; var pass: 1..2; (* the pass through the instructions *) pageline, (* lines printed on this listing page *) pagenumber: integer; (* the page in listing *) pass1errors:boolean; (* true if there were any errors in pass 1 *) pass2errors:boolean; (* true if there were any errors in pass 2 *) warnings: boolean; (* true for warnings in pass 2 *) error1list: array[ 0..max1errors] of boolean; (* pass 1 errors *) error2list: array[201..max2errors] of boolean; (* pass 2 errors *) numoferrors: integer; (* current number of errors per line *) errorline: elnrecord; (* the errors on the current line *) (* all of the words recognized by delila: *) wordlist: array[delilaword] of instword; numbers: set of '0'..'9'; (* parsed parts of the instructions *) chr: char; (* storage for characters read from file inst *) parsedword: instword; (* the latest word parsed from the instructions *) parsedlength: integer; (* the length of parsedword *) parsedlocation: integer; (* current location in parsedword *) word: delilaword; (* result of irword: the translation into type delila of an instruction word parsed by spaces *) inumber: integer; (* an integer, result of routine irnumber *) rnumber: real; (* a real number, result of routine irnumber *) keyname: name; (* a key name *) save: delilaword; (* to save the previous instruction word *) usernotes: lineptr; (* notes from the user to be inserted into the next header specified *) title: lineptr; (* the title of the book *) titleexists: boolean; (* true if a title was found in the first pass - a requirement *) instructioncount: integer; (* count number of instructions *) linecount: integer; (* count number of delila instruction lines *) lineposition: integer; (* keeps track of where chr is being read from a line in file inst. used by ichread and error routines only *) correct: boolean; (* set by many routines as output to show when the operation was successful *) eoinst: boolean; (* end of instruction- set to true by any ir routine that finds a ';' *) parentheses: integer; (* the number of open parentheses *) comment: boolean; (* true if we are inside a comment *) commentline: integer; (* the linecount where a comment was found *) quote: boolean; (* true if we are inside a quote *) linebreak: boolean; (* true when at end of a line in inst eof(inst) usually will not work since ichread does an automatic readln(inst) *) significant: boolean; (* a significant non-blank, non- comment letter in the instructions is now in chr *) numberword: packed array[1..numberlength] of char; (* the word 'number ', or the symbol '# '. this is the symbol written to the book to number items. used in procedure specify *) (* what items are currently specified: 0.5 = unspecified -0.5 = specified, not numbered integer = specified, numbered *) okspec,ckspec, mkspec,tkspec,gkspec,pkspec, rkspec,ekspec: real; (* interface between the library and the book. these variables temporarily store data from the library. *) ok: orgkey; ck: chrkey; mkoff: markerptr; mkon: markerptr; tk: trakey; gk: genkey; pk: pieceptr; (* a piece of DNA sequence *) rk: reckey; ek: enzymeptr; (* variables that determine a request for a piece of dna *) fromposition, (* the first base desired *) toposition, (* the last base desired *) (* the first base to provide from a circular piece of dna: *) cutposition: integer; dirwanted: direction; (* orientation of the dna *) (* pass 2 variables *) numofpositions: integer; (* number of position numbers on a line *) positionline: plnrecord; (* the dna positions on the current line *) numofnumbers: integer; (* the number of numbers on a line *) numofchanges: integer; (* the number of changes on a line, used to keep track of whether the mutations have been listed yet *) numberline: nlnrecord; (* the numbers on the current line *) { not implemented choice: delilaword; (* the choice of an . see also . *) } (* this is a set of those delilawords which represent structures in the library: *) structure: set of alldelila..eodelila; (* the previous from position for the *) previousfromposition: integer; sameusageisvalid: boolean; (* if true, use of 'same' is valid *) (* look for cases where both ends of a sequence have been reduced; warn the user in this case. *) reduced: boolean; (* amount to indent before writing information to listing *) indentlisting: integer; coordinateside: direction; (* which side of the coordinate system extra sequence should be added or deleted *) getcount: integer; (* count of the number of pieces gotten *) procedure startheader(var hea: header); (* initialize the values of a header *) begin (* starthea *) with hea do begin fulnam := nil; note := nil end end; (* starthea *) procedure initializedelila; (* initialize the variables of delila *) var i: integer; begin versioninbook:=false; (* version has not been written yet *) pagenumber:=0; (* first page is 1 *) (* set up conversion between letters and delilawords: *) (* 123456789 123456789 123456789 123456789 123456789 *) wordlist[alldelila]:='all '; wordlist[arrdelila]:='arrowlength '; wordlist[begdelila]:='beginning '; wordlist[chrdelila]:='chromosome '; wordlist[cledelila]:='cleave '; wordlist[comdelila]:='complement '; wordlist[condelila]:='continue '; wordlist[coodelila]:='coordinate '; wordlist[cutdelila]:='cut '; wordlist[defdelila]:='default '; wordlist[dirdelila]:='direction '; wordlist[doudelila]:='doubling '; wordlist[elsdelila]:='else '; wordlist[enddelila]:='ending '; wordlist[enzdelila]:='enzyme '; wordlist[evedelila]:='every '; wordlist[expdelila]:='expand '; wordlist[frodelila]:='from '; wordlist[gendelila]:='gene '; wordlist[getdelila]:='get '; wordlist[haldelila]:='halt '; wordlist[homdelila]:='homologous '; wordlist[ ifdelila]:='if '; wordlist[keydelila]:='key '; wordlist[mardelila]:='marker '; wordlist[moddelila]:='modify '; wordlist[namdelila]:='name '; wordlist[notdelila]:='note '; wordlist[numdelila]:='numbering '; wordlist[offdelila]:='off '; wordlist[ ondelila]:='on '; wordlist[orgdelila]:='organism '; wordlist[outdelila]:='out-of-range '; wordlist[piedelila]:='piece '; { wordlist[recdelila]:='recognition-class '; } wordlist[reddelila]:='reduce-range '; wordlist[samdelila]:='same '; wordlist[setdelila]:='set '; wordlist[sitdelila]:='site '; wordlist[sizdelila]:='size '; wordlist[thedelila]:='then '; wordlist[titdelila]:='title '; wordlist[ todelila]:='to '; wordlist[tradelila]:='transcript '; wordlist[witdelila]:='with '; wordlist[zerdelila]:='zero '; wordlist[nordelila]:='normal '; wordlist[ eodelila]:='end-of-delila-words '; (* no errors have been detected yet: *) pass1errors:=false; pass2errors:=false; for i:= 0 to max1errors do error1list[i]:=false; for i:=201 to max2errors do error2list[i]:=false; (* no special conditions yet: *) warnings:=false; title:=nil; titleexists := false; longname:=nil; usernotes:=nil; numberword:='# '; numbers:=['0','1','2','3','4','5','6','7','8','9']; structure:=[orgdelila, chrdelila, mardelila, tradelila, gendelila, piedelila {not implemented , recdelila, enzdelila} ]; getmarker(mkoff); mkon:=nil; new(pk); new(ek); (* start all lines at nil *) startheader(ok.hea); ok.mapunit:=nil; startheader(ck.hea); startheader(mkoff^.key.hea); mkoff^.key.phenotype:=nil; mkoff^.key.next:=nil; startheader(tk.hea); startheader(gk.hea); startheader(pk^.key.hea); startheader(rk.hea); startheader(ek^.key.hea); (* unspecify all the items. see procedure unspec *) okspec:=0.5; ckspec:=0.5; mkspec:=0.5; tkspec:=0.5; gkspec:=0.5; pkspec:=0.5; rkspec:=0.5; ekspec:=0.5; indentlisting:=2*instwidth + 4; (* extra space after each inst number *) coordinateside := plus; getcount := 0; end; (* initializedelila *) procedure pageheader; (* write delila page header *) begin pageline:=1; pagenumber:=succ(pagenumber); write(listing,' '); writedatetime(listing,datetime); write(listing,' delila ', version:4:2, ' pass ', pass:1, ' page ',pagenumber:1); writeln(listing); writeln(listing); pageline:= pageline+2; end; procedure startlistingline; (* start an instruction listing line *) begin if not eof(inst) then begin linecount:=succ(linecount); write(listing,' ', linecount:instwidth, ' ', instructioncount:instwidth, ' ') end end; procedure initreadinst; (* prepare to read the instructions starting from the top *) begin eoinst:=false; (* not at end of instruction *) parentheses:=0; (* no open ( yet *) comment:=false; (* not inside a comment *) quote:=false; (* not inside a quote *) chr:=' '; instructioncount:=1; linecount:=0; lineposition:=0; (* before the zeroth character *) numoferrors:=0; numofpositions:=0; numofnumbers:=0; numofchanges:=0; (* make sure that we are normal coordinate system by default *) def.coo := coornormal; pageheader; writeln(listing,' Parent Grand Parent', ' Library'); write (listing,'# Library Date/Time: '); write (listing,' Library Date/Time:'); write (listing,' Title:'); writeln(listing); pageline:= pageline+2; (* lll *) if copylibname(listing,lib1,1) then pageline:=succ(pageline); if copylibname(listing,lib2,2) then pageline:=succ(pageline); if copylibname(listing,lib3,3) then pageline:=succ(pageline); writeln(listing); writeln(listing,' inst statement'); writeln(listing,' line number '); pageline:= pageline+3; startlistingline; longnameexists := false; (* make sure that there is no previous name from the first pass: *) clearname(lastpiecename); (* make sure that there are no mutations from the first pass: *) mutations.number := 0; (* force *) (* clear the pieces *) pk^.dna := nil; libpie^.dna := nil; clearpiece(pk); clearpiece(libpie); (* initialize defaults *) with def do begin with key do begin note:=on; mar:=on; gen:=on; tra:=on end; with sit do begin expand:=on; modify:=off; cleave:=off end; defout:=rhalt; with num do begin sta:=on; (* turn OFF all numbering: *) for indnum:=orgnum to enznum do def.num.str[indnum]:=off; (* turn on only pieces initially [1997 Jan 10] *) def.num.str[pienum]:=on; item:=0; (* first value will be 1 *) end; doubling := off; arrowlength := 1.5; end; (* default values for variables for making marksdelila file *) insertlowerbits := -1.3; (* lowerbits for insertion symbol *) insertupperbits := -0.1; (* upperbits for insertion symbol *) deletelowerbits := -1.3; (* lowerbits for deletion symbol *) deleteupperbits := -0.1; (* upperbits for deletion symbol *) changelowerbits := +0.3; (* lowerbits for change symbol *) (* upperbits for change symbol: *) changeupperbits := def.arrowlength+changelowerbits; end; (* initreadinst *) function awarning(e: integer): boolean; (* this defines which errors are actually warnings *) begin awarning:= (e=206) or (e=208) or (e=209) or (e=210) or (e=212) or (e=213) or (e=215) or (e=216) or (e=219) or (e=220) or (e=221) or (e=3) end; procedure writepasserrors(var listing: text; pass: integer); (* write out what the error numbers marked in the listing correspond to in english *) var e: integer; (* index for errors *) i: integer; (* index to write piece key name *) warn1: boolean; (* warnings in pass 1 must be listed for pass 2... *) begin (* writepasserrors *) if (pass = 1) and (not pass1errors) then begin writeln(listing,'No syntax errors found.'); end; if (pass = 2) and (not pass2errors) then begin writeln(listing,'No extraction errors found.'); end; warn1 := ((pass = 2) and (error1list[3])); if (pass = 1) or warn1 then begin if pass1errors or warn1 then begin writeln(listing,'Key to error and warning numbers:'); for e:=0 to max1errors do begin if error1list[e] then begin if awarning(e) then write(listing,'WARNING! '); write(listing,e:4,': '); case e of 0: write(listing,'I do not know how to do this instruction yet, sorry'); 1: write(listing,'Instruction longer than expected (missing semicolon)'); 2: write(listing,'I can not recognize this word'); 3: write(listing,'Warning: this title was ignored'); 4: write(listing,'I can not recognize this integer'); 6: write(listing,'You are missing a specification (illegal tree traversal)'); 7: write(listing,'This word is not allowed here'); 8: write(listing,'A cut is only allowed with a piece'); 9: write(listing,'This word is too long for me (>',widinst:3,' chars)'); 10: begin writeln(listing,'Unclosed comment found at end of instructions.'); write (listing,' ':6, 'The comment started on line ',commentline:1,'.'); end; 11: write(listing,'Unclosed ( or ) in this instruction'); 12: write(listing,'Unclosed quote'); 13: write(listing,'No closing ;'); 14: write(listing,'This key name is too long (>',namelength:3,' chars)'); 15: write(listing,'The statement ended before I expected it to'); 16: write(listing,'Quote expected'); 17: write(listing,'A title is required'); 18: write(listing,'"from same" is not allowed; use it only as "to same"'); 19: begin writeln(listing,'mutation: Change must be identified by one of: acgtdi'); write (listing,' ':6, ' Instead, an illegal change character was', ' found'); end; 20: write(listing,'Old base must be: a, c, g, t'); 21: write(listing,'New base must be: a, c, g, t'); 22: write(listing,'Insertion bases must be one of a, c, g, or t.'); 23: write(listing,'A comma (,) is expected between coordinates for deletion.'); 24: write(listing,'No more than ',insertmax:1,' insertion bases allowed,', ' increase constant insertmax'); 25: write(listing,'A comma (,) is expected between coordinates for insertion.'); 26: write(listing,'A number was expected but no digits were found.'); 27: begin writeln(listing,'First coordinate number must not be larger than second:'); write(listing, ' reverse their order. See libdef for the reason.'); end; 28: write(listing,'Too many changes requested, increase constant changesetmax.'); 29: write(listing,'Extra dots not allowed.'); 30: write(listing,'Unidentified change command.'); 31: write(listing,'Numbers cannot have more than 1 sign (+ or -).'); 32: write(listing,'Insertion complement symbol must be before insertion sequence.'); 33: write(listing,'Non-printable ASCII character found. Eg: tabs not allowed.'); {zzz111} {EEE} end; writeln(listing) end end end end; { else (* if pass = 2 then *) begin } (* 2001 Mar 28 pass 2 errors must be written independently of pass 1 so that the warning about double titles does not block errors in pass 2 from being listed *) if pass = 2 then begin if pass2errors or warnings then begin writeln(listing,'Key to error or warning numbers:'); for e:=201 to max2errors do begin if error2list[e] then begin write(listing,e:4,': '); if awarning(e) then write(listing,'WARNING! '); case e of 201: write(listing,'I can not find this item in the library'); 202: write(listing,'This item was not previously specified'); 203: write(listing,'Out of range and default range = halt'); 204: write(listing,'Positions not consistent with requested direction'); 205: write(listing,'This thing was not on the previously specified piece'); 206: write(listing,'We do not know this limit'); 207: write(listing,'For cutting, the piece must be circular. it is linear'); 208: write(listing,'Out of range and default range = reduce'); 209: write(listing,'Out of range and default range = continue'); 210: begin write(listing,'Piece had two end reductions: length will be 1 base!'); write(output,'Piece '); with pk^.key.hea.keynam do begin for i:=1 to length do write(output,letters[i]); end; writeln(output,' had two end reductions: length will be 1 base!'); end; 211: write(listing,'The base you want to mutate is not what you said it is'); 212: write(listing,'mutation: inserted sequence alters coordinate system'); 213: write(listing,'mutation: deleted sequence alters coordinate system'); 214: begin writeln(listing,'mutation: the initial and final bases', ' are the same,'); write(listing, ' so you did not request any change!'); end; 215: write(listing,'mutation: Requested coordinate off piece', ' in 5'' direction'); 216: write(listing,'mutation: requested coordinate off piece', ' in 3'' direction'); { no longer valid: 217: write(listing,'mutation: The requested coordinates are out of order.'); } 218: write(listing,'mutation: The requested coordinates cannot be equal.'); 219: write(listing,'There will be no change to the sequence!'); 220: write(listing,'The entire sequence was deleted! NO PIECE WAS OUTPUT!'); 221: write(listing,'Deletion outside of sequence will have no effect!'); 222: write(listing,'Overlapping changes are not allowed.'); 223: write(listing,'Booksize would exceed ',maxbook:1, ' bases.'); end; writeln(listing) end end end end; writeln(listing); end; procedure bookhalt; (* put a 'halt' at the top of the book, and warn the user *) var bookcopy: text; procedure copyfile(var fromfile, tofile: text); (* copy one file into the other. no resets or rewrites are done *) var ch: char; begin while not eof(fromfile) do begin while not eoln(fromfile) do begin read (fromfile,ch); write( tofile,ch); end; readln (fromfile); writeln( tofile); end end; begin (* bookhalt *) if pass=2 then begin (* copy away the book *) reset(book); rewrite(bookcopy); copyfile(book,bookcopy) end; rewrite(book); writeln(book,'halt: error in Delila pass ',pass:1); if pass=2 then begin reset(bookcopy); copyfile(bookcopy,book) end end; procedure iwritenumber(c: char; n: integer); (* write a number n to listing, move the line position *) (* the number is proceeded by the character c *) var spots, (* positions moved *) absn, (* absolute value of n *) power: (* a power of 10 *) integer; begin if n>=0 then begin spots:=2; absn:=n end else begin spots:=3; absn:=-n end; power:=10; (* perform log function *) while absn>=power do begin power:=power*10; spots:=succ(spots) end; write(listing,c,n:(spots-1)); lineposition:=lineposition+spots end; procedure rwritenumber(c: char; n: real); (* write a real number n to listing, move the line position *) (* the number is proceeded by the character c *) var spots, (* positions moved *) power: (* a power of 10 *) integer; absn: real; (* absolute value of n *) begin if n>=0 then begin spots:=2; absn:=n end else begin spots:=3; absn:=-n end; power:=10; (* perform log function *) while absn>=power do begin power:=power*10; spots:=succ(spots) end; spots := spots+2; (* decimal and two digits *) write(listing,c,n:(spots-1):2); lineposition:=lineposition+spots end; (* error routines *********************************************************) procedure error(err:integer); (* enter an error in to 'errorline' recording the error number and where it occured on the line (ie. lineposition) since pass 2 implies a perfect pass 1, this routine is (usually) only active during pass 1. so in pass 2 this routine checks that the code is running correctly. error sets the appropriate passerror:=true correct is set *) begin if awarning(err) then writeln(output,'WARNING ',err:1) else writeln(output,'ERROR ',err:1); {crash;} {zzzfff} (* make warnings for some numbers, death for others *) correct:=false; if awarning(err) then begin warnings:=true; correct:=true; end; numoferrors:=succ(numoferrors); if numoferrors<=maxerrors then begin errorline.pos[numoferrors]:=lineposition; errorline.err[numoferrors]:=err; if pass=1 then begin pass1errors:=pass1errors or (not correct); error1list[err]:=true end else begin (* pass=2 *) pass2errors:=pass2errors or (not correct); error2list[err]:=true end end ;if debugging then writeln(debug,'error(',err:4,')'); end; (* begin module delila.plural *) procedure plural(var thefile: text; number: integer; blank: char); (* if the number is not 1, return an s. Otherwise return the blank character *) begin (* plural *) if number <> 1 then write(thefile,'s') else write(thefile,blank) { else write(thefile,'-') } end; (* plural *) (* end module delila.plural version = 2.09; (@ of rembla.p 1995 dec 21 *) procedure writeerrors; (* write out errors for a line from 'errorline'. clear errorline. indicate location of errors with a ^ *) var {EEPP} another: integer; (* another message is needed *) e: integer; (* counter of number of errors *) errornumber: integer; (* the ith error *) errorsonline: integer; (* number of errors found on a line *) warningsonline: integer; (* number of warnings found on a line *) begin if numoferrors <> 0 then begin { writeln(listing,'123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789'); } errorsonline := 0; warningsonline := 0; for e := 1 to numoferrors do begin if not awarning(errorline.err[e]) then errorsonline := errorsonline + 1 else warningsonline := warningsonline + 1; { writeln(listing,'e = ',e:1); writeln(listing,'errorsonline = ',errorsonline:1); writeln(listing,'warningsonline = ',warningsonline:1); } end; if errorsonline > 0 then begin write(listing,' ---error'); plural(listing,errorsonline,'-'); write(listing, '------') (* 16 characters *) end else begin write(listing,' ---warning'); plural(listing,warningsonline,'-'); write(listing, '----') (* 16 characters *) end; { write(listing,' ---error(s)----'); (* 16 characters *) } lineposition:=1; errornumber:=1; while (numoferrors>0) and (errornumber<=maxerrors) do begin if lineposition >= errorline.pos[errornumber] then begin iwritenumber('^',errorline.err[errornumber]); if errorline.err[errornumber] in [19,21,22,23,25,30,211,215,216,217,218,222] {EEE} then another := errorline.err[errornumber]; errornumber:=succ(errornumber); numoferrors:=pred(numoferrors) end else begin lineposition:=succ(lineposition); write(listing,'-') end end; if numoferrors>0 then write(listing,' and other errors'); numoferrors:=0; writeln(listing); { writeln(listing); writeln(listing,'lineposition = ',lineposition:1); write (listing,' '); writeln(listing,'123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789'); write (listing,' '); writeln(listing,' 1 2 3 4 5 6 7 8'); } pageline:=succ(pageline); (* special additional information for mutations *) if another = 19 then begin writeln(listing,' Change must be identified by one of: acgtdi'); writeln(listing,' Instead, the illegal change character "', mutischar, '" was found'); pageline:=succ(pageline); end; if another = 21 then begin writeln(listing,' '); writeln(listing,' Change bases must be one of a, c, g, or t.'); writeln(listing,' Instead, "',mutischar,'" was found'); pageline:=succ(pageline); end; if another = 22 then begin writeln(listing,' Insertion bases must be one of a, c, g, or t.'); writeln(listing,' Instead, "',mutischar,'" was found'); pageline:=succ(pageline); end; if another = 25 then begin writeln(listing,' A comma (,) must separate insertion parts.'); writeln(listing,' Instead, "',mutischar,'" was found'); pageline:=succ(pageline); end; if another = 23 then begin writeln(listing,' A comma (,) must separate deletion parts.'); writeln(listing,' Instead, "',mutischar,'" was found'); pageline:=succ(pageline); end; if another = 30 then begin write (listing,' Unidentified change command:'); writeln(listing,' "',mutischar,'"'); pageline:=succ(pageline); end; if another = 211 then begin writeln(listing,'On the positively oriented strand,', ' the old base at ',mutposition1:1, ' is NOT ',mutnotchar,'!', ' It is ',mutischar,'.'); pageline:=succ(pageline); end; if another = 215 then begin writeln(listing,'The requested coordinate', ' at ',mutposition1:1, ' is off the piece', ' in the 5'' direction'); pageline:=succ(pageline); end; if another = 216 then begin writeln(listing,'The requested coordinate', ' at ',mutposition1:1, ' is off the piece', ' in the 3'' direction'); pageline:=succ(pageline); end; if another = 217 then begin writeln(listing,'The requested coordinates', ' ',mutposition1:1, ' and ',mutposition2:1, ' are out of order.'); pageline:=succ(pageline); end; if another = 218 then begin writeln(listing,'The requested coordinates', ' "',mutposition1:1, ',',mutposition2:1, '" cannot be equal.'); pageline:=succ(pageline); end; {EEE} if another = 222 then begin write(listing,'Overlapping changes: '); wchangedata(listing,mutcd1); write(listing,' overlaps '); wchangedata(listing,mutcd2); writeln(listing); {zzzppp} pageline:=succ(pageline); end; end end; (* pass 2 number routines ***********************************************) procedure ivaluenumber(p: integer); (* put a number into the number line array like procedure error *) begin numofnumbers:=succ(numofnumbers); if numofnumbers<=maxnumbers then begin numberline.pos[numofnumbers]:=lineposition-1; numberline.num[numofnumbers]:=p end ;if debugging then writeln(debug,'ivaluenumber',p:5); end; procedure writenumbers; (* write the numbers of number line out like procedure writeerrors *) var i: integer; number: integer; (* the ith number *) begin if numofnumbers <> 0 then begin for i:=1 to indentlisting do write(listing,' '); lineposition:=0; number:=1; while (numofnumbers>0) and (number<=maxnumbers) do begin if lineposition >= numberline.pos[number] then begin iwritenumber('#',numberline.num[number]); number:=succ(number); numofnumbers:=pred(numofnumbers); end else begin lineposition:=succ(lineposition); write(listing,' ') end end; if numofnumbers>0 then write(listing,'and others (sorry)'); numofnumbers:=0; writeln(listing); pageline:=succ(pageline) end end; (**************************************************************************) (**************************************************************************) (**************************************************************************) {OOO} {OOO} {OOO} {OOO} {OOO} {OOO} (**************************************************************************) (* pass 2 position routines ***********************************************) (**************************************************************************) procedure ivalueposition(p: integer); (* put a position into the position line array like procedure error *) begin numofpositions:=succ(numofpositions); if numofpositions<=maxpositions then begin positionline.pos[numofpositions]:=lineposition-1; positionline.psn[numofpositions]:=p; positionline.rea[numofpositions]:=p; end ;if debugging then writeln(debug,'ivalueposition',p:5); end; procedure rvalueposition(r: real); (* put a real number into the position line array like procedure error *) begin numofpositions:=succ(numofpositions); if numofpositions<=maxpositions then begin positionline.pos[numofpositions]:=lineposition-1; (* flag for real number: *) positionline.psn[numofpositions]:=round(r)+10; positionline.rea[numofpositions]:=r end end; procedure writepositions; (* write the positions of position line out like procedure writeerrors *) var i: integer; number: integer; (* the ith number *) begin if numofpositions <> 0 then begin for i:=1 to indentlisting do write(listing,' '); lineposition:=1; number:=1; while (numofpositions>0) and (number<=maxpositions) do begin if lineposition >= positionline.pos[number] then begin if ( positionline.psn[number] <> round(positionline.rea[number])) then rwritenumber('^',positionline.rea[number]) else iwritenumber('^',positionline.psn[number]); number:=succ(number); numofpositions:=pred(numofpositions); end else begin lineposition:=succ(lineposition); write(listing,' ') end end; if numofpositions>0 then write(listing,'and others (sorry)'); numofpositions:=0; writeln(listing); pageline:=succ(pageline) end end; procedure writechanges; (* write out mutational changes about the instruction *) var i: integer; (* index for writing spaces *) begin if numofchanges > 0 then begin write(listing,'Mutation'); if mutations.number > 1 then write(listing,'s:') else write(listing,': '); (* +6 is normal spacing, -11 is "Mutations: " *) for i:=1 to ((instwidth*2)+6-12) do write(listing,' '); describechangeset(listing, mutations); writeln(listing); end; numofchanges:=0; end; procedure writelineinformation; (* write out information about the previous line *) begin writenumbers; writepositions; writechanges; end; (**************************************************************************) procedure readinstruction; (* read and process one instruction up to the next ';' *) (**************************************************************************) (* specification procedures ***********************************************) (* low level routines *****************************************************) var indnum: numberedstructure; (* an index for def.num.str *) function noder(var word: delilaword): node; (* convert a delila to a node *) begin case word of orgdelila: noder:=orgnode; chrdelila: noder:=chrnode; mardelila: noder:=marnode; tradelila: noder:=tranode; gendelila: noder:=gennode; piedelila: noder:=pienode; recdelila: noder:=recnode; enzdelila: noder:=enznode; end; end; procedure ichread; (* read a single character from the instruction file into chr *) (* this and error routines are the only routines which know where we are in the instruction file *) begin linebreak:=false; (* old loop: while (eoln(inst) and (not eof(inst))) do begin *) (* new loop version: *) if not eof(inst) then if eoln(inst) then repeat if debugging then writeln(listing,'ichread: LINE BREAK'); linebreak:=true; readln(inst); writeln(listing); pageline:=succ(pageline); if pass=2 then begin writelineinformation; end; writeerrors; if (pageline >= pagelength) and not eof(inst) then begin page(listing); pageheader end; if not eof(inst) then startlistingline; lineposition:=0; until eof(inst) or (not eoln(inst)); if eof(inst) then begin if quote then error(12); if not eoinst then error(13); eoinst:=true; (* there are no more instructions *) if pass=2 then begin writelineinformation; end; writeerrors; end else begin read(inst,chr); lineposition:=succ(lineposition); (*if debugging then writeln(debug,'ichread',lineposition:3,chr); *) write(listing,chr) end; (* check that the character is a normal ASCII and not a control character. This will block tab characters. *) (* | 0 NUL| 1 SOH| 2 STX| 3 ETX| 4 EOT| 5 ENQ| 6 ACK| 7 BEL | 8 BS | 9 HT | 10 NL | 11 VT | 12 NP | 13 CR | 14 SO | 15 SI | 16 DLE| 17 DC1| 18 DC2| 19 DC3| 20 DC4| 21 NAK| 22 SYN| 23 ETB | 24 CAN| 25 EM | 26 SUB| 27 ESC| 28 FS | 29 GS | 30 RS | 31 US | 32 SP | 33 ! | 34 " | 35 # | 36 $ | 37 % | 38 & | 39 ' | 40 ( | 41 ) | 42 * | 43 + | 44 , | 45 - | 46 . | 47 / | 48 0 | 49 1 | 50 2 | 51 3 | 52 4 | 53 5 | 54 6 | 55 7 | 56 8 | 57 9 | 58 : | 59 ; | 60 < | 61 = | 62 > | 63 ? | 64 @ | 65 A | 66 B | 67 C | 68 D | 69 E | 70 F | 71 G | 72 H | 73 I | 74 J | 75 K | 76 L | 77 M | 78 N | 79 O | 80 P | 81 Q | 82 R | 83 S | 84 T | 85 U | 86 V | 87 W | 88 X | 89 Y | 90 Z | 91 [ | 92 \ | 93 ] | 94 ^ | 95 _ | 96 ` | 97 a | 98 b | 99 c |100 d |101 e |102 f |103 g |104 h |105 i |106 j |107 k |108 l |109 m |110 n |111 o |112 p |113 q |114 r |115 s |116 t |117 u |118 v |119 w |120 x |121 y |122 z |123 { |124 | |125 } |126 ~ |127 DEL *) if (ord(chr) < 32) or (ord(chr) > 126) then begin writeln(output,'bad character found in inst file'); error(33); end; end; procedure findword; (* find the next significant item, nonblank, noncomment *) procedure findnonblank; (* find the next non blank character, skipping up to one comment *) procedure skipblanks; (* skip blank characters *) begin while (chr=' ') and not eof(inst) do ichread (* now at a non-blank character or end of the file *) end; begin (* findnonblank *) skipblanks; if not eof(inst) then begin if chr='(' then begin parentheses:=succ(parentheses); ichread; if chr='*' then begin (* a comment *) commentline := linecount; comment:=true; while (chr<>')') and not eof(inst) do begin while ((chr<>'*') and (not eof(inst))) do ichread; if eof(inst) then error(10) else ichread; (* get past the } *) end; comment:=false; if not eof(inst) then ichread; (* get past the ) *) parentheses:=pred(parentheses) end end else if chr='{' (* new comment type *) then begin commentline := linecount; while (chr<>'}') and not eof(inst) do ichread; if eof(inst) then error(10) else ichread; (* get past the } *) end else if chr=')' then begin parentheses:=pred(parentheses); ichread; (* get a move on.. *) skipblanks end else if chr = ';' then begin eoinst:=true; ichread (* get past the ; *) end else if chr<>' ' then significant:=true; end end; begin (* findword *) significant:=false; (* keep trying to find a significant non blank *) while (not significant) and (not eof(inst)) do findnonblank end; procedure showparsed(var debug: text; parsedword: instword; (* the latest word parsed from the instructions *) parsedlength: integer; (* the length of parsedword *) parsedlocation: integer (* current location in parsedword *) ); (* show the parsed word *) var index: integer; (* index to parsed word *) begin write(debug,'parsedword: "'); for index := 1 to parsedlength do write(debug,parsedword[index]); writeln(debug,'" <- parsedlength = ',parsedlength:1); write(debug,' '); for index := 1 to parsedlocation-1 do write(debug,' '); writeln(debug,'^ parsedlocation = ',parsedlocation:1); end; procedure parse(anumber : boolean); (* parse out a new word, the result is in parsedword the parser assumes that the first character of the word is in chr. If anumber is true then stop parsing when the next character is not a number *) var j: integer; (* index to clear the remainder of parsedword *) stopcharacter: boolean; (* characters to stop the parse *) stopcondition: boolean; (* condition to stop the parse *) begin for parsedlength:=1 to widinst do parsedword[parsedlength]:=' '; parsedlength:=1; parsedword[parsedlength]:=chr; {zzzfff} { if eoln(inst) then writeln(output,'parse: eoln'); writeln(output,'in parse: parsedword[',parsedlength:1,']=', parsedword[parsedlength]); } if not eoln(inst) then begin repeat (* record up to space, ';', end of line or end of file *) correct:=true; if parsedlength >= widinst then error(9); if correct then begin ichread; { 2000 Oct 17: the linebreak prevented correct parsing if the ";" was on the line following the with. That is, for: get from 4021131 -50 to same +30 direction + with a4021131c ; the 'c' was lost. if not eof(inst) and not linebreak then begin } if not eof(inst) then begin parsedlength:=succ(parsedlength); parsedword[parsedlength]:=chr { ; writeln(output,'parse chr = ',chr); } end end; stopcharacter := (chr=';') or (chr in [' ','(',')']) or ( (not (chr in numbers)) and (not (chr in ['.','+','-']) and anumber)); stopcondition := linebreak or eof(inst) or (not correct); until stopcharacter or stopcondition; if stopcharacter then parsedlength:=pred(parsedlength); (* bump parentheses *) if chr in ['(',')'] then begin case chr of '(': parentheses:=succ(parentheses); ')': parentheses:=pred(parentheses) end; ichread (* get past the parentheses *) end; if parsedlengthnil then begin (* find last line *) while workline^.next<>nil do begin if workline = workline^.next then begin writeln(output, 'PROGRAM ERROR:', ' An infinite end loop was found in irquote'); workline^.next := nil; end else workline:=workline^.next; end; getline(nextline); (* grab new line *) workline^.next:=nextline; (* string it in *) workline:=nextline; (* move pointer *) end else begin getline(workline); (* make line *) aline:=workline (* make aline know about this *) end; aquote:=chr; if chr in ['''','"'] then begin quote:=true; ichread; while (chr<>aquote) and correct do begin if (workline^.length=linelength) or linebreak then begin getline(nextline); workline^.next:=nextline; workline:=nextline end; workline^.length:=succ(workline^.length); workline^.letters[workline^.length]:=chr; ichread end; quote:=false; if correct then ichread (* get past the quote *) end else begin error(16); geteoinst end end end; procedure irtitle; (* get the book title from the instructions *) var extra: lineptr; begin (*if debugging then writeln(debug,'irtitle');*) irquote(title); if correct then begin (* now remove extra title: *) extra:=title^.next; while extra<>nil do clearline(extra); title^.next:=nil end end; procedure irlongname; (* get the book longname from the instructions *) var extra: lineptr; begin (*if debugging then writeln(debug,'irlongname');*) (* clear away any previous junk *) clearline(longname); irquote(longname); if correct then longnameexists := true else longnameexists := false; (* clear any extra lines *) extra := longname^.next; if extra <> nil then begin writeln(output,'WARNING: the extra letters of this long name', ' are being ignored:'); write (output,'"'); writeline(output, longname, false); writeln(output,'"'); writeln(output,'The extra letters are: '); write (output,'"'); writeline(output, extra, false); writeln(output,'"'); emptyline(extra); longname^.next := nil; end; end; procedure irword; (* this routine reads and translates an instruction word into the delila word called 'word'. it sets correct to true if the read was successful, otherwise sets it to false *) var matching: boolean; (* have we identified the parsed word? *) i: integer; (* index for parsedword and wordlist *) begin (*if debugging then writeln(debug,'irword');*) findword; if not eoinst then begin parse(false); (* result is in parsedword *) if correct then begin matching:=false; word:=alldelila; (* start at first delilaword *) (* now identify the parsed word in wordlist: *) while ((not matching) and (ord(word) < ord(eodelila))) (* scan for the worddelila *) do begin (*if debugging then writeln(debug,'ird1',ord(word):5); if debugging then writeln(debug,wordlist[word]); *) matching:=true; i:=0; while ((parsedword[i+1] <> ' ') and matching and (i < widinst-1)) (* scan this worddelila *) do begin i:=succ(i); (*if debugging then writeln(debug,'ird4',i:5,parsedword); *) matching:=(parsedword[i]=wordlist[word][i]) end; if matching then if i <= minword then if wordlist[word][i+1] <> ' ' then matching:=false; if not matching then word:=succ(word) end; correct:=matching end; if not correct then begin (* the end of the word is (probably) before this point *) lineposition := pred(lineposition); error(2); lineposition := succ(lineposition); geteoinst end end else error(15) end; procedure irkeyname; (* read a key name into the variable keyname. left justified, blanks trailing. *) var i: integer; begin if not eoinst then begin findword; parse(false); if correct then begin with keyname do begin length:=0; (* stuff the parsed word into the key name *) for i:=1 to namelength do begin letters[i]:=parsedword[i]; if parsedword[i] <> ' ' then length:=succ(length) end; if length=namelength *) if not correct then begin error(14); geteoinst end (*;if debugging then writeln(debug,'keyname number',letters,'number');*) end end end else error(15) end; procedure irnumber; (* read an integer and put it in the global inumber. If the number has a decimal place, the truncated version is returned in inumber and the complete number is returned in rnumber. check that the number has proper format. *) var i, (* index to parsedword *) sign, (* of the number *) start, (* reading the number from place in parsedword *) stop, (* to *) power, (* powers of 10 for this number *) increment: (* to inumber *) integer; areal: boolean; (* a real number was encountered *) procedure digitize; (* convert a characte to a digit *) begin case parsedword[i] of '0': increment:=0; '1': increment:=1; '2': increment:=2; '3': increment:=3; '4': increment:=4; '5': increment:=5; '6': increment:=6; '7': increment:=7; '8': increment:=8; '9': increment:=9 end; end; procedure signblank; (* allow blanks after the sign *) begin if parsedword[i+1]=' ' then begin findword; parse(true); start:=1 end end; begin inumber := -maxint; (* flag for reading error *) areal := false; if not eoinst then begin findword; parse(true); if correct then begin i:=1; if parsedword[i] = '+' then begin sign:=+1; start:=2; signblank end else if parsedword[i] = '-' then begin sign:=-1; start:=2; signblank end else begin sign:=+1; start:=1 end; i:=start; correct:=true; while (i < widinst) and (parsedword[i] <> ' ') and correct and (not areal) do begin if parsedword[i] = '.' then areal := true; correct:=(parsedword[i] in numbers) or areal; if not areal then i:=succ(i) end; if correct then begin (* build that number.. *) stop:=i-1; (* parsedword[i] is blank *) power:=1; (* start at one"s place *) inumber:=0; (* start sum at zero *) for i:=stop downto start do begin digitize; inumber:=inumber+power*increment; power:=power*10 end; inumber:=sign*inumber; rnumber := inumber; if areal then begin start := stop+2; power:=10; (* start at tens place *) i := start; while parsedword[i] <> ' ' do begin digitize; { writeln(output,'parsedword[',i,']=', parsedword[i]); writeln(output,'increment=', increment:1); writeln(output,'rnumber=', rnumber:10:5); } rnumber:=rnumber+increment/power; power:=power*10; i := succ(i); end; end; if areal then rvalueposition(rnumber) else ivalueposition(inumber); end end; if not correct then begin error(4); geteoinst end end else error(15) { ;writeln(output,'inumber=', inumber:1); } end; (* ************ MUTATION PROCEDURES ***************************************) (* To allow these routines to skip to the end of the instruction, they must be below geteoinst *) (* begin module delila.readchangeset *) (* routines for handling mutations: changeset *) procedure mutationparseerror(e: integer); (* show the error at the correct location on the previously parsed word *) var shift: integer; (* how much to shift the recording line position to account for the position of the error in the parsedword *) begin {EEE} shift := (parsedlength - parsedlocation) + 1; { showparsed; writeln(output,'mutationparseerror: shift = ',shift:1); } lineposition := lineposition - shift; error(e); lineposition := lineposition + shift; end; procedure grabcharacter(var c: char; var done: boolean); (* Get the next character in parsedword *) begin { writeln(output,'grabcharacter ---------'); } if parsedlocation >= parsedlength then begin done := true; c := ' '; end else begin parsedlocation := succ(parsedlocation); c := parsedword[parsedlocation]; done := false; end { else begin writeln(output,'refused grabcharacter, correct = ',correct); end; } { ;showparsed; writeln(output,'grabcharacter: parsedlocation = ',parsedlocation:1); writeln(output,'grabcharacter: ',c:1); writeln(output,'grabcharacter: correct = ',correct); } end; procedure grabnumber(var number: integer; var donereading, correct: boolean); (* get the next number in parsedword. This can't be done by the regular Delila routines because it is imbedded in the middle of the parsed word. If we are at the end of the parse string, done is true. If the number is correct, correct is true. *) var digit: integer; (* a digit in the number *) digits: integer; (* number of digits obtained *) donenumber: boolean; (* done reading the number. This is different from being done reading the parsed string. *) sign: integer; (* +1 is positive, -1 is negative *) signnumber: integer; (* Number of signs for this number *) procedure signstuff; (* handle stuff for + for - signs *) begin if digits > 0 then begin mutationparseerror(31); writeln(output,'Numbers cannot have more than 1 sign (+ or -).'); correct := false; end; signnumber := succ(signnumber); end; begin { debugging := true; ;writeln(output,'grabnumber BEGIN DDD correct =',correct); } if correct then begin donenumber := false; number := 0; digits := 0; sign := +1; signnumber := 0; while not donenumber do begin grabcharacter(mutischar, donereading); if donereading then donenumber := true; {zzz deal with 1-54} if not donenumber then begin if mutischar in numbers then begin digit := ord(mutischar) - ord('0'); digits := succ(digits); number := 10*number + digit; { writeln(output,'digit = ',digit:1); writeln(output,'digits = ',digits:1); } end else if mutischar = '-' then begin sign := -1; signstuff end else if mutischar = '+' then begin sign := +1; signstuff end else begin (* the last location parsed is one base earlier *) parsedlocation := pred(parsedlocation); donenumber := true; {DDD} end end end; if digits = 0 then begin mutationparseerror(26); correct := false; end; if signnumber > 1 then begin {EEE} mutationparseerror(31); correct := false; end; if sign < 0 then number := -number end { ;showparsed; ; writeln(output,'grabnumber: number = ',round(number):1); ;writeln(output,'grabnumber END DDD correct =',correct); ;debugging := false; } end; procedure readchanges(var c: changeset; var done: boolean); (* Read the base changes from f in the form 'g2343c'. If there is an error in reading, set correct to false. *) begin { writeln(output,'BEGIN readchanges ==========================='); } with c.data[c.number] do begin changetype := 'c'; { nextnonblank(f); read(f, baseold); nextnonblank(f); read(f, basecoo1); } grabcharacter(baseold, done); if not (baseold in ['a','c','g','t']) then begin writeln(output,'ERROR: old base usually should be a, c, g, t'); mutationparseerror(20); end; grabnumber(basecoo1, done, correct); { writeln(output,'basecoo1 = ',basecoo1:1); } if correct then begin basecoo2 := basecoo1; grabcharacter(basenew, done); {zzzfff} { if done then begin writeln(output,'zzzzfff'); end; } if not (basenew in ['a','c','g','t']) then begin writeln(output,'ERROR: new base should be a, c, g, or t'); mutischar := basenew; mutationparseerror(21); end; inserts := 0; end; end; { writeln(output,'baseold = ',baseold); writeln(output,'basenew = ',basenew); writeln(output,'basecoo1 = ',basecoo1:1); writeln(output,'correct: ',correct); writeln(output,'END readchanges ==========================='); } end; procedure checknumberorder(basecoo1, basecoo2: real); (* Make sure that the value of basecoo1 is lessthan or equal to basecoo2, following the Libdef *) begin if basecoo1 > basecoo2 then begin error(27); {EEE} correct := false end end; procedure readinsertion(var c: changeset; var done: boolean); (* Read the insertion from f in the form 'i449,450tt'. *) var doneinsert: boolean; (* true if we are done with the insert *) begin { writeln(output,'readinsertion =================================='); } with c.data[c.number] do begin { read(f, changetype); read(f, basecoo1); skipblanks(f); read(f, comma); read(f, basecoo2); } grabcharacter(changetype, done); if done then correct := false; if correct then begin grabnumber(basecoo1, done, correct); if done then correct := false; if correct then begin grabcharacter(mutischar, done); if done then correct := false; if mutischar <> ',' then begin writeln(output,' comma expected between', ' coordinates for insertion'); mutationparseerror(25); end; if correct then begin grabnumber(basecoo2,done, correct); if correct then checknumberorder(basecoo1, basecoo2); if correct then begin inserts := 0; doneinsert := false; insertcomplement := false; repeat grabcharacter(mutischar, done); { showparsed; } if not done then begin { writeln(output,'readinsertion 1'); writeln(output,'inserts =',inserts:1,' mutischar = "',mutischar,'"'); } if (mutischar = '.') then doneinsert := true else if mutischar = '~' then begin if inserts <> 0 then begin doneinsert := true; correct := false; mutationparseerror(32); {zzz111} end else insertcomplement := true; end else if mutischar in ['a','c','g','t'] then begin inserts := inserts + 1; { writeln(output,'readinsertion 2'); writeln(output,'inserts =',inserts:1,' mutischar = "',mutischar,'"'); showparsed; } insert[inserts] := mutischar; if inserts > insertmax then begin writeln(output,' no more than ',insertmax:1, ' insertion bases allowed,', ' increase constant insertmax'); mutationparseerror(24); doneinsert := true; correct := false end; end else begin (* bad character detected *) doneinsert := true; correct := false; mutationparseerror(22); {zzz111} end; if (parsedlocation >= parsedlength) then doneinsert := true end else begin doneinsert := true; end until doneinsert end end end end end; end; procedure readdeletion(var c: changeset; var done: boolean); (* Read the deletion from parsedword in the form 'M55114 d449,450'. *) begin { writeln(output,'readdeletion'); } with c.data[c.number] do begin { read(f, changetype); read(f, basecoo1); skipblanks(f); read(f, comma); read(f, basecoo2); ;writeln(output,'readdeletion BEGIN DDD correct =',correct); } inserts := 0; grabcharacter(changetype, done); if correct then begin grabnumber(basecoo1,done, correct); if correct then begin grabcharacter(mutischar, done); if correct then begin if mutischar <> ',' then begin writeln(output,'comma expected between coordinates for deletion'); mutationparseerror(23); done := true; correct := false end; { ;writeln(output,'readdeletion MID 1 DDD correct =',correct); } if correct then grabnumber(basecoo2, done, correct); { ;writeln(output,'readdeletion MID 2 DDD correct =',correct); } if correct then checknumberorder(basecoo1, basecoo2); end end end end { ;writeln(output,'readdeletion END DDD correct =',correct); } end; procedure readchangeset(var c: changeset); (* read in the changes for a sequence. *) var done: boolean; (* we found the end of the statement so we are done *) dotfound: boolean; (* we have found a separation dot *) lastcharacter: char; (* have we found a dot at the end? *) begin { writeln(output,'readchangeset'); } c.number := 0; done := false; dotfound := false; findword; parse(false); if correct then begin (* make sure there is no initial dot *) parsedlocation := 0; (* about to read *) grabcharacter(mutischar, done); if mutischar = '.' then begin writeln(output,'Extra dots not allowed'); mutationparseerror(29); correct := false end; (* start main reading loop *) parsedlocation := 0; (* about to read *) while not done do if correct then with c do begin { writeln(output,'LOOP in READCHANGESET --------------------------------------------'); } {EEE} grabcharacter(mutischar, done); if not done then begin if mutischar = '.' then begin if dotfound then begin writeln(output,'Extra dots not allowed'); mutationparseerror(29); done := true end else dotfound := true; end else if mutischar in ['a','c','g','t','d','i'] then begin number := succ(number); if number > changesetmax then begin write(output,'Too many changes requested,', ' increase constant changesetmax.'); mutationparseerror(28); correct := false; done := true; end; if correct then begin dotfound := false; parsedlocation := pred(parsedlocation); { writeln(marksdelila,'* LOOP in READCHANGESET ----- number = ', number:1,'-----------------'); } case mutischar of 'a','c','g','t': readchanges(c, done); 'd': readdeletion(c, done); 'i': readinsertion(c, done); end; end end else begin write (output,'Unidentified change command:'); writeln(output,' "',mutischar,'"'); mutationparseerror(30); done := true; correct := false end; lastcharacter := mutischar end; if not correct then done := true; end; if lastcharacter = '.' then begin (* no dots at the end either *) writeln(output,'Extra dots not allowed'); mutationparseerror(29); end; end; { writechangeset(output,c); ;writeln(output,'END OF READCHANGESET --------------------------------------------'); } end; (* end module delila.readchangeset *) (* last common version: dbmutate.readchangeset version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* ************************************************************************ *) (* ************************************************************************ *) (* ************************************************************************ *) (* routines for marking the lister map *) (* begin module sortchanges *) procedure sortchanges(unsorted: changeset; var sorted: changeset; phenotype: char); (* sort the changeset unsorted and put the result into sorted. Since marks need to be in increasing position order (as currently defined in lister) it is nice to sort the changes for each piece. * When phenotype = 'w', the sequence is considered wild-type. In this case insertions are a single position and are listed BEFORE deletions, since deletions will extend further to the 3'. * When phenotype = 'm', the sequence is considered mutant. In this case deletions are a single position and are listed BEFORE insertions, since insertions will extend further to the 3'. *) type position = integer; (* make quicksort be happy but still standard *) function lessthan(a, b: integer): boolean; (* see quicksort *) var aposition, bposition: integer; (* adjusted locations of a and b *) atype, btype: char; (* the types of a and b *) function adjust(x: integer): integer; (* Sorting does NOT use the ends of an insert, but rather the middle bases. This allows g1t.i1,4cc to be properly sorted with the g1t in front. Sorting is sufficiently rare that we'll just compute this here. *) var xposition: integer; (* adjusted locations of x *) begin xposition := sorted.data[x].internal1; if sorted.data[x].changetype = 'i' { then xposition := xposition + 1; } then with sorted.data[x] do begin {*U WILD external unsorted: i1,2tct.i5,6ggg.d6,8.d1,1 *U WILD external sorted: i1,2tct.d1,1.i5,6ggg.d6,8 should not have sorted, so if internal1 + 1 < internal2 } if internal1 + 1 <= internal2 then xposition := xposition + 1; end; adjust := xposition end; (* adjust *) begin {SSS} aposition := adjust(a); bposition := adjust(b); if (aposition = bposition) and (a <> b) then begin atype := sorted.data[a].changetype; btype := sorted.data[b].changetype; { Interesting chunk of code, but base changes can be next to insertions and deletions too. (* this can only occur when insertions are right next to deletions *) if not( ((atype='i') and (btype='d')) or ((atype='d') and (btype='i')) ) then begin writeln(output,'sortchanges: DELILA SORTING ERROR'); writeln(output,'Please contact toms@ncifcrf.gov!'); write(output,'while sorting: '); wchangedata(output,sorted.data[a]); write(output,' and '); wchangedata(output,sorted.data[b]); writeln(output); write(output,'in: '); writechangeset(output,unsorted); writeln(output); write(output,'On the '); case phenotype of 'w': write(output,'wild-type'); 'm': write(output,'mutant'); end; writeln(output,'sequence.'); writeln(output,'At position ',aposition:1,'=',bposition:1); crash end; } case phenotype of (* the doubleY wins always *) 'w': begin (* wild-type situation *) if atype = 'i' then lessthan := true else lessthan := false end; 'm': begin (* mutant situation *) if atype = 'd' then lessthan := true else lessthan := false end; end end else lessthan:=(aposition < bposition); { old: lessthan:=(sorted.data[a].internal1 < sorted.data[b].internal1) } end; (* lessthan *) procedure swap(a,b: integer); (* see quicksort *) var hold: changedata; begin hold:=sorted.data[a]; sorted.data[a]:=sorted.data[b]; sorted.data[b]:=hold (*;write(output,'a=',a:1,', b=',b:1); print (@ for testing *) end; (* swap *) (* begin module quicksort *) procedure quicksort(left, right: position); (* quick sort a list between positions left and right, into ascending order. a position is simply a scalar of the form 0..max. the array to be sorted is dimensioned 1..max. (the difference in the ranges is important to the correct operation of the sort...) two external routines are used: function lessthan(a, b: position): boolean is a generalized test for value-at-a < value-at-b. procedure swap(a, b: position) switches the items at positions a and b. since these routines are external, the procedure is general. this procedure taken from the book 'algorithms + data structures = programs' by niklaus wirth, prentice-hall, inc., englewood cliffs, n.j.(1976), pp. 76-82 *) var lower, upper: position; (* the positions looked at currently *) center: position; (* the rough center of the region being sorted *) begin lower := left; center := (left + right) div 2; upper := right; repeat while lessthan(lower, center) do lower := succ(lower); while lessthan(center, upper) do upper := pred(upper); if lower <= upper then begin (* keep track of the center through the map: *) if lower = center then center:=upper else if upper = center then center:=lower; swap(lower, upper); lower := succ(lower); upper := pred(upper) end until lower > upper; if left < upper then quicksort(left, upper); if lower < right then quicksort(lower, right) end; (* end module quicksort version = 4.21; (@ of prgmod.p 1997 October 22 *) begin sorted := unsorted; quicksort(1,sorted.number) end; (* end module sortchanges version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* begin module propagateonechange *) procedure propagateonechange(var changes: changeset; m: integer; pie: pieceptr); (* Propagate the one change at m through the rest of the changeset. *) var location: integer; (* location of a change *) shift: integer; (* amount to shift a change *) b1, b2: integer; (* internal coordinate values for the basecoo1 and basecoo2 of the changeset *) deletes: integer; (* amount to delete the sequence *) pielength: integer; (* current length of the piece *) procedure loop; (* propagate the change to all other all changes. The second part of the boolean test, for changes at internal+1 with inserts is required for the case: get from 1 to 6 with i1,2tct.i1,3; Without it, the i1,3 doesn't get propagated past the i1,2tct. *) var n: integer; (* counter for the changes, later places *) begin {writeln(output,'loop : location = ',location:1);} for n := 1 to changes.number do begin {write(output,'propagateonechange: TRY ');} {wchangedata(output,changes.data[n]);} {writeln(output);} if n <> m then if (location <= changes.data[n].internal1) or ((location <= changes.data[n].internal1+1) and (changes.data[n].changetype = 'i')) then with changes.data[n] do begin {write(output,'propagateonechange: shifting ');} {wchangedata(output,changes.data[n]);} {writeln(output);} {write (output,'loop PRE: internal1 = ',internal1:1);} {writeln(output, ' internal2 = ',internal2:1);} internal1 := internal1 + shift; internal2 := internal2 + shift; {write (output,'loop POS: internal1 = ',internal1:1);} {writeln(output, ' internal2 = ',internal2:1);} end end; end; begin (* propagateonechange *) { writeln(output,'propagateonechange: m = ',m:1,' --------------------'); write(output,'propagateonechange: '); wchangedata(output,changes.data[m]); writeln(output); } with changes do begin location := data[m].internal1; (* make location point to actual insertion start *) {write (output,'changetype = data[',m:1,'] = ',data[m].changetype);} {writeln(output,' location = ',location:1);} pielength := piecelength(pie); with data[m] do begin b1 := internal1; b2 := internal2; case changetype of 'c': ; (* base changes cause no downstream changes *) 'i': begin location := location +1; (* this code is from changesequence *) if b1 > pielength then b1 := pielength; if b2 > pielength then b2 := pielength + 1; if b1 < 0 then b1 := 0; if b2 < 1 then b2 := 1; (* original *) deletes := b2 - b1 - 1; shift := inserts - deletes; loop end; (* NOTE: i and d code are NOT identical and cannot be condensed! *) 'd': begin (* this code is from changesequence *) if b1 < 1 then b1 := 1; if b2 < 1 then b2 := 1; if b1 > pielength then b1 := pielength; if b2 > pielength then b2 := pielength; shift := b2 - b1 + 1; shift := -shift; loop end; end; end; end; end; (* end module propagateonechange *) (* begin module delila.writemarks *) procedure doubleYmark(var markspots: text; internal1, internal2: integer; pie: pieceptr; insertlowerbits,insertupperbits: real); (* Put a doubleY mark between internal1 and internal 2 in file markspots for piece pie. The difference between these two coordinates must be 1, and internal1 < internal2. *) var b1: integer; (* adjusted coordinate for internal1 *) b2: integer; (* adjusted coordinate for internal2 *) shiftY: integer; (* amount in bases to shift the doubleY *) begin writeln(markspots,'* doubleY: internal1 = ',internal1:1); writeln(markspots,'* doubleY: internal2 = ',internal2:1); if internal2-internal1 <> 1 then begin writeln(output,'doubleY: program error: ', 'internal2-internal1 <> 1'); crash end; if withininternal(pie, internal2) then begin (* If internal2 is inside the piece, then put the doubleY to the left of it. *) b1 := internal2; b2 := internal2; shiftY := 0; writeln(markspots,'* doubleY: case 1'); end else begin if withininternal(pie, internal1) then begin (* Ok, so internal 1 IS inside, so put the doubleY to the right of it by shifting 1 *) b1 := internal1; b2 := b1; shiftY := 1; {zzzYYY} writeln(markspots,'* doubleY: case 2a'); end else begin (* Neither are inside; this occurs when one deletes on the left edge. So reduce to the edge and put the mark to the left of that base (no shift). *) b1 := internal2; reduceposition(pie, b1); b2 := b1; shiftY := 0; writeln(markspots,'* doubleY: case 2b'); end; end; writeln(markspots,'U', ' ',inttopie(b1,pie):widbase, ' ',insertlowerbits:widbits:decbits, (* lowerbits *) ' ',inttopie(b2,pie):widbase, ' ',insertupperbits:widbits:decbits, (* upperbits *) ' ',shiftY:1, ' doubleY') end; procedure createmark(var markspots: text; b1, b2: integer; (* internal coordinates *) boxposition: char; changetype: char; pie: pieceptr; insertupperbits: real; (* upperbits for insertion symbol *) insertlowerbits: real; (* upperbits for insertion symbol *) deleteupperbits: real; (* upperbits for deletion symbol *) deletelowerbits: real (* upperbits for deletion symbol *) ); (* Create box marks for an insertion or deletion. Boxes are defined to be AROUND the given position. Internal coordinates are used to allow the routine to be called in a loop. The routine writes external (piece) coordinates out. No marks are made if either end of the object is off the piece. *) var lowerbits, upperbits: real; (* the vertical range of the mark *) piecebase1, piecebase2: integer; (* piece coordinates corresponding to the internal coordinates b1 and b2 *) procedure kind(c: char); (* show the kind of thing we are marking *) begin case c of 'd': write(markspots,'deletion'); 'i': write(markspots,'insertion'); 'f': write(markspots,'full'); 'l': write(markspots,'left'); 'm': write(markspots,'mid'); 'r': write(markspots,'right'); end; end; begin { writeln(markspots,'* createmark, changetype = ',changetype); } if withininternal(pie,b1) and withininternal(pie,b2) then begin case changetype of 'd': begin lowerbits := deletelowerbits; upperbits := deleteupperbits end; 'i': begin lowerbits := insertlowerbits; upperbits := insertupperbits end; end; piecebase1 := inttopie(b1,pie); piecebase2 := inttopie(b2,pie); write(markspots,'U', ' ',piecebase1:widbase, ' ',lowerbits:widbits:decbits, ' ',piecebase2:widbase, ' ',upperbits:widbits:decbits, ' ',0:widbits); (* no displacement *) write(markspots,' ('); kind(boxposition); kind(changetype); write(markspots,') '); kind(boxposition); kind(changetype); writeln(markspots); end end; procedure multimarks(var markspots: text; internal1, internal2: integer; changetype: char; pie: pieceptr; insertupperbits: real; (* upperbits for insertion symbol *) insertlowerbits: real; (* upperbits for insertion symbol *) deleteupperbits: real; (* upperbits for deletion symbol *) deletelowerbits: real (* upperbits for deletion symbol *) ); (* create multiple marks to file markspots to represent a deletion (changetype = 'd') or insertion (changetype = 'i'), for the range internal1 to internal2. The routine only uses internal coordinates. *) var partindex: integer; (* an index for the marker parts of insertion or deletion *) begin writeln(markspots,'* Multimarks, changetype = ',changetype); writeln(markspots,'* internal1 = ',internal1:1); writeln(markspots,'* internal2 = ',internal2:1); if internal1 = internal2 then createmark(markspots,internal1,internal2,'f',changetype,pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits) else begin createmark(markspots,internal1,internal1,'l',changetype,pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits); partindex := internal1+1; while partindex < internal2 do begin createmark(markspots,partindex,partindex,'m',changetype,pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits); { writeln(markspots,'* multimarks: partindex = ',partindex:1); } partindex := partindex + 1; end; createmark(markspots,internal2,internal2,'r',changetype,pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits); end end; (* multimarks *) procedure setinternal(var changes: changeset; pie: pieceptr); (* piece for these changes *) (* Determine the internal coordinates for each change in the changeset *) var hold: integer; (* holding variable to switch internal1 and internal2 *) n: integer; (* counter for the changes *) begin { writeln(marksdelila,'* setinternal=========================='); } (* Obtain internal coordinates for each changedata *) for n :=1 to changes.number do with changes.data[n] do begin (* Convert to internal coordinates, without wrapping *) internal1 := nwpietoint(basecoo1, pie); internal2 := nwpietoint(basecoo2, pie); {writeln(marksdelila,'* setinternal: internal1 = ',internal1:1);} {writeln(marksdelila,'* setinternal: internal2 = ',internal2:1);} (* handle reversed coordinate system *) if internal1 > internal2 then begin hold := internal2; internal2 := internal1; internal1 := hold; { writeln(marksdelila,'* setinternal: FLIP '); } end; { end; } { else writeln(marksdelila,'* setinternal: NO FLIP '); } (* move insertions to the end of their regions - simplifies all later code! *) if changetype = 'i' then begin (* handle cases off ends. Note that nwpietoint sets out of bounds values to the ends *) if (internal1 = 0) and (internal2 = 0) then begin internal1 := -1; end; if (internal1 > piecelength(pie)) and (internal2 > piecelength(pie)) then begin internal1 := internal1 - 1; (* ie, piecelength *) internal2 := internal1 + 1; end; (* If the insert locations are far apart, then mark the part deleted. Otherwise leave it alone. *) if abs(internal1-internal2) = 1 then begin insertasdeletion := false end else begin insertasdeletion := true end; {PPP} end; { writeln(marksdelila,'* setinternal: n = "',n:1,'"'); writeln(marksdelila,'* setinternal: changetype = "',changetype:1,'"'); writeln(marksdelila,'* setinternal: basecoo1 = ',basecoo1:1); writeln(marksdelila,'* setinternal: basecoo2 = ',basecoo2:1); writeln(marksdelila,'* setinternal: internal1 = ',internal1:1); writeln(marksdelila,'* setinternal: internal2 = ',internal2:1); writeln(marksdelila,'* adjusted external = ',inttopie(internal1,pie):1); writeln(marksdelila,'* adjusted external = ',inttopie(internal2,pie):1); writeln(marksdelila,'* insertasdeletion: ',insertasdeletion); } { bwpie(output,pie); } end; end; procedure skippiece(var marksdelila: text); (* skip to the next piece *) begin writeln(marksdelila); writeln(marksdelila,'piece #',def.num.item:1,' - skip ahead to next piece'); end; procedure writewildtypemarks(var markspots: text; changes: changeset; insertupperbits: real; (* upperbits for insertion symbol *) insertlowerbits: real; (* upperbits for insertion symbol *) deleteupperbits: real; (* upperbits for deletion symbol *) deletelowerbits: real; (* upperbits for deletion symbol *) changeupperbits: real; (* upperbits for change symbol *) changelowerbits: real; (* upperbits for change symbol *) pie: pieceptr; (* piece for these changes *) thenumber: integer); (* the piece number *) (* Write the marks to file markspots, at the locations defined. Writemarks works with two sequences, first is the wild-type sequence and second is the mutant sequence. Definition: This routine assumes that all previous actions have placed us onto the wild-type sequence. NOTE: the internal values changes get modified locally here but are not altered in the original copy. *) const shiftdown = 1; (* shift down the change mark arrow on wt sequence. {zzz} This should be removed by changing the definition. *) var chorient: integer; (* orientation of the change on a piece. This is needed to switch the order of deletion marks. *) n: integer; (* counter for the changes *) sorted: changeset; (* the changes sorted by the exact positions of the marks, must be done for both wild and mutant sequences *) thelength: integer; (* the length of pie *) unsorted: changeset; (* the exact positions of the marks, unsorted *) begin writeln(markspots); { writeln(markspots,'* writeWILDTYPEmarks'); } write(markspots,'* piece #',thenumber:1,' '); writechangeset(markspots,changes); writeln(markspots); if pie^.key.piedir = pie^.key.coodir then chorient := +1 else chorient := -1; { not needed displacement := chorient; } thelength := piecelength(pie); { writeln(markspots,'* chorient = ',chorient:1); writeln(markspots,'* thelength = ',thelength:1); } (* first do the wild type sequence: *) {zzzppp} {Othis should probably be removed - no propagation for wild type because the change coordinates refer to wild type } { OLD CODE: propagatechanges(changes,unsorted,true); (* wild type *) } (* No propagation for wild type because the change coordinates refer to wild type *) unsorted := changes; (* adjust the locations of the marks *) {zzz probably eliminate entire section, it's not used now} with unsorted do for n := 1 to number do with data[n] do begin { writeln(markspots,'* basecoo1 = ',basecoo1:1); writeln(markspots,'* basecoo2 = ',basecoo2:1); writeln(markspots,'* internal1 = ',internal1:1); writeln(markspots,'* internal2 = ',internal2:1); } case changetype of 'c': ; (* wait until the next loop *) 'i': begin { shifting of the inserts should be done at output, so I am blocking this segment for now. it is not clear how badly not doing this will affect sorting if chorient = +1 then begin internal1 := shiftbase(internal1,+1); internal2 := shiftbase(internal2,-1); end else begin internal1 := shiftbase(internal1,-1); internal2 := shiftbase(internal2,+1); end; writeln(markspots,'* chorient internal1 = ',internal1:1); writeln(markspots,'* chorient internal2 = ',internal2:1); } end; 'd': begin (* mark deletion as a red bar on the wild type sequence *) end; end; end; sortchanges(unsorted, sorted, 'w'); { writeln(markspots,'** MARKSPOTS *********************************************'); write(markspots,'*U WILD external '); write(markspots,'unsorted: '); writechangeset(markspots,unsorted); writeln(markspots); write(markspots,'*U WILD external '); write(markspots,'sorted: '); writechangeset(markspots,sorted); writeln(markspots); write(markspots,'*U WILD internal '); write(markspots,'unsorted: '); checkchangeset(markspots,unsorted); writeln(markspots); write(markspots,'*U WILD internal '); write(markspots,'sorted: '); checkchangeset(markspots,sorted); writeln(markspots); writeln(markspots,'********************************************************'); } (* print the marks out for wildtype sequence *) with sorted do for n := 1 to number do with data[n] do begin case changetype of 'c': begin writeln(markspots,'U', ' ',inttopie(internal1,pie):widbits, ' ',changeupperbits-shiftdown:widbits:decbits, ' ',inttopie(internal2,pie):widbits, ' ',changelowerbits-shiftdown:widbits:decbits, ' ',0.0:widbits:decbits, ' (', baseold, '->', basenew, ') changeworra ') {TTT this is old!!!} { writeln(markspots,'U', ' ',inttopie(internal1-displacement,pie):widbits:decbits, ' ',changeupperbits-shiftdown:widbits:decbits, ' ',inttopie(internal2-displacement,pie):widbits:decbits, ' ',changelowerbits-shiftdown:widbits:decbits, ' ',+displacement:widbits:decbits, ' (', baseold, '->', basenew, ') changeworra ') } end; 'i': begin (* mark insertion location as a black bar on the wild type sequence *) { this is probably important to keep - done below - is it right? (* handle the special cases of insertion at the ends of the sequence *) if (nwpietoint(trunc(basecoo1),pie) = 0) then begin basecoo1 := basecoo1 + 1.0; basecoo2 := basecoo2 + 1.0; displacement := -1.8; (* I don't know why this needs more *) end; if (nwpietoint(trunc(basecoo1),pie) > piecelength(pie)) then begin basecoo1 := basecoo1 - 1.0; basecoo2 := basecoo2 - 1.0; displacement := 0.0; end; } (* handle the special cases of insertion at the ends of the sequence *) if internal1 = 0 then begin { internal1 := internal1 + 1; internal2 := internal2 + 1; messes up } end; {NNN} if internal1 > thelength then begin { internal1 := thelength; internal2 := thelength; handled in setinternal } {note: I'm guessing that this is the right thing to do rather than just subtracting 1! } end; writeln(markspots,'* alive'); writeln(markspots,'* basecoo1 = ',basecoo1:1); writeln(markspots,'* basecoo2 = ',basecoo2:1); writeln(markspots,'* internal1 = ',internal1:1); writeln(markspots,'* internal2 = ',internal2:1); {YY} { writeln(markspots,'* insertasdeletion: ',insertasdeletion); } { simplify: if (1 = abs(internal1 - internal2)) and (not insertasdeletion) } if insertasdeletion then multimarks(markspots,internal1+1, internal2-1,'d',pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits) else doubleYmark(markspots, internal1, internal2, pie, insertlowerbits,insertupperbits) end; 'd': begin (* mark deletion as a red bar on the wild type sequence *) multimarks(markspots,internal1,internal2,'d',pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits); end; end; end; end; {WWWW 3} procedure writemutantmarks(var markspots: text; changes: changeset; insertupperbits: real; (* upperbits for insertion symbol *) insertlowerbits: real; (* upperbits for insertion symbol *) deleteupperbits: real; (* upperbits for deletion symbol *) deletelowerbits: real; (* upperbits for deletion symbol *) changeupperbits: real; (* upperbits for change symbol *) changelowerbits: real; (* upperbits for change symbol *) pie: pieceptr; (* piece for these changes *) thenumber: integer); (* the piece number *) (* Write the marks to file markspots, at the locations defined. Writemarks works with two sequences, first is the wild-type sequence and second is the mutant sequence. Definition: This routine assumes that all previous actions have placed us onto the wild-type sequence. NOTE: the internal values changes get modified locally here but are not altered in the original copy. *) var displacement: integer; (* number of bases to displace the mark backwards *) chorient: integer; (* orientation of the change on a piece. This is needed to switch the order of deletion marks. *) n: integer; (* counter for the changes *) sorted: changeset; (* the changes sorted by the exact positions of the marks, must be done for both wild and mutant sequences *) thelength: integer; (* the length of pie *) unsorted: changeset; (* the exact positions of the marks, unsorted *) begin writeln(markspots); { writeln(markspots,'* writeMUTANTmarks'); } write(markspots,'* piece #',thenumber:1,' '); writechangeset(markspots,changes); writeln(markspots); if pie^.key.piedir = pie^.key.coodir then chorient := +1 else chorient := -1; displacement := chorient; thelength := piecelength(pie); { writeln(markspots,'* chorient = ',chorient:1); writeln(markspots,'* thelength = ',thelength:1); } {original: write(markspots,'p - skip ahead to mutated piece'); writeln(markspots,' #',(thenumber+1):1); (* anticipate its number *) } {zzzppp} { propagatechanges(changes,unsorted,false); (* mutant *) } unsorted := changes; with unsorted do for n := 1 to number do propagateonechange(unsorted,n,pie); (* create marks for mutant sequence *) (* better title? set coordinates of marks for mutant sequence *) with unsorted do for n := 1 to number do with data[n] do begin { writeln(markspots,'* writeMUTANTmarks basecoo1 = ',basecoo1:1); writeln(markspots,'* writeMUTANTmarks basecoo2 = ',basecoo2:1); writeln(markspots,'* writeMUTANTmarks inserts = ',inserts:6); writeln(markspots,'* writeMUTANTmarks internal1 = ',internal1:1); writeln(markspots,'* writeMUTANTmarks internal2 = ',internal2:1); } case changetype of 'c': begin { NO!!??? internal1 := internal1 - displacement; (* still valid? *) internal2 := internal1; } end; 'i': begin { writeln(markspots,'* VOILA ADJUST internal1 = ',internal1:1); writeln(markspots,'* VOILA ADJUST internal2 = ',internal2:1); } if inserts = 0 then begin internal2 := internal1 + 1; { writeln(markspots,'* ADJUST a: ON THE NOSE'); } end else begin { writeln(markspots,'* ADJUST b: ON THE NOSE'); } internal1 := internal1 + 1; internal2 := internal1 + inserts - 1; end end; 'd': begin end; end; end; sortchanges(unsorted, sorted,'m'); { writeln(markspots,'*2***********************************************************'); write(markspots,'* MUTANT external '); write(markspots,'unsorted: '); writechangeset(markspots,unsorted); writeln(markspots); write(markspots,'* MUTANT external '); write(markspots,'sorted: '); writechangeset(markspots,sorted); writeln(markspots); write(markspots,'* MUTANT internal '); write(markspots,'unsorted: '); checkchangeset(markspots,unsorted); writeln(markspots); write(markspots,'* MUTANT internal '); write(markspots,'sorted: '); checkchangeset(markspots,sorted); writeln(markspots); writeln(markspots,'*3***********************************************************'); } (* mutant *) with sorted do for n := 1 to number do with data[n] do begin case changetype of 'c': begin writeln(markspots,'U', ' ',inttopie(internal1,pie):widbits, ' ',changeupperbits:widbits:decbits, ' ',inttopie(internal2,pie):widbits, ' ',changelowerbits:widbits:decbits, ' ',0.0:widbits:decbits, { ' ',displacement:widbits:decbits, } ' (', baseold, '->', basenew, ') change '); end; 'i': begin (* handle the special cases of insertion at the ends of the sequence *) if internal1 = 0 then begin internal1 := internal1 + 1; internal2 := internal2 + 1; end; if internal1 > thelength then begin (* writeln(output,'nwpietoint(trunc(basecoo1),pie) = ',nwpietoint(trunc(basecoo1),pie):1); writeln(output,'nwpietoint(trunc(basecoo2),pie) = ',nwpietoint(trunc(basecoo2),pie):1); writeln(output,'piecelength(pie) = ',piecelength(pie):1); *) internal1 := internal1 - 1; (* assumes internal1 = thelength *) internal2 := internal2 - 1; end; writeln(markspots,'* NOWH internal1 = ',internal1:1); writeln(markspots,'* NOWH internal2 = ',internal2:1); if inserts = 0 then doubleYmark(markspots, internal1, internal2, pie, insertlowerbits,insertupperbits) else multimarks(markspots, internal1, internal2,'i',pie, insertupperbits, insertlowerbits, deleteupperbits, deletelowerbits); writeln(markspots,'* VOILA internal1 = ',internal1:1); writeln(markspots,'* VOILA internal2 = ',internal2:1); end; 'd': begin (* mark the point between bases of the deletion *) (* slide back to point just before deletion: *) internal1 := internal1 - 1; (* make second base match first base: *) internal2 := internal1 + 1; writeln(markspots,'* writeMUTANTmarks DELETION internal1 = ',internal1:1); writeln(markspots,'* writeMUTANTmarks DELETION internal2 = ',internal2:1); doubleYmark(markspots, internal1, internal2, pie, insertlowerbits,insertupperbits) end; end; end; { write(markspots,'*N sorted: '); checkchangeset(markspots,sorted); writeln(markspots); write(markspots,'*N unsorted: '); checkchangeset(markspots,unsorted); writeln(markspots); write(markspots,'*N changes: '); checkchangeset(markspots,changes); writeln(markspots); write(markspots,'*N sorted: '); writechangeset(markspots,sorted); writeln(markspots); write(markspots,'*N unsorted: '); writechangeset(markspots,unsorted); writeln(markspots); write(markspots,'*N changes: '); writechangeset(markspots,changes); writeln(markspots); } { writeln(markspots,'p - skip ahead to next piece'); } end; {WW} {WWW} {WWW} (* originally from version = 1.89; (@ of dbmutate.p 1999 April 26 *) (* end module delila.writemarks *) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* routines for extracting pieces *) procedure dnaconcat(Adna: dnaptr; Bdna: dnaptr; var Cdna: dnaptr); (* concatenate Adna to Bdna, put in Cdna *) var aDNAptr: dnaptr; (* pointer into Adna *) bDNAptr: dnaptr; (* pointer into Bdna *) cpart: dnaptr; (* a new part for Cdna *) pb: integer; (* position in b *) pc: integer; (* position in c *) begin emptydna(Cdna); getdna(Cdna); cpart:=Cdna; (* copy a into c *) aDNAptr:=Adna; while aDNAptr<>nil do begin cpart^.part:=aDNAptr^.part; cpart^.length:=aDNAptr^.length; aDNAptr:=aDNAptr^.next; if aDNAptr<>nil then begin getdna(cpart^.next); cpart:=cpart^.next end else cpart^.next:=nil end; (* now copy b into c and compact *) pc := cpart^.length; pb := 1; bDNAptr:=Bdna; repeat if cpart^.length = dnamax then begin getdna(cpart^.next); cpart:=cpart^.next; cpart^.next := nil; cpart^.length := 1 end else begin cpart^.length := succ(cpart^.length); end; cpart^.part[cpart^.length] := bDNAptr^.part[pb]; pb := succ(pb); if pb > bDNAptr^.length then begin bDNAptr := bDNAptr^.next; pb := 1; end; until bDNAptr = nil; end; procedure dnacomplement(Adna: dnaptr; var Bdna: dnaptr); (* make the complement of Adna in Bdna *) var a: dnaptr; (* pointer into Adna *) bpart: dnaptr; (* a new part for Bdna *) p: integer; (* position in a *) done: boolean; (* done complementing? *) begin { writeln(output,'dnacomplement'); } emptydna(Bdna); getdna(Bdna); bpart := Bdna; a := Adna; done := false; while not done do begin bpart^.length := a^.length; for p := 1 to a^.length do bpart^.part[a^.length - p + 1] := complement(a^.part[p]); if a^.next <> nil then begin new(bpart); { writeln(output,'new(bpart)'); } bpart^.next := Bdna; Bdna := bpart; a := a^.next; end else done := true; { writeln(output,'(a^.next = nil) is ', (a^.next = nil)); } end; { writeln(output,'about to crash in dnacomplement'); crash; } end; procedure getdnasegment( big: dnaptr; var little: dnaptr; s, e: integer); (* create the dna segment from s (start) to e (end) of a big dna. Both s and e are INTERNAL coordinates, 1 to the length of a dna. *) var b: integer; (* position in a bDNAptr segment *) bDNAptr: dnaptr; (* pointer to big dna *) l: integer; (* position in an ldna segment *) ldna: dnaptr; (* pointer to little dna *) p: integer; (* position from s to e *) begin (* writeln(output, 'about to getdnasegment !!!!!!!!!!!!!!!!!!!!'); {zzzddd} *) if s > e then begin writeln(output,'getdnasegment: s < e is required'); halt end; bDNAptr := big; { writeln(output,'s = ',s:1); writeln(output,'e = ',e:1); writeln(output,'bDNAptr^.length = ',bDNAptr^.length:1); } (* get p to s-1 in big *) p := bDNAptr^.length; if p < s - 1 then begin while p < s - 1 do begin if bDNAptr^.next = nil then begin writeln(output,'getdnasegment: request off end of piece'); halt; end; bDNAptr := bDNAptr^.next; p := p + bDNAptr^.length; { ;writeln(output,'STEPPPPPPPP p = ',p:1); } end; b := bDNAptr^.length - (p-s) - 1; (* the point of start in this segment is p-s, but back up one position so that the loop can advance us into place *) end else b := s - 1; (* back one position so that the loop can advance us into place *) emptydna(little); getdna(little); ldna := little; l := 0; for p := s to e do begin (* prepare the little to receive *) l := succ(l); if l > dnamax then begin (* shift to next little segment *) (* writeln(output,'getdnasegment: NEW SEGMENT'); {zzzddd} *) ldna^.length := dnamax; getdna(ldna^.next); ldna := ldna^.next; l := 1; (* writeln(output,'ldna^.length = ',ldna^.length:1);{zzzddd} *) end; b := succ(b); {zzzooo} if b > bDNAptr^.length then begin (* shift to next big segment *) (* writeln(output,'b = ',b:1); *) if bDNAptr^.next = nil then begin writeln(output,'getdnasegment: request e (' ,e:1,') is beyond piece end'); crash; halt; end; { writeln(output,'next segment'); } b := 1; bDNAptr := bDNAptr^.next; end; (* transfer the base *) ldna^.part[l] := bDNAptr^.part[b]; { write (output,'getdnasegment '); write (output,basetochar(ldna^.part[l])); write (output,' p=',p:1); write (output,' l=',l:1); writeln(output,' b=',b:1); if basetochar(ldna^.part[l]) = 'X' then crash; zzzddd } end; ldna^.length := l; (* writeln(output,'l = ',l:1); writeln(output,' ', '123456789 123456789 123456789 123456789 123456789 123456789 123456789 1234567890'); writeln(output,' ', ' 1 2 3 4 5 6 7 8'); writeln(output, 'getdnasegment -----BIG------------'); {zzzddd} LIKEbwdna(output,big); {zzzddd} writeln(output, 'getdnasegment -----LITTLE---------'); {zzzddd} LIKEbwdna(output,little); {zzzddd} writeln(output, 'END to getdnasegment !!!!!!!!!!!!!!!!!!!!'); {zzzddd} *) end; (* getdnasegment *) procedure circledna( big: pieceptr; var little: pieceptr; s,e: integer; inverting: boolean); (* Extract a segment of a big dna, treated as a circle. Both s (start) and e (end) are internal coordinates that must be inside the big piece. There are 4 possible cases: 1<-s e<-n inverting 1 s->e n not inverting 1 e<-s n inverting 1->e s->n not inverting 1 and n are the ends of big. -> means the segment that is given 5' to 3'. *) const showcase = false; (* show the cases *) {zzzddd} var first: dnaptr; (* first dna bit to grab *) n: integer; (* length of big *) second: dnaptr; (* second dna bit to grab *) procedure invert(var d: dnaptr); (* complement the dna in d *) var i: dnaptr; begin i := d; d := nil; dnacomplement(i, d); (* 1999 April 27 forgetting the following cleardna caused a memory leak! *) cleardna(i); end; begin n := piecelength(big); {zzzddd} { writeln(output, 'about to circledna --------------------'); writeln(output,'circledna, n = ',n:1); writeln(output,'s = ',s:1); writeln(output,'e = ',e:1); writeln(output,'inverting = ',inverting); } if (s < e) then begin if inverting then begin (* 1<-s e<-n inverting *) if showcase then writeln(output,' 1<-s e<-n inverting'); first := nil; second := nil; getdnasegment(big^.dna,first,1,s); invert(first); getdnasegment(big^.dna,second,e,n); invert(second); dnaconcat(first,second,little^.dna); cleardna(first); cleardna(second); end else begin (* 1 s->e n not inverting *) if showcase then writeln(output, '1 s->e n not inverting'); getdnasegment(big^.dna,little^.dna,s,e) end end else if s > e then begin if inverting then begin (* 1 e<-s n inverting *) if showcase then writeln(output,' 1 e<-s n inverting'); getdnasegment(big^.dna,little^.dna,e,s); { writeln(output,'gotsegment little^.dna'); } invert(little^.dna) (* ;writeln(output,'inverted little^.dna'); ;crash;{zzzddd} *) end else begin (* 1->e s->n not inverting *) if showcase then writeln(output, '1->e s->n not inverting'); first := nil; second := nil; getdnasegment(big^.dna,first,s,n); getdnasegment(big^.dna,second,1,e); dnaconcat(first,second,little^.dna); cleardna(first); cleardna(second); end end else begin (* s = e, so get one base *) getdnasegment(big^.dna,little^.dna,s,s); if inverting then invert(little^.dna) end (* ; writeln(output, 'about to circledna -----BIG------------'); {zzzddd} LIKEbwdna(output, big^.dna); {zzzddd} writeln(output, 'about to circledna -----LITTLE---------'); {zzzddd} LIKEbwdna(output, little^.dna); {zzzddd} writeln(output, 'about to circledna -----DONE-----------'); {zzzddd} *) end; procedure compress(var pie: pieceptr); (* after deletion of bases, it is possible to have extra dna parts at the end of a sequence. These must be removed so that they don't get accidentally used. *) var p: integer; (* the last base of the dna part *) pielength: integer; (* current length of the piece *) previous: integer; (* length of the all previous segments *) workdna: dnaptr; (* current working dna segment *) begin { writeln(output,'compresssssssssssssssssssssssssssssssssssss'); bwpie(output,pie); } pielength := piecelength(pie); (* find end of used pie *) workdna:=pie^.dna; p := workdna^.length; previous := 0; while (pielength > p) and (workdna^.next <> nil) do begin { writeln(output,'NNNNNNNNNNNNNNNNNNNNnnooooooooooooooooooooo'); } previous := previous + workdna^.length; workdna := workdna^.next; p := p + workdna^.length; end; (* At this point, p is the total length including the workdna segment. p is larger than pielength. previous is the amount prior to this link. So we need to set the workdna length to pielength - previous. *) { writeln(output, 'pielength = ',pielength:1); writeln(output, 'p = ',p:1); writeln(output, 'pielength - previous = ',(pielength - previous):1); } workdna^.length := pielength - previous; if workdna^.next <> nil then emptydna(workdna^.next); { bwpie(output,pie); writeln(output,'END sssssssssssssssssssssssssssssssssssss'); } {zzzyyy} end; (* begin module fixpiececoordinate *) procedure fixpiececoordinate(var pie: pieceptr; excess: integer; coordinateside: direction); (* Fix the piece coordinates for insertions or deletions. Coordinateside is the end of the coordinate system that gets changed *) begin { writeln(output,'fixpiececoordinate: excess = ',excess:1); } with pie^.key do begin case coordinateside of minus: begin case piedir of minus: begin (* piedir minus coordinateside minus *) pieend := pieend - excess; coobeg := coobeg - excess; end; plus: begin (* piedir plus coordinateside minus *) piebeg := piebeg - excess; coobeg := coobeg - excess; end; end; end; plus: begin case piedir of minus: begin (* piedir minus coordinateside plus *) piebeg := piebeg + excess; cooend := cooend + excess; end; plus: begin (* piedir plus coordinateside plus *) pieend := pieend + excess; cooend := cooend + excess; end; end; end; end; end; end; (* end module fixpiececoordinate version = 7.52; {of delmod.p 2000 Jul 30} *) (* begin module book.putbase *) procedure putbase(b: base; position: integer; var pie: pieceptr; coordinateside: direction); (* put a base b into the nth position (internal coordinates) of the piece. Protection is made against positions outside the piece. NEW: * If the base is before coordinate 1, the program halts. * If the base is after the end of the sequence, extra space is made, and the coordinate system is changed on the coordinate side given. *) var excess: integer; (* the implied insertion size *) workdna: dnaptr; (* working dna segment *) p: integer; (* the last base of the dna part or current length *) pielength: integer; (* current length of the piece *) { z: integer; (* index to the piece for debuging stuff *) } begin pielength := piecelength(pie); { writeln(output,'putbase =================================='); write (output,'putbase: b = ',basetochar(b),', pielength=',pielength:1); writeln(output,', position=',position:1); } if position < 1 then begin writeln(output,'putbase: can not put bases', ' before the start of the piece (position < 1)'); writeln(output,'Program error! Please report it to', ' toms@ncifcrf.gov.'); halt; end; workdna:=pie^.dna; if position > pielength then begin (* add to end of piece *) { writeln(output,'putbase: INCREMENT WORKDNA'); } (* since position > pielength, excess is always positive *) excess := position - pielength; fixpiececoordinate(pie,excess,coordinateside); { writeln(output,'putbase: excess=',excess:1); } (* find the last segment *) p := workdna^.length; while workdna^.next <> nil do begin { writeln(output,'looping'); } workdna := workdna^.next; p := p + workdna^.length; end; { writeln(output,'USED SPACE is p=',p:1); } if workdna^.length + excess <= dnamax then begin (* fill into the available segment *) { writeln(output,'NO NEW SEGMENTS NEEDED'); } workdna^.length := workdna^.length + excess; workdna^.part[workdna^.length] := b end else begin (* add new segments *) { write(output,'NEW SEGMENTS NEEDED:'); write(output,' p =',p:1); write(output,' excess =',excess:1); writeln(output); } (* make the rest of the last segment useable by increasing its length up to dnamax, first increment p: *) p := p + (dnamax - workdna^.length); workdna^.length := dnamax; (* build enough segments to accommodate the new sequence *) while p < position do begin { writeln(output,' start adding segment, p = ',p:1,'--------------++++++++++'); } (* now make a new segment: *) p := p + dnamax; getdna(workdna^.next); workdna := workdna^.next; workdna^.length := dnamax; (* make it all available to fill *) { for z := 1 to dnamax do workdna^.part[z] := x; writeln(output,'----------add segment, p = ',p:1); } end; (* put the base in at the end *) workdna^.length := dnamax -(p-position); workdna^.part[workdna^.length] := b { ;writeln(output,'placement, workdna^.length = ',workdna^.length:1); showsegments(output,pie^.dna); } end; end { ;writeln(output,'putbase: final workdna^.length=',workdna^.length:1); } else begin (* insert into the current piece *) { writeln(output,'putbase: NO EXCESS'); } (* locate segment in which to put the base *) p := workdna^.length; while position > p do begin { writeln(output,'THE PLACE =============!! p = ',p:1); } if workdna^.next <> nil then workdna := workdna^.next; p := p + workdna^.length; (* only count the filled part *) end; { writeln(output,'FINAL =============== p = ',p:1); writeln(output,'FINAL =============== position=',position:1); writeln(output,'FINAL =============== (position > p)=', (position > p)); writeln(output,'FINAL =============== (workdna^.next <> nil)=', (workdna^.next <> nil)); ;writeln(output,'putbase: about to smash: workdna^.length=',workdna^.length:1); } workdna^.part[ workdna^.length - (p-position) ] := b; { ;writeln(output,'putbase: after smash: workdna^.length=',workdna^.length:1); } { writeln(output,'base ',b,' placed into ', workdna^.length - (p-position):1); showdnasegment(output,workdna, workdna^.length); writeln(output); } end; { showsegments(output,pie^.dna); ;writeln(output,'putbase: p=',p:1); ;writeln(output,'putbase: dnamax=',dnamax:1); ;writeln(output,'putbase: position=',position:1); ;writeln(output,'=========================================='); } end; (* end module book.putbase version = 7.52; {of delmod.p 2000 Jul 30} *) function getbasech(position: integer; pie: pieceptr):char; (* get a base from the piece at position, return as a character. Use internal coordinates. *) begin getbasech := basetochar(getbase(position, pie)) end; procedure putbasech(c: char; position: integer; pie: pieceptr; coordinateside: direction); (* put the character c as a base at the position in piece pie. Use internal coordinates. Change the coordinate system on the coordinate side given. *) begin putbase(chartobase(c), position, pie, coordinateside) end; (* begin module numbar.info *) (* numbar: routines that make it easy to write a bar of numbers out on the page. these are very convenient for making graphs. the function numberdigit takes a number as input, and the log (base ten) of the place value (called logplace) desired. a character for that spot is returned. example: numberdigit(246,2) = 2 the function numbersize determines the number of characters required for printing the number (including the sign). the procedure numberbar uses numberdigit and numbersize to write a bar of numbers into a file, with several spaces before. the range of numbers is specified, and the number of lines written is returned. example: numberbar(output,5,-20,20,lines); gives ----------- +++++++++++ 21111111111--------- +++++++++11111111112 09876543210987654321012345678901234567890 3 lines were written. *) (* end module numbar.info version = 4.47; (@ of prgmod.p 2000 Sep 9 *) (* begin module numberdigit *) function numberdigit(number, logplace:integer): char; (* return the digit at the place value ('logplace') position of number. example: numberdigit(13625, 3) = 3 numberdigit(13625, 4) = 1 2000 July 30 'myabsolute' replaced 'absolute', which is apparently a keyword for GPC. The name is kept to keep the code looking similar to its origin. *) var place: integer; (* the exponent of logplace *) count: integer; (* used to make place *) myabsolute: integer; (* the absolute value of number *) acharacter: char; (* the character to be returned *) procedure digit; (* extract a digit at the place position *) var tenplace: integer; (* ten times place *) z: integer; (* an intermediate value *) d: integer; (* the digit extracted *) begin (* digit *) tenplace:=10*place; z:=myabsolute-((myabsolute div tenplace)*tenplace); if place = 1 then d:=z else d:= z div place; case d of 0: acharacter:='0'; 1: acharacter:='1'; 2: acharacter:='2'; 3: acharacter:='3'; 4: acharacter:='4'; 5: acharacter:='5'; 6: acharacter:='6'; 7: acharacter:='7'; 8: acharacter:='8'; 9: acharacter:='9'; end end; (* digit *) procedure sign; (* put a negative sign out or a positive sign *) begin (* sign *) if number <0 then acharacter:='-' else acharacter:='+' end; (* sign *) begin (* numberdigit *) place:=1; for count:=1 to logplace do place:=10*place; if number=0 then begin if place=1 then acharacter:='0' else acharacter:=' ' end else begin myabsolute:=abs(number); if myabsolute < (place div 10) then acharacter:=' ' else if myabsolute >= place then digit else sign end; numberdigit:=acharacter end; (* numberdigit *) (* end module numberdigit version = 4.47; (@ of prgmod.p 2000 Sep 9 *) (* begin module numbersize *) function numbersize(n: integer):integer; (* calculate amount of space to be reserved for the integer n *) const ln10 = 2.30259; (* natural log of 10 - for conversion to log base 10 *) epsilon = 0.00001; (* a small number to correct log base 10 errors *) var size: integer; (* intermediate result *) begin (* numbersize *) if n = 0 then numbersize:=1 else begin size:=trunc(ln(abs(n))/ln10 + epsilon) + 1; (* the 1 is for the last digit *) (* the epsilon assures that we do not lose a place due to roundoff. eg, sometimes log base 10 of 10 would be 0.9999 instead of 1, and we would not do it right... note: this will fail for very large numbers on the order of 1/epsilon. *) if n < 0 then size := succ(size); (* account for minus sign *) numbersize := size; end end; (* numbersize *) (* end module numbersize version = 4.47; (@ of prgmod.p 2000 Sep 9 *) procedure namechangeset(var pie: pieceptr; changes: changeset); (* put the changeset into the piece name *) (* this could be upgraded to be more general like linechangeset *) var i: integer; (* index to insertion *) n: integer; (* index to changes *) procedure putin(c: char); (* put c into the name at s *) begin { writeln(output,'putin(',c,')'); } with pie^.key.hea.keynam do if length < namelength then begin length := succ(length); letters[length] := c end; end; procedure putnumber(r: real); var c: char; (* a character in the number *) n: integer; (* index to the characters in number *) i: integer; (* round of r *) begin i := round(r); for n := numbersize(i)-1 downto 0 do begin c := numberdigit(i,n); { writeln(output,'i=',i); writeln(output,'c=',c); writeln(output,'n=',n); } if c <> '+' then putin(c); { writename(output,pie^.key.hea.keynam); } end; end; begin with changes do for n := 1 to number do begin putin('.'); with data[n] do case changetype of 'c': begin putin(baseold); putnumber(basecoo1); putin(basenew); end; 'i': begin putin('i'); putnumber(basecoo1); putin(','); putnumber(basecoo2); {zzz111} if insertcomplement then putin('~'); for i := 1 to inserts do putin(insert[i]); end; 'd': begin putin('d'); putnumber(basecoo1); putin(','); putnumber(basecoo2); end; end; end; end; procedure putinline(var pie: pieceptr; c: char); (* put c onto the end of the full name of pie *) begin with pie^.key.hea.fulnam^ do if length < linelength then begin length := succ(length); letters[length] := c end; end; procedure putnumline(var pie: pieceptr; r: real); (* put a number onto the end of the full name of pie *) var c: char; (* a character in the number *) n: integer; (* index to the characters in number *) i: integer; (* round of r *) begin i := round(r); for n := numbersize(i)-1 downto 0 do begin c := numberdigit(i,n); { writeln(output,'i=',i); writeln(output,'c=',c); writeln(output,'n=',n); } if c <> '+' then putinline(pie, c); end; end; procedure linechangeset(var pie: pieceptr; changes: changeset); (* put the changeset into the piece fulnam *) var i: integer; (* index to insertion *) n: integer; (* index to changes *) begin with changes do for n := 1 to number do begin if n > 1 then begin putinline(pie,','); putinline(pie,' '); end; with data[n] do case changetype of 'c': begin putinline(pie,'a'); putinline(pie,'t'); putinline(pie,' '); putnumline(pie,basecoo1); putinline(pie,' '); putinline(pie,baseold); putinline(pie,'-'); putinline(pie,'>'); putinline(pie,basenew); end; 'i': begin putinline(pie,'i'); putinline(pie,'n'); putinline(pie,'s'); putinline(pie,'e'); putinline(pie,'r'); putinline(pie,'t'); putinline(pie,' '); if insertcomplement then putinline(pie,'~'); if inserts > 0 then for i := 1 to inserts do putinline(pie,insert[i]) else begin putinline(pie,'N'); (* NOTHING *) putinline(pie,'O'); putinline(pie,'T'); putinline(pie,'H'); putinline(pie,'I'); putinline(pie,'N'); putinline(pie,'G'); end; putinline(pie,' '); putinline(pie,'b'); putinline(pie,'e'); putinline(pie,'t'); putinline(pie,'w'); putinline(pie,'e'); putinline(pie,'e'); putinline(pie,'n'); putinline(pie,' '); putnumline(pie,basecoo1); putinline(pie,' '); putinline(pie,'a'); putinline(pie,'n'); putinline(pie,'d'); putinline(pie,' '); putnumline(pie,basecoo2); end; 'd': begin putinline(pie,'d'); putinline(pie,'e'); putinline(pie,'l'); putinline(pie,'e'); putinline(pie,'t'); putinline(pie,'e'); putinline(pie,' '); putnumline(pie,basecoo1); putinline(pie,' '); putinline(pie,'t'); putinline(pie,'o'); putinline(pie,' '); putnumline(pie,basecoo2); end; end; end; end; (* linechangeset *) function realbetween(a,b,c: real): boolean; (* is b between a and c? *) (* this is an inclusive between *) begin realbetween:=((a<=b) and (b<=c)) or ((c<=b) and (b<=a)) end; procedure checkoneoverlap(c: changeset; x: integer); (* check that one overlap in the changeset c does not overlap any other changes. *) var y: integer; (* another change *) x1, x2: real; (* internal coordinates for x *) y1, y2: real; (* internal coordinates for x *) procedure setrange(n: integer; var first, last: real); (* Determine the range (first to last) for then nth change, adjusting it so that insertions are for the region affected, just like deletions and changes. *) begin with c.data[n] do begin first := internal1; last := internal2; if changetype = 'i' then begin if first + 1 = last then begin first := first + 0.5; last := last - 0.5; end else begin first := first + 1; last := last - 1; end end end { ; writeln(output,'first = ',first:1); ; writeln(output,'last = ',last :1); } end; (* setrange *) begin { write(output,'checkoverlap: '); writechangeset(output,c); writeln(output); } { writeln(output,'checkoverlap: x = ',x:1); } setrange(x,x1,x2); for y := 1 to x-1 do begin { writeln(output,'checkoverlap: y=',y:1); } setrange(y,y1,y2); if realbetween(x1,y1,x2) or realbetween(x1,y2,x2) then begin {EEE} mutcd1 := c.data[x]; mutcd2 := c.data[y]; write(output,'In '); writechangeset(output,c); write(output,' '); wchangedata(output,c.data[x]); write(output,' overlaps '); wchangedata(output,c.data[y]); writeln(output); error(222); end {zzzppp} end end; procedure changesequence(changes: changeset; var thesequence: pieceptr; coordinateside: direction; var deleted: boolean); (* Make changes to the sequence. Setinternal must be called prior to this procedure to set up the internal variables used. Account for the rule in libdef that inserts have the direction as the coordinate system. Change the coordinate system on the coordinate side given. If the sequence was deleted by mutation(s), set deleted true. *) var b1, b2: integer; (* internal coordinate values for the basecoo1 and basecoo2 of the changeset *) deletes: integer; (* amount to delete the sequence *) i: integer; (* counter for inserts *) n: integer; (* counter for the changes *) shift: integer; (* amount to shift a portion of the sequence *) pielength: integer; (* current length of the piece *) begin { writeln(marksdelila,'*4**************************************************'); write(marksdelila,'* SSSSSSSSSSS changesequence: '); write(marksdelila,'* change: '); writechangeset(marksdelila,changes); writeln(marksdelila); } namechangeset(thesequence, changes); pielength := piecelength(thesequence); { checkoverlap(changes, thesequence); } deleted := false; (* First check that the c-type changes match the sequence. In the new scheme, this is done against the normal sequence. Since propagation is done in the next loop, thesequence changes, so it can't be done there. *) with changes do for n := 1 to number do with data[n] do begin checkoneoverlap(changes, n); b1 := internal1; b2 := internal2; case changetype of 'c': begin if baseold = basenew then begin write(output,'You wrote '); writechangeset(output,changes); writeln(output,' '); writeln(output,'but the initial and final bases are the same,'); writeln(output,'so you did not request any change!'); error(214); geteoinst; end else begin if thesequence^.key.coodir <> thesequence^.key.piedir then begin (* AaaHa! They are working on the complementary strand, so switch around and complement the bases!! *) { write (marksdelila,'* switching baseold and basenew:'); write (marksdelila,'* baseold = ',baseold); writeln(marksdelila,'* basenew = ',basenew); } baseold := basetochar(complement(chartobase(baseold))); basenew := basetochar(complement(chartobase(basenew))); { write (marksdelila,'* switched baseold and basenew to:'); write (marksdelila,'* baseold = ',baseold); writeln(marksdelila,'* basenew = ',basenew); } end; if (b1 < 1) then begin writeln(output,'A requested mutation coordinate is', ' off the piece in the 5'' direction'); mutposition1 := basecoo1; (* mutation position *) error(215); geteoinst; end else if (b1 > pielength) then begin writeln(output,'A mutation requested coordinate is', ' off the piece in the 3'' direction'); mutposition1 := basecoo1; (* mutation position *) error(216); geteoinst; end else if baseold <> getbasech(b1,thesequence) then begin (* it must be a wrong base *) mutposition1 := basecoo1; (* mutation position *) mutnotchar := baseold; (* what mutation is not *) mutischar := getbasech(b1,thesequence); (* what mutation is *) if thesequence^.key.coodir <> thesequence^.key.piedir then begin (* AaaHa! They are working on the complementary strand, so switch around and complement the bases!! *) { write (marksdelila,'* switching baseold and basenew:'); write (marksdelila,'* baseold = ',baseold); writeln(marksdelila,'* basenew = ',basenew); } baseold := chomplement(baseold); basenew := chomplement(basenew); mutischar := chomplement(mutischar); mutnotchar := chomplement(mutnotchar); { write (marksdelila,'* switched baseold and basenew to:'); write (marksdelila,'* baseold = ',baseold); writeln(marksdelila,'* basenew = ',basenew); } end; writeln(output,'On the positively oriented strand,', ' the old base at ',basecoo1:1, ' is NOT ',baseold,'!', ' It is ',mutischar,'.'); (* mutation: old base is not correctly identified *) error(211); geteoinst; end; end; end; 'i': ; (* no checks here at this time *) 'd': ; (* no checks here at this time *) end; end; (* make changes to thesequence and propagate them one by one in the changes changeset. *) with changes do for n := 1 to number do with data[n] do begin {zzzyyy} { writeln(output,'* changeloop n= ',n:1,'------------------------------'); writeln(marksdelila,'* changesequence: basecoo1 = ',basecoo1:1); writeln(marksdelila,'* changesequence: basecoo2 = ',basecoo2:1); writeln(marksdelila,'* changesequence: internal1 = ',internal1:1,' old'); writeln(marksdelila,'* changesequence: internal2 = ',internal2:1,' old'); writeln(marksdelila,'* changesequence: within(thesequence,basecoo1) = ', within(thesequence,basecoo1):1); writeln(marksdelila,'* changesequence: within(thesequence,basecoo2) = ', within(thesequence,basecoo2):1); } b1 := internal1; b2 := internal2; case changetype of 'c': begin if ((b1 < 1) or (b2 > pielength)) then error(219) else putbasech(basenew, b1, thesequence,coordinateside); end; 'i': begin if (b1 < 0) and (b2 > pielength) and (inserts = 0) then begin error(220); (* 2005 Sep 6 *) geteoinst; deleted := true; end else begin { writeln(output,'start of insert: b1= ',b1:1,' b2= ',b2:1); } (* make sure insertion is at the end of the piece *) if b1 > pielength then b1 := pielength; if b2 > pielength then b2 := pielength + 1; if b1 < 0 then b1 := 0; if b2 < 1 then b2 := 1; (* original *) { writeln(output,'after correction: b1= ',b1:1,' b2= ',b2:1); } if b1 = b2 then begin writeln(output,' The first base, ',basecoo1:1, ', must not equal ', ' the second base, ',basecoo2:1, ' for insertion'); mutposition1 := basecoo1; (* mutation position *) mutposition2 := basecoo2; (* mutation position *) error(218); geteoinst; end; deletes := b2 - b1 - 1; shift := inserts - deletes; pielength := pielength + shift; { writeln(marksdelila,'* changesequence: inserts= ',inserts:1); writeln(marksdelila,'* changesequence: shift= ',shift:1); writeln(marksdelila,'* changesequence: new pielength= ',pielength:1); } if shift > 0 then begin (* insert *) { writeln(marksdelila,'* changesequence: INSERTION'); } (* shift the rest of the sequence out of the way *) (* if insertion is at the 5' end, don't try to grab base off the end! *) { writeln(output,'old b2= ',b2:1); writeln(output,'shift= ',shift:1); } case thesequence^.key.piedir of plus: if b2-shift < 1 then b2 := b2 + shift; minus: if b2-shift < 1 then b2 := b2 + shift; end; { if thesequence^.key.piedir = plus then writeln(output,'PLUS') else writeln(output,'MINUS'); } { zzzyyy writeln(output,'basecoo1= ',basecoo1:4:1); writeln(output,'basecoo2= ',basecoo2:4:1); writeln(output,'new b2= ',b2:1); writeln(output,'b1= ',b1:1); writeln(output); writeln(output,'BEFORE:'); bwpie(output, thesequence); writeln(output,'-- pielength= ',pielength:1); (* the NEW piece length *) (* piecelength will be the actual length after the first insert *) } {zzzCCC} for i := pielength downto b2 do putbasech(getbasech(i-shift,thesequence), i,thesequence,coordinateside); { writeln(output,'AFTER:'); bwpie(output, thesequence); for i := pielength downto b2 do write(output,getbasech(i-shift,thesequence)); writeln(output); } {zzznnn} { writeln(output); writeln(output,'AFTER'); bwpie(output, thesequence); } { (* wrecking: *) for i := pielength downto b2 do putbasech('g', i,thesequence,coordinateside); for i := 2 to pielength do putbasech('g', i,thesequence,coordinateside); writeln(output); writeln(output,'AFTER WRECKING'); bwpie(output, thesequence); } end else if shift < 0 then begin (* delete *) { writeln(marksdelila,'* changesequence: insert type DELETION'); writeln(marksdelila,'* BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB'); writeln(marksdelila,'* shift < 0 in changesequence'); writeln(marksdelila,'* basecoo1= ',basecoo1:4:1); writeln(marksdelila,'* basecoo2= ',basecoo2:4:1); writeln(marksdelila,'* new b2= ',b2:1); writeln(marksdelila,'* b1= ',b1:1); writeln(marksdelila); } { writeln(output,'BEFORE:'); bwpie(output, thesequence); } {BBB} for i := b2+shift to pielength do begin putbasech(getbasech(i-shift,thesequence), i,thesequence,coordinateside); end; (* shift < 0 so to make deletion, take the positive: *) fixpiececoordinate(thesequence, shift,coordinateside); { writeln(output,'AFTER:'); bwpie(output, thesequence); writeln(output,'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'); } if shift < 0 then compress(thesequence); end else begin (* shift = 0 *) if inserts = 0 then begin {EEE} error(219); end; end; { writeln(output); if not (thesequence^.key.coodir = thesequence^.key.piedir) then write(output,'* NOT') else write(output,'* '); writeln(output,' (thesequence^.key.coodir = thesequence^.key.piedir)'); if not insertcomplement then write(output,'* NOT') else write(output,'* '); writeln(output,' insertcomplement'); } (* insert the new material *) (* account for the rule that the insert has the same orientation as the coordinate system. Let same := (thesequence^.key.coodir = thesequence^.key.piedir); Then there are 4 cases: insertcomplement false true same true Homologous Complement same false Complement Homologous This logic can be done by testing that the booleans are equal or not (thanks to Karen Lewis for this suggestion!): *) if not ((thesequence^.key.coodir = thesequence^.key.piedir) = insertcomplement) {zzz111} then begin (* HOMOLOGOUS case *) for i := 1 to inserts do putbasech(insert[i], b1 + i, thesequence,coordinateside); { writeln(output,'* changesequence: EQUAL => HOMOLOGOUS'); writeln(marksdelila,'* changesequence: EQUAL => HOMOLOGOUS'); } end else begin (* COMPLEMENT case *) for i := 1 to inserts (* the order of insertion doesn't matter. What counts here is that the bases be put in in reverse position. This solves the bug of 2000 May 25 where inserts were not reversed *) do putbasech( chomplement(insert[i]), b1 + (inserts - i + 1), (* reverse insertion order! *) thesequence,coordinateside); { writeln(output,'* changesequence: FLIP => COMPLEMENT'); writeln(marksdelila,'* changesequence: FLIP => COMPLEMENT'); } end { writeln(output); writeln(output,'piecelength(thesequence)=',piecelength(thesequence):1); writeln(output,'AFTER INSERTION'); bwpie(output, thesequence); } end end; 'd': begin (* delete inclusively from b1 to b2 *) {DDD} { writeln(output,'changesequence: delete: basecoo1 = ',basecoo1:1); writeln(output,'changesequence: delete: basecoo2 = ',basecoo2:1); writeln(output,'changesequence: delete: b1 = ',b1:1); writeln(output,'changesequence: delete: b2 = ',b2:1); } if ((b1 <= 1) and (b2 >= pielength)) then begin error(220); (* 2005 Sep 6 *) geteoinst; deleted := true; end else if ( (b1 > pielength) and (b2 > pielength)) or (((b1 < 1) and (b2 < 1))) then begin (* b1 and b2 were both off the same end, consider shifting the piece coordinates. *) error(221); { writeln(output,'changesequence: delete: b1 = ',b1:1); writeln(output,'changesequence: delete: b2 = ',b2:1); writeln(output,'changesequence: delete: shift = ',shift:1); } { Not implemented: see notes in shiftpiececoordinate Sss if (coordinateside = minus) and ((b1 >= pielength) and (b2 >= pielength)) then begin shiftpiececoordinate(thesequence, basecoo1, basecoo2, minus); end else if (coordinateside = plus) and ((b1 <= 1) and (b2 <= 1)) then begin shiftpiececoordinate(thesequence, basecoo1, basecoo2, plus); end (* otherwise don't do the shift *) } end else begin if b1 < 1 then b1 := 1; if b2 < 1 then b2 := 1; if b1 > pielength then b1 := pielength; if b2 > pielength then b2 := pielength; { writeln(output,'changesequence: delete: b1 = ',b1:1); writeln(output,'changesequence: delete: b2 = ',b2:1); } shift := b2 - b1 + 1; pielength := pielength - shift; for i := b1 to pielength do putbasech(getbasech(i+shift,thesequence),i, thesequence,coordinateside); (* account for length change *) (* shift > 0 so to make deletion, take the negative: *) fixpiececoordinate(thesequence, -shift,coordinateside); {zzzyyy} { writeln(output,'BEFORE compress:'); bwpie(output, thesequence); halt; } if shift > 0 then compress(thesequence); { writeln(output,'AFTER compress:'); bwpie(output, thesequence); } end end; end; {zzzppp} { writeln(output); writeln(output,'changesequence: about to propagateonechange, n = ',n:1); } propagateonechange(changes,n,thesequence); { write(output,'changesequence: changes = '); writechangeset(output,changes); writeln(output); } end; { write(marksdelila,'* EEEEEEEEEE, changesequence: '); writechangeset(marksdelila,changes); writeln(marksdelila); writeln(output,'END OF changesequence:'); bwpie(output, thesequence); } {zzzyyy} end; (* changesequence *) procedure extractpiece(libpie: pieceptr; var pie: pieceptr; cutposition: integer {; not implemented var markers: markerptr}); (* extract the piece pie from the libpie *) (* formerly: lrpiedna, library read piece dna *) (* at the onset, the reading point is just before the first base of the piece. there are three major cases for reading in the dna. case 1. a. this is the simplest. a linear piece in the library can only have a linear fragment in the book. action: skip to the start of the region of interest. read the dna. case 1. b. a circle of dna is stored linearly in the library (and labeled 'circular'). if one wants a fragment of a circle and the fragment lies completely within the stored dna, then one still can apply case 1. case 2. one wants all of a circular sequence. one just reads the whole dna, and breaks at the requested cut point. case 3. this is the most complicated. one wants a linear fragment of a circle, but the fragment happens to lie across the boundary cut by storing the sequence linearly in the library (acrossboundary). (in case 1b, it could not lie across the boundary since its parent is linear.) one must: read from the piece beginning base to the stop of the fragment. skip around to the start of the fragment. read up to the piece ending base. fuse the two fragments into the final product. linear in library --> linear request --> case 1 or circular in library --> circular request --> case 2 or linear request --> across boundary --> case 3 or not across boundary --> case 1 *) var reversing: boolean; (* if true reverse the dna *) acrossboundary: boolean; (* request is across the boundary of a circular dna *) (* the direction of the piece of dna being read is positive with respect to its coordinate system: *) startread, stopread: integer; (* range of bases to read in the library *) procedure setvariables; (* set up variables *) begin (* we will complement if the direction was changed *) reversing:=(pie^.key.piedir <> libpie^.key.piedir); (* figure out acrossboundary: *) if (pie^.key.coocon = circular) and (pie^.key.piecon = linear) then acrossboundary:=within(pie, libpie^.key.piebeg) and within(pie, libpie^.key.pieend) else acrossboundary:=false; (* acrossboundary may not be right yet, since we must handle the special case where we want a linearized form of the circle: the piece and library piece ends will be the same, so we want to do case 1, not case 3 *) if acrossboundary then if reversing then begin if (libpie^.key.piebeg=pie^.key.pieend) and (libpie^.key.pieend=pie^.key.piebeg) then acrossboundary:=false end else begin if (libpie^.key.piebeg=pie^.key.piebeg) and (libpie^.key.pieend=pie^.key.pieend) then acrossboundary:=false end; (* get first dna *) getdna(pie^.dna); { writeln(output,'setvariables:'); writeln(output,'acrossboundary is ', acrossboundary); writeln(output,'pie^.key.piecon is ', pie^.key.piecon); writeln(output,'pk^.key.piecon is ', pk^.key.piecon); writeln(output,'libpie^.key.piecon is ', libpie^.key.piecon); } end; procedure case1; begin if debugging then writeln(debug,'case1'); { writeln(output,'==== BEG case1'); } startread := pietoint(pie^.key.piebeg,libpie); stopread := pietoint(pie^.key.pieend,libpie); { writeln(output,'pie^.key.piebeg = ', pie^.key.piebeg:1,' TAG'); writeln(output,'pie^.key.pieend = ', pie^.key.pieend:1,' TAG'); writeln(output,'startread = ',startread:1); writeln(output,' stopread = ',stopread:1); writeln(output,'case1: reversing = ',reversing); } circledna(libpie, pie, startread, stopread, reversing); end; procedure case2; (* circular in library --> circular request --> case 2 *) begin if debugging then writeln(debug,'case2'); { writeln(output,'case2'); } if reversing then begin startread:=pietoint(cutposition,libpie); if startread < piecelength(libpie) then stopread := startread + 1 else stopread := 1; circledna(libpie, pie, startread, stopread, reversing); end else begin startread:=pietoint(cutposition,libpie); if startread > 1 then stopread := startread - 1 else stopread := piecelength(libpie); circledna(libpie, pie, startread, stopread, reversing); end; end; procedure case3; (* circular in library --> linear request --> across boundary --> case 3 *) begin if debugging then writeln(debug,'case3'); { writeln(output,'case3'); } startread:=pietoint(pie^.key.piebeg,libpie); stopread :=pietoint(pie^.key.pieend,libpie); circledna(libpie, pie, startread, stopread, reversing); end; begin (* extractpiece *) { writeln(output,'BEGIN extractpiece'); writeln(output,'pie^.key.piecon is ', pie^.key.piecon); writeln(output,'pk^.key.piecon is ', pk^.key.piecon); writeln(output,'libpie^.key.piecon is ', libpie^.key.piecon); } setvariables; if acrossboundary then case3 else if (pie^.key.piecon = linear) then case1 else case2; { writeln(output,'OUT extractpiece'); } end; (* extractpiece *) procedure gozero(var Apiece: pieceptr); (* Transform a piece to a zero coordinate system. *) begin if def.coo = coorzero then begin Apiece^.key.coocon := linear; Apiece^.key.coodir := plus; case Apiece^.key.piedir of plus: begin zeroBS := zerobase-zeroshift; Apiece^.key.coobeg := Apiece^.key.piebeg-zeroBS; Apiece^.key.cooend := Apiece^.key.pieend-zeroBS; Apiece^.key.piecon := linear; Apiece^.key.piedir := plus; Apiece^.key.piebeg := Apiece^.key.piebeg-zeroBS; Apiece^.key.pieend := Apiece^.key.pieend-zeroBS; end; minus: begin zeroBS := zerobase+zeroshift; Apiece^.key.coobeg := zeroBS-Apiece^.key.piebeg; Apiece^.key.cooend := zeroBS-Apiece^.key.pieend; Apiece^.key.piecon := linear; Apiece^.key.piedir := plus; Apiece^.key.piebeg := zeroBS-Apiece^.key.piebeg; Apiece^.key.pieend := zeroBS-Apiece^.key.pieend; end; end; end; end; (**************************************************************************) (**************************************************************************) procedure addnumber(var hea: header; n: integer); (* add a number n to the header *) const precision = 1e-7; (* numbers are rounded to this precision. There was a bug on a Sun 4, where the 1000th piece was labeled 000 because of truncation. *) var numbernote: lineptr; (* the number of this item, for insertion into the notes of the header *) signspace: integer; (* the number of spaces reserved for the sign of the number *) numdigits, (* number of digits in the number *) digit, (* a digit *) index: integer; charnum: char; noteindex: lineptr; (* used to scan the header note *) begin getline(numbernote); with numbernote^ do begin for index:=1 to numberlength do letters[index]:=numberword[index]; if n<0 then begin (* handle negative sign *) letters[numberlength+1]:='-'; signspace:=1; n:=-n end else signspace:=0; (* for positive n *) if n=0 then numdigits:=1 (* old form that caused the bug: else numdigits:=trunc(ln(n)/ln(10)+1); *) (* the new form that avoids the bug by causing rounding from point precision *) else numdigits:=trunc(ln(n)/ln(10)+1+precision); length:=numberlength+signspace+numdigits; for index:=length downto numberlength+signspace+1 do begin digit:= n mod 10; n:= n div 10; case digit of 0: charnum:='0'; 1: charnum:='1'; 2: charnum:='2'; 3: charnum:='3'; 4: charnum:='4'; 5: charnum:='5'; 6: charnum:='6'; 7: charnum:='7'; 8: charnum:='8'; 9: charnum:='9'; end; letters[index]:=charnum end; (* insert number before the other notes *) next:=hea.note; (* end of number points to header *) hea.note:=numbernote (* header starts with the number *) end; (* insert user notes below the other notes *) if usernotes <> nil then begin if hea.note<>nil then begin (* point to first note *) noteindex:=hea.note; (* look for last note *) while noteindex^.next<>nil do noteindex:=noteindex^.next; (* link user notes to end of last header note *) noteindex^.next:=usernotes end else hea.note:=usernotes; (* clear usernotes *) usernotes:=nil end end; (* specification *) procedure spec(var hea: header; numberable: state; var object: real); (* specify the object and maybe give it a number and the current usernotes *) begin (* if library notes are off, then destroy them *) if def.key.note = off then while hea.note<>nil do clearline(hea.note); if numberable=on then begin def.num.item:=succ(def.num.item); object:=def.num.item; ivaluenumber(trunc(object)) end else object:=-0.5; (* specified but not numbered *) if (def.num.sta = on) and (numberable = on) then begin (* insert a number into the object"s note *) addnumber(hea, def.num.item); end; end; procedure specified(object: real); (* was this object specified? *) begin correct:=object<>0.5 end; function numbered(object: real): boolean; (* determine if this object was numbered *) begin numbered:=(abs(object)<>0.5) end; procedure specpiece(ref: reference); (* specify the piece referred to if it is not currently specified. 2001 Mar 16: read in the LIBRARY piece, libpie *) begin specified(pkspec); if (not correct) or (ref.pienam.letters <> pk^.key.hea.keynam.letters) then lrpiece(ref.pienam, libpie); (* the piece will become specified, but we do not want to give it a number until we get to the request *) pkspec:=-0.5; end; procedure unspec(var object: real); (* unspecify this object *) begin object:=0.5 end; procedure unspecbelowck; (* unspecify the objects below ck *) begin unspec(mkspec); unspec(tkspec); unspec(gkspec); unspec(pkspec); end; (**************************************************************************) (**************************************************************************) procedure checksize(pk: pieceptr); (* check that the size of the piece pk added to previous pieces would not exceed maxbook. *) begin booksize := booksize + piecelength(pk); if booksize > maxbook then begin booksize := 0; (* count of bases *) getcount := 0; (* count of pieces *) error(223); geteoinst; end; end; procedure dopiece; (* called in pass 2 only *) (* put together a call for a piece of dna, using the interface to the library. these global variables must be provided: fromposition, toposition, directionwanted (as a plus or minus value), cutposition, mkon, and configuration (pk^.key.piecon) the above variables are not modified by dopiece. dopiece calculates piedir, piebeg and pieend. modify the sequence according to the changeset. finally, bwpie is called. *) var vdirwanted, (* value (+1 or -1) of direction wanted *) vdirimplied: (* value of direction implied by fromposition *) integer; (* and toposition *) libstart: integer; (* the library coordinate boundary closest to the piece begin *) deleted: boolean; (* the sequence was deleted - don't write it out *) { sortedmutations: changeset; (* the mutations but sorted for printing *) } function dirvalue(d: direction): integer; (* convert a direction to a +1 or a -1 value *) begin case d of plus: dirvalue:=+1; minus: dirvalue:=-1; dirhomologous: dirvalue:=0; dircomplement: dirvalue:=0; end end; function sign(thetimes: integer): integer; (* find the sign of the times *) begin if thetimes >= 0 then sign:=+1 else sign:=-1 end; function BooleanXOR(a,b: boolean):boolean; (* exclusive or function *) begin BooleanXOR:=(a and (not b)) or (b and (not a)) end; begin (* dopiece *) { writeln(output,'dopiece-IN'); writeln(output,'pk^.key.piecon is ', pk^.key.piecon); writeln(output,'libpie^.key.piecon is ', libpie^.key.piecon); } with pk^.key do begin (* we first convert directionwanted into the vdirwanted value. only plus, minus values are honored. *) { if dirwanted in [plus, minus] then vdirwanted:=dirvalue(dirwanted) else vdirwanted:=dirvalue(libpie^.key.piedir); } (* 1999 mar 17:the reason for that restriction is lost in the mists of time! *) if dirwanted in [plus, minus] then vdirwanted:=dirvalue(dirwanted) else if dirwanted = dirhomologous then vdirwanted := +dirvalue(libpie^.key.piedir) else vdirwanted := -dirvalue(libpie^.key.piedir); (* we now calculate the direction implied by the from and to positions: vdirimplied *) if toposition = fromposition then begin vdirimplied:=vdirwanted; (* new as of 1999 mar 17!! *) end else begin case coocon of linear: vdirimplied:=sign(toposition-fromposition); circular: case libpie^.key.piecon of circular: vdirimplied:=vdirwanted; (* can not check it *) linear: if within(libpie,coobeg) and within(libpie,cooend) then begin (* coord. boundary may be in the way *) (* take care of piece direction *) case libpie^.key.piedir of plus: libstart:=libpie^.key.cooend; minus: libstart:=libpie^.key.coobeg; end; (* which side of the boundary are the positions on? if (t and not f) or (f and not t) then opposite *) if BooleanXOR(between(libpie^.key.piebeg,toposition,libstart), between(libpie^.key.piebeg,fromposition,libstart)) (* opposite sides: *) then vdirimplied:=-sign(toposition-fromposition) (* same sides: *) else vdirimplied:=+sign(toposition-fromposition) end else vdirimplied:=sign(toposition-fromposition) end end; case dirwanted of plus, minus: if vdirwanted <> vdirimplied then error(204); dirhomologous: vdirwanted:= vdirimplied; dircomplement: if libpie^.key.piecon=circular then vdirwanted:=-vdirimplied else error(204) end; end; if correct then begin piebeg:=fromposition; pieend:=toposition; case vdirwanted of +1: piedir:=plus; -1: piedir:=minus; end; if debugging then begin writeln(debug,'piebeg:',piebeg:5,' pieend:',pieend:5,' cut:',cutposition:5); writeln(debug,'ord(piecon)=',ord(piecon)); writeln(debug,'ord(piedir)=',ord(piedir)); end; { writeln(output,'piebeg:',piebeg:5,' pieend:',pieend:5,' cut:',cutposition:5); writeln(output,'ord(piecon)=',ord(piecon)); writeln(output,'ord(piedir)=',ord(piedir)); } (* start out with the library full name *) copyline(libpie^.key.hea.fulnam, pk^.key.hea.fulnam); (* manipulate the pk copy of libpie *) extractpiece(libpie,pk,cutposition{,mkon not implemented}); (* convert pk to zero based coordinate system if needed *) gozero(pk); {zzzwww}{with ... get} deleted := false; if mutations.number > 0 then begin (* make any mutations desired *) (* adjust the coordinate system if needed so that writemarks can use it. Changesequence also changes the sequence. *) setinternal(mutations, pk); (* are we doubling? *) if def.doubling = on then begin { write(output,'def.num.str[pienum] is '); if def.num.str[pienum] = on then writeln(output,'on') else writeln(output,'off'); } (* we first write out the wt sequence: *) checksize(pk); (* maxbook *) bwpie(book, pk); getcount := succ(getcount); (* then increment the piece number: *) spec(pk^.key.hea,def.num.str[pienum],pkspec); end; if def.doubling = on then begin writewildtypemarks(marksdelila, mutations, insertupperbits, (* upperbits for insertion symbol *) insertlowerbits, (* lowerbits for insertion symbol *) deleteupperbits, (* upperbits for deletion symbol *) deletelowerbits, (* lowerbits for deletion symbol *) changeupperbits, (* upperbits for change symbol *) changelowerbits, (* lowerbits for change symbol *) pk, (* the piece *) def.num.item); (* the piece number *) (* skip piece corresponding to the piece. Note that at this point, since there are mutations, withused must be true so it is not necessary to test for it here. *) skippiece(marksdelila); end; changesequence(mutations, pk, coordinateside, deleted); (* writeln(output,'AFTER changesequence:');{zzzyyy} bwpie(output, pk);{zzzyyy} *) (* always write mutant marks. This allows multiple changes for one wild type sequence. *) if not deleted then writemutantmarks(marksdelila, mutations, insertupperbits, (* upperbits for insertion symbol *) insertlowerbits, (* lowerbits for insertion symbol *) deleteupperbits, (* upperbits for deletion symbol *) deletelowerbits, (* lowerbits for deletion symbol *) changeupperbits, (* upperbits for change symbol *) changelowerbits, (* lowerbits for change symbol *) pk, (* the piece *) def.num.item); (* the piece number *) clearline(pk^.key.hea.fulnam); (* remove old *) getline(pk^.key.hea.fulnam); (* remove old *) linechangeset(pk, mutations); end; (* if mutation(s) did not delete the entire piece, write it out *) if not deleted then begin (* Put the long name into the piece from the name command. This will override the fulnam given if there were mutations. *) if longnameexists then begin copyline(longname, pk^.key.hea.fulnam); longnameexists := false; (* use it only once *) end; (* skippiece corresponding to the piece *) if withused then skippiece(marksdelila); (* output the result to the book *) checksize(pk); (* maxbook *) bwpie(book, pk); getcount := succ(getcount); end else begin (* skippiece corresponding to the piece *) (* (Note: when deleted is true, withused must be true, but check anyway ... *) if withused then skippiece(marksdelila); end end end ;if debugging then writeln(debug,'dopiece-out') end; (* end of dopiece *) (**************************************************************************) (**************************************************************************) (**************************************************************************) (* instruction reading routines *******************************************) procedure irdirection; (* *) (* read a direction into directionwanted *) begin if not eoinst then begin findword; if correct then begin if chr in ['+','-'] then begin { writeln(output,'First chr: "', chr,'"'); } parse(false); { writeln(output,'after parse: parsedword[1]: "', parsedword[1],'"'); writeln(output,'after parse: parsedword[2]: "', parsedword[2],'"'); } if correct then begin if parsedword[2]=' ' then case parsedword[1] of '+': dirwanted:=plus; '-': dirwanted:=minus; ' ': begin writeln(output, 'irdirection: Delila parse fault'); writeln(output,' chr: "', chr,'"'); writeln(output,' parsedword[1]: "', parsedword[1],'"'); writeln(output,' parsedword[2]: "', parsedword[2],'"'); {zzzfff} halt; end; end else error(2) { else begin writeln(output,'parsedword[1]: "', parsedword[1],'"'); writeln(output,'parsedword[2]: "', parsedword[2],'"'); writeln(listing,'boop'); halt; error(2); end } end end else begin irword; if correct then begin if word in [comdelila,homdelila] then case word of comdelila: dirwanted:=dircomplement; homdelila: dirwanted:=dirhomologous end else if word in [mardelila,tradelila,gendelila,piedelila] then begin if pass=2 then begin case word of mardelila: begin specified(mkspec); if correct then dirwanted:= mkoff^.key.ref.refdir else error(202) end; tradelila: begin specified(tkspec); if correct then dirwanted:=tk.ref.refdir end; gendelila: begin specified(gkspec); if correct then dirwanted:=gk.ref.refdir end; piedelila: begin specified(pkspec); if correct then dirwanted:=pk^.key.piedir end end end end else begin error(7); geteoinst end end end end end else error(15) end; procedure irwith; (* *) (* read markers into mkon *) (* check that: 1) each marker refers to the currently specified piece (warning) 2) each marker is within the piece (warning message) 3) there are no overlapping markers since one can not make them recombine (fatal error message) then string those markers that do not violate 1 and 2 into mkon in the order of the piece. the with instruction may not need parenthesis, see the section on numbering default. see also respecification for use of parenthesis *) begin { writeln(output,'IN irwith'); } (* 2000 Oct 18: Set withused on during the first pass so that if there are statements without a 'with', they will trigger stepping the piece forward even if they are before the withs that come later in the inst. This keeps the marks properly aligned with the sequences. (Prior to this date the test was "(not withused) and (pass = 2)".) *) if not withused then begin rewrite(marksdelila); marksautomate(marksdelila); withused := true; end; readchangeset(mutations); { writechangeset(output, mutations); writeln(output,'after writechangeset'); } numofchanges := mutations.number; end; procedure irdirwit; (* *) (* sets directionwanted and mkon *) begin if (word in [dirdelila,witdelila]) then begin if (word=witdelila) then irwith else if (word=dirdelila) then begin irdirection; if correct then begin if not eoinst then findword; if not eoinst then begin if word<>witdelila then irword; if correct then if word=witdelila then irwith else error(7) end end end else error(7) end end; (* higher level procedures ************************************************) procedure outofrange(pie: pieceptr; var p: integer); (* p is out of range of the piece pie *) begin with pie^.key do begin case def.defout of rreduce: begin (* force it into the piece *) error(208); reduceposition(pie,p); ivalueposition(p); (* If the position is the same as before ... *) if p = previousfromposition then begin (* if this is the second reduce for the instruction *) if reduced then begin (* this is a reduction from one side and is likely be a mistake *) error(210); end else reduced := true; end; correct:=true; end; rcontinue: begin (* ignore it with a warning *) error(209); correct:=false; geteoinst (* skip rest of statement *) end; rhalt: begin (* die horribly *) error(203); correct:=false; geteoinst (* skip rest of statement *) end end end end; procedure known(var theposition: integer); (* is the position known? the convention is that a reference outside the coordinate system is unknown *) begin if not between(libpie^.key.coobeg, theposition, libpie^.key.cooend) then begin error(206); outofrange(libpie,theposition) end end; procedure okposition(var p: integer); (* is the position ok? (within the range of the piece). 2007 Dec 06: if not within AND circular piece always reduce *) begin if not within(libpie, p) then begin if libpie^.key.piecon = circular {qqq} then reduceposition(libpie,p) (* new as of 2007 Dec 06 *) else outofrange(libpie, p); end; end; procedure posit(var theposition: integer); (* read a *) { not implemented var (* variables for *) remkoff: markerptr; retk: trakey; regk: genkey; } procedure relative; (* read a number *) begin findword; if correct and (chr in ['+','-']) then begin zerobase := theposition; (* pick up the default coordinate base *) {writeln(output,'proc.relative: zerobase = ',zerobase:1);} irnumber; theposition:=theposition+inumber; if correct then ivalueposition(theposition) end end; procedure limitref(ref: reference); (* handle s of references *) begin with ref do begin if pass=2 then if ref.pienam.letters <> libpie^.key.hea.keynam.letters then error(205); if correct then begin irword; if correct then begin if word in [begdelila, enddelila] then begin if (pass=2) then begin case word of begdelila: theposition:=refbeg; enddelila: theposition:=refend end; known(theposition) end end else error(7) end end end end; procedure limitmkoff(mkoff: markerptr); (* get limits for marker *) begin if numbered(mkspec) then ivaluenumber(trunc(mkspec)); limitref(mkoff^.key.ref) end; procedure limittk(tk: trakey); (* get limits for transcript *) begin if numbered(tkspec) then ivaluenumber(trunc(tkspec)); limitref(tk.ref) end; procedure limitgk(gk: genkey); (* get limits for gene *) begin if numbered(gkspec) then ivaluenumber(trunc(gkspec)); limitref(gk.ref) end; procedure limitpk(pk: pieceptr); (* get limits for piece *) begin ivaluenumber(trunc(pkspec)); irword; if correct then begin if word in [begdelila,enddelila] then begin if (pass=2) then begin with pk^.key do begin case word of begdelila: theposition:=piebeg; enddelila: theposition:=pieend; end end end end else error(7) end end; procedure limitco(pk: pieceptr); (* get limits for coordinate system *) begin ivaluenumber(trunc(pkspec)); irword; if correct then begin if word in [begdelila,enddelila] then begin if (pass=2) then begin with pk^.key do begin case word of begdelila: theposition:=coobeg; enddelila: theposition:=cooend; end end end end else error(7) end end; begin (* *) findword; if (chr in numbers) or (chr in ['+','-']) then begin (* *) irnumber; (* does the ivalueposition *) if correct then begin theposition:=inumber; previousfromposition := theposition; relative end end else if chr = 's' then begin (* should be a 'same' *) irword; if correct then begin if word = samdelila then begin if sameusageisvalid then begin theposition := previousfromposition; relative end else begin error(18); geteoinst end end end end else begin (* *) if parentheses=1 then begin (* *) irword; { not implemented if correct then begin if word in [mardelila,tradelila,gendelila] then begin save:=word; irkeyname; if correct then begin case save of mardelila: begin remkoff := nil; getmarker(remkoff); startheader(remkoff^.key.hea); with remkoff^.key do begin phenotype := nil; next := nil end; if pass=2 then lrmarker(keyname,remkoff) else itemfound:=true; if itemfound then limitref(remkoff^.key.ref) else error(201); clearheader(remkoff^.key.hea); clearmarker(remkoff) end; tradelila: begin startheader(retk.hea); if pass=2 then lrtrakey(keyname,retk) else itemfound:=true; if itemfound then limitref(retk.ref) else error(201); clearheader(retk.hea) end; gendelila: begin startheader(regk.hea); if pass=2 then lrgenkey(keyname,regk) else itemfound:=true; if itemfound then limitref(regk.ref) else error(201); clearheader(regk.hea) end end end end else error(7) end; } if correct then begin ivalueposition(theposition); (* finish the last value *) relative end end else begin (* or coordinate *) irword; if correct then begin if word in [mardelila,tradelila,gendelila,piedelila, coodelila] then begin case word of (* *) mardelila: begin if pass=2 then specified(mkspec) else correct:=true; if correct then limitmkoff(mkoff) else error(202) end; tradelila: begin if pass=2 then specified(tkspec) else correct:=true; if correct then limittk(tk) else error(202) end; gendelila: begin if pass=2 then specified(gkspec) else correct:=true; if correct then limitgk(gk) else error(202) end; piedelila: begin if pass=2 then specified(pkspec) else correct:=true; if correct then limitpk(pk) else error(202) end; (* coordinates *) coodelila: begin if pass=2 then specified(pkspec) else correct:=true; if correct then limitco(pk) else error(202) end end end else error(7) end; if correct then begin ivalueposition(theposition); (* finish the last value *) relative end end end; if (pass = 2) and correct then okposition(theposition) end; procedure posits; (* read in the from, to. the frodelila has already been read *) begin if debugging then writeln(debug,'positions-in'); previousfromposition := 0; sameusageisvalid := false; posit(fromposition); if correct then begin irword; sameusageisvalid := true; if correct then if (word=todelila) then posit(toposition) else error(7) end ;if debugging then writeln(debug, 'positions-out',fromposition:5,toposition:5); end; procedure absdirection(var wanted, reference: direction); (* calculate the absolute direction wanted (plus or minus) relative to a reference direction *) begin case wanted of plus:; minus:; dircomplement: case reference of plus: wanted:=minus; minus: wanted:=plus end; dirhomologous: wanted:=reference end end; procedure allref(ref: reference); (* does the all from a reference *) begin with ref do begin (* decide absolute direction *) absdirection(dirwanted,refdir); (* now decide which base is first *) if dirwanted = refdir then begin fromposition:=refbeg; toposition:= refend end else begin fromposition:=refend; toposition:= refbeg end; ivalueposition(fromposition); known(fromposition); if correct then begin ivalueposition(toposition); known(toposition); if correct then begin (* there are no circular refered objects *) pk^.key.piecon:=linear; dopiece end end end end; procedure ircondition; (* read the condition of an if statement, return the value in the variable choice thedelila = then (true) should be done elsdelila = else (false) should be done *) begin (* dummy *) end; procedure numberstructure; (* allow the structure indicated by word to be numbered see the section on numbering default *) var indnum: numberedstructure; (* an index for def.num.str *) begin if word in structure then case word of orgdelila: def.num.str[orgnum]:=on; chrdelila: def.num.str[chrnum]:=on; mardelila: def.num.str[marnum]:=on; tradelila: def.num.str[tranum]:=on; gendelila: def.num.str[gennum]:=on; piedelila: def.num.str[pienum]:=on; recdelila: def.num.str[recnum]:=on; enzdelila: def.num.str[enznum]:=on; end else if word = alldelila (* turn all on *) then for indnum:=orgnum to enznum do def.num.str[indnum]:=on else begin error(7); geteoinst end end; (* the core of delila follows *) (**************************************************************************) (**************************************************************************) (**************************************************************************) begin (* begin of readinstruction *) (* read an instruction from the *) (* the conception of this part of the code is due to paul morrissett *) eoinst:=true; (* still at end of last instruction *) reduced:=false; (* piece ends have not been reduced yet *) {zzzppp} mutations.number := 0; (* force it to be clear of mutations *) findword; (* move to first word *) if not eof(inst) then begin eoinst:=false; (* now we are not at the end of an instruction.... *) irword; if correct then begin if (word=titdelila) then begin (* a *) if instructioncount=1 then begin irtitle; titleexists := true; if pass = 2 then clearline(title) (* save space *) end else begin if pass = 1 then error(3); geteoinst end end else if (word in structure) then begin (* a <specification> *) if tvrschecks(noder(word)) then begin irkeyname; if correct then begin if word in [mardelila,recdelila,enzdelila] then error(0); (* not implemented *) if (pass=2) then begin case word of orgdelila: begin lrorgkey(keyname,ok); if itemfound then begin spec(ok.hea,def.num.str[orgnum],okspec); tvrsbook(orgnode); bworgkey(book,ok); (* make lower nodes unspecified *) unspec(ckspec); unspecbelowck end else error(201) end; chrdelila: begin lrchrkey(keyname,ck); if itemfound then begin spec(ck.hea,def.num.str[chrnum],ckspec); tvrsbook(chrnode); bwchrkey(book,ck); unspecbelowck end else error(201) end; { not implemented mardelila: begin lrmarker(keyname,mkoff); if itemfound then begin spec(mkoff^.key.hea,def.num.str[marnum],mkspec); if def.key.mar=on then bwmarkers(mkoff) else clearmarker(mkoff); specpiece(mkoff^.key.ref) end else error(201) end; tradelila: begin lrtrakey(keyname,tk); if itemfound then begin spec(tk.hea,def.num.str[tranum],tkspec); if def.key.tra = on then bwtrakey(tk); specpiece(tk.ref) end else error(201) end; } gendelila: begin lrgenkey(keyname,gk); if itemfound then begin spec(gk.hea,def.num.str[gennum],gkspec); if def.key.gen = on then bwgen(book, gk); specpiece(gk.ref); end else error(201) end; piedelila: begin lrpiece(keyname,libpie); if itemfound then pkspec:=-0.5 else error(201) (* do not write piece here, but rather, at the request *) end; { not implemented recdelila: begin lrreckey(keyname,rk); if itemfound then begin spec(rk.hea,def.num.str[recnum],rkspec); bwreckey(rk); unspec(ekspec) end else error(201) end; enzdelila: begin lrenzyme(keyname,ek); if itemfound then begin spec(ek^.key.hea,def.num.str[enznum],ekspec); bwenzyme(ek) end else error(201) end; } end; (* case *) end;(* pass *) end (* correct keyname *) end (* traverse check *) else begin error(6); geteoinst end end {zzzwww} else if (word=getdelila) then begin (* a <request> *) (* give the piece a number *) if pass=2 then begin specified(pkspec); (* was the piece (fake) specified? *) if correct then begin { writeln(output,'spec(pk^.key.hea,def.num.str[pienum],pkspec) ****'); writeln(output,'def.num.str[pienum] = ',def.num.str[pienum]:1); writeln(output,'pkspec = ',pkspec:1); } (* put the number into the pk now so that it gets listed properly on the listing *) copyline(libpie^.key.hea.note, pk^.key.hea.note); spec(pk^.key.hea,def.num.str[pienum],pkspec); {libpkset} { writeln(output,'getdelila @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'); write(output,'libpie^.key.hea.keynam: "'); writename(output,libpie^.key.hea.keynam); writeln(output,'"'); write(output,'pk^.key.hea.keynam: "'); writename(output,pk^.key.hea.keynam); writeln(output,'"'); } (* we better remove any previous pieces or we will run out of memory pretty darn quick!! *) clearpiece(pk); (* since we are starting to get the piece, fill it out as much as we can *) copyname(libpie^.key.hea.keynam, pk^.key.hea.keynam); {zzzDDD} (* make sure the number is in the piece: *) addnumber(pk^.key.hea, def.num.item); pk^.key.mapbeg := libpie^.key.mapbeg; pk^.key.coocon := libpie^.key.coocon; pk^.key.coodir := libpie^.key.coodir; pk^.key.coobeg := libpie^.key.coobeg; pk^.key.cooend := libpie^.key.cooend; pk^.key.piecon := libpie^.key.piecon; pk^.key.piedir := libpie^.key.piedir; pk^.key.piebeg := libpie^.key.piebeg; pk^.key.pieend := libpie^.key.pieend; { writeln(output,'getdelila: pk^.key.coocon ', pk^.key.coocon); writeln(output,'getdelila: libpie^.key.coocon ', libpie^.key.coocon); writeln(output,'getdelila: pk^.key.piecon is ', pk^.key.piecon); writeln(output,'getdelila: libpie^.key.piecon is ', libpie^.key.piecon); } {zzzsss} end else begin error(202); geteoinst end end; if correct then irword; (* <range> *) if correct then begin (* set up defaults *) dirwanted:=dirhomologous; cutposition:=pk^.key.coobeg; (* default cut position *) while mkon<>nil do clearmarker(mkon); (* figure out the kind of range *) if (word=frodelila) then begin posits; if correct then begin if not eoinst then begin irword; if correct then irdirwit end; if correct and (pass=2) then begin pk^.key.piecon:=linear; (* positions forces linearity *) dopiece; geteoinst; end end else geteoinst end else if (word=alldelila) then begin (* get <all> *) irword; if correct then begin (* do cut, direction and with *) if (word in structure) then begin save:=word; if not eoinst then findword; if not eoinst then begin (* parse: cut, direction, with *) irword; if correct then begin if word in [cutdelila,dirdelila,witdelila] then begin if word = cutdelila then begin if save = piedelila then begin if pass = 2 then if libpie^.key.piecon <> circular then error(207); if correct then begin posit(cutposition); (* <cut> *) if correct and not eoinst then begin irword; if correct then irdirwit end end end else begin error(8); geteoinst end end else irdirwit end else error(7) end end; (* now do the all *) if correct then if save in [orgdelila,chrdelila,mardelila, recdelila,enzdelila] then error(0); if correct and (pass=2) then begin case save of orgdelila: ; chrdelila: ; mardelila: begin specified(mkspec); if correct then allref(mkoff^.key.ref) else error(202) end; tradelila: begin specified(tkspec); if correct then allref(tk.ref) else error(202) end; gendelila: begin specified(gkspec); if correct then allref(gk.ref) else error(202) end; piedelila: begin (* we checked the specification after the get *) with pk^.key do begin absdirection(dirwanted,libpie^.key.piedir); if (piecon=linear) then begin if dirwanted = libpie^.key.piedir then begin fromposition:=libpie^.key.piebeg; toposition:= libpie^.key.pieend; end else begin fromposition:=libpie^.key.pieend; toposition:= libpie^.key.piebeg; end end else begin (* piecon = circular *) fromposition:=cutposition; (* figure out toposition. there are 6 cases to consider for this base: reversing cut at (coordinate position): begin end other yes succ(cut) begin succ(cut) no end pred(cut) pred(cut) *) (* if reversing . . . *) if dirwanted<>libpie^.key.piedir then if cutposition = cooend then toposition:=coobeg else toposition:=succ(cutposition) else if cutposition = coobeg then toposition:=cooend else toposition:=pred(cutposition); end; ivalueposition(fromposition); ivalueposition(toposition); {zzzsss} { writeln(output,'allcall'); writeln(output,'pk^.key.piecon is ', pk^.key.piecon); writeln(output,'libpie^.key.piecon is ', libpie^.key.piecon); writeln(output,'PASS 2 at get from/to dopiece call <<<<<<<<<<<<<<<<<<<<<'); } dopiece; geteoinst; {zzzthe following is probably irrelevant now and so geteoinst can be removed. do so when life is calm... however, another place that calls dopiece has geteoinst, so ????} (* Previously (prior to 1999 March 21) geteoinst was done before dopiece. But to allow the changeset reading routines to violate the stupid ichread mechanism, I removed the geteoinst from irwith. This allows the listing to have the #s and mutations done neatly. *) end end; recdelila: ; enzdelila: ; end end end end end else if (word=evedelila) then begin (* get <every> *) irword; if correct then begin if (word in structure) then begin error(0) (* dummy. for the most part this could be a search for the type and then a copy library to book (except gen and tra, and non-printing due to def.key.xxx=off *) end else error(7); end end else error(7); end end else if (word=ifdelila) then begin (* an <if> statement *) (* this is a recursive call.. readinst is used. *) ircondition; if correct then begin irword; if correct then begin if word = thedelila then begin (* then *) error(0) (* dummy *) end else error(7) end end end else if (word=defdelila) or (word=setdelila) (* <SET> *) then begin (* a <default> *) irword; if correct then begin if (word=keydelila) then begin (* a <key default> *) irword; if correct then begin if (word in [notdelila,mardelila,gendelila,tradelila]) then begin save:=word; irword; if correct then begin if word in [ondelila,offdelila] then begin if (pass=2) then begin case save of notdelila:if (word=ondelila) then def.key.note:=on else def.key.note:=off; mardelila:if (word=ondelila) then def.key.mar:=on else def.key.mar:=off; gendelila:if (word=ondelila) then def.key.gen:=on else def.key.gen:=off; tradelila:if (word=ondelila) then def.key.tra:=on else def.key.tra:=off; end end end else error(7); end end else error(7); end end else if (word=sitdelila) then begin (* <recognition site default> *) irword; if correct then begin if (word in [expdelila,moddelila,cledelila]) then begin save:=word; irword; if correct then begin if word in [ondelila,offdelila] then begin if eoinst then begin if (pass=2) then case save of expdelila: if (word = ondelila) then def.sit.expand := on else def.sit.expand := off; moddelila: if (word = ondelila) then def.sit.modify := on else def.sit.modify := off; cledelila: if (word = ondelila) then def.sit.cleave := on else def.sit.cleave := off end end end else error(7); end end else error(7); end end else if (word=outdelila) then begin (* a <range default> *) irword; if correct then begin if (word in [reddelila,condelila,haldelila])then begin if (pass=2) then case word of reddelila: def.defout := rreduce; condelila: def.defout := rcontinue; haldelila: def.defout := rhalt end end else error(7); end end (* new code 1995 January 24 *) else if (word=coodelila) then begin (* a <coordinate default> *) {writeln(output,'parsing<coordinate default>');} findword; if (chr in numbers) or (chr in ['+','-']) then begin {writeln(output,'parsing<coordinate default>: presume it is a number');} (* new code 1995 Nov 13 *) irnumber; (* does the ivalueposition, but so what? *) if correct then begin {writeln(output,'the number read is ',inumber:1);} def.coo := coorzero; zeroshift := inumber end end {zzzlll} else begin irword; if correct then begin {writeln(output,'parsing<coordinate default>: word read ok');} if (word in [nordelila,zerdelila]) then begin {writeln(output,'parsing<coordinate default>: word', ' was nordelila or zerdelila');} if (pass=2) then {writeln(output,'parsing<coordinate default>: pass=2');} case word of (* coordinatetype = (coornormal,coorzero); *) nordelila: def.coo := coornormal; zerdelila: def.coo := coorzero; end end end end end (* new code 1999 March 20 *) else if (word=doudelila) then begin (* a <doubling default> *) {writeln(output,'parsing<doubling default>');} irword; if correct then begin if word in [ondelila,offdelila] then begin if word=ondelila then def.doubling := on else def.doubling := off end else error(7); end end {zzzlll} (* new code 1999 March 20 *) else if (word=arrdelila) then begin (* a <arrow default> *) {writeln(output,'parsing<arrowlength default>');} findword; if (chr in numbers) or (chr in ['+','-']) then begin irnumber; if correct then begin { writeln(output,'the number read is ',rnumber:20:15); } (* force arrows to be downwards. This also prevents a postscript bomb if the arrow length were 0 (when changeupperbits=changelowerbits) *) if rnumber < 0.1 then rnumber := 0.1; def.arrowlength := rnumber; changeupperbits := def.arrowlength+changelowerbits; end end end else if (word=numdelila) then begin (* <numbering default> *) findword; (* get first character in chr *) if not eoinst then begin if (chr in numbers) or (chr in ['+','-']) then begin (* <signed number> *) irnumber; (* the extra findword will not hurt *) if correct then if pass = 2 then def.num.item:=inumber-1 end else begin (* state or structure *) irword; if correct then begin if word in [ondelila,offdelila] (* <state> *) then case word of ondelila: def.num.sta:=on; offdelila: def.num.sta:=off; end else begin (* <structure set> *) (* clear all numbering *) {iii} for indnum:=orgnum to enznum do def.num.str[indnum]:=off; numberstructure; (* the first one *) while not eoinst do begin (* all the others *) findword; if not eoinst then begin irword; if correct then numberstructure end end end end end end else error(15) end else error(7) end end else if (word = notdelila) then begin (* <note insertion> *) irquote(usernotes); if pass=1 then (* remove these extra notes *) while usernotes<>nil do clearline(usernotes) end else if (word = namdelila) then begin (* <name> *) irlongname; if pass=1 (* remove these extra notes *) then emptyline(longname); end else error(7); end else geteoinst; (* first instruction word not correct *) end; instructioncount:=succ(instructioncount); if not eoinst then (* catch all cases of extra stuff *) begin findword; (* give um one chance *) if not eoinst then begin error(1); geteoinst end end; if parentheses<>0 then begin error(11); parentheses:=0 (* prevents snowballing *) end; if chr=';' then ichread; (* get past the ; *) end; (* end of readinstruction *) begin (* librarian *) initializedelila; (* pass 1: check syntax and grammer *) pass:=1; writeln(output,'Pass 1'); initreadinst; while (not eof(inst)) do readinstruction; if not titleexists then error(17); if not tvrschecks(libnode) then begin error(6); writeerrors; end; writeln(listing); writepasserrors(listing, pass); writepasserrors(output, pass); (* if there were no errors in the first pass, then do pass 2 *) if not pass1errors then begin writeln(output,'Pass 2'); (* pass 2: read library, create book *) pass:=2; reset(inst); (* go back to start of instruction file *) bwbookheader(title); (* the first pass picked up the title *) page(listing); initreadinst; while ((not eof(inst)) and (not pass2errors)) do readinstruction; writeln(listing); (* finish last line *) writelineinformation; if pass2errors then begin writeerrors; (* for that line *) bookhalt end else begin (* close up shop *) tvrslibrary(libnode); tvrsbook(libnode) end; writeln(listing); writepasserrors(listing, pass); (* picks up warnings also *) writepasserrors(output, pass); (* picks up warnings also *) end else bookhalt; if pass2errors or pass1errors then writeln(output,'Error(s) in pass ', pass:1, ': see the listing file for details'); if warnings then writeln(output, 'There are warnings: see the listing'); write(output,getcount:1,' piece'); if getcount <> 1 then write(output,'s'); writeln(output,' extracted'); write(output,booksize:1,' base'); (* compared to maxbook *) if booksize <> 1 then write(output,'s'); writeln(output,' extracted'); end; (* librarian *) (* main library calls *****************************************************) begin (* delila *) debugging := false; (* for debugging purposes *) initlibrarian; librarian; 1: end. (* delila *)