protpars.c

一个神经网络工具箱

#include "phylip.h"
#include "seq.h"

/* version 3.6. (c) Copyright 1993-2002 by the University of Washington.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#define maxtrees        100   /* maximum number of tied trees stored */

typedef enum {
  universal, ciliate, mito, vertmito, flymito, yeastmito
} codetype;

/* nodes will form a binary tree */

typedef struct gseq {
  seqptr seq;
  struct gseq *next;
} gseq;

#ifndef OLDC
/* function prototypes */
void   protgnu(gseq **);
void   protchuck(gseq *);
void   code(void);
void   setup(void);
void   getoptions(void);
void   protalloctree(void);
void   allocrest(void);
void   doinit(void);
void   protinputdata(void);

void   protmakevalues(void);
void   doinput(void);
void   protfillin(node *, node *, node *);
void   protpreorder(node *);
void   protadd(node *, node *, node *);
void   protre_move(node **, node **);
void   evaluate(node *);
void   protpostorder(node *);
void   protreroot(node *);
void   protsavetraverse(node *, long *, boolean *);

void   protsavetree(long *, boolean *);
void   tryadd(node *, node **, node **);
void   addpreorder(node *, node *, node *);
void   tryrearr(node *, boolean *);
void   repreorder(node *, boolean *);
void   rearrange(node **);
void   protgetch(Char *);
void   protaddelement(node **, long *, long *, boolean *);
void   prottreeread(void);
void   protancestset(long *, long *, long *, long *, long *);
                
void   prothyprint(long , long , boolean *, node *, boolean *, boolean *);
void   prothyptrav(node *, sitearray *, long, long, long *, boolean *,
                sitearray);
void   prothypstates(long *);
void   describe(void);
void   maketree(void);
void   reallocnode(node* p);
void   reallocchars(void);
/* function prototypes */
#endif



Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], outtreename[FNMLNGTH], weightfilename[FNMLNGTH];
node *root;
long chars, col, msets, ith, njumble, jumb;
/*   chars = number of sites in actual sequences */
long inseed, inseed0;
boolean jumble, usertree, weights, thresh, trout, progress, stepbox,
    justwts, ancseq, mulsets, firstset;
codetype whichcode;
long fullset, fulldel;
pointarray treenode;   /* pointers to all nodes in tree */
double threshold;
steptr threshwt;
longer seed;
long *enterorder;
sitearray translate[(long)quest - (long)ala + 1];
aas trans[4][4][4];
long **fsteps;
bestelm *bestrees;
boolean dummy;
gseq *garbage;
node *temp, *temp1;
Char ch;
aas tmpa;
char *progname;

/* Local variables for maketree, propagated globally for c version: */
long minwhich;
double like, bestyet, bestlike, minsteps, bstlike2;
boolean lastrearr, recompute;
node *there;
double nsteps[maxuser];
long *place;
boolean *names;


void protgnu(gseq **p)
{
  /* this and the following are do-it-yourself garbage collectors.
     Make a new node or pull one off the garbage list */
  if (garbage != NULL) {
    *p = garbage;
    free((*p)->seq);
    (*p)->seq = (seqptr)Malloc(chars*sizeof(sitearray));
    garbage = garbage->next;
  } else {
    *p = (gseq *)Malloc(sizeof(gseq));
    (*p)->seq = (seqptr)Malloc(chars*sizeof(sitearray));
  }
  (*p)->next = NULL;
}  /* protgnu */


void protchuck(gseq *p)
{
  /* collect garbage on p -- put it on front of garbage list */
  p->next = garbage;
  garbage = p;
}  /* protchuck */


