protpars.c

一个神经网络工具箱

    case his:
      ch = 'H';
      break;

    case ileu:
      ch = 'I';
      break;

    case lys:
      ch = 'K';
      break;

    case leu:
      ch = 'L';
      break;

    case met:
      ch = 'M';
      break;

    case asn:
      ch = 'N';
      break;

    case pro:
      ch = 'P';
      break;

    case gln:
      ch = 'Q';
      break;

    case arg:
      ch = 'R';
      break;

    case ser:
      ch = 'S';
      break;

    case ser1:
      ch = 'S';
      break;

    case ser2:
      ch = 'S';
      break;

    case thr:
      ch = 'T';
      break;

    case trp:
      ch = 'W';
      break;

    case tyr:
      ch = 'Y';
      break;

    case val:
      ch = 'V';
      break;

    case glx:
      ch = 'Z';
      break;

    case del:
      ch = '-';
      break;

    case stop:
      ch = '*';
      break;

    case unk:
      ch = 'X';
      break;

    case quest:
      ch = '?';
      break;
    }
    if (!(*bottom))
      dot = (r->siteset[i] [0] == treenode[r->back->index - 1]->siteset[i][0]
             || ((r->siteset[i][0] &
                  (~((1L << ((long)ser1)) | (1L << ((long)ser2)) |
                                         (1L << ((long)ser))))) == 0 &&
                (treenode[r->back->index - 1]->siteset[i] [0] &
                  (~((1L << ((long)ser1)) | (1L << ((long)ser2)) |
                                         (1L << ((long)ser))))) == 0));
    else
      dot = false;
    if (dot)
      putc('.', outfile);
    else
      putc(ch, outfile);
    if ((i + 1) % 10 == 0)
      putc(' ', outfile);
  }
  putc('\n', outfile);
}  /* prothyprint */


void prothyptrav(node *r, sitearray *hypset, long b1, long b2, long *k,
                        boolean *bottom, sitearray nothing)
{
  boolean maybe, nonzero;
  long i;
  aas aa;
  long anc = 0, hset;
  gseq *ancset, *temparray;

  protgnu(&ancset);
  protgnu(&temparray);
  maybe = false;
  nonzero = false;
  for (i = b1 - 1; i < b2; i++) {
    if (!r->tip) {
      protancestset(hypset[i], r->next->back->siteset[i],
                r->next->next->back->siteset[i], temparray->seq[i], k);
      memcpy(r->siteset[i], temparray->seq[i], sizeof(sitearray));
    }
    if (!(*bottom))
      anc = treenode[r->back->index - 1]->siteset[i][0];
    if (!r->tip) {
      hset = r->siteset[i][0];
      r->seq[i] = quest;
      for (aa = ala; (long)aa <= (long)stop; aa = (aas)((long)aa + 1)) {
        if (hset == 1L << ((long)aa))
          r->seq[i] = aa;
      }
      if (hset == ((1L << ((long)asn)) | (1L << ((long)asp))))
        r->seq[i] = asx;
      if (hset == ((1L << ((long)gln)) | (1L << ((long)gly))))
        r->seq[i] = glx;
      if (hset == ((1L << ((long)ser1)) | (1L << ((long)ser2))))
        r->seq[i] = ser;
      if (hset == fullset)
        r->seq[i] = unk;
    }
    nonzero = (nonzero || (r->siteset[i][0] & anc) == 0);
    maybe = (maybe || r->siteset[i][0] != anc);
  }
  prothyprint(b1, b2,bottom,r,&nonzero,&maybe);
  *bottom = false;
  if (!r->tip) {
    memcpy(temparray->seq, r->next->back->siteset, chars*sizeof(sitearray));
    for (i = b1 - 1; i < b2; i++)
      protancestset(hypset[i], r->next->next->back->siteset[i], nothing,
                ancset->seq[i], k);
    prothyptrav(r->next->back, ancset->seq, b1, b2,k,bottom,nothing );
    for (i = b1 - 1; i < b2; i++)
      protancestset(hypset[i], temparray->seq[i], nothing, ancset->seq[i],k);
    prothyptrav(r->next->next->back, ancset->seq, b1, b2, k,bottom,nothing);
  }
  protchuck(temparray);
  protchuck(ancset);
}  /* prothyptrav */


