source: trunk/GDE/PHYLIP/dnainvar.c

Last change on this file was 2175, checked in by westram, 22 years ago

upgrade to PHYLIP 3.6
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.0 KB
Line 
1
2#include "phylip.h"
3#include "seq.h"
4
5/* version 3.6. (c) Copyright 1993-2002 by the University of Washington.
6   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
7   Permission is granted to copy and use this program provided no fee is
8   charged for it and provided that this copyright notice is not removed. */
9
10#define maxsp           4   /* maximum number of species -- must be 4 */
11
12typedef enum {
13  xx, yy, zz, ww
14} simbol;
15
16#ifndef OLDC
17/* function prototypes */
18void getoptions(void);
19void allocrest(void);
20void doinit(void);
21void dnainvar_sitecombine(void);
22void makeweights(void);
23void doinput(void);
24void prntpatterns(void);
25void makesymmetries(void);
26
27void prntsymbol(simbol);
28void prntsymmetries(void);
29void tabulate(long,long,long,long,double *,double *,double *,double *);
30void dnainvar_writename(long);
31void writetree(long, long, long, long);
32void exacttest(long, long);
33void invariants(void);
34void makeinv(void);
35void reallocsites(void);
36/* function prototypes */
37#endif
38
39
40extern sequence y;
41
42Char infilename[FNMLNGTH], outfilename[FNMLNGTH], weightfilename[FNMLNGTH];
43long sites, msets, ith;
44boolean weights, progress, prntpat, printinv, mulsets, firstset, justwts;
45steptr aliasweight;
46
47long f[(long)ww - (long)xx + 1][(long)ww - (long)xx + 1]
48       [(long)ww - (long)xx + 1]; /* made global from being local to makeinv */
49
50void getoptions()
51{
52  /* interactively set options */
53  long loopcount, loopcount2;
54  boolean done;
55  Char ch, ch2;
56
57  fprintf(outfile, "\nNucleic acid sequence Invariants ");
58  fprintf(outfile, "method, version %s\n\n",VERSION);
59  putchar('\n');
60  printdata = false;
61  weights = false;
62  dotdiff = true;
63  progress = true;
64  prntpat = true;
65  printinv = true;
66  interleaved = true;
67  loopcount = 0;
68  do {
69    cleerhome();
70    printf("\nNucleic acid sequence Invariants ");
71    printf("method, version %s\n\n",VERSION);
72    printf("Settings for this run:\n");
73    printf("  W                       Sites weighted?  %s\n",
74           (weights ? "Yes" : "No"));
75    printf("  M           Analyze multiple data sets?");
76    if (mulsets)
77      printf("  Yes, %2ld %s\n", msets,
78              (justwts ? "sets of weights" : "data sets"));
79    else
80      printf("  No\n");
81    printf("  I          Input sequences interleaved?");
82    if (interleaved)
83      printf("  Yes\n");
84    else
85      printf("  No, sequential\n");
86    printf("  0   Terminal type (IBM PC, ANSI, none)?");
87    if (ibmpc)
88      printf("  IBM PC\n");
89    if (ansi)
90      printf("  ANSI\n");
91    if (!(ibmpc || ansi))
92      printf("  (none)\n");
93    printf("  1    Print out the data at start of run");
94    if (printdata)
95      printf("  Yes\n");
96    else
97      printf("  No\n");
98    if (printdata)
99      printf("  .  Use dot-differencing to display them  %s\n",
100           dotdiff ? "Yes" : "No");
101    printf("  2  Print indications of progress of run  %s\n",
102           (progress ? "Yes" : "No"));
103    printf("  3      Print out the counts of patterns");
104    if (prntpat)
105      printf("  Yes\n");
106    else
107      printf("  No\n");
108    printf("  4              Print out the invariants");
109    if (printinv)
110      printf("  Yes\n");
111    else
112        printf("  No\n");
113    if(weights && justwts){
114        printf(
115              "WARNING:  W option and Multiple Weights options are both on.  ");
116        printf(
117         "The W menu option is unnecessary and has no additional effect. \n");
118    }
119    printf("\n  Y to accept these or type the letter for one to change\n");
120#ifdef WIN32
121    phyFillScreenColor();
122#endif
123    scanf("%c%*[^\n]", &ch);
124    getchar();
125    uppercase(&ch);
126    done = (ch == 'Y');
127    if (!done) {
128      if (strchr("WMI01.234",ch) != NULL) {
129        switch (ch) {
130
131        case 'W':
132          weights = !weights;
133          break;
134
135        case 'M':
136          mulsets = !mulsets;
137          if (mulsets){
138            printf("Multiple data sets or multiple weights?");
139            loopcount2 = 0;
140            do {
141              printf(" (type D or W)\n");
142#ifdef WIN32
143              phyFillScreenColor();
144#endif
145              scanf("%c%*[^\n]", &ch2);
146              getchar();
147              if (ch2 == '\n')
148                ch2 = ' ';
149                uppercase(&ch2);
150                countup(&loopcount2, 10);
151            } while ((ch2 != 'W') && (ch2 != 'D'));
152            justwts = (ch2 == 'W');
153            if (justwts)
154              justweights(&msets);
155            else
156              initdatasets(&msets);
157            }
158            break;
159
160        case 'I':
161          interleaved = !interleaved;
162          break;
163
164        case '0':
165          initterminal(&ibmpc, &ansi);
166          break;
167
168        case '1':
169          printdata = !printdata;
170          break;
171
172        case '.':
173          dotdiff = !dotdiff;
174          break;
175
176        case '2':
177          progress = !progress;
178              break;
179       
180        case '3':
181          prntpat = !prntpat;
182          break;
183
184        case '4':
185          printinv = !printinv;
186          break;
187        }
188      } else
189        printf("Not a possible option!\n");
190    }
191    countup(&loopcount, 100);
192  } while (!done);
193}  /* getoptions */
194
195void reallocsites(void) 
196{
197  long i;
198
199  for (i=0;  i < spp; i++) {
200    free(y[i]);
201    y[i] = (Char *)Malloc(sites*sizeof(Char));
202  }
203  free(weight);
204  free(alias);
205  free(aliasweight);
206  weight       = (steptr)Malloc(sites * sizeof(long));
207  alias        = (steptr)Malloc(sites * sizeof(long));
208  aliasweight  = (steptr)Malloc(sites * sizeof(long));
209}
210
211
212void allocrest()
213{
214  long i;
215
216  y       = (Char **)Malloc(spp*sizeof(Char *));
217  for (i = 0; i < spp; i++)
218    y[i] = (Char *)Malloc(sites*sizeof(Char));
219  nayme        = (naym *)Malloc(maxsp * sizeof(naym));
220  weight       = (steptr)Malloc(sites * sizeof(long));
221  alias        = (steptr)Malloc(sites * sizeof(long));
222  aliasweight  = (steptr)Malloc(sites * sizeof(long));
223}
224
225
226void doinit()
227{
228  /* initializes variables */
229
230  inputnumbers(&spp, &sites, &nonodes, 1);
231  if (spp > maxsp){
232    printf("TOO MANY SPECIES: only 4 allowed\n");
233    exxit(-1);}
234  getoptions();
235  if (printdata)
236    fprintf(outfile, "%2ld species, %3ld  sites\n", spp, sites);
237  allocrest();
238}  /* doinit*/
239
240
241void dnainvar_sitecombine()
242{
243  /* combine sites that have identical patterns */
244  long i, j, k;
245  boolean tied;
246
247  i = 1;
248  while (i < sites) {
249    j = i + 1;
250    tied = true;
251    while (j <= sites && tied) {
252      k = 1;
253      while (k <= spp && tied) {
254        tied = (tied &&
255            y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
256        k++;
257      }
258      if (tied && aliasweight[j - 1] > 0) {
259        aliasweight[i - 1] += aliasweight[j - 1];
260        aliasweight[j - 1] = 0;
261      }
262      j++;
263    }
264    i = j - 1;
265  }
266}  /* dnainvar_sitecombine */
267
268
269void makeweights()
270{
271  /* make up weights vector to avoid duplicate computations */
272  long i;
273
274  for (i = 1; i <= sites; i++) {
275    alias[i - 1] = i;
276    aliasweight[i - 1] = weight[i - 1];
277  }
278  sitesort(sites, aliasweight);
279  dnainvar_sitecombine();
280  sitescrunch2(sites, 1, 2, aliasweight);
281  for (i = 1; i <= sites; i++) {
282    weight[i - 1] = aliasweight[i - 1];
283    if (weight[i - 1] > 0)
284      endsite = i;
285  }
286}  /* makeweights */
287
288
289void doinput()
290{
291  /* reads the input data */
292  long i;
293
294  if (justwts) {
295    if (firstset)
296      inputdata(sites);
297    for (i = 0; i < sites; i++)
298      weight[i] = 1;
299    inputweights(sites, weight, &weights);
300    if (justwts) {
301      fprintf(outfile, "\n\nWeights set # %ld:\n\n", ith);
302      if (progress)
303        printf("\nWeights set # %ld:\n\n", ith);
304    }
305    if (printdata)
306      printweights(outfile, 0, sites, weight, "Sites");
307  } else {
308    if (!firstset){
309      samenumsp(&sites, ith);
310      reallocsites();
311    }
312    inputdata(sites);
313    for (i = 0; i < sites; i++)
314      weight[i] = 1;
315    if (weights) {
316      inputweights(sites, weight, &weights);
317      if (printdata)
318        printweights(outfile, 0, sites, weight, "Sites");
319    }
320  }
321
322  makeweights();
323}  /* doinput */
324
325
326
327void prntpatterns()
328{
329  /* print out patterns */
330  long i, j;
331
332  fprintf(outfile, "\n   Pattern");
333  if (prntpat)
334    fprintf(outfile, "   Number of times");
335  fprintf(outfile, "\n\n");
336  for (i = 0; i < endsite; i++) {
337    fprintf(outfile, "     ");
338    for (j = 0; j < spp; j++)
339      putc(y[j][alias[i] - 1], outfile);
340    if (prntpat)
341      fprintf(outfile, "  %8ld", weight[i]);
342    putc('\n', outfile);
343  }
344  putc('\n', outfile);
345}  /* prntpatterns */
346
347
348void makesymmetries()
349{
350  /* get frequencies of symmetrized patterns */
351  long i, j;
352  boolean drop, usedz;
353  Char ch, ch1, zchar;
354  simbol s1, s2, s3;
355  simbol t[maxsp - 1];
356
357  for (s1 = xx; (long)s1 <= (long)ww; s1 = (simbol)((long)s1 + 1)) {
358    for (s2 = xx; (long)s2 <= (long)ww; s2 = (simbol)((long)s2 + 1)) {
359      for (s3 = xx; (long)s3 <= (long)ww; s3 = (simbol)((long)s3 + 1))
360        f[(long)s1 - (long)xx][(long)s2 - (long)xx]
361          [(long)s3 - (long)xx] = 0;
362    }
363  }
364  for (i = 0; i < endsite; i++) {
365    drop = false;
366    for (j = 0; j < spp; j++) {
367      ch = y[j][alias[i] - 1];
368      drop = (drop ||
369              (ch != 'A' && ch != 'C' && ch != 'G' && ch != 'T' && ch != 'U'));
370    }
371    ch1 = y[0][alias[i] - 1];
372    if (!drop) {
373      usedz = false;
374      zchar = ' ';
375      for (j = 2; j <= spp; j++) {
376        ch = y[j - 1][alias[i] - 1];
377        if (ch == ch1)
378          t[j - 2] = xx;
379        else if ((ch1 == 'A' && ch == 'G') || (ch1 == 'G' && ch == 'A') ||
380                 (ch1 == 'C' && (ch == 'T' || ch == 'U')) ||
381                 ((ch1 == 'T' || ch1 == 'U') && ch == 'C'))
382          t[j - 2] = yy;
383        else if (!usedz) {
384          t[j - 2] = zz;
385          usedz = true;
386          zchar = ch;
387        } else if (usedz && ch == zchar)
388          t[j - 2] = zz;
389        else if (usedz && ch != zchar)
390          t[j - 2] = ww;
391      }
392      f[(long)t[0] - (long)xx][(long)t[1] - (long)xx]
393        [(long)t[2] - (long)xx] += weight[i];
394    }
395  }
396}  /* makesymmetries */
397
398
399void prntsymbol(simbol s)
400{
401  /* print 1, 2, 3, 4 as appropriate */
402  switch (s) {
403
404  case xx:
405    putc('1', outfile);
406    break;
407
408  case yy:
409    putc('2', outfile);
410    break;
411
412  case zz:
413    putc('3', outfile);
414    break;
415
416  case ww:
417    putc('4', outfile);
418    break;
419  }
420}  /* prntsymbol */
421
422
423void prntsymmetries()
424{
425  /* print out symmetrized pattern numbers */
426  simbol s1, s2, s3;
427
428  fprintf(outfile, "\nSymmetrized patterns (1, 2 = the two purines  ");
429  fprintf(outfile, "and  3, 4 = the two pyrimidines\n");
430  fprintf(outfile, "                  or  1, 2 = the two pyrimidines  ");
431  fprintf(outfile, "and  3, 4 = the two purines)\n\n");
432  for (s1 = xx; (long)s1 <= (long)ww; s1 = (simbol)((long)s1 + 1)) {
433    for (s2 = xx; (long)s2 <= (long)ww; s2 = (simbol)((long)s2 + 1)) {
434      for (s3 = xx; (long)s3 <= (long)ww; s3 = (simbol)((long)s3 + 1)) {
435        if (f[(long)s1 - (long)xx][(long)s2 - (long)xx]
436            [(long)s3 - (long)xx] > 0) {
437          fprintf(outfile, "     1");
438          prntsymbol(s1);
439          prntsymbol(s2);
440          prntsymbol(s3);
441          if (prntpat)
442            fprintf(outfile, "   %7ld",
443                    f[(long)s1 - (long)xx][(long)s2 - (long)xx]
444                    [(long)s3 - (long)xx]);
445          putc('\n', outfile);
446        }
447      }
448    }
449  }
450}  /* prntsymmetries */
451
452
453void tabulate(long mm, long nn, long pp, long qq, double *mr,
454                        double *nr, double *pr, double *qr)
455{
456  /* make quadratic invariant, table, chi-square */
457  long total;
458  double k, TEMP;
459
460  fprintf(outfile, "\n   Contingency Table\n\n");
461  fprintf(outfile, "%7ld%6ld\n", mm, nn);
462  fprintf(outfile, "%7ld%6ld\n\n", pp, qq);
463  *mr = (long)(mm);
464  *nr = (long)(nn);
465  *pr = (long)pp;
466  *qr = (long)qq;
467  total = mm + nn + pp + qq;
468  if (printinv)
469    fprintf(outfile, "   Quadratic invariant = %15.1f\n\n",
470            (*nr) * (*pr) - (*mr) * (*qr));
471  fprintf(outfile, "   Chi-square = ");
472  TEMP = (*mr) * (*qr) - (*nr) * (*pr);
473  k = total * (TEMP * TEMP) / (((*mr) + (*nr)) * ((*mr) + (*pr)) *
474                               ((*nr) + (*qr)) * ((*pr) + (*qr)));
475  fprintf(outfile, "%10.5f", k);
476  if ((*mr) * (*qr) > (*nr) * (*pr) && k > 2.71)
477    fprintf(outfile, " (P < 0.05)\n");
478  else
479    fprintf(outfile, " (not significant)\n");
480  fprintf(outfile, "\n\n");
481}  /* tabulate */
482
483
484void dnainvar_writename(long m)
485{
486  /* write out a species name */
487  long i, n;
488
489  n = nmlngth;
490  while (nayme[m - 1][n - 1] == ' ')
491    n--;
492  if (n == 0)
493    n = 1;
494  for (i = 0; i < n; i++)
495    putc(nayme[m - 1][i], outfile);
496}  /* dnainvar_writename */
497
498
499void writetree(long i, long j, long k, long l)
500{
501  /* write out tree topology ((i,j),(k,l)) using names */
502  fprintf(outfile, "((");
503  dnainvar_writename(i);
504  putc(',', outfile);
505  dnainvar_writename(j);
506  fprintf(outfile, "),(");
507  dnainvar_writename(k);
508  putc(',', outfile);
509  dnainvar_writename(l);
510  fprintf(outfile, "))\n");
511}  /* writetree */
512
513
514void exacttest(long m, long n)
515{
516  /* exact binomial test that m <= n */
517  long i;
518  double p, sum;
519
520  p = 1.0;
521  for (i = 1; i <= m + n; i++)
522    p /= 2.0;
523  sum = p;
524  for (i = 1; i <= n; i++) {
525    p = p * (m + n - i + 1) / i;
526    sum += p;
527  }
528  fprintf(outfile, "      %7.4f", sum);
529  if (sum <= 0.05)
530    fprintf(outfile, "              yes\n");
531  else
532    fprintf(outfile, "               no\n");
533}  /* exacttest */
534
535
536void invariants()
537{
538  /* compute invariants */
539  long  m, n, p, q;
540  double L1, L2, L3;
541  double mr,nr,pr,qr;
542
543  fprintf(outfile, "\nTree topologies (unrooted): \n\n");
544  fprintf(outfile, "    I:  ");
545  writetree(1, 2, 3, 4);
546  fprintf(outfile, "   II:  ");
547  writetree(1, 3, 2, 4);
548  fprintf(outfile, "  III:  ");
549  writetree(1, 4, 2, 3);
550  fprintf(outfile, "\n\nLake's linear invariants\n");
551  fprintf(outfile,
552    " (these are expected to be zero for the two incorrect tree topologies.\n");
553  fprintf(outfile,
554          "  This is tested by testing the equality of the two parts\n");
555  fprintf(outfile,
556          "  of each expression using a one-sided exact binomial test.\n");
557  fprintf(outfile,
558    "  The null hypothesis is that the first part is no larger than the second.)\n\n");
559  fprintf(outfile, " Tree                           ");
560  fprintf(outfile, "  Exact test P value    Significant?\n\n");
561  m = f[(long)yy - (long)xx][(long)zz - (long)xx]
562      [(long)ww - (long)xx] + f[0][(long)zz - (long)xx]
563      [(long)zz - (long)xx];
564  n = f[(long)yy - (long)xx][(long)zz - (long)xx]
565      [(long)zz - (long)xx] + f[0][(long)zz - (long)xx]
566      [(long)ww - (long)xx];
567  fprintf(outfile, "   I  %5ld    - %5ld   = %5ld", m, n, m - n);
568  exacttest(m, n);
569  m = f[(long)zz - (long)xx][(long)yy - (long)xx]
570      [(long)ww - (long)xx] + f[(long)zz - (long)xx][0]
571      [(long)zz - (long)xx];
572  n = f[(long)zz - (long)xx][(long)yy - (long)xx]
573      [(long)zz - (long)xx] + f[(long)zz - (long)xx][0]
574      [(long)ww - (long)xx];
575  fprintf(outfile, "   II %5ld    - %5ld   = %5ld", m, n, m - n);
576  exacttest(m, n);
577  m = f[(long)zz - (long)xx][(long)ww - (long)xx]
578      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
579      [(long)zz - (long)xx][0];
580  n = f[(long)zz - (long)xx][(long)zz - (long)xx]
581      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
582      [(long)ww - (long)xx][0];
583  fprintf(outfile, "   III%5ld    - %5ld   = %5ld", m, n, m - n);
584  exacttest(m, n);
585  fprintf(outfile, "\n\nCavender's quadratic invariants (type L)");
586  fprintf(outfile, " using purines vs. pyrimidines\n");
587  fprintf(outfile,
588          " (these are expected to be zero, and thus have a nonsignificant\n");
589  fprintf(outfile, "  chi-square, for the correct tree topology)\n");
590  fprintf(outfile, "They will be misled if there are substantially\n");
591  fprintf(outfile, "different evolutionary rate between sites, or\n");
592  fprintf(outfile, "different purine:pyrimidine ratios from 1:1.\n\n");
593  fprintf(outfile, "  Tree I:\n");
594  m = f[0][0][0] + f[0][(long)yy - (long)xx]
595      [(long)yy - (long)xx] + f[0][(long)zz - (long)xx]
596      [(long)zz - (long)xx];
597  n = f[0][0][(long)yy - (long)xx] + f[0][0]
598      [(long)zz - (long)xx] + f[0][(long)yy - (long)xx][0] + f[0]
599      [(long)yy - (long)xx][(long)zz - (long)xx] + f[0]
600      [(long)zz - (long)xx][0] + f[0][(long)zz - (long)xx]
601      [(long)yy - (long)xx] + f[0][(long)zz - (long)xx]
602      [(long)ww - (long)xx];
603  p = f[(long)yy - (long)xx][0][0] + f[(long)yy - (long)xx]
604      [(long)yy - (long)xx]
605      [(long)yy - (long)xx] + f[(long)yy - (long)xx]
606      [(long)zz - (long)xx]
607      [(long)zz - (long)xx] + f[(long)zz - (long)xx][0]
608      [0] + f[(long)zz - (long)xx][(long)yy - (long)xx]
609      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
610      [(long)zz - (long)xx]
611      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
612      [(long)ww - (long)xx][(long)ww - (long)xx];
613  q = f[(long)yy - (long)xx][0][(long)yy - (long)xx] +
614      f[(long)yy - (long)xx][0][(long)zz - (long)xx] +
615      f[(long)yy - (long)xx][(long)yy - (long)xx][0] +
616      f[(long)yy - (long)xx][(long)yy - (long)xx][(long)zz - (long)xx] +
617      f[(long)yy - (long)xx][(long)zz - (long)xx][0] +
618      f[(long)yy - (long)xx][(long)zz - (long)xx][(long)yy - (long)xx] +
619      f[(long)yy - (long)xx][(long)zz - (long)xx][(long)ww - (long)xx] +
620      f[(long)zz - (long)xx][0][(long)yy - (long)xx] +
621      f[(long)zz - (long)xx][0][(long)zz - (long)xx] +
622      f[(long)zz - (long)xx][0][(long)ww - (long)xx] +
623      f[(long)zz - (long)xx][(long)yy - (long)xx][0] +
624      f[(long)zz - (long)xx][(long)yy - (long)xx][(long)zz - (long)xx] +
625      f[(long)zz - (long)xx][(long)yy - (long)xx][(long)ww - (long)xx] +
626      f[(long)zz - (long)xx][(long)zz - (long)xx][0] +
627      f[(long)zz - (long)xx][(long)zz - (long)xx][(long)yy - (long)xx] +
628      f[(long)zz - (long)xx][(long)zz - (long)xx][(long)ww - (long)xx] +
629      f[(long)zz - (long)xx][(long)ww - (long)xx][0] +
630      f[(long)zz - (long)xx][(long)ww - (long)xx][(long)yy - (long)xx] +
631      f[(long)zz - (long)xx][(long)ww - (long)xx][(long)zz - (long)xx];
632
633  nr = n;
634  pr = p;
635  mr = m;
636  qr = q;
637  L1 = nr * pr - mr * qr;
638  tabulate(m, n, p, q, &mr,&nr,&pr,&qr);
639  fprintf(outfile, "  Tree II:\n");
640  m = f[0][0][0] + f[(long)yy - (long)xx][0]
641      [(long)yy - (long)xx] + f[(long)zz - (long)xx][0]
642      [(long)zz - (long)xx];
643  n = f[0][0][(long)yy - (long)xx] + f[0][0]
644      [(long)zz - (long)xx] + f[(long)yy - (long)xx][0]
645      [0] + f[(long)yy - (long)xx][0]
646      [(long)zz - (long)xx] + f[(long)zz - (long)xx][0]
647      [0] + f[(long)zz - (long)xx][0]
648      [(long)yy - (long)xx] + f[(long)zz - (long)xx][0]
649      [(long)ww - (long)xx];
650  p = f[0][(long)yy - (long)xx][0] + f[(long)yy - (long)xx]
651      [(long)yy - (long)xx]
652      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
653      [(long)yy - (long)xx][(long)zz - (long)xx] + f[0]
654      [(long)zz - (long)xx][0] + f[(long)yy - (long)xx]
655      [(long)zz - (long)xx]
656      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
657      [(long)zz - (long)xx]
658      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
659      [(long)ww - (long)xx][(long)zz - (long)xx];
660  q = f[0][(long)yy - (long)xx][(long)yy - (long)xx] + f[0]
661      [(long)yy - (long)xx][(long)zz - (long)xx] +
662      f[(long)yy - (long)xx][(long)yy - (long)xx][0] +
663      f[(long)yy - (long)xx][(long)yy - (long)xx][(long)zz - (long)xx] +
664      f[(long)zz - (long)xx][(long)yy - (long)xx][0] +
665      f[(long)zz - (long)xx][(long)yy - (long)xx][(long)yy - (long)xx] +
666      f[(long)zz - (long)xx][(long)yy - (long)xx][(long)ww - (long)xx] +
667      f[0][(long)zz - (long)xx][(long)yy - (long)xx] + f[0]
668      [(long)zz - (long)xx][(long)zz - (long)xx] + f[0]
669      [(long)zz - (long)xx][(long)ww - (long)xx] +
670      f[(long)yy - (long)xx][(long)zz - (long)xx][0] +
671      f[(long)yy - (long)xx][(long)zz - (long)xx][(long)zz - (long)xx] +
672      f[(long)yy - (long)xx][(long)zz - (long)xx][(long)ww - (long)xx] +
673      f[(long)zz - (long)xx][(long)zz - (long)xx][0] +
674      f[(long)zz - (long)xx][(long)zz - (long)xx][(long)yy - (long)xx] +
675      f[(long)zz - (long)xx][(long)zz - (long)xx][(long)ww - (long)xx] +
676      f[(long)zz - (long)xx][(long)ww - (long)xx][0] +
677      f[(long)zz - (long)xx][(long)ww - (long)xx][(long)yy - (long)xx] +
678      f[(long)zz - (long)xx][(long)ww - (long)xx][(long)ww - (long)xx];
679  nr = n;
680  pr = p;
681  mr = m;
682  qr = q;
683  L2 = nr * pr - mr * qr;
684  tabulate(m, n, p, q, &mr,&nr,&pr,&qr);
685  fprintf(outfile, "  Tree III:\n");
686  m = f[0][0][0] + f[(long)yy - (long)xx][(long)yy - (long)xx]
687      [0] + f[(long)zz - (long)xx][(long)zz - (long)xx][0];
688  n = f[(long)yy - (long)xx][0][0] + f[(long)zz - (long)xx][0]
689      [0] + f[0][(long)yy - (long)xx][0] + f[(long)zz - (long)xx]
690      [(long)yy - (long)xx][0] + f[0][(long)zz - (long)xx]
691      [0] + f[(long)yy - (long)xx][(long)zz - (long)xx]
692      [0] + f[(long)zz - (long)xx][(long)ww - (long)xx][0];
693  p = f[0][0][(long)yy - (long)xx] + f[(long)yy - (long)xx]
694      [(long)yy - (long)xx]
695      [(long)yy - (long)xx] + f[(long)zz - (long)xx]
696      [(long)zz - (long)xx][(long)yy - (long)xx] + f[0][0]
697      [(long)zz - (long)xx] + f[(long)yy - (long)xx]
698      [(long)yy - (long)xx]
699      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
700      [(long)zz - (long)xx]
701      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
702      [(long)zz - (long)xx][(long)ww - (long)xx];
703  q = f[(long)yy - (long)xx][0][(long)yy - (long)xx] + f[(long)zz - (long)xx]
704      [0][(long)yy - (long)xx] + f[0][(long)yy - (long)xx][(long)yy - (long)xx] +
705      f[(long)zz - (long)xx][(long)yy - (long)xx][(long)yy - (long)xx] +
706      f[0][(long)zz - (long)xx]
707      [(long)yy - (long)xx] + f[(long)yy - (long)xx][(long)zz - (long)xx]
708      [(long)yy - (long)xx] + f[(long)zz - (long)xx][(long)ww - (long)xx]
709      [(long)yy - (long)xx] + f[(long)yy - (long)xx][0]
710      [(long)zz - (long)xx] + f[(long)zz - (long)xx][0]
711      [(long)zz - (long)xx] + f[0][(long)zz - (long)xx]
712      [(long)ww - (long)xx] + f[0][(long)yy - (long)xx]
713      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
714      [(long)yy - (long)xx]
715      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
716      [(long)yy - (long)xx][(long)ww - (long)xx] + f[0]
717      [(long)zz - (long)xx]
718      [(long)zz - (long)xx] + f[(long)yy - (long)xx]
719      [(long)zz - (long)xx]
720      [(long)zz - (long)xx] + f[(long)zz - (long)xx]
721      [(long)ww - (long)xx]
722      [(long)ww - (long)xx] + f[(long)zz - (long)xx][0]
723      [(long)ww - (long)xx] + f[(long)yy - (long)xx]
724      [(long)zz - (long)xx][(long)ww - (long)xx] +
725      f[(long)zz - (long)xx][(long)ww - (long)xx][(long)zz - (long)xx];
726  nr = n;
727  pr = p;
728  mr = m;
729  qr = q;
730  L3 = nr * pr - mr * qr;
731  tabulate(m, n, p, q, &mr,&nr,&pr,&qr);
732  fprintf(outfile, "\n\nCavender's quadratic invariants (type K)");
733  fprintf(outfile, " using purines vs. pyrimidines\n");
734  fprintf(outfile,
735          " (these are expected to be zero for the correct tree topology)\n");
736  fprintf(outfile, "They will be misled if there are substantially\n");
737  fprintf(outfile, "different evolutionary rate between sites, or\n");
738  fprintf(outfile, "different purine:pyrimidine ratios from 1:1.\n");
739  fprintf(outfile, "No statistical test is done on them here.\n\n");
740  fprintf(outfile, "  Tree I:   %15.1f\n", L2 - L3);
741  fprintf(outfile, "  Tree II:  %15.1f\n", L3 - L1);
742  fprintf(outfile, "  Tree III: %15.1f\n\n", L1 - L2);
743}  /* invariants */
744
745
746void makeinv()
747{
748  /* print out patterns and compute invariants */
749
750  prntpatterns();
751  makesymmetries();
752  prntsymmetries();
753  if (printinv)
754    invariants();
755}  /* makeinv */
756
757
758int main(int argc, Char *argv[])
759{  /* DNA Invariants */
760#ifdef MAC
761  argc = 1;                /* macsetup("Dnainvar","");                */
762  argv[0] = "Dnainvar";         
763#endif
764  init(argc,argv);
765  openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
766  openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
767
768  ibmpc = IBMCRT;
769  ansi = ANSICRT;
770  mulsets = false;
771  firstset = true;
772  msets = 1;
773  doinit();
774  if (weights || justwts)
775    openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
776  for (ith = 1; ith <= msets; ith++) {
777    doinput();
778    if (ith == 1)
779      firstset = false;
780    if (msets > 1 && !justwts) {
781      if (progress)
782        printf("\nData set # %ld:\n",ith);
783      fprintf(outfile, "Data set # %ld:\n\n",ith);
784    }
785    makeinv();
786  }
787  if (progress) {
788    putchar('\n');
789    printf("Output written to output file \"%s\"\n", outfilename);
790    putchar('\n');
791  }
792  FClose(outfile);
793  FClose(infile);
794#ifdef MAC
795  fixmacfile(outfilename);
796#endif
797  printf("Done.\n\n");
798#ifdef WIN32
799  phyRestoreConsoleAttributes();
800#endif
801  return 0;
802}  /* DNA Invariants */
Note: See TracBrowser for help on using the repository browser.