program compan(cmp,anal,companp,output);
(* compan: chi-squared analysis of a composition.
   by Gary Stormo
module libraries needed: delman, delmods, auxmods  *)

label 1;  (* end of program *)

const
(* begin module version *)
version = 3.25;  (* of compan.p 2001 Mar 26
2001 Mar 26, 3.25: spelling correction
previous change 1994 Sep 5
origin before 1985 apr 18 *)
(* end module version *)

(* begin module describe.compan *)
(*
name
      compan: composition analysis.

synopsis
      compan(cmp: in, anal: out, companp: in, output: out)

files
      cmp: the input composition, which is the output of program comp;
      anal: the output analysis of this program;
      companp: for parameters; should contain a single integer which specifies
         the maximum level for which the composition is analyzed.  the
         maximum allowed level is 4, or the maximum level for which the
         composition was determined.
      output: for user messages;

description
      calculates chi squared from a composition using:
         1) assumption of equal frequencies to calculate mono, di, tri
            and tet expecteds;
         2) mono frequencies to calculate di, tri  and tet expecteds;
         3) di frequencies to calculate tri and tet expecteds;
         4) tri frequencies to calculate tet expecteds;
      the partial chi squared values are given for each oligo.
      the 'information' content of the composition is also calculated,
      using the standard information theory definition:
         information = -sum(frequency * log(frequency)),
      where the sum is over each oligonucleotide of a given length
      and the log is taken to the base 2.  this gives the information
      in bits.

see also
      comp.p

author
      gary stormo

bugs
      the program cannot do calculations for compositions larger than 4
*)
(* end module describe.compan *)

      defcompmax = 4;  (* the limit, and default value, for analysis *)

type  base = (a,c,g,t);
(* begin module comp.type *)
            (* the composition is stored in a tree of these nodes *)
      compnodeptr = ^compnode; (* points to a node of the tree *)
      compnode = record    (* a node of the composition tree *)
         count : integer; (* the number of oligos for this node *)
         son : array[base] of compnodeptr; (* the pointers to 'descendants'
                                           of this node of the tree *)
      end;

            (* spiders are used to make the composition tree *)
      spiderptr = ^spider; (* points to a 'spider' *)
      spider = record  (* a spider climbs the composition tree, its path
                        determined by the sequences, and increments 'count' at
                        all the nodes it passes, thereby determining the
                        composition *)
         depth : integer; (* the level of the node now at *)
         place : compnodeptr; (* a pointer to the current node *)
         next : spiderptr; (* the next spider in the collection *)
      end;

            (* the total number of composition entries at a given level
               is stored in the linked list of type comptotal *)
      compptr = ^comptotal;
      comptotal = record
         count: integer;  (* the number at a given level *)
         next: compptr;  (* pointer to the next level totals *)
      end;

           (* 'path' is used in printing the tree *)
      pathptr = ^path; (* pointer into a path on the tree *)
      path = record  (* the path of bases to get to a particular node *)
         bas : base;
         next : pathptr;
      end;
