seq.c

一个神经网络工具箱

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

/* version 3.6. (c) Copyright 1993-2000 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. */

long nonodes, endsite, outgrno, nextree, which;
boolean interleaved, printdata, outgropt, treeprint, dotdiff, transvp;
steptr weight, category, alias, location, ally;
sequence y;

void free_all_x_in_array (long nonodes, pointarray treenode)
{
  /* used in dnaml & dnamlk */
  long i, j, k;
  node *p;

  /* Zero thru spp are tips, */
  for (i = 0; i < spp; i++) {
    for (j = 0; j < endsite; j++)
      free(treenode[i]->x[j]);
    free(treenode[i]->x);
  }

  /* The rest are rings (i.e. triads) */
  for (i = spp; i < nonodes; i++) {
    if (treenode[i] != NULL) {
      p = treenode[i];
      for (j = 1; j <= 3; j++) {
        for (k = 0; k < endsite; k++)
          free(p->x[k]);
        free(p->x);
        p = p->next;
      }
    }
  }
}  /* free_all_x_in_array        */


void free_all_x2_in_array (long nonodes, pointarray treenode)
{
  /* used in restml */
  long i, j;
  node *p;

  /* Zero thru spp are tips */
  for (i = 0; i < spp; i++)
    free(treenode[i]->x2);

  /* The rest are rings (i.e. triads) */
  for (i = spp; i < nonodes; i++) {
    p = treenode[i];
    for (j = 1; j <= 3; j++) {
      free(p->x2);
      p = p->next;
    }
  }
}  /* free_all_x2_in_array        */

void alloctemp(node **temp, long *zeros, long endsite)
{
  /*used in dnacomp and dnapenny */
  *temp = (node *)Malloc(sizeof(node));
  (*temp)->numsteps = (steptr)Malloc(endsite*sizeof(long));
  (*temp)->base = (baseptr)Malloc(endsite*sizeof(long));
  (*temp)->numnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
  memcpy((*temp)->base, zeros, endsite*sizeof(long));
  memcpy((*temp)->numsteps, zeros, endsite*sizeof(long));
  zeronumnuc(*temp, endsite);
}  /* alloctemp */


void freetemp(node **temp)
{
  /* used in dnacomp, dnapars, & dnapenny */
  free((*temp)->numsteps);
  free((*temp)->base);
  free((*temp)->numnuc);
  free(*temp);
}  /* freetemp */


void freetree2 (pointarray treenode, long nonodes)
{
  /* The natural complement to alloctree2.  Free all elements of all
  the rings (normally triads) in treenode */
  long i;
  node *p, *q;

  /* The first spp elements are just nodes, not rings */
  for (i = 0; i < spp; i++)
    free (treenode[i]);

  /* The rest are rings */
  for (i = spp; i < nonodes; i++) {
    p = treenode[i]->next;
    while (p != treenode[i]) {
      q = p->next;
      free (p);
      p = q;
    }
    /* p should now point to treenode[i], which has yet to be freed */
    free (p);
  }
  free (treenode);
}  /* freetree2 */


void inputdata(long chars)
{
  /* input the names and sequences for each species */
  /* used by dnacomp, dnadist, dnainvar, dnaml, dnamlk, dnapars, & dnapenny */
  long i, j, k, l, basesread, basesnew=0;
  Char charstate;
  boolean allread, done;

  if (printdata)
    headings(chars, "Sequences", "---------");
  basesread = 0;
  allread = false;
  while (!(allread)) {
    if (eoln(infile))
      scan_eoln(infile);
    i = 1;
    while (i <= spp) {
      if ((interleaved && basesread == 0) || !interleaved)
        initname(i-1);
      j = (interleaved) ? basesread : 0;
      done = false;
      while (!done && !eoff(infile)) {
        if (interleaved)
          done = true;
        while (j < chars && !(eoln(infile) || eoff(infile))) {
          charstate = gettc(infile);
          if (charstate == '\n')
            charstate = ' ';
          if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
            continue;
          uppercase(&charstate);
          if ((strchr("ABCDGHKMNRSTUVWXY?O-",charstate)) == NULL){
            printf("ERROR: bad base: %c at site %5ld of species %3ld\n",
                   charstate, j+1, i);
            if (charstate == '.') {
              printf("       Periods (.) may not be used as gap characters.\n");
              printf("       The correct gap character is (-)\n");
            }
            exxit(-1);
          }
          j++;
          y[i - 1][j - 1] = charstate;
        }
        if (interleaved)
          continue;
        if (j < chars) 
          scan_eoln(infile);
        else if (j == chars)
          done = true;
      }
      if (interleaved && i == 1)
        basesnew = j;

      scan_eoln(infile);
    
      if ((interleaved && j != basesnew) ||
          (!interleaved && j != chars)) {
        printf("\nERROR: sequences out of alignment at position %ld", j+1);
        printf(" of species %ld\n\n", i);
        exxit(-1);
      }
      i++;
    }
    
    if (interleaved) {
      basesread = basesnew;
      allread = (basesread == chars);
    } else
      allread = (i > spp);
  }
  if (!printdata)
    return;
  for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
    for (j = 1; j <= spp; j++) {
      for (k = 0; k < nmlngth; k++)
        putc(nayme[j - 1][k], outfile);
      fprintf(outfile, "   ");
      l = i * 60;
      if (l > chars)
        l = chars;
      for (k = (i - 1) * 60 + 1; k <= l; k++) {
        if (dotdiff && (j > 1 && y[j - 1][k - 1] == y[0][k - 1]))
          charstate = '.';
        else
          charstate = y[j - 1][k - 1];
        putc(charstate, outfile);
        if (k % 10 == 0 && k % 60 != 0)
          putc(' ', outfile);
      }
      putc('\n', outfile);
    }
    putc('\n', outfile);
  }
  putc('\n', outfile);
}  /* inputdata */