void code()
{
  /* make up table of the code 1 = u, 2 = c, 3 = a, 4 = g */
  trans[0][0][0] = phe;
  trans[0][0][1] = phe;
  trans[0][0][2] = leu;
  trans[0][0][3] = leu;
  trans[0][1][0] = ser;
  trans[0][1][1] = ser1;
  trans[0][1][2] = ser1;
  trans[0][1][3] = ser1;
  trans[0][2][0] = tyr;
  trans[0][2][1] = tyr;
  trans[0][2][2] = stop;
  trans[0][2][3] = stop;
  trans[0][3][0] = cys;
  trans[0][3][1] = cys;
  trans[0][3][2] = stop;
  trans[0][3][3] = trp;
  trans[1][0][0] = leu;
  trans[1][0][1] = leu;
  trans[1][0][2] = leu;
  trans[1][0][3] = leu;
  trans[1][1][0] = pro;
  trans[1][1][1] = pro;
  trans[1][1][2] = pro;
  trans[1][1][3] = pro;
  trans[1][2][0] = his;
  trans[1][2][1] = his;
  trans[1][2][2] = gln;
  trans[1][2][3] = gln;
  trans[1][3][0] = arg;
  trans[1][3][1] = arg;
  trans[1][3][2] = arg;
  trans[1][3][3] = arg;
  trans[2][0][0] = ileu;
  trans[2][0][1] = ileu;
  trans[2][0][2] = ileu;
  trans[2][0][3] = met;
  trans[2][1][0] = thr;
  trans[2][1][1] = thr;
  trans[2][1][2] = thr;
  trans[2][1][3] = thr;
  trans[2][2][0] = asn;
  trans[2][2][1] = asn;
  trans[2][2][2] = lys;
  trans[2][2][3] = lys;
  trans[2][3][0] = ser2;
  trans[2][3][1] = ser2;
  trans[2][3][2] = arg;
  trans[2][3][3] = arg;
  trans[3][0][0] = val;
  trans[3][0][1] = val;
  trans[3][0][2] = val;
  trans[3][0][3] = val;
  trans[3][1][0] = ala;
  trans[3][1][1] = ala;
  trans[3][1][2] = ala;
  trans[3][1][3] = ala;
  trans[3][2][0] = asp;
  trans[3][2][1] = asp;
  trans[3][2][2] = glu;
  trans[3][2][3] = glu;
  trans[3][3][0] = gly;
  trans[3][3][1] = gly;
  trans[3][3][2] = gly;
  trans[3][3][3] = gly;
  if (whichcode == mito)
    trans[0][3][2] = trp;
  if (whichcode == vertmito) {
    trans[0][3][2] = trp;
    trans[2][3][2] = stop;
    trans[2][3][3] = stop;
    trans[2][0][2] = met;
  }
  if (whichcode == flymito) {
    trans[0][3][2] = trp;
    trans[2][0][2] = met;
    trans[2][3][2] = ser2;
  }
  if (whichcode == yeastmito) {
    trans[0][3][2] = trp;
    trans[1][0][2] = thr;
    trans[2][0][2] = met;
  }
} /* code */