(* end module comp.type version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

var   cmp,anal,companp : text;
      root: compnodeptr;  (* of the composition tree *)
      monocomptotal: compptr; (* beginning of comp totals list *)
      compmax,  (* maximum length oligo for which comp was determined *)
      readmax: integer;  (* maximum amount of the comp used *)

(* begin module package.primitive *)
(* ************************************************************************ *)
(* 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.51 85 apr 17 tds/gds' *)

(* begin module unlimitln *)
procedure unlimitln(var afile: text);
(* this procedure removes a stupid system dependent limit on the number of
lines that one can write to a file.  you may remove it from the code if
your system does not want or need this.  suggested method: place comments
around the contents of the procedure. *)
begin
(*    linelimit(afile, maxint); (@ set 'infinite' lines allowed for afile *)
end;
(* end module unlimitln version = 'delmod 6.51 85 apr 17 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.51 85 apr 17 tds/gds' *)

(* begin module copylines *)
function copylines(var fin, fout: text; n: integer): integer;
(* copy n lines of file fin to file fout.
the actual number of lines copied is returned. *)
var
      index: integer; (* the current line number *)
begin (* copylines *)
      index := 0;
      while (not eof(fin)) and (index < n) do begin
         copyaline(fin, fout);
         index := succ(index)
      end;

      copylines := index
end; (* copylines *)
(* end module copylines version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.primitive version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module basetochar *)
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;
(* end module basetochar version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

(* begin module comp.readcomp *)

(* begin module skipoligo *)
procedure skipoligo(var thefile: text);
(* this procedure is used to skip over an oligonucleotide string.
   it reads until it comes to a base, and then continues to read until
   it comes to a blank. no checking is done for non-base characters
   in between. *)
var ch : char;
begin
      repeat read(thefile,ch)
      until (ch in ['a','c','g','t']);
      repeat read(thefile,ch)
      until (ch = ' ');
end;
(* end module skipoligo version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

procedure readcomp(var comp: text; var compmax: integer; readmax: integer;
                   var root: compnodeptr; var monocomptotal: compptr);
(* this procedure requires modules: comp.type, skipoligo, halt;
   this procedure reads from file 'comp' a composition and puts it into
   the tree pointed to by the 'root' pointer.  'compmax' is the depth
   of the composition tree which is stored.  it is the minimum of
   'readmax', the requested depth, and 'detcomp', the depth for which
   the input file composition was determined.  'monocomptotal' points
   to the beginning (for the monos) of a linked list which gives the
   totals for each level of the composition. *)
type  (* for efficient reading  of the information into the tree, each
         node that has a successor is stored in a queue as well as put
         into the tree.  the variables of the type list are the queue
         elements. *)
    listptr = ^list;
    list = record
      item: compnodeptr;  (* points to the composition tree node *)
      next: listptr;  (* points to the next list item in the queue *)
    end;

var  freeitem,  (* the list of unused list pointers *)
     listitem, (* an item in the queue *)
     first,  (* the first item in the queue *)
     last: listptr;  (* the last item in the queue *)
     comptot,  (* the comp total for a given level *)
     newcomptot: compptr;  (* for adding to the string of comptot"s *)
     detcomp,  (* the level to which the input composition was determined *)
     level,  (* of the composition being read, i.e., monos, dis, ... *)
     number: integer;  (* read from the 'comp' file *)
     ch: char;  (* for reading from 'comp' *)
     ba: base;  (* an index *)

(* getitem and clearitem provide efficient use of linked list storage by
   keeping a list of unused pointers that can be allocated instead of
   always creating 'new' ones  *)
procedure getitem(var l: listptr);
(* obtain a listitem from the free list or by making a new one *)
begin
      if freeitem<>nil
         then begin
            l:=freeitem;
            freeitem:=freeitem^.next
         end
         else new(l);
      l^.next:=nil
end;


procedure clearitem(var l: listptr);
(* return a listitem to the free list *)
var   lptr: listptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freeitem;
         freeitem:=lptr
      end
end;

begin  (* readcomp *)
      reset(comp);
      if eof(comp) then begin
         writeln(output,' error: no composition file provided');
         halt;
      end;
      readln(comp); (* skip the program identification *)
      readln(comp); (* skip the book identification *)
      readln(comp,detcomp); (* obtain the determined composition *)
      readln(comp); (* skip the blank line *)

            (* determine the level of composition to be stored *)
      if (readmax < 1) then begin
         writeln(output);
         writeln(output,' warning: 0 or negative oligo length requested');
         writeln(output,' composition used is depth ',detcomp:1);
         writeln(output);
         compmax := detcomp
      end
      else if (readmax > detcomp) then begin
         writeln(output);
         writeln(output,' warning: requested composition oligo',
                        ' length (',readmax:1,')');
         writeln(output,' is larger than the determined composition',
                        ' oligo length (',detcomp:1,').');
         writeln(output,' composition used is to depth ',detcomp:1);
         writeln(output);
         compmax := detcomp
      end
      else compmax := readmax;

      new(root);
      new(monocomptotal);
      new(comptot);
      new(first);
      new(last);
      first^.item := root;
      first^.next := nil;
      last := first;
      new(freeitem);
      freeitem^.next := nil;

         (* read in the total number of bases from the composition *)
      readln(comp); (* skip the * *)
      readln(comp); (* skip the information line *)
      readln(comp,number);
      root^.count := number;

      repeat
         read(comp,ch); (* skip the space *)
         read(comp,ch); (* the ch determines what to do next *)
         if (ch = '*') then begin  (* determine the level about to be read *)
            readln(comp);
            readln(comp,level);
            if (level = 1) then monocomptotal := comptot
            else begin
               new(newcomptot);
               comptot^.next := newcomptot;
               comptot := newcomptot;
            end;
            comptot^.count := 0;
            comptot^.next := nil
         end
         else begin  (* read the composition values *)
            for ba := a to t do begin
               skipoligo(comp);
               read(comp,number);
               if (number <> 0) then begin
                  new(first^.item^.son[ba]);
                  first^.item^.son[ba]^.count := number;
                  comptot^.count := comptot^.count + number;
                  getitem(listitem);
                  last^.next := listitem;
                  last := listitem;
                  last^.next := nil;
                  last^.item := first^.item^.son[ba]
               end
               else first^.item^.son[ba] := nil;
            end;
            clearitem(first);
            readln(comp);
         end;
      until level = compmax;
         (* for the last level of compositions to be read we do not
            need to store the nodes in the queue because their successors
            will not be read in. *)
      repeat
         for ba := a to t do begin
            skipoligo(comp);
            read(comp,number);
            if (number <> 0) then begin
               new(first^.item^.son[ba]);
               first^.item^.son[ba]^.count := number;
               comptot^.count := comptot^.count + number;
            end
            else first^.item^.son[ba] := nil;
         end;
         clearitem(first);
         readln(comp);
      until (first = nil);
end;  (* readcomp *)

function getcount(root: compnodeptr; start: pathptr): integer;
(* this function follows the tree from the node 'root' (which
   may be the root of a subtree) along the path initiated with
   'start', and returns the count of the resulting node.  if the
   path through the tree ever hits a null node in the process
   the count is returned to be zero. *)
var place: compnodeptr;  (* a place in the composition tree *)
    point: pathptr;  (* a point in the path through the tree *)
begin
      point:= start;
      place := root;
      while ((place <> nil) and (point <> nil)) do begin
         place := place^.son[point^.bas];
         point := point^.next;
      end;
      if (place = nil) then getcount := 0
      else getcount := place^.count;
end;

(* end module comp.readcomp version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

procedure header;
(* write the header to the anal file *)
begin
      reset(cmp);
      writeln(anal,' compan version ',version:4:2,
            '; chi-squared analysis of composition');
      if copylines(cmp,anal,2)<>2 then begin
         writeln(output, ' composition file is too short');
         halt
      end;

      writeln(anal);
      writeln(anal,' ',compmax:1,'-long oligos are the longest analyzed');
      write  (anal,' (+) means observed > expected;');
      writeln(anal,' (-) means observed < expected');
      writeln(anal);
end;

procedure chisq(observed, expected: real; var chi: real; var sign: char);
(* given observed and expected values, calculate chisquared and return
   as variable chi; set sign to '+' if observed is >= expected, else '-' *)
var diff : real; (* difference between observed and expected *)

begin
      diff := observed - expected;
      if (diff >= 0) then sign := '+'
      else sign := '-';
      if (expected > 0) then chi := diff * diff / expected
      else chi := 0;
end;

procedure calc0; (* calculate expecteds assuming equal frequencies *)
var   b1,b2,b3,b4 : base;
      comptot,  (* total number of oligos of a certain length *)
      observed: integer;  (* number of oligos of a certain composition *)
      expected,  (* number of oligos of a certain composition *)
      chi,  (* the partial chisquared for each entry *)
      totchi, (* the total chisquared for each section *)
      frequency,  (* of each oligo *)
      inf: real; (* the 'information' at each level of composition, as defined
         by information theory: inf = -sum(frequency(log(frequency)))  *)
      start,  (* of the path through the tree *)
      point: pathptr;  (* on the path through the tree *)
      sign : char; (* '+' or '-', as observed is more or less than expected *)

begin
      writeln(anal,' **** chi squareds assuming equal frequencies ****');
      writeln(anal);

          (* chi squareds for monos *)
      totchi := 0;
      inf := 0;
      comptot := monocomptotal^.count;
      expected := comptot/4;
      writeln(anal,' partials');
      new(start);
      start^.next := nil;
      for b1 := a to t do begin
         start^.bas := b1;
         observed := getcount(root,start);
         chisq(observed,expected,chi,sign);
         write(anal,' ':5,basetochar(b1),chi:8:2,' (',sign,');');
         totchi := totchi + chi;
         if (observed > 0) then begin
            frequency := observed / comptot;
            inf := inf - frequency * ln(frequency) / ln(2);
         end;
      end;
      writeln(anal);
      writeln(anal,' * total chisquared = ',totchi:8:2);
      writeln(anal,' * information = ',inf:5:3,' bits per base');
      writeln(anal);

         (* chi squareds for dis *)
      if (compmax >= 2) then begin
         totchi := 0;
         inf := 0;
         comptot := monocomptotal^.next^.count;
         expected := comptot/16;
         writeln(anal,' partials');
         new(point);
         point^.next := nil;
         start^.next := point;
         for b1 := a to t do begin
            start^.bas := b1;
            for b2 := a to t do begin
               point^.bas := b2;
               observed := getcount(root,start);
               chisq(observed,expected,chi,sign);
               write(anal,' ':4,basetochar(b1),basetochar(b2),
                  chi:8:2,' (',sign,');');
               totchi := totchi + chi;
               if (observed > 0) then begin
                  frequency := observed / comptot;
                  inf := inf - frequency * ln(frequency) / ln(2);
               end;
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
         writeln(anal,' * information = ',inf/2:5:3,' bits per base');
         writeln(anal);
      end;

         (* chi squareds for tris *)
      if (compmax >= 3) then begin
         totchi := 0;
         inf := 0;
         comptot := monocomptotal^.next^.next^.count;
         expected := comptot/64;
         writeln(anal,' partials');
         new(point);
         point^.next := nil;
         start^.next^.next := point;
         for b1 := a to t do begin
            start^.bas := b1;
            for b2 := a to t do begin
               start^.next^.bas := b2;
               for b3 := a to t do begin
                  point^.bas := b3;
                  observed := getcount(root,start);
                  chisq(observed,expected,chi,sign);
                  write(anal,' ':3,basetochar(b1),basetochar(b2),basetochar(b3),
                     chi:8:2,' (',sign,');');
                  totchi := totchi + chi;
                  if (observed > 0) then begin
                     frequency := observed / comptot;
                     inf := inf - frequency * ln(frequency) / ln(2);
                  end;
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
         writeln(anal,' * information = ',inf/3:5:3,' bits per base');
         writeln(anal);
      end;

         (* chi squareds for tets *)
      if (compmax >= 4) then begin
         totchi := 0;
         inf := 0;
         comptot := monocomptotal^.next^.next^.next^.count;
         expected := comptot/256;
         writeln(anal,' partials');
         new(point);
         point^.next := nil;
         start^.next^.next^.next := point;
         for b1 := a to t do begin
            start^.bas := b1;
            for b2 := a to t do begin
               start^.next^.bas := b2;
               for b3 := a to t do begin
                  start^.next^.next^.bas := b3;
                  for b4 := a to t do begin
                     point^.bas := b4;
                     observed := getcount(root,start);
                     chisq(observed,expected,chi,sign);
                     write(anal,' ':2,basetochar(b1),basetochar(b2),
                        basetochar(b3),basetochar(b4),chi:8:2,' (',sign,');');
                     totchi := totchi + chi;
                     if (observed > 0) then begin
                        frequency := observed / comptot;
                        inf := inf - frequency * ln(frequency) / ln(2);
                     end;
                  end;
                  writeln(anal);
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
         writeln(anal,' * information = ',inf/4:5:3,' bits per base');
         writeln(anal);
      end;
end;

procedure calc1;  (* calculate chi squareds using mono frequencies *)
var b1,b2,b3,b4 : base;
      observed,  (* number of oligos of a certain composition *)
      expected,  (* number of oligos of a certain composition *)
                 (* using monos, the expected is the frequency of each
                    mono in the oligo times the total number
                    of oligos of a particular length *)
      ratio,  (* for calculating expecteds *)
      n1,n2,n3,n4, (* for storing the results of 'getcount' calls *)
      chi,  (* partial chisquareds *)
      totchi : real; (* total chisquareds *)
      sign : char; (* '+' or '-', as observed is more or less than expected *)
      point1,  (* a point on a path through the composition tree *)
      point2,               (*  '  *)
      point3,               (*  '  *)
      point4: pathptr;      (*  '  *)

begin
      page(anal);
      writeln(anal,' **** chi squareds using mono frequencies ****');
      writeln(anal);

         (* chi squareds for dis *)
      totchi := 0;
      ratio := ln(monocomptotal^.next^.count) - 2*ln(monocomptotal^.count);
      writeln(anal,' partials');
      new(point1);
      new(point2);
      point2^.next := nil;
      for b1 := a to t do begin
         point1^.bas := b1;
         point1^.next := nil;
         n1 := getcount(root,point1);
         for b2 := a to t do begin
            point2^.bas := b2;
            n2 := getcount(root,point2);
            if ((n1 > 0) and (n2 > 0)) then
               expected := exp(ln(n1) + ln(n2) + ratio)
            else expected := 0;
            point1^.next := point2;
            observed := getcount(root,point1);
            chisq(observed,expected,chi,sign);
            write(anal,' ':4,basetochar(b1),basetochar(b2),
               chi:8:2,' (',sign,');');
            totchi := totchi + chi;
         end;
         writeln(anal);
      end;
      writeln(anal,' * total chisquared = ',totchi:8:2);
      writeln(anal);

         (* chi squareds for tris *)
      if (compmax >= 3) then begin
         totchi := 0;
         ratio := ln(monocomptotal^.next^.next^.count) -
             3 * ln(monocomptotal^.count);
         writeln(anal,' partials');
         new(point3);
         point3^.next := nil;
         for b1 := a to t do begin
            point1^.bas := b1;
            point1^.next := nil;
            n1 := getcount(root,point1);
            for b2 := a to t do begin
               point2^.bas := b2;
               point2 ^.next := nil;
               n2 := getcount(root,point2);
               point1^.next := point2;
               for b3 := a to t do begin
                  point3^.bas := b3;
                  n3 := getcount(root,point3);
                  if ((n1 > 0) and (n2 > 0) and (n3 > 0)) then
                     expected := exp(ln(n1) + ln(n2) + ln(n3) + ratio)
                  else expected := 0;
                  point2^.next := point3;
                  observed := getcount(root,point1);
                  chisq(observed,expected,chi,sign);
                  write(anal,' ':3,basetochar(b1),basetochar(b2),basetochar(b3),
                     chi:8:2,' (',sign,');');
                  totchi := totchi + chi;
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
         writeln(anal);
      end;

         (* chi squareds for tets *)
      if (compmax >= 4) then begin
         totchi := 0;
         ratio := ln(monocomptotal^.next^.next^.next^.count) -
            4 * ln(monocomptotal^.count);
         writeln(anal,' partials');
         new(point4);
         point4^.next := nil;
         for b1 := a to t do begin
            point1^.bas := b1;
            point1^.next := nil;
            n1 := getcount(root,point1);
            for b2 := a to t do begin
               point2^.bas := b2;
               point2^.next := nil;
               n2 := getcount(root,point2);
               point1^.next := point2;
               for b3 := a to t do begin
                  point3^.bas := b3;
                  point3^.next := nil;
                  n3 := getcount(root,point3);
                  point2^.next := point3;
                  for b4 := a to t do begin
                     point4^.bas := b4;
                     n4 := getcount(root,point4);
                     if ((n1 > 0) and (n2 > 0) and (n3 > 0) and (n4 > 0)) then
                        expected := exp(ln(n1)+ ln(n2)+ ln(n3)+ ln(n4)+ ratio)
                     else expected := 0;
                     point3^.next := point4;
                     observed := getcount(root,point1);
                     chisq(observed,expected,chi,sign);
                     write(anal,' ':3,basetochar(b1),basetochar(b2),
                        basetochar(b3),basetochar(b4),chi:8:2,' (',sign,');');
                     totchi := totchi + chi;
                  end;
                  writeln(anal);
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
         writeln(anal);
      end;
end;

procedure calc2;  (* calculate chi squareds using mono and di frequencies *)
var   b1,b2,b3,b4 : base;
      observed,  (* number of a particular oligo *)
      expected,  (* number of a particular oligo, see sections for details *)
      ratio,  (* for calculating expecteds *)
      ratio2,  (* for calculating expecteds *)
      ratio3,  (* for calculating expecteds *)
      n1,n2,n3, (* for storing values from 'getcount' calls *)
      chi,  (* partial chisquareds *)
      totchi : real; (* total chisquareds *)
      sign : char; (* '+' or '-', as observed is more or less than expected *)
      point1,  (* on the path through the tree *)
      point2,            (*   '   *)
      point3,            (*   '   *)
      point4: pathptr;   (*   '   *)

begin
      page(anal);
      writeln(anal,' **** chi squareds using di frequencies ****');
      writeln(anal);

         (* chi squareds for tris *)
         (* calculating expected tris from dis uses this formula:
            expected[b1,b2,b3] = #[b1,b2]/2totals * #[b2,b3]/2totals
               * 1totals/#[b2] * 3totals.  *)
      totchi := 0;
         (* ratio = ln( 1totals * 3totals / (2totals)**2) <= 0  *)
      ratio := ln(monocomptotal^.count) +
         ln(monocomptotal^.next^.next^.count) -
         2 * ln(monocomptotal^.next^.count);
      writeln(anal,' partials');

      new(point1);
      new(point2);
      new(point3);
      point3^.next := nil;
      for b1 := a to t do begin
         point1^.bas := b1;
         for b2 := a to t do begin
            point2^.bas := b2;
            point1^.next := point2;
            point2^.next := nil;
            n1 := getcount(root,point1);
            n2 := getcount(root,point2);
            if (n1 > 0) then
               ratio2 := ratio + ln(n1) - ln(n2)
                  (* ratio2 <= 0 unless n1 = 0 *)
            else ratio2 := 1;
            for b3 := a to t do begin
               point3^.bas := b3;
               point2^.next := point3;
               n2 := getcount(root,point2);
               if ((n2 > 0) and (ratio2 < 1)) then
                  expected := exp(ln(n2) + ratio2)
               else expected := 0;
               observed := getcount(root,point1);
               chisq(observed,expected,chi,sign);
               write(anal,' ':3,basetochar(b1),basetochar(b2),basetochar(b3),
                  chi:8:2,' (',sign,');');
               totchi := totchi + chi;
            end;
            writeln(anal);
         end;
         writeln(anal);
      end;
      writeln(anal,' * total chisquared = ',totchi:8:2);
      writeln(anal);

         (* chi squareds for tets *)
         (* expected tets using dis is calculated by this formula:
            expected[b1,b2,b3,b4] = #[b1,b2]/2totals * #[b2,b3]/2totals *
               #[b3,b4]/2totals * 1totals/#[b2] * 1totals/#[b3] * 4totals *)

      if (compmax >= 4) then begin
         totchi := 0;
           (* ratio = ln( 1totals * 1totals * 4totals / (2totals)**3) <= 0  *)
         ratio := 2 * ln(monocomptotal^.count) +
            ln(monocomptotal^.next^.next^.next^.count) -
            3 * ln(monocomptotal^.next^.count);
         writeln(anal,' partials');

         new(point4);
         point4^.next := nil;
         for b1 := a to t do begin
            point1^.bas := b1;
            for b2 := a to t do begin
               point2^.bas := b2;
               point1^.next := point2;
               point2^.next := nil;
               n1 := getcount(root,point1);
               n2 := getcount(root,point2);
               if (n1 > 0) then
                  ratio2 := ratio + ln(n1) - ln(n2)
                     (* ratio2 <= 0 unless n1 = 0 *)
               else ratio2 := 1;
               for b3 := a to t do begin
                  point3^.bas := b3;
                  point2^.next := point3;
                  point3^.next := nil;
                  n2 := getcount(root,point2);
                  n3 := getcount(root,point3);
                  if ((n2 > 0) and (ratio2 < 1)) then
                     ratio3 := ratio2 + ln(n2) -  ln(n3)
                     (* ratio3 <= 0 unles n2 = 0 or n1 = 0 *)
                  else ratio3 := 1;
                  for b4 := a to t do begin
                     point4^.bas := b4;
                     point3^.next := point4;
                     n3 := getcount(root,point3);
                     if ((n3 > 0) and (ratio3 < 1)) then
                        expected := exp(ratio3 + ln(n3))
                     else expected := 0;
                     observed := getcount(root,point1);
                     chisq(observed,expected,chi,sign);
                     write(anal,' ':3,basetochar(b1),basetochar(b2),
                        basetochar(b3),basetochar(b4),chi:8:2,' (',sign,');');
                     totchi := totchi + chi;
                  end;
                  writeln(anal);
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal,' * total chisquared = ',totchi:8:2);
      end;
end;


procedure calc3; (* calculate chi squareds using tri frequencies *)
var b1,b2,b3,b4 : base;
      observed,  (* number of a particular oligo *)
      expected, (* number of a particular oligo *)
                (* expected tets are calculated from tris by this formula:
                      expected[b1,b2,b3,b4] = #[b1,b2,b3]/3totals *
                      #[b2,b3,b4]/3totals * 2totals/#[b2,b3] * 4totals  *)
      ratio,  (*   = 2totals * 4totals / (3totals)**2  *)
      ratio2,  (*   = ratio * #[b1,b2,b3] / #[b2,b3]  *)
      n1,n2, (* for storing values from 'getcount' calls *)
      chi,  (* partial chisquareds *)
      totchi : real;  (* total chisquareds *)
      sign : char; (* '+' or '-', as observed is more or less than expected *)
      point1,  (* on the path throught the tree *)
      point2,               (*   '   *)
      point3,               (*   '   *)
      point4: pathptr;      (*   '   *)

begin
      page(anal);
      writeln(anal,' **** chi squareds using tri frequenies ****');
      writeln(anal);

                        (* chi squareds for tets *)
      totchi := 0;
      ratio := ln(monocomptotal^.next^.count) +
         ln(monocomptotal^.next^.next^.next^.count) -
         2 * ln(monocomptotal^.next^.next^.count);
      writeln(anal,' partials');

      new(point1);
      new(point2);
      point1^.next := point2;
      new(point3);
      point2^.next := point3;
      new(point4);
      point4^.next := nil;
      for b1 := a to t do begin
         point1^.bas := b1;
         for b2 := a to t do begin
            point2^.bas := b2;
            for b3 := a to t do begin
               point3^.bas := b3;
               point3^.next := nil;
               n1 := getcount(root,point1);
               n2 := getcount(root,point2);
               if (n1 > 0) then
                  ratio2 := ratio + ln(n1) -   ln(n2)
                  (* ratio2 <= 0 unless n1 = 0 *)
               else ratio2 := 1;
               for b4 := a to t do begin
                  point4^.bas := b4;
                  point3^.next := point4;
                  n2 := getcount(root,point2);
                  if ((n2 > 0) and (ratio2 < 1)) then
                     expected :=  exp(ratio2 + ln(n2))
                  else expected := 0;
                  observed := getcount(root,point1);
                  chisq(observed,expected,chi,sign);
                  write(anal,' ':3,basetochar(b1),basetochar(b2),
                     basetochar(b3),basetochar(b4),chi:8:2,' (',sign,');');
                  totchi := totchi + chi;
               end;
               writeln(anal);
            end;
            writeln(anal);
         end;
         writeln(anal);
      end;
      writeln(anal,' * total chisquared = ',totchi:8:2);
end;

begin   (**** compan ****)
      writeln(output,'compan ',version:4:2);
      rewrite(anal);
      reset(cmp);
      if eof(cmp) then begin
         writeln(output,' no composition file provided');
         halt;
      end;
      reset(companp);
      if eof(companp) then readmax := defcompmax else readln(companp,readmax);
      if (readmax > defcompmax) then readmax := defcompmax;
      readcomp(cmp,compmax,readmax,root,monocomptotal);
      header;
      calc0;
      if (compmax >= 2) then calc1;
      if (compmax >= 3) then calc2;
      if (compmax >= 4) then calc3;
1: end.