void prothypstates(long *k)
{
  /* fill in and describe states at interior nodes */
  boolean bottom;
  sitearray nothing;
  long i, n;
  seqptr hypset;

  fprintf(outfile, "\nFrom    To     Any Steps?    State at upper node\n");
  fprintf(outfile, "                             ");
  fprintf(outfile, "( . means same as in the node below it on tree)\n\n");
  memcpy(nothing, translate[(long)quest - (long)ala], sizeof(sitearray));
  hypset = (seqptr)Malloc(chars*sizeof(sitearray));
  for (i = 0; i < (chars); i++)
    memcpy(hypset[i], nothing, sizeof(sitearray));
  bottom = true;
  for (i = 1; i <= ((chars - 1) / 40 + 1); i++) {
    putc('\n', outfile);
    n = i * 40;
    if (n > chars)
      n = chars;
    bottom = true;
    prothyptrav(root, hypset, i * 40 - 39, n, k,&bottom,nothing);
  }
  free(hypset);
}  /* prothypstates */


void describe()
{
  /* prints ancestors, steps and table of numbers of steps in
     each site */
  long i,j,k;

  if (treeprint)
    fprintf(outfile, "\nrequires a total of %10.3f\n", like / -10);
  if (stepbox) {
    putc('\n', outfile);
    if (weights)
      fprintf(outfile, "weighted ");
    fprintf(outfile, "steps in each position:\n");
    fprintf(outfile, "      ");
    for (i = 0; i <= 9; i++)
      fprintf(outfile, "%4ld", i);
    fprintf(outfile, "\n     *-----------------------------------------\n");
    for (i = 0; i <= (chars / 10); i++) {
      fprintf(outfile, "%5ld", i * 10);
      putc('!', outfile);
      for (j = 0; j <= 9; j++) {
        k = i * 10 + j;
        if (k == 0 || k > chars)
          fprintf(outfile, "    ");
        else
          fprintf(outfile, "%4ld", root->numsteps[k - 1] / 10);
      }
      putc('\n', outfile);
    }
  }
  if (ancseq) {
    prothypstates(&k);
    putc('\n', outfile);
  }
  putc('\n', outfile);
  if (trout) {
    col = 0;
    treeout(root, nextree, &col, root);
  }
}  /* describe */


void maketree()
{
  /* constructs a binary tree from the pointers in treenode.
     adds each node at location which yields highest "likelihood"
     then rearranges the tree for greatest "likelihood" */
  long i, j, numtrees;
  double gotlike;
  node *item, *nufork, *dummy;

  if (!usertree) {
    for (i = 1; i <= (spp); i++)
      enterorder[i - 1] = i;
    if (jumble)
      randumize(seed, enterorder);
    root = treenode[enterorder[0] - 1];
    recompute = true;
    protadd(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
        treenode[spp]);
    if (progress) {
      printf("\nAdding species:\n");
      writename(0, 2, enterorder);
    }
    lastrearr = false;
    for (i = 3; i <= (spp); i++) {
      bestyet = -30.0*spp*chars;
      there = root;
      item = treenode[enterorder[i - 1] - 1];
      nufork = treenode[spp + i - 2];
      addpreorder(root, item, nufork);
      protadd(there, item, nufork);
      like = bestyet;
      rearrange(&root);
      if (progress)
        writename(i - 1, 1, enterorder);
      lastrearr = (i == spp);
      if (lastrearr) {
        if (progress) {
          printf("\nDoing global rearrangements\n");
          printf("  !");
          for (j = 1; j <= nonodes; j++)
            putchar('-');
          printf("!\n");
        }
        bestlike = bestyet;
        if (jumb == 1) {
          bstlike2 = bestlike;
          nextree = 1;
        }
        do {
          if (progress)
            printf("   ");
          gotlike = bestlike;
          for (j = 0; j < (nonodes); j++) {
            bestyet = -30.0*spp*chars;
            item = treenode[j];
            if (item != root) {
              nufork = treenode[treenode[j]->back->index - 1];
              protre_move(&item, &nufork);
              there = root;
              addpreorder(root, item, nufork);
              protadd(there, item, nufork);
            }
            if (progress) {
              putchar('.');
              fflush(stdout);
            }
          }
          if (progress)
            putchar('\n');
        } while (bestlike > gotlike);
      }
    }
    if (progress)
      putchar('\n');
    for (i = spp - 1; i >= 1; i--)
      protre_move(&treenode[i], &dummy);
    if (jumb == njumble) {
      if (treeprint) {
        putc('\n', outfile);
        if (nextree == 2)
          fprintf(outfile, "One most parsimonious tree found:\n");
        else
          fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
      }
      if (nextree > maxtrees + 1) {
        if (treeprint)
          fprintf(outfile, "here are the first%4ld of them\n", (long)maxtrees);
        nextree = maxtrees + 1;
      }
      if (treeprint)
        putc('\n', outfile);
      recompute = false;
      for (i = 0; i <= (nextree - 2); i++) {
        root = treenode[0];
        protadd(treenode[0], treenode[1], treenode[spp]);
        for (j = 3; j <= (spp); j++)
          protadd(treenode[bestrees[i].btree[j - 1] - 1], treenode[j - 1],
              treenode[spp + j - 2]);
        protreroot(treenode[outgrno - 1]);
        protpostorder(root);
        evaluate(root);
        printree(root, 1.0);
        describe();
        for (j = 1; j < (spp); j++)
          protre_move(&treenode[j], &dummy);
      }
    }
  } else {
    openfile(&intree,INTREE,"input tree file", "r",progname,intreename);
    numtrees = countsemic(&intree);
    if (treeprint) {
      fprintf(outfile, "User-defined tree");
      if (numtrees > 1)
        putc('s', outfile);
      fprintf(outfile, ":\n\n\n\n");
    }
    which = 1;
    while (which <= numtrees) {
      prottreeread();
      if (outgropt)
        protreroot(treenode[outgrno - 1]);
      protpostorder(root);
      evaluate(root);
      printree(root, 1.0);
      describe();
      which++;
    }
    FClose(intree);
    putc('\n', outfile);
    if (numtrees > 1 && chars > 1 )
      standev(chars, numtrees, minwhich, minsteps, nsteps, fsteps, seed);
  }
  if (jumb == njumble && progress) {
    printf("Output written to file \"%s\"\n\n", outfilename);
    if (trout)
      printf("Trees also written onto file \"%s\"\n\n", outtreename);
  }
}  /* maketree */