void setup()
{
  /* set up set table to get aasets from aas */
  aas a, b;
  long i, j, k, l, s;

  for (a = ala; (long)a <= (long)stop; a = (aas)((long)a + 1)) {
    translate[(long)a - (long)ala][0] = 1L << ((long)a);
    translate[(long)a - (long)ala][1] = 1L << ((long)a);
  }
  for (i = 0; i <= 3; i++) {
    for (j = 0; j <= 3; j++) {
      for (k = 0; k <= 3; k++) {
        for (l = 0; l <= 3; l++) {
          translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[l][j][k]);
          translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[i][l][k]);
          translate[(long)trans[i][j][k]][1] |= (1L << (long)trans[i][j][l]);
        }
      }
    }
  }
  translate[(long)del - (long)ala][1] = 1L << ((long)del);
  fulldel = (1L << ((long)stop + 1)) - (1L << ((long)ala));
  fullset = fulldel & (~(1L << ((long)del)));
  translate[(long)asx - (long)ala][0]
    = (1L << ((long)asn)) | (1L << ((long)asp));
  translate[(long)glx - (long)ala][0]
    = (1L << ((long)gln)) | (1L << ((long)glu));
  translate[(long)ser - (long)ala][0]
    = (1L << ((long)ser1)) | (1L << ((long)ser2));
  translate[(long)unk - (long)ala][0] = fullset;
  translate[(long)quest - (long)ala][0] = fulldel;
  translate[(long)asx - (long)ala][1] = translate[(long)asn - (long)ala][1]
                                       | translate[(long)asp - (long)ala][1];
  translate[(long)glx - (long)ala][1] = translate[(long)gln - (long)ala][1]
                                       | translate[(long)glu - (long)ala][1];
  translate[(long)ser - (long)ala][1] = translate[(long)ser1 - (long)ala][1]
                                       | translate[(long)ser2 - (long)ala][1];
  translate[(long)unk - (long)ala][1] = fullset;
  translate[(long)quest - (long)ala][1] = fulldel;
  for (a = ala; (long)a <= (long)quest; a = (aas)((long)a + 1)) {
    s = 0;
    for (b = ala; (long)b <= (long)stop; b = (aas)((long)b + 1)) {
      if (((1L << ((long)b)) & translate[(long)a - (long)ala][1]) != 0)
        s |= translate[(long)b - (long)ala][1];
    }
    translate[(long)a - (long)ala][2] = s;
  }
}  /* setup */