void alloctree(pointarray *treenode, long nonodes, boolean usertree)
{
  /* allocate treenode dynamically */
  /* used in dnapars, dnacomp, dnapenny & dnamove */
  long i, j;
  node *p, *q;

  *treenode = (pointarray)Malloc(nonodes*sizeof(node *));
  for (i = 0; i < spp; i++) {
    (*treenode)[i] = (node *)Malloc(sizeof(node));
    (*treenode)[i]->tip = true;
    (*treenode)[i]->index = i+1;
    (*treenode)[i]->iter = true;
    (*treenode)[i]->branchnum = i+1;
    (*treenode)[i]->initialized = true;
  }
  if (!usertree)
    for (i = spp; i < nonodes; i++) {
      q = NULL;
      for (j = 1; j <= 3; j++) {
        p = (node *)Malloc(sizeof(node));
        p->tip = false;
        p->index = i+1;
        p->iter = true;
        p->branchnum = i+1;
        p->initialized = false;
        p->next = q;
        q = p;
      }
      p->next->next->next = p;
      (*treenode)[i] = p;
    }
} /* alloctree */


void allocx(long nonodes, long rcategs, pointarray treenode, boolean usertree)
{
  /* allocate x dynamically */
  /* used in dnaml & dnamlk */
  long i, j, k;
  node *p;

  for (i = 0; i < spp; i++){
    treenode[i]->x = (phenotype)Malloc(endsite*sizeof(ratelike));
    for (j = 0; j < endsite; j++)
      treenode[i]->x[j]  = (ratelike)Malloc(rcategs*sizeof(sitelike));
  }
  if (!usertree) {
    for (i = spp; i < nonodes; i++) {
      p = treenode[i];
      for (j = 1; j <= 3; j++) {
        p->x = (phenotype)Malloc(endsite*sizeof(ratelike));
        for (k = 0; k < endsite; k++)
          p->x[k] = (ratelike)Malloc(rcategs*sizeof(sitelike));
        p = p->next;
      }
    }
  }
}  /* allocx */


void prot_allocx(long nonodes, long rcategs, pointarray treenode, 
                        boolean usertree)
{
  /* allocate x dynamically */
  /* used in proml            */
  long i, j, k;
  node *p;

  for (i = 0; i < spp; i++){
    treenode[i]->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
    for (j = 0; j < endsite; j++)
      treenode[i]->protx[j]  = (pratelike)Malloc(rcategs*sizeof(psitelike));
  }  
  if (!usertree) {
    for (i = spp; i < nonodes; i++) {
      p = treenode[i];
      for (j = 1; j <= 3; j++) {
        p->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
        for (k = 0; k < endsite; k++)
          p->protx[k] = (pratelike)Malloc(rcategs*sizeof(psitelike));
        p = p->next;
      } 
    }  
  } 
}  /* prot_allocx */


void allocx2(long nonodes, long endsite, long sitelength, pointarray treenode,
   boolean usertree)
{
  /* allocate x2 dynamically */
  /* used in restml */
  long i, j, k, l;
  node *p;

  for (i = 0; i < spp; i++)
    treenode[i]->x2 = (phenotype2)Malloc(endsite*sizeof(sitelike2));
  if (!usertree) {
    for (i = spp; i < nonodes; i++) {
      p = treenode[i];
      for (j = 1; j <= 3; j++) {
        p->x2 = (phenotype2)Malloc(endsite*sizeof(sitelike2));
        for (k = 0; k < endsite; k++)
          for (l = 0; l < sitelength; l++)
             p->x2[k][l] = 1.0;
        p = p->next;
      }
    }
  }
}  /* allocx2 */


void setuptree(pointarray treenode, long nonodes, boolean usertree)
{
  /* initialize treenodes */
  long i;
  node *p;

  for (i = 1; i <= nonodes; i++) {
    if (i <= spp || !usertree) {
      treenode[i-1]->back = NULL;
      treenode[i-1]->tip = (i <= spp);
      treenode[i-1]->index = i;
      treenode[i-1]->numdesc = 0;
      treenode[i-1]->iter = true;
      treenode[i-1]->initialized = true;
      treenode[i-1]->tyme =  0.0;
    }
  }
  if (!usertree) {
    for (i = spp + 1; i <= nonodes; i++) {
      p = treenode[i-1]->next;
      while (p != treenode[i-1]) {
        p->back = NULL;
        p->tip = false;
        p->index = i;
        p->numdesc = 0;
        p->iter = true;
        p->initialized = false;
        p->tyme = 0.0;
        p = p->next;
      }
    }
  }
} /* setuptree */