void usage(void)
{
    fprintf(stderr,"Usage is %s [options]\n", progname);
    fprintf(stderr,"Options\n");
	fprintf(stderr,"  -j<number>  Number of times to jumble (0 = Not randomize input order)\n");
	fprintf(stderr,"  -s<number>  Seed\n");
    fprintf(stderr,"  -o<number>  Number of outgroup species\n");
	fprintf(stderr,"  -t<number>  Parsimony threshold (0 = Not use Threshold parsimony)\n");
    fprintf(stderr,"  -c<number>  Genetic code\n");
	exit (8);
}


int main(int argc, Char *argv[])
{  /* Protein parsimony by uphill search */
	long tempint=0;
#ifdef MAC
  argc = 1;         /* macsetup("Protpars","");                */
  argv[0] = "Protpars";
#endif
  init(argc,argv);
  progname = argv[0];


  printf("---------------------------------------\n");
  printf("  MBETOOLBOX_PROTPAR\n");
  printf("  The modified version of PROTPAR\n");
  printf("  (in phylip 3.6a3) for MEBToolbox\n");
  printf("---------------------------------------\n\n");


  if ((argc == 1) || (strcmp(argv[1],"-help") == 0)) {usage();};


  openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
  openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);

  ibmpc = IBMCRT;
  ansi = ANSICRT;
  garbage = NULL;
  mulsets = false;
  msets = 1;
  firstset = true;
  code();
  setup();
  doinit();

      while ((argc > 1) && (argv[1][0] == '-')) {
        switch (argv[1][1]) {
            case 'j':
                tempint = atoi(&argv[1][2]);
                if(tempint==0){
					jumble = false;
				}else{
					jumble = true;	njumble=tempint;
				};
                break;
            case 's':
                inseed0 = atoi(&argv[1][2]);
                inseed = atoi(&argv[1][2]);
                break;
			case 'o':
                tempint = atoi(&argv[1][2]);
                if(tempint==0){
					outgropt = false;
				    // fprintf(stderr,"  -o<number>  Number of outgroup species\n");
				}else{
					outgropt = true;	outgrno=tempint;
				    // fprintf(stderr,"  -o%d Number of outgroup species\n",outgrno );
				};
                break;
            case 't':
				tempint = atoi(&argv[1][2]);
                if(tempint==0){
					thresh = false;
				}else{
					thresh = true;	threshold=tempint;
				};
                break;
            case 'c':
                whichcode = atoi(&argv[1][2]);
                break;
			default:
                fprintf(stderr,"Bad option %s\n", argv[1]);
				usage();
        }
        /*
         * move the argument list up one
         * move the count down one
         */
        ++argv;
        --argc;
    }


  if (weights || justwts)
    openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
  if (trout)
    openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
  for (ith = 1; ith <= msets; ith++) {
    doinput();
    if (ith == 1)
      firstset = false;
    if (msets > 1 && !justwts) {
      fprintf(outfile, "Data set # %ld:\n\n",ith);
      if (progress)
        printf("Data set # %ld:\n\n",ith);
    }
    for (jumb = 1; jumb <= njumble; jumb++)
      maketree();
  }
  FClose(infile);
  FClose(outfile);
  FClose(outtree);
#ifdef MAC
  fixmacfile(outfilename);
  fixmacfile(outtreename);
#endif
  return 0;
}  /* Protein parsimony by uphill search */