void getoptions()
{
  /* interactively set options */
  long loopcount, loopcount2;
  Char ch, ch2;

  fprintf(outfile, "\nProtein parsimony algorithm, version %s\n\n",VERSION);
  putchar('\n');
  jumble = false;
  njumble = 1;
  outgrno = 1;
  outgropt = false;
  thresh = false;
  trout = true;
  usertree = false;
  weights = false;

  whichcode = universal;
  printdata = false;
  progress = true;
  treeprint = true;
  stepbox = false;
  ancseq = false;
  interleaved = false;

/*
  loopcount = 0;
  for (;;) {
    cleerhome();
    printf("\nProtein parsimony algorithm, version %s\n\n",VERSION);
    printf("Setting for this run:\n");
    printf("  U                 Search for best tree?  %s\n",
           (usertree ? "No, use user trees in input file" : "Yes"));
    if (!usertree) {
      printf("  J   Randomize input order of sequences?");
      if (jumble)
        printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
      else
        printf("  No. Use input order\n");
    }
    printf("  O                        Outgroup root?");
    if (outgropt)
      printf("  Yes, at sequence number%3ld\n", outgrno);
    else
      printf("  No, use as outgroup species%3ld\n", outgrno);
    printf("  T              Use Threshold parsimony?");
    if (thresh)
      printf("  Yes, count steps up to%4.1f per site\n", threshold);
    else
      printf("  No, use ordinary parsimony\n");
    printf("  C               Use which genetic code?  %s\n",
      (whichcode == universal) ? "Universal"                  :
      (whichcode == ciliate)   ? "Ciliate"                    :
      (whichcode == mito)      ? "Universal mitochondrial"    :
      (whichcode == vertmito)  ? "Vertebrate mitochondrial"   :
      (whichcode == flymito)   ? "Fly mitochondrial"          :
      (whichcode == yeastmito) ? "Yeast mitochondrial"        : "");
    printf("  W                       Sites weighted?  %s\n",
           (weights ? "Yes" : "No"));
    printf("  M           Analyze multiple data sets?");
    if (mulsets)
        printf("  Yes, %2ld %s\n", msets,
               (justwts ? "sets of weights" : "data sets"));
    else
      printf("  No\n");
    printf("  I          Input sequences interleaved?  %s\n",
           (interleaved ? "Yes" : "No"));
    printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
           (ibmpc ? "IBM PC" : ansi ? "ANSI" : "(none)"));
    printf("  1    Print out the data at start of run  %s\n",
           (printdata ? "Yes" : "No"));
    printf("  2  Print indications of progress of run  %s\n",
           (progress ? "Yes" : "No"));
    printf("  3                        Print out tree  %s\n",
           (treeprint ? "Yes" : "No"));
    printf("  4          Print out steps in each site  %s\n",
           (stepbox ? "Yes" : "No"));
    printf("  5  Print sequences at all nodes of tree  %s\n",
           (ancseq ? "Yes" : "No"));
    printf("  6       Write out trees onto tree file?  %s\n",
           (trout ? "Yes" : "No"));
    if(weights && justwts){
        printf(
         "WARNING:  W option and Multiple Weights options are both on.  ");
        printf(
         "The W menu option is unnecessary and has no additional effect. \n");
    }
    printf(
   "\nAre these settings correct? (type Y or the letter for one to change)\n");
    scanf("%c%*[^\n]", &ch);
    getchar();
    uppercase(&ch);
    if (ch == 'Y')
      break;
    if (strchr("WCJOTUMI1234560",ch) != NULL){
      switch (ch) {
        
      case 'J':
        jumble = !jumble;
        if (jumble)
          initjumble(&inseed, &inseed0, seed, &njumble);
        else njumble = 1;
        break;
      
      case 'W':
        weights = !weights;
        break;
  
      case 'O':
        outgropt = !outgropt;
        if (outgropt)
          initoutgroup(&outgrno, spp);
        else outgrno = 1;
        break;
        
      case 'T':
        thresh = !thresh;
        if (thresh)
          initthreshold(&threshold);
        break;

      case 'C':
        printf("\nWhich genetic code?\n");
        printf(" type         for\n\n");
        printf("   U           Universal\n");
        printf("   M           Mitochondrial\n");
        printf("   V           Vertebrate mitochondrial\n");
        printf("   F           Fly mitochondrial\n");
        printf("   Y           Yeast mitochondrial\n\n");
        loopcount2 = 0;
        do {
          printf("type U, M, V, F, or Y\n");
          scanf("%c%*[^\n]", &ch);
          getchar();
          if (ch == '\n')
            ch = ' ';
          uppercase(&ch);
          countup(&loopcount2, 10);
        } while (ch != 'U' && ch != 'M' && ch != 'V'
                  && ch != 'F' && ch != 'Y');
        switch (ch) {

        case 'U':
          whichcode = universal;
          break;

        case 'M':
          whichcode = mito;
          break;

        case 'V':
          whichcode = vertmito;
          break;

        case 'F':
          whichcode = flymito;
          break;

        case 'Y':
          whichcode = yeastmito;
          break;
        }
        break;

      case 'M':
        mulsets = !mulsets;
        if (mulsets){
            printf("Multiple data sets or multiple weights?");
          loopcount2 = 0;
          do {
            printf(" (type D or W)\n");
#ifdef WIN32
            phyFillScreenColor();
#endif
            scanf("%c%*[^\n]", &ch2);
            getchar();
            if (ch2 == '\n')
              ch2 = ' ';
            uppercase(&ch2);
            countup(&loopcount2, 10);
          } while ((ch2 != 'W') && (ch2 != 'D'));
          justwts = (ch2 == 'W');
          if (justwts)
            justweights(&msets);
          else
            initdatasets(&msets);
          if (!jumble) {
            jumble = true;
            initjumble(&inseed, &inseed0, seed, &njumble);
          }
        }
        break;
        
      case 'I':
        interleaved = !interleaved;
        break;
        
      case 'U':
        usertree = !usertree;
        break;
        
      case '0':
        initterminal(&ibmpc, &ansi);
        break;
        
      case '1':
        printdata = !printdata;
        break;
        
      case '2':
        progress = !progress;
        break;
        
      case '3':
        treeprint = !treeprint;
        break;
        
      case '4':
        stepbox = !stepbox;
        break;
        
      case '5':
        ancseq = !ancseq;
        break;