program hist(inst,book,hst,histp,output);
(* makes a table (3-d histogram) of occurence of oligonucleotides
      in a book of sequences aligned by their instructions.
      gary stormo
modules from: delman, prgmods, delmods *)

label 1;  (* end of program *)

const
(* begin module version *)
version = 4.27; (* of hist.p 1994 sep 6
origin: before 1986 Nov 14 *)
(* end module version *)

(* begin module describe.hist *)
(*
name
      hist: make a histogram of aligned sequences.

synopsis
      hist(inst: in, book: in, hst: out, histp: in, output: out);

files
      inst: the instructions which generated the book, used to
            align the sequences; the format of the instructions must be
            'get from # -x to # +y;' - the alignment is done on the base #.
            if this file is empty the sequences are aligned by their 5
            prime ends.
      book: the sequences, not longer than 'dnamax' (see technical note);
      hst: the histogram table, giving the position occurence
            of all oligonucleotides from oligomin to oligomax  in length
            (see file histp).  note - if the histogram is wider than can
            be printed on a page, use the program split to print hst.
      histp: to set the length of oligonucleotides searched for; contains two
            integers, one per line, the first specifying oligomin
            (the length of the shortest oligonucleotide which is searched for),
            the second oligomax (the longest oligonucleotide searched for,
            which cannot be greater than 4);
            note: if oligomin is zero, then the bases are counted.
      output: for error messages.

description
      makes a histogram of the occurences of oligonucleotides at positions
      relative to some aligned base.  this is done for all oligonucleotides
      with lengths from  'oligomin' to 'oligomax', set by file histp.

see also
      split.p, alist.p, histan.p

author
      gary stormo

bugs
      none known

technical note
      the constant 'dnamax'  from the module book.const may be too
      large for efficient use of this program.  if you do not expect to do
      histograms on aligned sequences of greater than, say, 120 nucleotides
      you can go into the source program (hists) and change 'dnamax' before
      compiling.
*)
(* end module describe.hist *)

(* begin module book.const *)
(* constants needed for book manipulations *)

      dnamax = 3000;  (* length of dna arrays *)
      namelength = 20; (* maximum key name length *)
      linelength = 80; (* maximum line readable in book *)

(* end module book.const version = 'delmod 6.61 91 Mar 19 tds/gds' *)