void setuptree2(tree a)
{
  /* initialize a tree */
  /* used in dnaml, dnamlk, & restml */

  a.likelihood = -999999.0;
  a.start = a.nodep[0]->back;
  a.root = NULL;
} /* setuptree2 */


void alloctip(node *p, long *zeros)
{ /* allocate a tip node */
  /* used by dnacomp, dnapars, & dnapenny */

  p->numsteps = (steptr)Malloc(endsite*sizeof(long));
  p->oldnumsteps = (steptr)Malloc(endsite*sizeof(long));
  p->base = (baseptr)Malloc(endsite*sizeof(long));
  p->oldbase = (baseptr)Malloc(endsite*sizeof(long));
  memcpy(p->base, zeros, endsite*sizeof(long));
  memcpy(p->numsteps, zeros, endsite*sizeof(long));
  memcpy(p->oldbase, zeros, endsite*sizeof(long));
  memcpy(p->oldnumsteps, zeros, endsite*sizeof(long));
}  /* alloctip */


void alloctrans(transptr *trans, long nonodes, long sitelength)
{
  /* used by restml */
  long i, j;

  *trans = (transptr)Malloc(nonodes*sizeof(transmatrix));
  for (i = 0; i < nonodes; ++i){
    (*trans)[i] = (transmatrix)Malloc((sitelength + 1) * sizeof(double *));
    for (j = 0;j < sitelength + 1; ++j)
      (*trans)[i][j] = (double *)Malloc((sitelength + 1) * sizeof(double));
  }
}  /* alloctrans */


void getbasefreqs(double freqa, double freqc, double freqg, double freqt,
            double *freqr, double *freqy, double *freqar, double *freqcy,
            double *freqgr, double *freqty, double *ttratio, double *xi,
            double *xv, double *fracchange, boolean freqsfrom,
            boolean printdata)
{
  /* used by dnadist, dnaml, & dnamlk */
  double aa, bb;

  if (printdata) {
    putc('\n', outfile);
    if (freqsfrom)
      fprintf(outfile, "Empirical ");
    fprintf(outfile, "Base Frequencies:\n\n");
    fprintf(outfile, "   A    %10.5f\n", freqa);
    fprintf(outfile, "   C    %10.5f\n", freqc);
    fprintf(outfile, "   G    %10.5f\n", freqg);
    fprintf(outfile, "  T(U)  %10.5f\n", freqt);
  }
  *freqr = freqa + freqg;
  *freqy = freqc + freqt;
  *freqar = freqa / *freqr;
  *freqcy = freqc / *freqy;
  *freqgr = freqg / *freqr;
  *freqty = freqt / *freqy;
  aa = *ttratio * (*freqr) * (*freqy) - freqa * freqg - freqc * freqt;
  bb = freqa * (*freqgr) + freqc * (*freqty);
  *xi = aa / (aa + bb);
  *xv = 1.0 - *xi;
  if (*xi < 0.0) {
    printf("\n WARNING: This transition/transversion ratio\n");
    printf(" is impossible with these base frequencies!\n");
    *xi = 0.0;
    *xv = 1.0;
    (*ttratio) = (freqa*freqg+freqc*freqt)/((*freqr)*(*freqy));
    printf(" Transition/transversion parameter reset\n");
    printf("  so transition/transversion ratio is %10.6f\n\n", (*ttratio));
  }
  if (freqa <= 0.0)
    freqa = 0.000001;
  if (freqc <= 0.0)
    freqc = 0.000001;
  if (freqg <= 0.0)
    freqg = 0.000001;
  if (freqt <= 0.0)
    freqt = 0.000001;
  *fracchange = (*xi) * (2 * freqa * (*freqgr) + 2 * freqc * (*freqty)) +
      (*xv) * (1.0 - freqa * freqa - freqc * freqc - freqg * freqg
      - freqt * freqt);
}  /* getbasefreqs */


void empiricalfreqs(double *freqa, double *freqc, double *freqg,
                        double *freqt, steptr weight, pointarray treenode)
{
  /* Get empirical base frequencies from the data */
  /* used in dnaml & dnamlk */
  long i, j, k;
  double sum, suma, sumc, sumg, sumt, w;

  *freqa = 0.25;
  *freqc = 0.25;
  *freqg = 0.25;
  *freqt = 0.25;
  for (k = 1; k <= 8; k++) {
    suma = 0.0;
    sumc = 0.0;
    sumg = 0.0;
    sumt = 0.0;
    for (i = 0; i < spp; i++) {
      for (j = 0; j < endsite; j++) {
        w = weight[j];