forester/archive/RIO/others/phylip_mod/src/protpars.c

   1
   2 #include "phylip.h"
   3 #include "seq.h"
   4
   5 /* version 3.6. (c) Copyright 1993-2004 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 maxtrees        100   /* maximum number of tied trees stored */
  11
  12 typedef enum {
  13   universal, ciliate, mito, vertmito, flymito, yeastmito
  14 } codetype;
  15
  16 /* nodes will form a binary tree */
  17
  18 typedef struct gseq {
  19   seqptr seq;
  20   struct gseq *next;
  21 } gseq;
  22
  23 #ifndef OLDC
  24 /* function prototypes */
  25 void   protgnu(gseq **);
  26 void   protchuck(gseq *);
  27 void   code(void);
  28 void   setup(void);
  29 void   getoptions(void);
  30 void   protalloctree(void);
  31 void   allocrest(void);
  32 void   doinit(void);
  33 void   protinputdata(void);
  34
  35 void   protmakevalues(void);
  36 void   doinput(void);
  37 void   protfillin(node *, node *, node *);
  38 void   protpreorder(node *);
  39 void   protadd(node *, node *, node *);
  40 void   protre_move(node **, node **);
  41 void   evaluate(node *);
  42 void   protpostorder(node *);
  43 void   protreroot(node *);
  44 void   protsavetraverse(node *, long *, boolean *);
  45
  46 void   protsavetree(long *, boolean *);
  47 void   tryadd(node *, node **, node **);
  48 void   addpreorder(node *, node *, node *);
  49 void   tryrearr(node *, boolean *);
  50 void   repreorder(node *, boolean *);
  51 void   rearrange(node **);
  52 void   protgetch(Char *);
  53 void   protaddelement(node **, long *, long *, boolean *);
  54 void   prottreeread(void);
  55 void   protancestset(long *, long *, long *, long *, long *);
  56
  57 void   prothyprint(long , long , boolean *, node *, boolean *, boolean *);
  58 void   prothyptrav(node *, sitearray *, long, long, long *, boolean *,
  59                 sitearray);
  60 void   prothypstates(long *);
  61 void   describe(void);
  62 void   maketree(void);
  63 void   reallocnode(node* p);
  64 void   reallocchars(void);
  65 /* function prototypes */
  66 #endif
  67
  68
  69
  70 Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], outtreename[FNMLNGTH], weightfilename[FNMLNGTH];
  71 node *root;
  72 long chars, col, msets, ith, njumble, jumb;
  73 /*   chars = number of sites in actual sequences */
  74 long inseed, inseed0;
  75 boolean jumble, usertree, weights, thresh, trout, progress, stepbox,
  76     justwts, ancseq, mulsets, firstset;
  77 codetype whichcode;
  78 long fullset, fulldel;
  79 pointarray treenode;   /* pointers to all nodes in tree */
  80 double threshold;
  81 steptr threshwt;
  82 longer seed;
  83 long *enterorder;
  84 sitearray translate[(long)quest - (long)ala + 1];
  85 aas trans[4][4][4];
  86 long **fsteps;
  87 bestelm *bestrees;
  88 boolean dummy;
  89 gseq *garbage;
  90 node *temp, *temp1;
  91 Char ch;
  92 aas tmpa;
  93 char *progname;
  94
  95 /* Local variables for maketree, propagated globally for c version: */
  96 long minwhich;
  97 double like, bestyet, bestlike, minsteps, bstlike2;
  98 boolean lastrearr, recompute;
  99 node *there;
 100 double nsteps[maxuser];
 101 long *place;
 102 boolean *names;
 103
 104
 105 void protgnu(gseq **p)
 106 {
 107   /* this and the following are do-it-yourself garbage collectors.
 108      Make a new node or pull one off the garbage list */
 109   if (garbage != NULL) {
 110     *p = garbage;
 111     free((*p)->seq);
 112     (*p)->seq = (seqptr)Malloc(chars*sizeof(sitearray));
 113     garbage = garbage->next;
 114   } else {
 115     *p = (gseq *)Malloc(sizeof(gseq));
 116     (*p)->seq = (seqptr)Malloc(chars*sizeof(sitearray));
 117   }
 118   (*p)->next = NULL;
 119 }  /* protgnu */
 120
 121
 122 void protchuck(gseq *p)
 123 {
 124   /* collect garbage on p -- put it on front of garbage list */
 125   p->next = garbage;
 126   garbage = p;
 127 }  /* protchuck */
 128
 129
 130 void code()
 131 {
 132   /* make up table of the code 1 = u, 2 = c, 3 = a, 4 = g */
 133   trans[0][0][0] = phe;
 134   trans[0][0][1] = phe;
 135   trans[0][0][2] = leu;
 136   trans[0][0][3] = leu;
 137   trans[0][1][0] = ser1;
 138   trans[0][1][1] = ser1;
 139   trans[0][1][2] = ser1;
 140   trans[0][1][3] = ser1;
 141   trans[0][2][0] = tyr;
 142   trans[0][2][1] = tyr;
 143   trans[0][2][2] = stop;
 144   trans[0][2][3] = stop;
 145   trans[0][3][0] = cys;
 146   trans[0][3][1] = cys;
 147   trans[0][3][2] = stop;
 148   trans[0][3][3] = trp;
 149   trans[1][0][0] = leu;
 150   trans[1][0][1] = leu;
 151   trans[1][0][2] = leu;
 152   trans[1][0][3] = leu;
 153   trans[1][1][0] = pro;
 154   trans[1][1][1] = pro;
 155   trans[1][1][2] = pro;
 156   trans[1][1][3] = pro;
 157   trans[1][2][0] = his;
 158   trans[1][2][1] = his;
 159   trans[1][2][2] = gln;
 160   trans[1][2][3] = gln;
 161   trans[1][3][0] = arg;
 162   trans[1][3][1] = arg;
 163   trans[1][3][2] = arg;
 164   trans[1][3][3] = arg;
 165   trans[2][0][0] = ileu;
 166   trans[2][0][1] = ileu;
 167   trans[2][0][2] = ileu;
 168   trans[2][0][3] = met;
 169   trans[2][1][0] = thr;
 170   trans[2][1][1] = thr;
 171   trans[2][1][2] = thr;
 172   trans[2][1][3] = thr;
 173   trans[2][2][0] = asn;
 174   trans[2][2][1] = asn;
 175   trans[2][2][2] = lys;
 176   trans[2][2][3] = lys;
 177   trans[2][3][0] = ser2;
 178   trans[2][3][1] = ser2;
 179   trans[2][3][2] = arg;
 180   trans[2][3][3] = arg;
 181   trans[3][0][0] = val;
 182   trans[3][0][1] = val;
 183   trans[3][0][2] = val;
 184   trans[3][0][3] = val;
 185   trans[3][1][0] = ala;
 186   trans[3][1][1] = ala;
 187   trans[3][1][2] = ala;
 188   trans[3][1][3] = ala;
 189   trans[3][2][0] = asp;
 190   trans[3][2][1] = asp;
 191   trans[3][2][2] = glu;
 192   trans[3][2][3] = glu;
 193   trans[3][3][0] = gly;
 194   trans[3][3][1] = gly;
 195   trans[3][3][2] = gly;
 196   trans[3][3][3] = gly;
 197   if (whichcode == mito)
 198     trans[0][3][2] = trp;
 199   if (whichcode == vertmito) {
 200     trans[0][3][2] = trp;
 201     trans[2][3][2] = stop;
 202     trans[2][3][3] = stop;
 203     trans[2][0][2] = met;
 204   }
 205   if (whichcode == flymito) {
 206     trans[0][3][2] = trp;
 207     trans[2][0][2] = met;
 208     trans[2][3][2] = ser2;
 209   }
 210   if (whichcode == yeastmito) {
 211     trans[0][3][2] = trp;
 212     trans[1][0][2] = thr;
 213     trans[2][0][2] = met;
 214   }
 215 } /* code */
 216
 217
 218 void setup()
 219 {
 220   /* set up set table to get aasets from aas */
 221   aas a, b;
 222   long i, j, k, l, s;
 223
 224   for (a = ala; (long)a <= (long)stop; a = (aas)((long)a + 1)) {
 225     translate[(long)a - (long)ala][0] = 1L << ((long)a);
 226     translate[(long)a - (long)ala][1] = 1L << ((long)a);
 227   }
 228   for (i = 0; i <= 3; i++) {
 229     for (j = 0; j <= 3; j++) {
 230       for (k = 0; k <= 3; k++) {
 231         for (l = 0; l <= 3; l++) {
 232           translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[l][j][k]);
 233           translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[i][l][k]);
 234           translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[i][j][l]);
 235         }
 236       }
 237     }
 238   }
 239   translate[(long)del - (long)ala][1] = 1L << ((long)del);
 240   fulldel = (1L << ((long)stop + 1)) - (1L << ((long)ala));
 241   fullset = fulldel & (~(1L << ((long)del)));
 242   translate[(long)asx - (long)ala][0]
 243     = (1L << ((long)asn)) | (1L << ((long)asp));
 244   translate[(long)glx - (long)ala][0]
 245     = (1L << ((long)gln)) | (1L << ((long)glu));
 246   translate[(long)ser - (long)ala][0]
 247     = (1L << ((long)ser1)) | (1L << ((long)ser2));
 248   translate[(long)unk - (long)ala][0] = fullset;
 249   translate[(long)quest - (long)ala][0] = fulldel;
 250   translate[(long)asx - (long)ala][1] = translate[(long)asn - (long)ala][1]
 251                                        | translate[(long)asp - (long)ala][1];
 252   translate[(long)glx - (long)ala][1] = translate[(long)gln - (long)ala][1]
 253                                        | translate[(long)glu - (long)ala][1];
 254   translate[(long)ser - (long)ala][1] = translate[(long)ser1 - (long)ala][1]
 255                                        | translate[(long)ser2 - (long)ala][1];
 256   translate[(long)unk - (long)ala][1] = fullset;
 257   translate[(long)quest - (long)ala][1] = fulldel;
 258   for (a = ala; (long)a <= (long)quest; a = (aas)((long)a + 1)) {
 259     s = 0;
 260     for (b = ala; (long)b <= (long)stop; b = (aas)((long)b + 1)) {
 261       if (((1L << ((long)b)) & translate[(long)a - (long)ala][1]) != 0)
 262         s |= translate[(long)b - (long)ala][1];
 263     }
 264     translate[(long)a - (long)ala][2] = s;
 265   }
 266 }  /* setup */
 267
 268
 269 void getoptions()
 270 {
 271   /* interactively set options */
 272   long loopcount, loopcount2;
 273   Char ch, ch2;
 274
 275   fprintf(outfile, "\nProtein parsimony algorithm, version %s\n\n",VERSION);
 276   putchar('\n');
 277   jumble = false;
 278   njumble = 1;
 279   outgrno = 1;
 280   outgropt = false;
 281   thresh = false;
 282   trout = true;
 283   usertree = false;
 284   weights = false;
 285   whichcode = universal;
 286   printdata = false;
 287   progress = true;
 288   treeprint = true;
 289   stepbox = false;
 290   ancseq = false;
 291   dotdiff = true;
 292   interleaved = true;
 293   loopcount = 0;
 294   for (;;) {
 295     cleerhome();
 296     printf("\nProtein parsimony algorithm, version %s\n\n",VERSION);
 297     printf("Setting for this run:\n");
 298     printf("  U                 Search for best tree?  %s\n",
 299            (usertree ? "No, use user trees in input file" : "Yes"));
 300     if (!usertree) {
 301       printf("  J   Randomize input order of sequences?");
 302       if (jumble)
 303         printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
 304       else
 305         printf("  No. Use input order\n");
 306     }
 307     printf("  O                        Outgroup root?");
 308     if (outgropt)
 309       printf("  Yes, at sequence number%3ld\n", outgrno);
 310     else
 311       printf("  No, use as outgroup species%3ld\n", outgrno);
 312     printf("  T              Use Threshold parsimony?");
 313     if (thresh)
 314       printf("  Yes, count steps up to%4.1f per site\n", threshold);
 315     else
 316       printf("  No, use ordinary parsimony\n");
 317     printf("  C               Use which genetic code?  %s\n",
 318       (whichcode == universal) ? "Universal"                  :
 319       (whichcode == ciliate)   ? "Ciliate"                    :
 320       (whichcode == mito)      ? "Universal mitochondrial"    :
 321       (whichcode == vertmito)  ? "Vertebrate mitochondrial"   :
 322       (whichcode == flymito)   ? "Fly mitochondrial"          :
 323       (whichcode == yeastmito) ? "Yeast mitochondrial"        : "");
 324     printf("  W                       Sites weighted?  %s\n",
 325            (weights ? "Yes" : "No"));
 326     printf("  M           Analyze multiple data sets?");
 327     if (mulsets)
 328         printf("  Yes, %2ld %s\n", msets,
 329                (justwts ? "sets of weights" : "data sets"));
 330     else
 331       printf("  No\n");
 332     printf("  I          Input sequences interleaved?  %s\n",
 333            (interleaved ? "Yes" : "No"));
 334     printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
 335            (ibmpc ? "IBM PC" : ansi ? "ANSI" : "(none)"));
 336     printf("  1    Print out the data at start of run  %s\n",
 337            (printdata ? "Yes" : "No"));
 338     printf("  2  Print indications of progress of run  %s\n",
 339            (progress ? "Yes" : "No"));
 340     printf("  3                        Print out tree  %s\n",
 341            (treeprint ? "Yes" : "No"));
 342     printf("  4          Print out steps in each site  %s\n",
 343            (stepbox ? "Yes" : "No"));
 344     printf("  5  Print sequences at all nodes of tree  %s\n",
 345            (ancseq ? "Yes" : "No"));
 346     if (ancseq || printdata)
 347       printf("  .  Use dot-differencing to display them  %s\n",
 348               dotdiff ? "Yes" : "No");
 349     printf("  6       Write out trees onto tree file?  %s\n",
 350            (trout ? "Yes" : "No"));
 351     if(weights && justwts){
 352         printf(
 353          "WARNING:  W option and Multiple Weights options are both on.  ");
 354         printf(
 355          "The W menu option is unnecessary and has no additional effect. \n");
 356     }
 357     printf(
 358    "\nAre these settings correct? (type Y or the letter for one to change)\n");
 359     scanf("%c%*[^\n]", &ch);
 360     getchar();
 361     uppercase(&ch);
 362     if (ch == 'Y')
 363       break;
 364     if (strchr("WCJOTUMI12345.60",ch) != NULL) {
 365       switch (ch) {
 366
 367       case 'J':
 368         jumble = !jumble;
 369         if (jumble)
 370           initjumble(&inseed, &inseed0, seed, &njumble);
 371         else njumble = 1;
 372         break;
 373
 374       case 'W':
 375         weights = !weights;
 376         break;
 377
 378       case 'O':
 379         outgropt = !outgropt;
 380         if (outgropt)
 381           initoutgroup(&outgrno, spp);
 382         else outgrno = 1;
 383         break;
 384
 385       case 'T':
 386         thresh = !thresh;
 387         if (thresh)
 388           initthreshold(&threshold);
 389         break;
 390
 391       case 'C':
 392         printf("\nWhich genetic code?\n");
 393         printf(" type         for\n\n");
 394         printf("   U           Universal\n");
 395         printf("   M           Mitochondrial\n");
 396         printf("   V           Vertebrate mitochondrial\n");
 397         printf("   F           Fly mitochondrial\n");
 398         printf("   Y           Yeast mitochondrial\n\n");
 399         loopcount2 = 0;
 400         do {
 401           printf("type U, M, V, F, or Y\n");
 402           scanf("%c%*[^\n]", &ch);
 403           getchar();
 404           if (ch == '\n')
 405             ch = ' ';
 406           uppercase(&ch);
 407           countup(&loopcount2, 10);
 408         } while (ch != 'U' && ch != 'M' && ch != 'V'
 409                   && ch != 'F' && ch != 'Y');
 410         switch (ch) {
 411
 412         case 'U':
 413           whichcode = universal;
 414           break;
 415
 416         case 'M':
 417           whichcode = mito;
 418           break;
 419
 420         case 'V':
 421           whichcode = vertmito;
 422           break;
 423
 424         case 'F':
 425           whichcode = flymito;
 426           break;
 427
 428         case 'Y':
 429           whichcode = yeastmito;
 430           break;
 431         }
 432         break;
 433
 434       case 'M':
 435         mulsets = !mulsets;
 436         if (mulsets){
 437             printf("Multiple data sets or multiple weights?");
 438           loopcount2 = 0;
 439           do {
 440             printf(" (type D or W)\n");
 441 #ifdef WIN32
 442             phyFillScreenColor();
 443 #endif
 444             scanf("%c%*[^\n]", &ch2);
 445             getchar();
 446             if (ch2 == '\n')
 447               ch2 = ' ';
 448             uppercase(&ch2);
 449             countup(&loopcount2, 10);
 450           } while ((ch2 != 'W') && (ch2 != 'D'));
 451           justwts = (ch2 == 'W');
 452           if (justwts)
 453             justweights(&msets);
 454           else
 455             initdatasets(&msets);
 456           if (!jumble) {
 457             jumble = true;
 458             initjumble(&inseed, &inseed0, seed, &njumble);
 459           }
 460         }
 461         break;
 462
 463       case 'I':
 464         interleaved = !interleaved;
 465         break;
 466
 467       case 'U':
 468         usertree = !usertree;
 469         break;
 470
 471       case '0':
 472         initterminal(&ibmpc, &ansi);
 473         break;
 474
 475       case '1':
 476         printdata = !printdata;
 477         break;
 478
 479       case '2':
 480         progress = !progress;
 481         break;
 482
 483       case '3':
 484         treeprint = !treeprint;
 485         break;
 486
 487       case '4':
 488         stepbox = !stepbox;
 489         break;
 490
 491       case '5':
 492         ancseq = !ancseq;
 493         break;
 494
 495       case '.':
 496         dotdiff = !dotdiff;
 497         break;
 498
 499       case '6':
 500         trout = !trout;
 501         break;
 502       }
 503     } else
 504         printf("Not a possible option!\n");
 505     countup(&loopcount, 100);
 506   }
 507 }  /* getoptions */
 508
 509
 510 void protalloctree()
 511 { /* allocate treenode dynamically */
 512   long i, j;
 513   node *p, *q;
 514
 515   treenode = (pointarray)Malloc(nonodes*sizeof(node *));
 516   for (i = 0; i < (spp); i++) {
 517     treenode[i] = (node *)Malloc(sizeof(node));
 518     treenode[i]->numsteps = (steptr)Malloc(chars*sizeof(long));
 519     treenode[i]->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 520     treenode[i]->seq = (aas *)Malloc(chars*sizeof(aas));
 521   }
 522   for (i = spp; i < (nonodes); i++) {
 523     q = NULL;
 524     for (j = 1; j <= 3; j++) {
 525       p = (node *)Malloc(sizeof(node));
 526       p->numsteps = (steptr)Malloc(chars*sizeof(long));
 527       p->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 528       p->seq = (aas *)Malloc(chars*sizeof(aas));
 529       p->next = q;
 530       q = p;
 531     }
 532     p->next->next->next = p;
 533     treenode[i] = p;
 534   }
 535 }  /* protalloctree */
 536
 537
 538 void reallocnode(node* p)
 539 {
 540   free(p->numsteps);
 541   free(p->siteset);
 542   free(p->seq);
 543   p->numsteps = (steptr)Malloc(chars*sizeof(long));
 544   p->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 545   p->seq = (aas *)Malloc(chars*sizeof(aas));
 546 }
 547
 548
 549 void reallocchars(void)
 550 { /* reallocates variables that are dependand on the number of chars
 551    * do we need to reallocate the garbage list too? */
 552   long i;
 553   node *p;
 554
 555   if (usertree)
 556     for (i = 0; i < maxuser; i++) {
 557       free(fsteps[i]);
 558       fsteps[i] = (long *)Malloc(chars*sizeof(long));
 559     }
 560
 561   for (i = 0; i < nonodes; i++) {
 562     reallocnode(treenode[i]);
 563     if (i >= spp) {
 564       p=treenode[i]->next;
 565       while (p != treenode[i])  {
 566         reallocnode(p);
 567         p = p->next;
 568       }
 569     }
 570   }
 571
 572   free(weight);
 573   free(threshwt);
 574   free(temp->numsteps);
 575   free(temp->siteset);
 576   free(temp->seq);
 577   free(temp1->numsteps);
 578   free(temp1->siteset);
 579   free(temp1->seq);
 580
 581   weight = (steptr)Malloc(chars*sizeof(long));
 582   threshwt = (steptr)Malloc(chars*sizeof(long));
 583   temp->numsteps = (steptr)Malloc(chars*sizeof(long));
 584   temp->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 585   temp->seq = (aas *)Malloc(chars*sizeof(aas));
 586   temp1->numsteps = (steptr)Malloc(chars*sizeof(long));
 587   temp1->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 588   temp1->seq = (aas *)Malloc(chars*sizeof(aas));
 589 }
 590
 591
 592 void allocrest()
 593 { /* allocate remaining global arrays and variables dynamically */
 594   long i;
 595
 596   if (usertree) {
 597     fsteps = (long **)Malloc(maxuser*sizeof(long *));
 598     for (i = 0; i < maxuser; i++)
 599       fsteps[i] = (long *)Malloc(chars*sizeof(long));
 600   }
 601   bestrees = (bestelm *)Malloc(maxtrees*sizeof(bestelm));
 602   for (i = 1; i <= maxtrees; i++)
 603     bestrees[i - 1].btree = (long *)Malloc(spp*sizeof(long));
 604   nayme = (naym *)Malloc(spp*sizeof(naym));
 605   enterorder = (long *)Malloc(spp*sizeof(long));
 606   place = (long *)Malloc(nonodes*sizeof(long));
 607   weight = (steptr)Malloc(chars*sizeof(long));
 608   threshwt = (steptr)Malloc(chars*sizeof(long));
 609   temp = (node *)Malloc(sizeof(node));
 610   temp->numsteps = (steptr)Malloc(chars*sizeof(long));
 611   temp->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 612   temp->seq = (aas *)Malloc(chars*sizeof(aas));
 613   temp1 = (node *)Malloc(sizeof(node));
 614   temp1->numsteps = (steptr)Malloc(chars*sizeof(long));
 615   temp1->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
 616   temp1->seq = (aas *)Malloc(chars*sizeof(aas));
 617 }  /* allocrest */
 618
 619
 620 void doinit()
 621 {
 622   /* initializes variables */
 623
 624   inputnumbers(&spp, &chars, &nonodes, 1);
 625   getoptions();
 626   if (printdata)
 627     fprintf(outfile, "%2ld species, %3ld  sites\n\n", spp, chars);
 628   protalloctree();
 629   allocrest();
 630 }  /* doinit*/
 631
 632
 633 void protinputdata()
 634 {
 635   /* input the names and sequences for each species */
 636   long i, j, k, l, aasread, aasnew = 0;
 637   Char charstate;
 638   boolean allread, done;
 639   aas aa;   /* temporary amino acid for input */
 640
 641   if (printdata)
 642     headings(chars, "Sequences", "---------");
 643   aasread = 0;
 644   allread = false;
 645   while (!(allread)) {
 646     /* eat white space -- if the separator line has spaces on it*/
 647     do {
 648       charstate = gettc(infile);
 649     } while (charstate == ' ' || charstate == '\t');
 650     ungetc(charstate, infile);
 651     if (eoln(infile)) {
 652       scan_eoln(infile);
 653     }
 654     i = 1;
 655     while (i <= spp) {
 656       if ((interleaved && aasread == 0) || !interleaved)
 657         initname(i - 1);
 658       j = interleaved ? aasread : 0;
 659       done = false;
 660       while (!done && !eoff(infile)) {
 661         if (interleaved)
 662           done = true;
 663         while (j < chars && !(eoln(infile) || eoff(infile))) {
 664           charstate = gettc(infile);
 665           if (charstate == '\n' || charstate == '\t')
 666             charstate = ' ';
 667           if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
 668             continue;
 669           uppercase(&charstate);
 670           if ((!isalpha(charstate) && charstate != '?' &&
 671                charstate != '-' && charstate != '*') || charstate == 'J' ||
 672               charstate == 'O' || charstate == 'U') {
 673             printf("WARNING -- BAD AMINO ACID:%c",charstate);
 674             printf(" AT POSITION%5ld OF SPECIES %3ld\n",j,i);
 675             exxit(-1);
 676           }
 677           j++;
 678           aa = (charstate == 'A') ?  ala :
 679                (charstate == 'B') ?  asx :
 680                (charstate == 'C') ?  cys :
 681                (charstate == 'D') ?  asp :
 682                (charstate == 'E') ?  glu :
 683                (charstate == 'F') ?  phe :
 684                (charstate == 'G') ?  gly : aa;
 685           aa = (charstate == 'H') ?  his :
 686                (charstate == 'I') ? ileu :
 687                (charstate == 'K') ?  lys :
 688                (charstate == 'L') ?  leu :
 689                (charstate == 'M') ?  met :
 690                (charstate == 'N') ?  asn :
 691                (charstate == 'P') ?  pro :
 692                (charstate == 'Q') ?  gln :
 693                (charstate == 'R') ?  arg : aa;
 694           aa = (charstate == 'S') ?  ser :
 695                (charstate == 'T') ?  thr :
 696                (charstate == 'V') ?  val :
 697                (charstate == 'W') ?  trp :
 698                (charstate == 'X') ?  unk :
 699                (charstate == 'Y') ?  tyr :
 700                (charstate == 'Z') ?  glx :
 701                (charstate == '*') ? stop :
 702                (charstate == '?') ? quest:
 703                (charstate == '-') ? del  :  aa;
 704
 705           treenode[i - 1]->seq[j - 1] = aa;
 706           memcpy(treenode[i - 1]->siteset[j - 1],
 707                  translate[(long)aa - (long)ala], sizeof(sitearray));
 708         }
 709         if (interleaved)
 710           continue;
 711         if (j < chars)
 712           scan_eoln(infile);
 713         else if (j == chars)
 714           done = true;
 715       }
 716       if (interleaved && i == 1)
 717         aasnew = j;
 718       scan_eoln(infile);
 719       if ((interleaved && j != aasnew) || ((!interleaved) && j != chars)){
 720         printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
 721         exxit(-1);}
 722       i++;
 723     }
 724     if (interleaved) {
 725       aasread = aasnew;
 726       allread = (aasread == chars);
 727     } else
 728       allread = (i > spp);
 729   }
 730   if (printdata) {
 731     for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
 732       for (j = 1; j <= (spp); j++) {
 733         for (k = 0; k < nmlngth; k++)
 734           putc(nayme[j - 1][k], outfile);
 735         fprintf(outfile, "   ");
 736         l = i * 60;
 737         if (l > chars)
 738           l = chars;
 739         for (k = (i - 1) * 60 + 1; k <= l; k++) {
 740           if (j > 1 && treenode[j - 1]->seq[k - 1] == treenode[0]->seq[k - 1])
 741             charstate = '.';
 742           else {
 743             tmpa = treenode[j-1]->seq[k-1];
 744               charstate =  (tmpa == ala) ? 'A' :
 745                            (tmpa == asx) ? 'B' :
 746                            (tmpa == cys) ? 'C' :
 747                            (tmpa == asp) ? 'D' :
 748                            (tmpa == glu) ? 'E' :
 749                            (tmpa == phe) ? 'F' :
 750                            (tmpa == gly) ? 'G' :
 751                            (tmpa == his) ? 'H' :
 752                            (tmpa ==ileu) ? 'I' :
 753                            (tmpa == lys) ? 'K' :
 754                            (tmpa == leu) ? 'L' : charstate;
 755               charstate =  (tmpa == met) ? 'M' :
 756                            (tmpa == asn) ? 'N' :
 757                            (tmpa == pro) ? 'P' :
 758                            (tmpa == gln) ? 'Q' :
 759                            (tmpa == arg) ? 'R' :
 760                            (tmpa == ser) ? 'S' :
 761                            (tmpa ==ser1) ? 'S' :
 762                            (tmpa ==ser2) ? 'S' : charstate;
 763               charstate =  (tmpa == thr) ? 'T' :
 764                            (tmpa == val) ? 'V' :
 765                            (tmpa == trp) ? 'W' :
 766                            (tmpa == unk) ? 'X' :
 767                            (tmpa == tyr) ? 'Y' :
 768                            (tmpa == glx) ? 'Z' :
 769                            (tmpa == del) ? '-' :
 770                            (tmpa ==stop) ? '*' :
 771                            (tmpa==quest) ? '?' : charstate;
 772         }
 773           putc(charstate, outfile);
 774           if (k % 10 == 0 && k % 60 != 0)
 775             putc(' ', outfile);
 776         }
 777         putc('\n', outfile);
 778       }
 779       putc('\n', outfile);
 780     }
 781     putc('\n', outfile);
 782   }
 783   putc('\n', outfile);
 784 }  /* protinputdata */
 785
 786
 787 void protmakevalues()
 788 {
 789   /* set up fractional likelihoods at tips */
 790   long i, j;
 791   node *p;
 792
 793   for (i = 1; i <= nonodes; i++) {
 794     treenode[i - 1]->back = NULL;
 795     treenode[i - 1]->tip = (i <= spp);
 796     treenode[i - 1]->index = i;
 797     for (j = 0; j < (chars); j++)
 798       treenode[i - 1]->numsteps[j] = 0;
 799     if (i > spp) {
 800       p = treenode[i - 1]->next;
 801       while (p != treenode[i - 1]) {
 802         p->back = NULL;
 803         p->tip = false;
 804         p->index = i;
 805         for (j = 0; j < (chars); j++)
 806           p->numsteps[j] = 0;
 807         p = p->next;
 808       }
 809     }
 810   }
 811 }  /* protmakevalues */
 812
 813
 814 void doinput()
 815 {
 816   /* reads the input data */
 817   long i;
 818
 819   if (justwts) {
 820     if (firstset)
 821       protinputdata();
 822     for (i = 0; i < chars; i++)
 823       weight[i] = 1;
 824     inputweights(chars, weight, &weights);
 825     if (justwts) {
 826       fprintf(outfile, "\n\nWeights set # %ld:\n\n", ith);
 827       if (progress)
 828         printf("\nWeights set # %ld:\n\n", ith);
 829     }
 830     if (printdata)
 831       printweights(outfile, 0, chars, weight, "Sites");
 832   } else {
 833     if (!firstset){
 834       samenumsp(&chars, ith);
 835       reallocchars();
 836     }
 837     for (i = 0; i < chars; i++)
 838       weight[i] = 1;
 839     if (weights) {
 840       inputweights(chars, weight, &weights);
 841     }
 842     if (weights)
 843       printweights(outfile, 0, chars, weight, "Sites");
 844     protinputdata();
 845   }
 846   if(!thresh)
 847     threshold = spp * 3.0;
 848   for(i = 0 ; i < (chars) ; i++){
 849     weight[i]*=10;
 850     threshwt[i] = (long)(threshold * weight[i] + 0.5);
 851   }
 852
 853   protmakevalues();
 854 }  /* doinput */
 855
 856
 857 void protfillin(node *p, node *left, node *rt)
 858 {
 859   /* sets up for each node in the tree the aa set for site m
 860      at that point and counts the changes.  The program
 861      spends much of its time in this function */
 862   boolean counted, done;
 863   aas aa;
 864   long s = 0;
 865   sitearray ls, rs, qs;
 866   long i, j, m, n;
 867
 868   for (m = 0; m < chars; m++) {
 869     if (left != NULL)
 870       memcpy(ls, left->siteset[m], sizeof(sitearray));
 871     if (rt != NULL)
 872       memcpy(rs, rt->siteset[m], sizeof(sitearray));
 873     if (left == NULL) {
 874       n = rt->numsteps[m];
 875       memcpy(qs, rs, sizeof(sitearray));
 876     }
 877     else if (rt == NULL) {
 878       n = left->numsteps[m];
 879       memcpy(qs, ls, sizeof(sitearray));
 880     }
 881     else {
 882       n = left->numsteps[m] + rt->numsteps[m];
 883       if ((ls[0] == rs[0]) && (ls[1] == rs[1]) && (ls[2] == rs[2])) {
 884         qs[0] = ls[0];
 885         qs[1] = ls[1];
 886         qs[2] = ls[2];
 887       }
 888       else {
 889         counted = false;
 890         for (i = 0; (!counted) && (i <= 3); i++) {
 891           switch (i) {
 892
 893             case 0:
 894               s = ls[0] & rs[0];
 895               break;
 896
 897             case 1:
 898               s = (ls[0] & rs[1]) | (ls[1] & rs[0]);
 899               break;
 900
 901             case 2:
 902               s = (ls[0] & rs[2]) | (ls[1] & rs[1]) | (ls[2] & rs[0]);
 903               break;
 904
 905             case 3:
 906               s = ls[0] | (ls[1] & rs[2]) | (ls[2] & rs[1]) | rs[0];
 907               break;
 908
 909           }
 910           if (s != 0) {
 911             qs[0] = s;
 912             counted = true;
 913           } else
 914               n += weight[m];
 915         }
 916         switch (i) {
 917           case 1:
 918             qs[1] = qs[0] | (ls[0] & rs[1]) | (ls[1] & rs[0]);
 919             qs[2] = qs[1] | (ls[0] & rs[2]) | (ls[1] & rs[1]) | (ls[2] & rs[0]);
 920             break;
 921           case 2:
 922             qs[1] = qs[0] | (ls[0] & rs[2]) | (ls[1] & rs[1]) | (ls[2] & rs[0]);
 923             qs[2] = qs[1] | ls[0] | (ls[1] & rs[2]) | (ls[2] & rs[1]) | rs[0];
 924             break;
 925           case 3:
 926             qs[1] = qs[0] | ls[0] | (ls[1] & rs[2]) | (ls[2] & rs[1]) | rs[0];
 927             qs[2] = qs[1] | ls[1] | (ls[2] & rs[2]) | rs[1];
 928             break;
 929           case 4:
 930             qs[1] = qs[0] | ls[1] | (ls[2] & rs[2]) | rs[1];
 931             qs[2] = qs[1] | ls[2] | rs[2];
 932             break;
 933         }
 934         for (aa = ala; (long)aa <= (long)stop; aa = (aas)((long)aa + 1)) {
 935           done = false;
 936           for (i = 0; (!done) && (i <= 1); i++) {
 937             if (((1L << ((long)aa)) & qs[i]) != 0) {
 938               for (j = i+1; j <= 2; j++)
 939                 qs[j] |= translate[(long)aa - (long)ala][j-i];
 940               done = true;
 941             }
 942           }
 943         }
 944       }
 945     }
 946     p->numsteps[m] = n;
 947     memcpy(p->siteset[m], qs, sizeof(sitearray));
 948   }
 949 }  /* protfillin */
 950
 951
 952 void protpreorder(node *p)
 953 {
 954   /* recompute number of steps in preorder taking both ancestoral and
 955      descendent steps into account */
 956   if (p != NULL && !p->tip) {
 957     protfillin (p->next, p->next->next->back, p->back);
 958     protfillin (p->next->next, p->back, p->next->back);
 959     protpreorder (p->next->back);
 960     protpreorder (p->next->next->back);
 961   }
 962 } /* protpreorder */
 963
 964
 965 void protadd(node *below, node *newtip, node *newfork)
 966 {
 967   /* inserts the nodes newfork and its left descendant, newtip,
 968      to the tree.  below becomes newfork's right descendant */
 969
 970   if (below != treenode[below->index - 1])
 971     below = treenode[below->index - 1];
 972   if (below->back != NULL)
 973     below->back->back = newfork;
 974   newfork->back = below->back;
 975   below->back = newfork->next->next;
 976   newfork->next->next->back = below;
 977   newfork->next->back = newtip;
 978   newtip->back = newfork->next;
 979   if (root == below)
 980     root = newfork;
 981   root->back = NULL;
 982
 983   if (recompute) {
 984     protfillin (newfork, newfork->next->back, newfork->next->next->back);
 985     protpreorder(newfork);
 986     if (newfork != root)
 987       protpreorder(newfork->back);
 988   }
 989 }  /* protadd */
 990
 991
 992 void protre_move(node **item, node **fork)
 993 {
 994   /* removes nodes item and its ancestor, fork, from the tree.
 995      the new descendant of fork's ancestor is made to be
 996      fork's second descendant (other than item).  Also
 997      returns pointers to the deleted nodes, item and fork */
 998   node *p, *q, *other;
 999
1000   if ((*item)->back == NULL) {
1001     *fork = NULL;
1002     return;
1003   }
1004   *fork = treenode[(*item)->back->index - 1];
1005   if ((*item) == (*fork)->next->back)
1006     other = (*fork)->next->next->back;
1007   else other = (*fork)->next->back;
1008   if (root == *fork)
1009     root = other;
1010   p = (*item)->back->next->back;
1011   q = (*item)->back->next->next->back;
1012   if (p != NULL) p->back = q;
1013   if (q != NULL) q->back = p;
1014   (*fork)->back = NULL;
1015   p = (*fork)->next;
1016   do {
1017     p->back = NULL;
1018     p = p->next;
1019   } while (p != (*fork));
1020   (*item)->back = NULL;
1021   if (recompute) {
1022     protpreorder(other);
1023     if (other != root) protpreorder(other->back);
1024   }
1025 }  /* protre_move */
1026
1027
1028 void evaluate(node *r)
1029 {
1030   /* determines the number of steps needed for a tree. this is the
1031      minimum number of steps needed to evolve sequences on this tree */
1032   long i, steps, term;
1033   double sum;
1034
1035   sum = 0.0;
1036   for (i = 0; i < (chars); i++) {
1037     steps = r->numsteps[i];
1038     if (steps <= threshwt[i])
1039       term = steps;
1040     else
1041       term = threshwt[i];
1042     sum += term;
1043     if (usertree && which <= maxuser)
1044       fsteps[which - 1][i] = term;
1045   }
1046   if (usertree && which <= maxuser) {
1047     nsteps[which - 1] = sum;
1048     if (which == 1) {
1049       minwhich = 1;
1050       minsteps = sum;
1051     } else if (sum < minsteps) {
1052       minwhich = which;
1053       minsteps = sum;
1054     }
1055   }
1056   like = -sum;
1057 }  /* evaluate */
1058
1059
1060 void protpostorder(node *p)
1061 {
1062   /* traverses a binary tree, calling PROCEDURE fillin at a
1063      node's descendants before calling fillin at the node */
1064   if (p->tip)
1065     return;
1066   protpostorder(p->next->back);
1067   protpostorder(p->next->next->back);
1068   protfillin(p, p->next->back, p->next->next->back);
1069 }  /* protpostorder */
1070
1071
1072 void protreroot(node *outgroup)
1073 {
1074   /* reorients tree, putting outgroup in desired position. */
1075   node *p, *q;
1076
1077   if (outgroup->back->index == root->index)
1078     return;
1079   p = root->next;
1080   q = root->next->next;
1081   p->back->back = q->back;
1082   q->back->back = p->back;
1083   p->back = outgroup;
1084   q->back = outgroup->back;
1085   outgroup->back->back = q;
1086   outgroup->back = p;
1087 }  /* protreroot */
1088
1089
1090 void protsavetraverse(node *p, long *pos, boolean *found)
1091 {
1092   /* sets BOOLEANs that indicate which way is down */
1093   p->bottom = true;
1094   if (p->tip)
1095     return;
1096   p->next->bottom = false;
1097   protsavetraverse(p->next->back, pos,found);
1098   p->next->next->bottom = false;
1099   protsavetraverse(p->next->next->back, pos,found);
1100 }  /* protsavetraverse */
1101
1102
1103 void protsavetree(long *pos, boolean *found)
1104 {
1105   /* record in place where each species has to be
1106      added to reconstruct this tree */
1107   long i, j;
1108   node *p;
1109   boolean done;
1110
1111   protreroot(treenode[outgrno - 1]);
1112   protsavetraverse(root, pos,found);
1113   for (i = 0; i < (nonodes); i++)
1114     place[i] = 0;
1115   place[root->index - 1] = 1;
1116   for (i = 1; i <= (spp); i++) {
1117     p = treenode[i - 1];
1118     while (place[p->index - 1] == 0) {
1119       place[p->index - 1] = i;
1120       while (!p->bottom)
1121         p = p->next;
1122       p = p->back;
1123     }
1124     if (i > 1) {
1125       place[i - 1] = place[p->index - 1];
1126       j = place[p->index - 1];
1127       done = false;
1128       while (!done) {
1129         place[p->index - 1] = spp + i - 1;
1130         while (!p->bottom)
1131           p = p->next;
1132         p = p->back;
1133         done = (p == NULL);
1134         if (!done)
1135           done = (place[p->index - 1] != j);
1136       }
1137     }
1138   }
1139 }  /* protsavetree */
1140
1141
1142 void tryadd(node *p, node **item, node **nufork)
1143 {
1144   /* temporarily adds one fork and one tip to the tree.
1145      if the location where they are added yields greater
1146      "likelihood" than other locations tested up to that
1147      time, then keeps that location as there */
1148   long pos;
1149   boolean found;
1150   node *rute, *q;
1151
1152   if (p == root)
1153     protfillin(temp, *item, p);
1154   else {
1155     protfillin(temp1, *item, p);
1156     protfillin(temp, temp1, p->back);
1157   }
1158   evaluate(temp);
1159   if (lastrearr) {
1160     if (like < bestlike) {
1161       if ((*item) == (*nufork)->next->next->back) {
1162         q = (*nufork)->next;
1163         (*nufork)->next = (*nufork)->next->next;
1164         (*nufork)->next->next = q;
1165         q->next = (*nufork);
1166       }
1167     }
1168     else if (like >= bstlike2) {
1169       recompute = false;
1170       protadd(p, (*item), (*nufork));
1171       rute = root->next->back;
1172       protsavetree(&pos,&found);
1173       protreroot(rute);
1174       if (like > bstlike2) {
1175         bestlike = bstlike2 = like;
1176         pos = 1;
1177         nextree = 1;
1178         addtree(pos, &nextree, dummy, place, bestrees);
1179       } else {
1180         pos = 0;
1181         findtree(&found, &pos, nextree, place, bestrees);
1182         if (!found) {
1183           if (nextree <= maxtrees)
1184             addtree(pos, &nextree, dummy, place, bestrees);
1185         }
1186       }
1187       protre_move (item, nufork);
1188       recompute = true;
1189     }
1190   }
1191   if (like >= bestyet) {
1192     bestyet = like;
1193     there = p;
1194   }
1195 }  /* tryadd */
1196
1197
1198 void addpreorder(node *p, node *item, node *nufork)
1199 {
1200   /* traverses a binary tree, calling PROCEDURE tryadd
1201      at a node before calling tryadd at its descendants */
1202
1203   if (p == NULL)
1204     return;
1205   tryadd(p, &item,&nufork);
1206   if (!p->tip) {
1207     addpreorder(p->next->back, item, nufork);
1208     addpreorder(p->next->next->back, item, nufork);
1209   }
1210 }  /* addpreorder */
1211
1212
1213 void tryrearr(node *p, boolean *success)
1214 {
1215   /* evaluates one rearrangement of the tree.
1216      if the new tree has greater "likelihood" than the old
1217      one sets success := TRUE and keeps the new tree.
1218      otherwise, restores the old tree */
1219   node *frombelow, *whereto, *forknode, *q;
1220   double oldlike;
1221
1222   if (p->back == NULL)
1223     return;
1224   forknode = treenode[p->back->index - 1];
1225   if (forknode->back == NULL)
1226     return;
1227   oldlike = bestyet;
1228   if (p->back->next->next == forknode)
1229     frombelow = forknode->next->next->back;
1230   else
1231     frombelow = forknode->next->back;
1232   whereto = treenode[forknode->back->index - 1];
1233   if (whereto->next->back == forknode)
1234     q = whereto->next->next->back;
1235   else
1236     q = whereto->next->back;
1237   protfillin(temp1, frombelow, q);
1238   protfillin(temp, temp1, p);
1239   protfillin(temp1, temp, whereto->back);
1240   evaluate(temp1);
1241   if (like <= oldlike) {
1242     if (p == forknode->next->next->back) {
1243       q = forknode->next;
1244       forknode->next = forknode->next->next;
1245       forknode->next->next = q;
1246       q->next = forknode;
1247     }
1248   }
1249   else {
1250     recompute = false;
1251     protre_move(&p, &forknode);
1252     protfillin(whereto, whereto->next->back, whereto->next->next->back);
1253     recompute = true;
1254     protadd(whereto, p, forknode);
1255     *success = true;
1256     bestyet = like;
1257   }
1258 }  /* tryrearr */
1259
1260
1261 void repreorder(node *p, boolean *success)
1262 {
1263   /* traverses a binary tree, calling PROCEDURE tryrearr
1264      at a node before calling tryrearr at its descendants */
1265   if (p == NULL)
1266     return;
1267   tryrearr(p,success);
1268   if (!p->tip) {
1269     repreorder(p->next->back,success);
1270     repreorder(p->next->next->back,success);
1271   }
1272 }  /* repreorder */
1273
1274
1275 void rearrange(node **r)
1276 {
1277   /* traverses the tree (preorder), finding any local
1278      rearrangement which decreases the number of steps.
1279      if traversal succeeds in increasing the tree's
1280      "likelihood", PROCEDURE rearrange runs traversal again */
1281   boolean success = true;
1282   while (success) {
1283     success = false;
1284     repreorder(*r, &success);
1285   }
1286 }  /* rearrange */
1287
1288
1289 void protgetch(Char *c)
1290 {
1291   /* get next nonblank character */
1292   do {
1293     if (eoln(intree))
1294       scan_eoln(intree);
1295     *c = gettc(intree);
1296     if (*c == '\n' || *c == '\t')
1297       *c = ' ';
1298   } while (!(*c != ' ' || eoff(intree)));
1299 }  /* protgetch */
1300
1301
1302 void protaddelement(node **p,long *nextnode,long *lparens,boolean *names)
1303 {
1304   /* recursive procedure adds nodes to user-defined tree */
1305   node *q;
1306   long i, n;
1307   boolean found;
1308   Char str[nmlngth];
1309
1310   protgetch(&ch);
1311
1312   if (ch == '(' ) {
1313     if ((*lparens) >= spp - 1) {
1314       printf("\nERROR IN USER TREE: TOO MANY LEFT PARENTHESES\n");
1315       exxit(-1);
1316     }
1317     (*nextnode)++;
1318     (*lparens)++;
1319     q = treenode[(*nextnode) - 1];
1320     protaddelement(&q->next->back, nextnode,lparens,names);
1321     q->next->back->back = q->next;
1322     findch(',', &ch, which);
1323     protaddelement(&q->next->next->back, nextnode,lparens,names);
1324     q->next->next->back->back = q->next->next;
1325     findch(')', &ch, which);
1326     *p = q;
1327     return;
1328   }
1329   for (i = 0; i < nmlngth; i++)
1330     str[i] = ' ';
1331   n = 1;
1332   do {
1333     if (ch == '_')
1334       ch = ' ';
1335     str[n - 1] = ch;
1336     if (eoln(intree))
1337       scan_eoln(intree);
1338     ch = gettc(intree);
1339     n++;
1340   } while (ch != ',' && ch != ')' && ch != ':' && n <= nmlngth);
1341   n = 1;
1342   do {
1343     found = true;
1344     for (i = 0; i < nmlngth; i++)
1345       found = (found && ((str[i] == nayme[n - 1][i]) ||
1346                          ((nayme[n - 1][i] == '_') && (str[i] == ' '))));
1347     if (found) {
1348       if (names[n - 1] == false) {
1349         *p = treenode[n - 1];
1350         names[n - 1] = true;
1351       } else {
1352         printf("\nERROR IN USER TREE: DUPLICATE NAME FOUND -- ");
1353         for (i = 0; i < nmlngth; i++)
1354           putchar(nayme[n - 1][i]);
1355         putchar('\n');
1356         exxit(-1);
1357       }
1358     } else
1359       n++;
1360   } while (!(n > spp || found));
1361   if (n <= spp)
1362     return;
1363   printf("CANNOT FIND SPECIES: ");
1364   for (i = 0; i < nmlngth; i++)
1365     putchar(str[i]);
1366   putchar('\n');
1367 }  /* protaddelement */
1368
1369
1370 void prottreeread()
1371 {
1372   /* read in user-defined tree and set it up */
1373   long nextnode, lparens, i;
1374
1375   root = treenode[spp];
1376   nextnode = spp;
1377   root->back = NULL;
1378   names = (boolean *)Malloc(spp*sizeof(boolean));
1379   for (i = 0; i < (spp); i++)
1380     names[i] = false;
1381   lparens = 0;
1382   protaddelement(&root, &nextnode,&lparens,names);
1383   if (ch == '[') {
1384     do
1385       ch = gettc(intree);
1386     while (ch != ']');
1387     ch = gettc(intree);
1388   }
1389   findch(';', &ch, which);
1390   scan_eoln(intree);
1391   free(names);
1392 }  /* prottreeread */
1393
1394
1395 void protancestset(long *a, long *b, long *c, long *d, long *k)
1396 {
1397   /* sets up the aa set array. */
1398   aas aa;
1399   long s, sa, sb;
1400   long i, j, m, n;
1401   boolean counted;
1402
1403   counted = false;
1404   *k = 0;
1405   for (i = 0; i <= 5; i++) {
1406     if (*k < 3) {
1407       s = 0;
1408       if (i > 3)
1409         n = i - 3;
1410       else
1411         n = 0;
1412       for (j = n; j <= (i - n); j++) {
1413         if (j < 3)
1414           sa = a[j];
1415         else
1416           sa = fullset;
1417         for (m = n; m <= (i - j - n); m++) {
1418           if (m < 3)
1419             sb = sa & b[m];
1420           else
1421             sb = sa;
1422           if (i - j - m < 3)
1423             sb &= c[i - j - m];
1424           s |= sb;
1425         }
1426       }
1427       if (counted || s != 0) {
1428         d[*k] = s;
1429         (*k)++;
1430         counted = true;
1431       }
1432     }
1433   }
1434   for (i = 0; i <= 1; i++) {
1435     for (aa = ala; (long)aa <= (long)stop; aa = (aas)((long)aa + 1)) {
1436       if (((1L << ((long)aa)) & d[i]) != 0) {
1437         for (j = i + 1; j <= 2; j++)
1438           d[j] |= translate[(long)aa - (long)ala][j - i];
1439       }
1440     }
1441   }
1442 }  /* protancestset */
1443
1444
1445 void prothyprint(long b1, long b2, boolean *bottom, node *r,
1446                         boolean *nonzero, boolean *maybe)
1447 {
1448   /* print out states in sites b1 through b2 at node */
1449   long i;
1450   boolean dot;
1451   Char ch = 0;
1452   aas aa;
1453
1454   if (*bottom) {
1455     if (!outgropt)
1456       fprintf(outfile, "      ");
1457     else
1458       fprintf(outfile, "root  ");
1459   } else
1460     fprintf(outfile, "%3ld   ", r->back->index - spp);
1461   if (r->tip) {
1462     for (i = 0; i < nmlngth; i++)
1463       putc(nayme[r->index - 1][i], outfile);
1464   } else
1465     fprintf(outfile, "%4ld      ", r->index - spp);
1466   if (*bottom)
1467     fprintf(outfile, "          ");
1468   else if (*nonzero)
1469     fprintf(outfile, "   yes    ");
1470   else if (*maybe)
1471     fprintf(outfile, "  maybe   ");
1472   else
1473     fprintf(outfile, "   no     ");
1474   for (i = b1 - 1; i < b2; i++) {
1475     aa = r->seq[i];
1476     switch (aa) {
1477
1478     case ala:
1479       ch = 'A';
1480       break;
1481
1482     case asx:
1483       ch = 'B';
1484       break;
1485
1486     case cys:
1487       ch = 'C';
1488       break;
1489
1490     case asp:
1491       ch = 'D';
1492       break;
1493
1494     case glu:
1495       ch = 'E';
1496       break;
1497
1498     case phe:
1499       ch = 'F';
1500       break;
1501
1502     case gly:
1503       ch = 'G';
1504       break;
1505
1506     case his:
1507       ch = 'H';
1508       break;
1509
1510     case ileu:
1511       ch = 'I';
1512       break;
1513
1514     case lys:
1515       ch = 'K';
1516       break;
1517
1518     case leu:
1519       ch = 'L';
1520       break;
1521
1522     case met:
1523       ch = 'M';
1524       break;
1525
1526     case asn:
1527       ch = 'N';
1528       break;
1529
1530     case pro:
1531       ch = 'P';
1532       break;
1533
1534     case gln:
1535       ch = 'Q';
1536       break;
1537
1538     case arg:
1539       ch = 'R';
1540       break;
1541
1542     case ser:
1543       ch = 'S';
1544       break;
1545
1546     case ser1:
1547       ch = 'S';
1548       break;
1549
1550     case ser2:
1551       ch = 'S';
1552       break;
1553
1554     case thr:
1555       ch = 'T';
1556       break;
1557
1558     case trp:
1559       ch = 'W';
1560       break;
1561
1562     case tyr:
1563       ch = 'Y';
1564       break;
1565
1566     case val:
1567       ch = 'V';
1568       break;
1569
1570     case glx:
1571       ch = 'Z';
1572       break;
1573
1574     case del:
1575       ch = '-';
1576       break;
1577
1578     case stop:
1579       ch = '*';
1580       break;
1581
1582     case unk:
1583       ch = 'X';
1584       break;
1585
1586     case quest:
1587       ch = '?';
1588       break;
1589     }
1590     if (!(*bottom) && dotdiff)
1591       dot = (r->siteset[i] [0] == treenode[r->back->index - 1]->siteset[i][0]
1592              || ((r->siteset[i][0] &
1593                   (~((1L << ((long)ser1)) | (1L << ((long)ser2)) |
1594                                          (1L << ((long)ser))))) == 0 &&
1595                 (treenode[r->back->index - 1]->siteset[i] [0] &
1596                   (~((1L << ((long)ser1)) | (1L << ((long)ser2)) |
1597                                          (1L << ((long)ser))))) == 0));
1598     else
1599       dot = false;
1600     if (dot)
1601       putc('.', outfile);
1602     else
1603       putc(ch, outfile);
1604     if ((i + 1) % 10 == 0)
1605       putc(' ', outfile);
1606   }
1607   putc('\n', outfile);
1608 }  /* prothyprint */
1609
1610
1611 void prothyptrav(node *r, sitearray *hypset, long b1, long b2, long *k,
1612                         boolean *bottom, sitearray nothing)
1613 {
1614   boolean maybe, nonzero;
1615   long i;
1616   aas aa;
1617   long anc = 0, hset;
1618   gseq *ancset, *temparray;
1619
1620   protgnu(&ancset);
1621   protgnu(&temparray);
1622   maybe = false;
1623   nonzero = false;
1624   for (i = b1 - 1; i < b2; i++) {
1625     if (!r->tip) {
1626       protancestset(hypset[i], r->next->back->siteset[i],
1627                 r->next->next->back->siteset[i], temparray->seq[i], k);
1628       memcpy(r->siteset[i], temparray->seq[i], sizeof(sitearray));
1629     }
1630     if (!(*bottom))
1631       anc = treenode[r->back->index - 1]->siteset[i][0];
1632     if (!r->tip) {
1633       hset = r->siteset[i][0];
1634       r->seq[i] = quest;
1635       for (aa = ala; (long)aa <= (long)stop; aa = (aas)((long)aa + 1)) {
1636         if (hset == 1L << ((long)aa))
1637           r->seq[i] = aa;
1638       }
1639       if (hset == ((1L << ((long)asn)) | (1L << ((long)asp))))
1640         r->seq[i] = asx;
1641       if (hset == ((1L << ((long)gln)) | (1L << ((long)gly))))
1642         r->seq[i] = glx;
1643       if (hset == ((1L << ((long)ser1)) | (1L << ((long)ser2))))
1644         r->seq[i] = ser;
1645       if (hset == fullset)
1646         r->seq[i] = unk;
1647     }
1648     nonzero = (nonzero || (r->siteset[i][0] & anc) == 0);
1649     maybe = (maybe || r->siteset[i][0] != anc);
1650   }
1651   prothyprint(b1, b2,bottom,r,&nonzero,&maybe);
1652   *bottom = false;
1653   if (!r->tip) {
1654     memcpy(temparray->seq, r->next->back->siteset, chars*sizeof(sitearray));
1655     for (i = b1 - 1; i < b2; i++)
1656       protancestset(hypset[i], r->next->next->back->siteset[i], nothing,
1657                 ancset->seq[i], k);
1658     prothyptrav(r->next->back, ancset->seq, b1, b2,k,bottom,nothing );
1659     for (i = b1 - 1; i < b2; i++)
1660       protancestset(hypset[i], temparray->seq[i], nothing, ancset->seq[i],k);
1661     prothyptrav(r->next->next->back, ancset->seq, b1, b2, k,bottom,nothing);
1662   }
1663   protchuck(temparray);
1664   protchuck(ancset);
1665 }  /* prothyptrav */
1666
1667
1668 void prothypstates(long *k)
1669 {
1670   /* fill in and describe states at interior nodes */
1671   boolean bottom;
1672   sitearray nothing;
1673   long i, n;
1674   seqptr hypset;
1675
1676   fprintf(outfile, "\nFrom    To     Any Steps?    State at upper node\n");
1677   fprintf(outfile, "                             ");
1678   fprintf(outfile, "( . means same as in the node below it on tree)\n\n");
1679   memcpy(nothing, translate[(long)quest - (long)ala], sizeof(sitearray));
1680   hypset = (seqptr)Malloc(chars*sizeof(sitearray));
1681   for (i = 0; i < (chars); i++)
1682     memcpy(hypset[i], nothing, sizeof(sitearray));
1683   bottom = true;
1684   for (i = 1; i <= ((chars - 1) / 40 + 1); i++) {
1685     putc('\n', outfile);
1686     n = i * 40;
1687     if (n > chars)
1688       n = chars;
1689     bottom = true;
1690     prothyptrav(root, hypset, i * 40 - 39, n, k,&bottom,nothing);
1691   }
1692   free(hypset);
1693 }  /* prothypstates */
1694
1695
1696 void describe()
1697 {
1698   /* prints ancestors, steps and table of numbers of steps in
1699      each site */
1700   long i,j,k;
1701
1702   if (treeprint)
1703     fprintf(outfile, "\nrequires a total of %10.3f\n", like / -10);
1704   if (stepbox) {
1705     putc('\n', outfile);
1706     if (weights)
1707       fprintf(outfile, "weighted ");
1708     fprintf(outfile, "steps in each position:\n");
1709     fprintf(outfile, "      ");
1710     for (i = 0; i <= 9; i++)
1711       fprintf(outfile, "%4ld", i);
1712     fprintf(outfile, "\n     *-----------------------------------------\n");
1713     for (i = 0; i <= (chars / 10); i++) {
1714       fprintf(outfile, "%5ld", i * 10);
1715       putc('!', outfile);
1716       for (j = 0; j <= 9; j++) {
1717         k = i * 10 + j;
1718         if (k == 0 || k > chars)
1719           fprintf(outfile, "    ");
1720         else
1721           fprintf(outfile, "%4ld", root->numsteps[k - 1] / 10);
1722       }
1723       putc('\n', outfile);
1724     }
1725   }
1726   if (ancseq) {
1727     prothypstates(&k);
1728     putc('\n', outfile);
1729   }
1730   putc('\n', outfile);
1731   if (trout) {
1732     col = 0;
1733     treeout(root, nextree, &col, root);
1734   }
1735 }  /* describe */
1736
1737
1738 void maketree()
1739 {
1740   /* constructs a binary tree from the pointers in treenode.
1741      adds each node at location which yields highest "likelihood"
1742      then rearranges the tree for greatest "likelihood" */
1743   long i, j, numtrees;
1744   double gotlike;
1745   node *item, *nufork, *dummy;
1746
1747   if (!usertree) {
1748     for (i = 1; i <= (spp); i++)
1749       enterorder[i - 1] = i;
1750     if (jumble)
1751       randumize(seed, enterorder);
1752     root = treenode[enterorder[0] - 1];
1753     recompute = true;
1754     protadd(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
1755         treenode[spp]);
1756     if (progress) {
1757       printf("\nAdding species:\n");
1758       writename(0, 2, enterorder);
1759     }
1760     lastrearr = false;
1761     for (i = 3; i <= (spp); i++) {
1762       bestyet = -30.0*spp*chars;
1763       there = root;
1764       item = treenode[enterorder[i - 1] - 1];
1765       nufork = treenode[spp + i - 2];
1766       addpreorder(root, item, nufork);
1767       protadd(there, item, nufork);
1768       like = bestyet;
1769       rearrange(&root);
1770       if (progress)
1771         writename(i - 1, 1, enterorder);
1772       lastrearr = (i == spp);
1773       if (lastrearr) {
1774         if (progress) {
1775           printf("\nDoing global rearrangements\n");
1776           printf("  !");
1777           for (j = 1; j <= nonodes; j++)
1778             if ( j % (( nonodes / 72 ) + 1 ) == 0 )
1779               putchar('-');
1780           printf("!\n");
1781         }
1782         bestlike = bestyet;
1783         if (jumb == 1) {
1784           bstlike2 = bestlike = -30.0*spp*chars;
1785           nextree = 1;
1786         }
1787         do {
1788           if (progress)
1789             printf("   ");
1790           gotlike = bestlike;
1791           for (j = 0; j < (nonodes); j++) {
1792             bestyet = -30.0*spp*chars;
1793             item = treenode[j];
1794             if (item != root) {
1795               nufork = treenode[treenode[j]->back->index - 1];
1796               protre_move(&item, &nufork);
1797               there = root;
1798               addpreorder(root, item, nufork);
1799               protadd(there, item, nufork);
1800             }
1801             if (progress) {
1802               if ( j % (( nonodes / 72 ) + 1 ) == 0 )
1803                 putchar('.');
1804               fflush(stdout);
1805             }
1806           }
1807           if (progress)
1808             putchar('\n');
1809         } while (bestlike > gotlike);
1810       }
1811     }
1812     if (progress)
1813       putchar('\n');
1814     for (i = spp - 1; i >= 1; i--)
1815       protre_move(&treenode[i], &dummy);
1816     if (jumb == njumble) {
1817       if (treeprint) {
1818         putc('\n', outfile);
1819         if (nextree == 2)
1820           fprintf(outfile, "One most parsimonious tree found:\n");
1821         else
1822           fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
1823       }
1824       if (nextree > maxtrees + 1) {
1825         if (treeprint)
1826           fprintf(outfile, "here are the first%4ld of them\n", (long)maxtrees);
1827         nextree = maxtrees + 1;
1828       }
1829       if (treeprint)
1830         putc('\n', outfile);
1831       recompute = false;
1832       for (i = 0; i <= (nextree - 2); i++) {
1833         root = treenode[0];
1834         protadd(treenode[0], treenode[1], treenode[spp]);
1835         for (j = 3; j <= (spp); j++)
1836           protadd(treenode[bestrees[i].btree[j - 1] - 1], treenode[j - 1],
1837               treenode[spp + j - 2]);
1838         protreroot(treenode[outgrno - 1]);
1839         protpostorder(root);
1840         evaluate(root);
1841         printree(root, 1.0);
1842         describe();
1843         for (j = 1; j < (spp); j++)
1844           protre_move(&treenode[j], &dummy);
1845       }
1846     }
1847   } else {
1848     openfile(&intree,INTREE,"input tree file", "r",progname,intreename);
1849     numtrees = countsemic(&intree);
1850     if (treeprint) {
1851       fprintf(outfile, "User-defined tree");
1852       if (numtrees > 1)
1853         putc('s', outfile);
1854       fprintf(outfile, ":\n\n\n\n");
1855     }
1856     which = 1;
1857     while (which <= numtrees) {
1858       prottreeread();
1859       if (outgropt)
1860         protreroot(treenode[outgrno - 1]);
1861       protpostorder(root);
1862       evaluate(root);
1863       printree(root, 1.0);
1864       describe();
1865       which++;
1866     }
1867     printf("\n");
1868     FClose(intree);
1869     putc('\n', outfile);
1870     if (numtrees > 1 && chars > 1 )
1871       standev(chars, numtrees, minwhich, minsteps, nsteps, fsteps, seed);
1872   }
1873   if (jumb == njumble && progress) {
1874     printf("Output written to file \"%s\"\n\n", outfilename);
1875     if (trout)
1876       printf("Trees also written onto file \"%s\"\n\n", outtreename);
1877   }
1878 }  /* maketree */
1879
1880
1881 int main(int argc, Char *argv[])
1882 {  /* Protein parsimony by uphill search */
1883 #ifdef MAC
1884   argc = 1;         /* macsetup("Protpars","");                */
1885   argv[0] = "Protpars";
1886 #endif
1887   init(argc,argv);
1888   progname = argv[0];
1889   openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
1890   openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
1891
1892   ibmpc = IBMCRT;
1893   ansi = ANSICRT;
1894   garbage = NULL;
1895   mulsets = false;
1896   msets = 1;
1897   firstset = true;
1898   code();
1899   setup();
1900   doinit();
1901   if (weights || justwts)
1902     openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
1903   if (trout)
1904     openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
1905   for (ith = 1; ith <= msets; ith++) {
1906     doinput();
1907     if (ith == 1)
1908       firstset = false;
1909     if (msets > 1 && !justwts) {
1910       fprintf(outfile, "Data set # %ld:\n\n",ith);
1911       if (progress)
1912         printf("Data set # %ld:\n\n",ith);
1913     }
1914     for (jumb = 1; jumb <= njumble; jumb++)
1915       maketree();
1916   }
1917   FClose(infile);
1918   FClose(outfile);
1919   FClose(outtree);
1920 #ifdef MAC
1921   fixmacfile(outfilename);
1922   fixmacfile(outtreename);
1923 #endif
1924   return 0;
1925 }  /* Protein parsimony by uphill search */