program embed(inst, book, mkvseqs, ranbook, embedp, embedbk, output); (* embed: embed an aligned set of DNA sequences into random sequences Elaine Bucheimer NCI-Frederick Bldg 469. Room 144 P.O. Box B Frederick, MD 21702-1201 (301) 846-5581 (-5532 for messages) bucheime@ncifcrf.gov http://www.ccrnp.ncifcrf.gov/~bucheime/ National Cancer Institute Laboratory of Experimental and Computational Biology Now maintained by: Dr. Thomas D. Schneider National Institutes of Health National Cancer Institute Center for Cancer Research Nanobiology Program Molecular Information Theory Group Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu http://www.ccrnp.ncifcrf.gov/~toms/ *) label 1; (* end of program *) const (* begin module version *) version = 1.32; (* of embed.p 2011 Oct 14 2011 Oct 14, 1.32: increase dnamax to 2048; clean up 2011 Oct 14, 1.31: crash bug identified 2011 Oct 03, 1.30: report warnings; crash bug 2011 Oct 03, 1.29: save version 2007 Nov 15, 1.28: improve documentation 2005 Sep 27, 1.27: coo coordinates are NOT SET BY PIECE! 2005 Sep 27, 1.26: crash fixed 2000 Nov 30, 1.25: bug record (TDS) 1998 Apr 1, 1.24: coordinates are set by new piece (TDS) 1997 Aug 18: last changes by Elaine origin 1997 July 1 *) updateversion = 1.02; (* defines lowest acceptable current parameter file *) (* end module version *) (* begin module describe.embed *) (* name embed: embed an aligned set of DNA sequences into random sequences synopsis embed(inst: in, book: in, mkvseqs, in: ranbook, in: embedp: in, embedbk: out, output: out) files inst: delila instructions of the form 'get from 56 -5 to 56 +10;' book: the book generated by delila using inst mkvseqs: random sequence output from the markov program ranbook: book made from random sequences using makebk program; either mkvseqs or ranbook must be contain sequence. If both contain sequence, then mkvseqs will be used as the source for random sequences. embedp: parameters to control the program. The file must contain the following parameters, one per line: parameterversion: The version number of the program. This allows the user to be warned if an old parameter file is used. alignmenttype: The type of alignment to use. f: first base, i: inst, b: book alignment 'b' is to be used when 'default coordinate zero;' is used in the inst file, resulting in a book whose coordinates do not match the inst coordinates. 'i' is to be used when the book contains a normal coordinate system corresponding to the inst file. 'f' simply aligns by the first base in the book. See alist.p for more details on alignmenttype. InFrom, InTo: the from-to range of the input sequences to be used. OutFrom, OutTo: the from-to range of the sequences to output. This includes the Infrom range AND the random sequences. embedbk: book created by the program. Contains the sequences embedded within random sequences to the specified range. output: messages to the user description Embed embeds a given set of aligned sequences into random sequences having a specified range. If there is an incomplete sequence in the region to be embedded, it is filled in with random sequences as well. This allows one to destroy a pattern in the aligned sequences, so that the sequences can be realigned to find other patterns nearby. The parameters OutFrom, InFrom, InTo, OutTo in embedp set the range to do the embedding. In order for the program to function correctly, the following must be true: OutFrom <= InFrom <= InTo <= InFrom. The sequence from InFrom to InTo is not changed, and random sequence is filled in around it from OutFrom on the left to OutTo on the right. See example below. If the orginal sequence is longer than the range OutFrom to OutTo then the book will contain the embedded sequence with orginal sequence on either side of the random sequence. The program stores the random sequence as a string and then uses it base by base until there is no more in the string. Then it reads another string of random sequence. In this way, none of the random sequence is "thrown away". If the program finds the end of mkvseqs or ranbook before it has embedded all the sequences, it gives a message that it is out of random sequence and halts. Why doesn't the program reuse the random sequence? This is not a good idea because the embedded sequences are designed to be fed into malign, and malign would pick up on this reused sequence and find unnatural sequence conservation. Aligned sequences can be viewed with the alist program. The random sequences are generated by the markov program. They can be read from either mkvseqs or ranbook. mkvseqs is directly generated from markov to a given composition and length. Ranbook can be made using the makebk program. If both files are present, mkvseqs is used. The output of this program is designed to be fed into the malign program for multiple alignment. examples With the following parameters from embedp the sequence would be embedded as shown below. -10 10 InFrom, InTo: range of input sequences to be used -30 30 OutFrom, OutTo: range of the sequences to output original: -----|-------------------<---------0--------->-------------------|----- -30 -10 +10 +30 OutFrom InFrom InTo OutTo embedded: ********************<---------0--------->******************** -30 random -10 original +10 random +30 sequence sequence sequence Note that if there is any sequence in the original alignment outside the range OutFrom to OutTo, it will be copied to the embedbk. documentation see also alist.p, markov.p, makebk.p, malign.p author Elaine Bucheimer bugs The program cannot handle sequences longer than dnamax. This is a fixable bug. A possible future addition to the program would be to allow the user to specify if they want the old sequence hanging around or if the sequence should be chopped outside the OutFrom and OutTo coordinates. It appears that the 'i' option does not embed correctly. The resulting book does not have the advertised coordinates. A temporary solution is to use the f option with appropriate ranges. technical notes *) (* end module describe.embed *) (* begin module book.const *) (* constants needed for book manipulations *) dnamax = 2048; (* length of dna arrays *) namelength = 100; (* maximum key name length *) linelength = 200; (* maximum line readable in book *) (* end module book.const version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module embed.const *) maxstring = 80000; (* the maximum string *) verbose = false; (* turn on to get blow by blow details *) (* end module embed.const version = 5.70; (@ of alist.p 1997 April 22 *) (* begin module filler.const *) fillermax = 50; (* the size of the filler array for a string *) (* end module filler.const version = 5.70; (@ of alist.p 1997 April 22 *) type (* begin module interact.type *) (* begin module string.type *) stringptr = ^string; (* pointer to a string *) string = record (* a string of characters *) letters: array[1..maxstring] of char; (* the letters in the string *) length: integer; (* the number of characters in the string *) current: integer; (* the letter we are working on *) next: stringptr; (* the next string in a series *) end; (* end module string.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *) (* end module interact.type version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module filler.type *) (* the following is an array used to fill a string. it is convenient to have it much shorter than the maxstring, so that it is easy to fill the string using procedure fillstring. the user must declare the value of constant fillermax. *) filler = packed array[1..fillermax] of char; (* end module filler.type version = 5.70; (@ of alist.p 1997 April 22 *) (* begin module trigger.type *) trigger = record (* an object to be searched for *) seek: string; (* the characters looked for *) state: integer; (* how close to triggering we are *) skip: boolean; (* trigger not found- skip the line *) found: boolean (* the trigger was found *) end; (* end module trigger.type version = 5.70; (@ of alist.p 1997 April 22 *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) var (* files used by this program *) book, (* book created by delila using the inst file *) embedbk, (* book created with embedded sequences *) embedp, (* parameter file for this program *) inst, (* delila instructions *) mkvseqs, (* random sequences generated by markov *) ranbook: text; (* book of random sequences *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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 = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module package.trigger *) (* ************************************************************************ *) (* begin module interact.clearstring *) procedure clearstring(var ribbon: string); (* empty the string *) var index: integer; (* to the ribbon *) begin (* clearstring *) with ribbon do begin for index := 1 to maxstring do letters[index] := ' '; length := 0; current := 0; end end; (* clearstring *) (* end module interact.clearstring version = 4.18; (@ of prgmod.p 1996 September 12 *) (* begin module interact.writestring *) procedure writestring(var tofile: text; var s: string); (* write the string s to file tofile, no writeln *) var i: integer; (* index to s *) begin (* writestring *) with s do for i := 1 to length do write(tofile, letters[i]) end; (* writestring *) (* end module interact.writestring version = 4.18; (@ of prgmod.p 1996 September 12 *) (* begin module filler.fillstring *) procedure fillstring(var s: string; a: filler); (* this procedure makes it reasonably easy to fill the string s with characters. one calls the procedure as: *) (* 1 2 3 4 5 *) (* 12345678901234567890123456789012345678901234567890 *) (* fillstring(s, 'this-is-the-string '); the two comments make it easy to line the characters up. also, for this example, it was assumed that the length of filler as defined by the constant fillermax was 50. *) var length: integer; (* of the string without trailing blanks *) index: integer; (* of s *) begin (* fillstring *) clearstring(s); length := fillermax; while (length > 1) and (a[length] = ' ') do length := pred(length); if (length = 1) and (a[length] = ' ') then begin writeln(output, 'fillstring: the string is empty'); halt end; for index := 1 to length do s.letters[index] := a[index]; s.length := length; s.current := 1 end; (* fillstring *) (* end module filler.fillstring version = 4.18; (@ of prgmod.p 1996 September 12 *) (* begin module filler.filltrigger *) procedure filltrigger(var t: trigger; a: filler); (* fill the trigger t *) begin (* filltrigger *) fillstring(t.seek,a) end; (* fillstring *) (* end module filler.filltrigger version = 4.18; (@ of prgmod.p 1996 September 12 *) (* begin module trigger.proc *) (* this module allows one to scan a series of characters, as from an array or a file, and to "trigger" or detect a simple string in the series. the advantage of the trigger is that several triggers can "observe" a stream of characters at once, each looking for a different thing. some other modules required: interact.const, interact.type *) procedure resettrigger(var t: trigger); (* reset the trigger to ground state *) begin (* resettrigger *) with t do begin state := 0; skip := false; found := false end end; (* resettrigger *) procedure testfortrigger(ch: char; var t: trigger); (* look at the character ch. if it is part of the trigger (at the current trigger state), then the trigger state goes higher. if it is not part of the trigger then the trigger state is reset, skip is true and one should skip onward to find the trigger. if the trigger is found, found is true. 1996 Sep 12: Bug found! In the case of a trigger "ab", the program used to miss it for situations like "aab". This was because at the first a it would step up. Then it would see the second a and recognize that was not part of ab. It would fail to realize that it could be the start of a new one. The code now accounts for that possibility. *) begin (* testfortrigger *) with t do begin state := succ(state); { writestring(list,seek); writeln(list,'testfortrigger seek.letters[',state:1,']:', seek.letters[state],' ch:',ch); } if seek.letters[state] = ch then begin skip := false; if state = seek.length then found := true else found := false end else begin (* it failed. But wait! It could be the beginning of a NEW trigger string! *) if seek.letters[1] = ch then begin state := 1; skip := false; found := false end else begin (* reset trigger *) state := 0; skip := true; found := false end end end end; (* testfortrigger *) (* end module trigger.proc version = 4.18; (@ of prgmod.p 1996 September 12 *) (* begin module skipblanks *) procedure skipblanks(var thefile: text); (* skip over blanks until a non-blank, or end of line, is found *) begin while (thefile^ = ' ') and not eoln(thefile) do get(thefile); end; procedure skipnonblanks(var thefile: text); (* skip over nonblanks until a blank, or end of line, is found *) begin while (thefile^ <> ' ') and not eoln(thefile) do get(thefile); end; procedure skipcolumn(var thefile: text); (* skip over a data column *) begin skipblanks(thefile); skipnonblanks(thefile) end; (* end module skipblanks version = 4.18; (@ of prgmod.p 1996 September 12 *) (* ************************************************************************ *) (* end module package.trigger version = 5.70; (@ of alist.p 1997 April 22 *) (* begin module package.align *) (* ************************************************************************ *) (* 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); (* clear the dna strutures to the free list *) 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 *) (* Note: the dirhomologous and dircomplement are treated as plus and minus directions, which MIGHT NOT BE RIGHT! *) var i: integer; (* an intermediate value *) begin with pie^.key do begin case piedir of dirhomologous, plus: if p>=piebeg then i:=p-piebeg+1 else i:=(p-coobeg)+(cooend-piebeg)+2; dircomplement, 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 *) (* Note: the dirhomologous and dircomplement are treated as plus and minus directions, which MIGHT NOT BE RIGHT! *) var p: integer; (* an intermediate value *) begin with pie^.key do begin case piedir of dirhomologous, plus: begin p:=piebeg+(i-1); if p>cooend then if coocon=circular then p:=p-(cooend-coobeg+1) end; dircomplement, 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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 ',linelength:1,' characters'); writeln(output,'* Only ',linelength:1,' characters read from book'); writeln(output,'***********************************************'); end; l^.length:=i; l^.next:=nil; readln(thefile); theline := succ(theline) end; (* end module book.brline version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module book.brdirect *) procedure brdirect(var thefile: text; var theline: integer; var direct: direction); (* read a direction *) var ch: char; begin skipstar(thefile); readln(thefile,ch); theline := succ(theline); if ch='+' then direct:=plus else direct:=minus end; (* end module book.brdirect version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module book.brconfig *) procedure brconfig(var thefile: text; var theline: integer; var config: configuration); (* read a configuration *) var ch: char; begin skipstar(thefile); readln(thefile,ch); theline := succ(theline); if ch='l' then config:=linear else config:=circular end; (* end module book.brconfig version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module book.brnotenumber *) procedure brnotenumber(var thefile: text; var theline: integer; var note: lineptr); (* book note reading to obtain the number of the object. the procedure returns the value of the number as a global. (this is not such a good practice, but we are stuck with it for now.) *) begin (* brnotenumber *) note:=nil; numbered := false; number := 0; (* force number to zero if there is no number at all *) (* the next character is n or * depending on whether there are notes *) if thefile^ = 'n' then begin readln(thefile); theline := succ(theline); if thefile^ <> '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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* ************************************************************************ *) (* end module package.brpiece version = 7.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* ************************************************************************ *) (* end module package.getpiece version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module findblank *) procedure findblank(var afile: text); (* read a file to find the next blank character *) var ch: char; begin repeat read(afile,ch) until ch = ' ' end; (* end module findblank version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module findnonblank *) procedure findnonblank(var afile: text; var ch: char); (* find the next non blank character in a file, return it in ch. *) begin ch:=' '; while (not eof(afile)) and (ch = ' ') do begin read(afile,ch); if eoln(afile) then readln(afile) end end; (* end module findnonblank version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module align.align *) procedure align(var inst, book: text; var theline: integer; var pie: pieceptr; var length, alignedbase: integer); (* documentation on align is in module info.align and delman.use.aligned.books. 1996 Sep 12: The routine now uses the trigger functions found in prgmod. The bug in the oldalign routine (that it misses the end of comments that end in a series of asterisks) has been fixed. It now checks that the piece corresponds to the book. *) const maximumrange = 10000; (* if the alignment point is more than this distance from the piece ends, the program halts in an attempt to catch the alignment bug... 1991 Jan 11 It appears that the rewrite of the code has removed the bug, but the check will be kept. *) semicolon = ';'; (* end of delila instruction *) var ch: char; (* a character in inst *) p: integer; (* index to a piece name *) p1: integer; (* another index to a piece name *) done: boolean; (* done finding an aligning get *) thebase: integer; (* the base read in *) indefault: boolean; (* true when within a default statement. These can contain the word 'piece', which must be ignored. *) gettrigger: trigger; (* trigger to find 'get' *) defaulttrigger: trigger; (* trigger to find 'default' *) nametrigger: trigger; (* trigger to find 'name' *) piecetrigger: trigger; (* trigger to find 'piece' *) settrigger: trigger; (* trigger to find 'set' *) begincomment: trigger; (* trigger to find '(-*' (ignore the dash!) *) endcomment: trigger; (* trigger to find '*-)' (ignore the dash!) *) begincurly: trigger; (* trigger to find comments: '{' *) endcurly: trigger; (* trigger to find comments: '}' *) quote1trigger: trigger; (* trigger to find single quote ' *) quote2trigger: trigger; (* trigger to find double quote " *) { procedure rd(var f: text; var ch: char); (* read ch from f allowing inspection of the result *) begin read(f,ch); write(output,ch); write(list,ch); write(output,'<',ch,'>'); end; procedure rdln(var f: text); (* readln f allowing inspection of the result *) begin readln(f); writeln(output); writeln(list); end; } procedure skipcomment(var f: text); (* skip an entire comment *) var comment: boolean; (* true means we are inside a comment *) begin (* skip to end of comment *) resettrigger(endcomment); comment := true; while comment do begin if eof(f) then begin writeln(output,'A comment does not end!'); halt end; if eoln(f) then readln(f) { rdln(f) } else begin {write(output,'<'); rd(f,ch); write(output,'>');} read(f,ch); testfortrigger(ch, endcomment); if endcomment.found then comment := false; end end end; procedure skipcurly(var f: text); (* skip an entire comment made by {}*) var comment: boolean; (* true means we are inside a comment *) begin (* skip to end of comment *) resettrigger(endcurly); comment := true; while comment do begin if eof(f) then begin writeln(output,'A comment does not end!'); halt end; if eoln(f) then readln(f) { rdln(f) } else begin {write(output,'<'); rd(f,ch); write(output,'>');} read(f,ch); testfortrigger(ch, endcurly); if endcurly.found then comment := false; end end end; procedure skipquote(quote: trigger); (* skip an entire quote of either the ' or " persuasion *) var kind: char; (* the kind of quote, ' or " *) begin kind := quote.seek.letters[1]; {writeln(output,'skipquote ',kind);} repeat findnonblank(inst,ch); (* get to the quote *) until (ch = kind) or eof(inst); if ch <> kind then begin writeln(output,'end of quote starting with ',kind,' not found'); halt; end; end; begin filltrigger(defaulttrigger,'default'); filltrigger(gettrigger,'get '); filltrigger(nametrigger,'name '); filltrigger(piecetrigger,'piece '); filltrigger(settrigger,'set '); filltrigger(begincomment,'(* '); filltrigger(endcomment,'*) '); filltrigger(begincurly,'{ '); filltrigger(endcurly,'} '); filltrigger(quote1trigger,''' '); filltrigger(quote2trigger,'" '); resettrigger(defaulttrigger); resettrigger(gettrigger); resettrigger(nametrigger); resettrigger(piecetrigger); resettrigger(settrigger); resettrigger(begincomment); resettrigger(begincurly); resettrigger(quote1trigger); resettrigger(quote2trigger); indefault := false; if not eof(book) then begin (* if there is still more to the book ... *) getpiece(book,theline,pie); (* read in the piece *) if not eof(book) then begin (* if we found a piece ... *) length:=pietoint(pie^.key.pieend,pie); (* calculate piece length *) (* now find in inst the next occurance of 'get' *) done := false; while not done do begin if eof(inst) then begin (* no instructions? *) alignedbase := 1; (* simply align by the first base *) done := true end else begin if eoln(inst) then readln(inst) {then rdln(inst)} else begin {rd(inst,ch);} read(inst,ch); testfortrigger(ch, begincomment); testfortrigger(ch, begincurly); if begincomment.found or begincurly.found then begin if ch = '*' then begin skipcomment(inst); resettrigger(begincomment); end else begin resettrigger(begincurly); skipcurly(inst); end end else begin (* we are not inside a comment *) testfortrigger(ch, gettrigger); if gettrigger.found then begin findnonblank(inst,ch); (* get to "from" *) findblank(inst); (* get past "from" *) read(inst,thebase); (* read in the alignedbase *) {writeln(output);writeln(output,'thebase = ',thebase:1);} alignedbase:=pietoint(thebase,pie); {writeln(output,'alignedbase=',alignedbase:1);} done := true end; testfortrigger(ch, quote1trigger); if quote1trigger.found then begin skipquote(quote1trigger); end; testfortrigger(ch, quote2trigger); if quote2trigger.found then begin skipquote(quote2trigger); end; testfortrigger(ch, defaulttrigger); if defaulttrigger.found then begin indefault := true; resettrigger(defaulttrigger) end; if ch = semicolon then indefault := false; testfortrigger(ch, settrigger); if settrigger.found then begin indefault := true; resettrigger(settrigger) end; if ch = semicolon then indefault := false; (* check that piece names are correct *) testfortrigger(ch, piecetrigger); if not indefault then if piecetrigger.found then begin skipblanks(inst); (* get to name *) with pie^.key.hea.keynam do begin for p := 1 to length do begin read(inst,ch); if letters[p] <> ch then begin writeln(output,'The piece name in the book: '); writeln(output,letters:length); writeln(output,'does not match', ' the inst file piece name:'); (* write the letters that matched: *) for p1 := 1 to p-1 do write(output,letters[p1]); (* write the offending letter: *) write(output, ch); (* get the rest of the name and show it: *) done := eoln(inst); while not done do begin done := eoln(inst); if not done then begin read(inst,ch); if (ch = ' ') or (ch = ';') then done := true; if not done then write(output,ch); end; end; writeln(output); (* mark the first letter that does not match: *) for p1 := 1 to p-1 do write(output,' '); write(output,'^'); writeln(output); halt end; end end; end; end end end end; if (alignedbase <= -maximumrange) or (alignedbase > length + maximumrange) then begin writeln(output,' In procedure align:'); writeln(output,' read in base was ',thebase:1); writeln(output,' in internal coordinates: ',alignedbase:1); writeln(output,' maximum range was ',maximumrange:1); writeln(output,' piece length was ',length:1); with pie^.key.hea.keynam do writeln(output,' piece name: ',letters:length); writeln(output,' piece number: ',number:1); writeln(output,' aligned base is too far away... see the code'); halt end end end end; (* end module align.align version = 7.70; {of delmod.p 2005 Sep 15} *) (* ************************************************************************ *) (* end module package.align version = 7.70; {of delmod.p 2005 Sep 15} *) (* 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 dirhomologous, (* handle case, may not be right *) plus: writeln(book,'+'); dircomplement, (* handle case, may not be right *) 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.70; {of delmod.p 2005 Sep 15} *) (* 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(thefile,org); orgopen := true; end; (* end module book.bworg version = 7.70; {of delmod.p 2005 Sep 15} *) (* 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(thefile,chr); chropen := true; end; (* end module book.bwchr version = 7.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (* 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.70; {of delmod.p 2005 Sep 15} *) (****************************************************************************) (* end module package.bwrite version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module copyaline *) procedure copyaline(var fin, fout: text); (* copy a line from file fin to file fout *) begin (* copyaline *) while not eoln(fin) do begin fout^ := fin^; put(fout); get(fin) end; readln(fin); writeln(fout); end; (* copyaline *) (* end module copyaline version = 7.70; {of delmod.p 2005 Sep 15} *) (* begin module embed.readparameters *) procedure readparameters(var alignmenttype: char; var InFrom, InTo, OutFrom, OutTo: integer); (* reads the parameters from embedp starting with alignmenttype *) begin readln(embedp, alignmenttype); readln(embedp, InFrom, InTo); readln(embedp, OutFrom, OutTo); if (OutFrom > InFrom) or (InFrom > InTo) or (InTo > OutTo) then begin writeln(output, 'Error in ordering of parameters in embedp'); writeln(output, 'check that OutFrom <= InFrom <= InTo <= OutTo'); halt; end; write(output, 'embedding sequence range '); write(output, InFrom:1, ' to ', InTo:1); writeln(output, ' into ', OutFrom:1, ' to ', OutTo:1); writeln(output); end; (* end module embed.readparameters *) (* begin module embed.readrandom *) procedure readrandom (var rand: string; Useranbook: boolean; OutTo, OutFrom :integer; var ranbook: text; var theline: integer); (* reads the random sequences from mkvseqs ands stores them as a string called rand, with length OutTo-OutFrom+1 *) var index: integer; (* index to place in array rand.letters *) ch:char; (* character in mkvseqs *) begin clearstring(rand); index:= 0; if Useranbook=true then begin ch:=getto(ranbook,theline, ['d']); if not eof(ranbook) then begin read(ranbook,ch); (* skipstar *) while ch = '*' do begin read(ranbook,ch); (* skip blank *) repeat read(ranbook,ch); if (ch in ['a','c','g','t']) and (index < maxstring) then begin index:=succ(index); rand.letters[index]:=ch; end; until eoln(ranbook); readln(ranbook); (* go to next line *) theline := succ(theline); read(ranbook,ch); (* ch is either '*' or 'd' *) end; readln(ranbook); theline := succ(theline); end; readln(ranbook); end; if Useranbook = false then begin while (index apiece^.dna^.length - alignedbase then begin if apiece^.dna^.length - alignedbase < InFrom then begin adjustmentb := InTo - InFrom +1; write(output, 'warning: no sequence between '); writeln(output, InFrom:1, ' and ', InTo:1, ' for this piece'); end else begin adjustmentb := alignedbase + InTo - apiece^.dna^.length {-1}; if verbose then begin write(output, 'warning: incomplete sequence between '); writeln(output, InFrom:1, ' and ', InTo:1, ' for this piece'); writeln(output, 'adjusting on right end of sequence') end; end; end else adjustmentb := 0; (* determines how much the sequence needs to be moved to the right to make room for random sequence on left and adjusts alignedbase and length to match this *) if (alignedbase + OutFrom) < 1 then begin differencefrom:= abs(alignedbase + OutFrom) + 1 ; if apiece^.dna^.length + differencefrom > dnamax then begin writeln(output,'differencefrom = ',differencefrom:1); writeln(output,'apiece^.dna^.length = ', apiece^.dna^.length:1); writeln(output,'Procedure fillrandom: dnamax (= ',dnamax:1,')', ' would be exceeded: piece too big!'); halt; end; apiece^.dna^.length:=apiece^.dna^.length + differencefrom; if apiece^.dna^.length + differencefrom > dnamax then begin {2011 Oct 14 Apparently Elaine was adding differencefrom TWICE for some reason! I'll leave the code alone for now ...} writeln(output,'dnamax = ',dnamax:1); writeln(output,'differencefrom = ',differencefrom:1); writeln(output,'apiece^.dna^.length = ', apiece^.dna^.length:1); writeln(output,'differencefrom = ',differencefrom:1); writeln(output,'apiece^.dna^.length + 2* differencefrom = ', apiece^.dna^.length+2*differencefrom :1); writeln(output,'That exceeds dnamax; increase dnamax.'); {zq} writeln(output,'Procedure fillrandom:'); halt; end; for index := apiece^.dna^.length downto 1 do begin apiece^.dna^.part[differencefrom + index] := apiece^.dna^.part[index]; end; if apiece^.key.piedir = plus then apiece^.key.piebeg := apiece^.key.piebeg - differencefrom else apiece^.key.piebeg := apiece^.key.piebeg + differencefrom; alignedbase := alignedbase + differencefrom; end; (* determines how much extra space needs to be added to the original sequence to make room for random sequence on the right and adjusts length to match this *) if (alignedbase + OutTo) > apiece^.dna^.length then begin differenceto :=(alignedbase + OutTo) - apiece^.dna^.length; apiece^.dna^.length:=apiece^.dna^.length + differenceto; if apiece^.key.piedir = plus then apiece^.key.pieend := apiece^.key.pieend + differenceto else apiece^.key.pieend := apiece^.key.pieend - differenceto; (* make coordinate system match piece: (TDS 1998 April 1) *) (* WOAH!! 2005 Sep 27, 1.27: coo coordinates are NOT SET BY PIECE! *) (* with apiece^.key do begin {zzz} {uuu} NO NO NO NO NO! The coo must be whatever it was! If a piece is complementary, this gives a weird negative coordinate system ... coocon := piecon; coodir := piedir; coobeg := piebeg; cooend := pieend; end; *) end; (* fill left side of sequence with random sequence *) for index := (alignedbase + OutFrom) to (alignedbase + InFrom -1 + adjustmenta) do begin rand.current:=succ(rand.current); if rand.current > rand.length then begin if Useranbook then begin (* using ranbook *) if ranbook^<>'p' then outofsequence else begin readrandom(rand, Useranbook, OutTo, OutFrom, ranbook, thelineran); rand.current:=succ(rand.current); end; end else begin (* not using ranbook *) if eof(mkvseqs) then outofsequence else begin readrandom(rand, Useranbook, OutTo, OutFrom, ranbook, thelineran); rand.current:=succ(rand.current); end; end; end; apiece^.dna^.part[index]:= chartobase(rand.letters[rand.current]); end; (* fill right side of sequence with random sequence *) for index := (alignedbase + InTo +1 - adjustmentb) to (alignedbase + OutTo) do begin rand.current:=succ(rand.current); if rand.current > rand.length then begin if Useranbook then begin (* using ranbook *) if ranbook^<>'p' then outofsequence else begin readrandom(rand, Useranbook, OutTo, OutFrom, ranbook, thelineran); rand.current:=succ(rand.current); end; end else begin (* not using ranbook *) if eof(mkvseqs) then outofsequence else begin readrandom(rand, Useranbook, OutTo, OutFrom, ranbook, thelineran); rand.current:=succ(rand.current); end; end; end; apiece^.dna^.part[index]:= chartobase(rand.letters[rand.current]); end; end; (* end module embed.fillrandom *) (* begin module embed.copybook *) procedure copybook; (* Copy all lines outside of pieces such as organism and chromosome lines. Literally: copy up to but not including the next piece. *) var ch: char; (* character in book *) index, (* index to position along line in book *) s: integer; (* index to top string *) top: string; (* a line of the book not beginning with "p" *) done: boolean; (* done reading the book? *) begin done := false; while not done do begin if eof(book) then done := true; if not done then begin if book^ = 'p' then begin done := true; end else copyaline(book,embedbk); end; end; end; (* end module embed.copybook *) (* begin module embed.themain *) procedure themain(var inst, book, embedp, mkvseqs, ranbook, embedbk: text); (* the main procedure of the program *) var parameterversion: real; (* parameter version number *) s: integer; (* index to name string for piece *) (* variables used by the align routines: *) alignmenttype: char; (* 'f' means alignment by first internal coordinate base, 'b' means alignment by Book, 'i' means alignment by Instructions *) apiece: pieceptr; (* original piece designations, includes original sequence *) length, alignedbase: integer; piececount: integer; (* count of pieces *) (* variables used for embedding routines *) OutFrom, InFrom, InTo, OutTo: integer; (* parameters controlling range to do the embedding *) rand: string; (* string of random sequence *) Useranbook: boolean; (* true if ranbook is being used for the random sequences, false if mkvseqs is used *) theline: integer; (* the current line of the book *) thelineran: integer; (* the current line of the random book, ranbook *) begin writeln(output,'embed ',version:4:2); reset(embedp); readln(embedp, parameterversion); if parameterversion < updateversion then begin writeln(output, 'You have an old parameter file!'); halt; end; reset(mkvseqs); reset(ranbook); if not eof(mkvseqs) then begin Useranbook := false; writeln(output, 'using mkvseqs for random sequence'); end; if eof(mkvseqs) and not eof(ranbook) then begin Useranbook := true; writeln(output, 'using ranbook for random sequence'); end; if eof(mkvseqs) and eof(ranbook) then begin writeln(output, 'Neither mkvseqs or ranbook exists!'); writeln(output, 'Create a source for random sequences.'); halt; end; reset(book); rewrite(embedbk); readparameters(alignmenttype, InFrom, InTo, OutFrom, OutTo); (* get an initial random string to use in the embedding *) readrandom(rand, Useranbook, OutTo, OutFrom, ranbook, thelineran); reset(inst); reset(book); new(apiece); piececount := 0; while not eof(book) do begin copybook; (* copy up to the next piece or end of book *) if not eof(book) then begin (* read in a piece *) case alignmenttype of 'i': begin align(inst,book,theline,apiece,length,alignedbase); end; 'b': begin getpiece(book,theline,apiece); (* read in the piece *) alignedbase := pietoint(0, apiece); end; 'f': begin getpiece(book,theline,apiece); alignedbase:=1; end; end; (* tell user which piece we are working on *) piececount := succ(piececount); write(output,'piece ',piececount:1,' '); for s:= 1 to apiece^.key.hea.keynam.length do write(output, apiece^.key.hea.keynam.letters[s]); writeln(output); (* replace appropriate parts of the piece with random material *) fillrandom(alignedbase, InFrom, InTo, OutFrom, OutTo, apiece, rand, Useranbook,ranbook, thelineran); (* write the piece out *) bwpie(embedbk, apiece); end; end; writeln(output,'embed is done'); end; (* end module embed.themain *) begin themain(inst, book, embedp, mkvseqs, ranbook, embedbk); 1: end.