type
(* 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;


      (* base types *)
      base = (a,c,g,t);
      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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

      alpieceptr = ^alpiece;
      alpiece = record   (* an aligned piece *)
         length : integer;
         alignedbase : integer;
         dna : dnaptr;
         next : alpieceptr
      end;

var
(* 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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

      hst,  (* the output histogram *)
      histp,  (*the parameter file *)
      book,inst : text;
      oligomin, (* shortest oligos  used to scan dna *)
      oligomax, (* longest  oligos  used to  scan dna  *)
      oligolength : integer;   (*  oligomin..oligomax *)
      seqnumb : integer;  (* number of pieces in book *)
      maxaligned,  (* maximum bases to (and including) aligned base *)
      maxafter : integer; (* maximum bases after aligned base *)
      alpie : alpieceptr;
      first : alpieceptr;
      apiece :  pieceptr; (* piece  used in align routine  *)

(* 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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* 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 *)
begin (* numbersize *)
      if n = 0
         then numbersize:=1
         else numbersize:=trunc(ln(abs(n))/ln10 + epsilon) + 2;
(* 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. *)
(* the 2 is for the sign and last digit *)
end; (* numbersize *)
(* end module numbersize version = 'prgmod 3.96  85 mar 18 tds'; *)

(* 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 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 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<coobeg
                  then if coocon=circular
                     then p:=p+(cooend-coobeg+1)
            end
         end;
         inttopie:=p
      end
end;

function piecelength(pie: pieceptr): integer;
(* return the length of the dna in pie *)
begin
      piecelength:=pietoint(pie^.key.pieend,pie)
end;

(* end module book.basis version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.getto *)
function getto(var thefile: text; ch: chset): char;
(* search the file for a character in the first line which is a
   member of the set ch. *)
var achar: char;
begin
      achar:=' ';
      while (not(achar in ch)) and (not eof(thefile))
         do readln(thefile,achar);
      if (achar in ch) then getto:=achar
                       else getto:=' '
end;
(* end module book.getto version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.skipstar *)
procedure skipstar(var thefile:text);
(* skip start of line (or star = '*'). *)
begin (* skipstar *)
      if thefile^ <> '*' 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; (* skipstar *)
(* end module book.skipstar version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brreanum *)
procedure brreanum(var thefile: text; var reanum: real);
(* read a real number from the file *)
begin
      skipstar(thefile);
      readln(thefile,reanum);
end;
(* end module book.brreanum version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brnumber *)
procedure brnumber(var thefile: text; var num: integer);
(* read a number from the file *)
begin
      skipstar(thefile);
      readln(thefile,num)
end;
(* end module book.brnumber version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brname *)
procedure brname(var thefile: text; 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<namelength
            then for i:=length+1 to namelength do letters[i]:=' '
      end;
      readln(thefile)
end; (* brname *)
(* end module book.brname version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brline *)
procedure brline(var thefile: text; var l: lineptr);
(* read a line from the file *)
var
      i,j: integer;
      acharacter: char;
begin
      skipstar(thefile);
      i:=0;
      while (not eoln(thefile)) do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      if i<l^.length then for j:=i+1 to l^.length do l^.letters[j]:=' ';
      l^.length:=i;
      l^.next:=nil;
      readln(thefile)
end;
(* end module book.brline version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brdirect *)
procedure brdirect(var thefile: text; var direct: direction);
(* read a direction *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='+' then direct:=plus
                else direct:=minus
end;
(* end module book.brdirect version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brconfig *)
procedure brconfig(var thefile: text; var config: configuration);
(* read a configuration *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='l' then config:=linear
                else config:=circular
end;
(* end module book.brconfig version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brnotenumber *)
procedure brnotenumber(var thefile: text; 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);
         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)
            until thefile^ = 'n';
            readln(thefile)
         end
         else readln(thefile)
      end
end; (* brnotenumber *)
(* end module book.brnotenumber version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brnote *)
procedure brnote(var thefile: text; 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);
         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,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)
         end
         else readln(thefile)
      end
end; (* brnote *)
(* end module book.brnote version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brheader *)
procedure brheader(var thefile: text; var hea: header);
(* read the header of a key. *)
begin
      with hea do begin
         (* read key name *)
         brname(thefile,keynam);

         (* read full name *)
         getline(fulnam);
         brline(thefile,fulnam);

         (* read note key *)
         if readnumber then brnotenumber(thefile,note)
                       else brnote(thefile,note)
      end
end;
(* end module book.brheader version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brpiekey *)
procedure brpiekey(var thefile: text; var pie: piekey);
(* read piece key *)
begin
      with pie do begin
         brheader(thefile,hea);
         brreanum(thefile,mapbeg);
         brconfig(thefile,coocon);
         brdirect(thefile,coodir);
         brnumber(thefile,coobeg);
         brnumber(thefile,cooend);
         brconfig(thefile,piecon);
         brdirect(thefile,piedir);
         brnumber(thefile,piebeg);
         brnumber(thefile,pieend);
      end
end;
(* end module book.brpiekey version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brdna *)
procedure brdna(var thefile: text; var dna: dnaptr);
(* read in dna from thefile *)
(* 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,['d']);
      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 *)
         read(thefile,ch); (* ch is either '*' or 'd' *)
      end;
      readln(thefile)
end;
(* end module book.brdna version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brpiece *)
procedure brpiece(var thefile: text; var pie: pieceptr);
(* read in a piece *)
begin
      brpiekey(thefile,pie^.key);
      if numbered or (not skipunnum)
         then brdna(thefile,pie^.dna)
end;
(* end module book.brpiece version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.brinit *)
procedure brinit(var book: text);
(* 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;
end; (* brinit *)
(* end module book.brinit version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* ************************************************************************ *)
(* end module package.brpiece version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module book.getpiece *)
procedure getpiece(var thefile: text; var pie: pieceptr);
(* move to and read in the next piece in the book *)
var   ch: char;
begin
      ch:=getto(thefile,['p']); (* get to the next p(iece) in the book *)
      if ch<>' ' then begin
         brpiece(thefile,pie);
         ch:=getto(thefile,['p']); (* read past closing p *)
      end
end;
(* end module book.getpiece version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* ************************************************************************ *)
(* end module package.getpiece version = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* 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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* 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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* 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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

(* begin module align.align *)
procedure align(var inst, book: text;
                 var pie:  pieceptr;
                 var length, alignedbase: integer);
(* documentation on align is in module info.align and
delman.use.aligned.books *)
const maximumrange = 500; (* 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. *)
var
   ch: char; (* a character in inst *)
   comment: boolean; (* true means we are inside a comment *)
   done: boolean; (* done finding an aligning get *)
   thebase: integer; (* the base read in *)
begin
   if not eof(book) then begin (* if there is still more to the book ... *)
      getpiece(book,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 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 inst^ = '(' then begin (* skip comment *)
                  get(inst);
                  if inst^ = '*' then comment := true;
    {if comment then write(output,'COMMENT: (');}
                  while comment do begin
                     if eof(inst) then begin
                        writeln(output,' in procedure align:');
                        writeln(output,' an instruction comment does not end!');
                        halt
                     end;
    {write(output,inst^);}
                     get(inst);
                     if inst^ = '*' then begin
                        get(inst);
    {if inst^ = ')' then writeln(output,'*)');}
                        if inst^ = ')' then comment := false
                     end;
                  end;
               end;
               if inst^ = 'g' then begin
                  get(inst);
                  if inst^ = 'e' then begin
                  get(inst);
                     if inst^ = 't' then begin
                        get(inst);
                        if inst^ = ' ' then begin
                           findnonblank(inst,ch); (* get to "from" *)
                           findblank(inst); (* get past "from" *)
                           read(inst,thebase); (* read in the alignedbase *)
    {writeln(output,'thebase=',thebase:1);}
                           alignedbase:=pietoint(thebase,pie);
                           done := true
                        end;
                     end;
                  end;
               end;
               get(inst); (* move along now *)
            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 = 'delmod 6.61 91 Mar 19 tds/gds' *)

procedure histogramwidth; (* this determines maxaligned and
               maxafter, which determine the histogram width *)
begin
      if maxaligned < alpie^.alignedbase then
         maxaligned := alpie^.alignedbase;
      if maxafter < (alpie^.length - alpie^.alignedbase) then
         maxafter := (alpie^.length - alpie^.alignedbase);
      if ((maxaligned + maxafter) > dnamax) then
      begin
         write(output,'dnamax too small');
         halt
      end;
end;


function getalbase(position: integer; pie: alpieceptr):base;
(* get a base from the nth position (internal coordinates) of the
   piece.  no protection is made against positions outside the piece *)
var
      workdna: dnaptr;
      p: integer; (* the last base of the dna part *)
begin
      workdna:=pie^.dna;
      p:=dnamax;
      while position>p do begin
         p:=p+dnamax;
         workdna:=workdna^.next
      end;
      getalbase:=workdna^.part[position-(p-dnamax)]
end;

procedure makehist;
(* this writes the relevant information to hst and then calls
   getnextletter, which creates the oligos that are searched for *)

var  histwidth, (* width of histogram *)
     histfield,  (* size of field for printing histogram values *)
     i : integer;
     oligoarray : array [1..4] of base; (* oligos used in scanning, max is 4 *)
     oligolength : integer;   (*  oligomin..oligomax, and an index *)

procedure getnextletter(k : integer);
(* fills oligoarray with bases up to length k, then calls scan to
   search for the oligonucleotide in oligoarray *)
var ba : base;


procedure scan;
(* scans the sequence for  the occurence of the oligonucleotide in oligoarray
  and when found increaes the tablearray  at that position by 1 *)
var i,j,k,l : integer;
    match : boolean;
    tablearray : array[1..dnamax] of integer;

procedure  tableout;
(* write to hst the histogram row for the oligonucleotide just
 searched for *)
var  i : integer;

begin
      for i := 0  to (oligomax - oligolength) do
         write(hst,' ');
      for i := (oligomax - oligolength + 1)  to oligomax do
         write(hst,basetochar(oligoarray[i]));
      write(hst,' ');
      for i := 1  to histwidth do write(hst,' ',tablearray[i]:histfield);
      writeln(hst);
         (* write a blank line after every fourth line in the histogram
            or after the count of each column *)
      if (oligoarray[oligomax] = t) or (oligolength = 0)
         then writeln(hst)
end; (* tableout *)

begin (* scan *)
      for i := 1 to histwidth  do tablearray[i] := 0;
      alpie := first;
      for j := 1 to seqnumb do
      begin
         for l := 1 to (alpie^.length - (oligolength - 1)) do
         begin
            k := 0;
            match := true;
            while ((k < oligolength) and match) do
               if getalbase(k+l,alpie) =  oligoarray[oligomax-oligolength+k+1]
                  then k := k+1
               else match := false;
            if match then tablearray[maxaligned - alpie^.alignedbase + l] :=
                succ(tablearray[maxaligned - alpie^.alignedbase + l]);
         end;
         alpie := alpie^.next;
      end;
      tableout;
end; (* scan *)

begin (* getnextletter *)
      if k = 0 then scan
      else for ba :=  a to t do begin
         oligoarray[oligomax - (k-1)] := ba;
         getnextletter(k-1)
      end;
end; (* getnextletter *)

begin (* makehist *)
            (* the width of the histogram is *)
      histwidth := maxafter + maxaligned;
            (* the field size for the histogram values is enough to print the
            largest possible value, i.e., the number of digits in seqnumb,
            or the number of digits in the position list at the top *)
      histfield := numbersize(seqnumb);
      if histfield < numbersize(1-maxaligned)
         then histfield := numbersize(1-maxaligned);
      if histfield < numbersize(maxafter)
         then histfield := numbersize(maxafter);

            (* write to hst the parameters *)
      writeln(hst,oligomin:4,' is the shortest oligo searched for');
      writeln(hst,oligomax:4,' is the longest oligo searched for');
      writeln(hst,seqnumb:4,' is the number of sequences in the book');
      writeln(hst,histwidth:4,' is the width (in bases) of the histogram');
      writeln(hst);
            (* write position line to hst *)
      for i := 0 to (oligomax + 1) do write(hst,' ');
      for i := (1-maxaligned) to maxafter do write(hst,' ',i:histfield);
      writeln(hst);
               (* generate a divider line before the histogram *)
      for i := 0 to (oligomax + 1) do write(hst,' ');
      for i := 1 to (histwidth * (histfield + 1)) do write(hst,'-');
      writeln(hst);
            (* generate the oligos and search for them *)
      for oligolength := oligomin to oligomax do getnextletter(oligolength);
end; (* makehist *)

procedure setparam; (* get the parameters from histp, or use defaults *)
begin
      if not eof(histp) then
      begin
         readln(histp,oligomin);
         if eof(histp) then begin
            writeln(output,'there must be a second parameter in histp');
            halt
         end;
         readln(histp,oligomax);
                  (* check for valid parameters *)
         if (oligomin > oligomax) then begin
            writeln(output,'oligomin cannot be greater than oligomax');
            halt
         end;
         if ((oligomin > 4) or (oligomax > 4)) then begin
            writeln(output,'oligomin and oligomax must be no larger than 4');
            halt
         end;
         if oligomin < 0 then begin
            writeln(output,'oligomin must be positive or zero');
            halt
         end

      end
      else begin
         oligomin := 0;
         oligomax := 2;
      end
end; (* setparam *)

begin  (* hist ************* *)
      writeln(output,'hist ',version:4:2);

      brinit(book);
      reset(inst);
      reset(histp);
      setparam;
      rewrite(hst);
      writeln(hst,'* hist ',version:4:2,' histogram of ');
      copyaline(book,hst);
      new(apiece);

      maxaligned :=  0;
      maxafter := 0;
      seqnumb  := 0;
      first := nil;

      while not eof(book) do
      begin
         new(alpie);
         align(inst,book,apiece,alpie^.length,alpie^.alignedbase);
         if not eof(book) then
         begin
            seqnumb  := succ(seqnumb);
            alpie^.dna := apiece^.dna;
            apiece^.dna := nil;
            alpie^.next  := first;
            first  := alpie;
            clearpiece(apiece);
            histogramwidth;
         end;
      end;

      makehist;  (* this makes each of the oligos
               from oligomin to oligomax, scans  the book for each,  and
               writes the table to hist *)
1: end.