scaleswn.c

序列对齐 Compare a protein sequence to a protein sequence database or a DNA sequence to a DNA sequenc

  double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0;
  double s, l, es;
  int i;

  for(i = start ; i < stop ; i++) {
    s = (double)sptr[i].score;
    l = (double)sptr[i].n1;
    es = exp(-lambda * s);
    sum2 += l * es;
    sum1 += l * s * es;
    sum3 += l * s * s * es;
  }

  sum1 += sumlenL*cenL*exp(-lambda*cenL) - sumlenH*cenH*exp(-lambda*cenH);
  sum2 += sumlenL*exp(-lambda * cenL) -  sumlenH*exp(-lambda * cenH);
  sum3 += sumlenL*cenL*cenL * exp(-lambda * cenL) -
    sumlenH*cenH*cenH * exp(-lambda * cenH);
  return ((sum1 * sum1) / (sum2 * sum2)) - (sum3 / sum2)
    - (1.0 / (lambda * lambda));
}

double mle2_func(double *params,
		 double *consts,
		 struct stat_str *values,
		 int n, int start, int stop);

void simplex(double *fitparams,
	     double *lambda,
	     int nparam,
	     double (*minfunc) (double *tryparams, double *consts, 
				struct stat_str *data, int ndata,
				int start, int stop),
	     double *consts,
	     void *data,
	     int ndata, int start, int stop
	     );

int
mle_cen2(struct stat_str *sptr, int n, int M, double fc,
	double *a0, double *a1, double *a2, double *b1) {

  double params[4], lambdas[4], consts[9];
  double avglenL, avglenH, avgcompL, avgcompH, cenL, cenH;
  int start, stop;
  int i, nf;

  nf = (fc/2.0) * n;
  start = nf;
  stop = n - nf;

  st_sort(sptr,n);

  /* choose arithmetic or geometic mean for compositions by appropriate commenting */

  if (nf > 0) {
    avglenL = avglenH = 0.0;
    avgcompL = avgcompH = 0.0;
    /* avgcompL = avgcompH = 1.0 */
    for (i=0; i<start; i++) {
      avglenL += (double)sptr[i].n1;
      avgcompL += (double)sptr[i].comp;
      /* avgcompL *= (double) sptr[i].comp; */
    }
    avglenL /= (double) start;
    avgcompL /= (double) start;
    /* avgcompL = pow(avgcompL, 1.0/(double) start); */
  
    for (i=stop; i<n; i++) {
      avglenH += (double)sptr[i].n1;
      avgcompH += (double)sptr[i].comp;
      /* avgcompH *= (double) sptr[i].comp; */
    }
    avglenH /= (double) (n - stop);
    avgcompH /= (double) (n - stop);
    /* avgcompL = pow(avgcompL, 1.0/(double) (n - stop)); */

    cenL = (double)sptr[start].score;
    cenH = (double)sptr[stop].score;
    if (cenL >= cenH) return -1;
  }
  else {
    avglenL = avglenH = cenL = cenH = 0.0;
    avgcompL = avgcompH = 1.0;
  }

  params[0] = 10.0;
  params[1] = -10.0;
  params[2] = 1.0;
  params[3] = 1.0;

  lambdas[0] = 1.0;
  lambdas[1] = 0.5;
  lambdas[2] = 0.1;
  lambdas[3] = 0.01;

  consts[0] = M;
  consts[1] = (double) start;
  consts[2] = (double) stop;
  consts[3] = cenL;
  consts[4] = cenH;
  consts[5] = avglenL;
  consts[6] = avglenH;
  consts[7] = avgcompL;
  consts[8] = avgcompH;

  simplex(params, lambdas, 4,
	  (double (*) (double *, double *, struct stat_str *, int, int, int) )mle2_func,
	  consts, sptr, n, start, stop);

  *a0 = params[0];
  *a1 = params[1];
  *a2 = params[2];
  *b1 = params[3];

  return 0;
}

double mle2_func(double *params,
		 double *consts,
		 struct stat_str *values,
		 int n, int start, int stop
		 ) {

  double a0, a1, a2, b1, M;
  double score, length, comp;
  double cenL, cenH, avglenL, avglenH, avgcompL, avgcompH;
  double L, y;

  int i;

  a0 = params[0];
  a1 = params[1];
  a2 = params[2];
  b1 = params[3];

  M = consts[0];
  /*
  start = (int) consts[1];
  stop = (int) consts[2];
  */
  cenL = consts[3];
  cenH = consts[4];
  avglenL = consts[5];
  avglenH = consts[6];
  avgcompL = consts[7];
  avgcompH = consts[8];

  L = 0;
  y = 0;

  if (start > 0) {
    y = -(cenL - (a0 + a1*avgcompL +a2*avgcompL*log(M*avglenL)))/(b1*avgcompL);
    L += (double) start * exp(y);
  }

  for(i = start ; i < stop ; i++) {
    score = (double) values[i].score;
    length = (double) values[i].n1;
    comp = (double) values[i].comp;

    y = - (score - (a0 + a1*comp + a2 * comp * log(M*length))) / (b1*comp);

    L += -y + exp(y) + log(b1 * comp);
  }

  if (stop < n) {
    y = -(cenH -(a0 + a1*avgcompH + a2*avgcompH*log(M*avglenH)))/(b1*avgcompH);
    L -= (double) (n - stop) * exp(y);
  }
  return L;
}

/* Begin Nelder-Mead simplex code: */

double evalfunc(double **param,
		double *vals,
		double *psums,
		double *ptry,
		int nparam,
		double (*minfunc) (double *params, double *consts,
				   struct stat_str *data, int ndata,
				   int start, int stop),
		double *consts,
		void *data,
		int ndata, int start, int stop,
		int ihi,
		double factor);

void simplex(double *fitparams,
	     double *lambda,
	     int nparam,
	     double (*minfunc) (double *tryparams, double *consts,
				struct stat_str *data, int ndata, 
				int start, int stop),
	     double *consts,
	     void *data,
	     int ndata,
	     int start,
	     int stop
	     )
{

  int i, j, ilo, ihi, inhi;
  double rtol, sum, tmp, ysave, ytry;
  double *psum, *vals, *ptry, **param;


  psum = (double *) calloc(nparam, sizeof(double));
  ptry = (double *) calloc(nparam, sizeof(double));

  vals = (double *) calloc(nparam + 1, sizeof(double));

  param = (double **) calloc(nparam + 1, sizeof(double *));
  param[0] = (double *) calloc((nparam + 1) * nparam, sizeof(double));
  for( i = 1 ; i < (nparam + 1) ; i++ ) {
    param[i] = param[0] + i * nparam;
  }

  /* Get our N+1 initial parameter values for the simplex */

  for( i = 0 ; i < nparam ; i++ ) {
    param[0][i] = fitparams[i];
  }

  for( i = 1 ; i < (nparam + 1) ; i++ ) {
    for( j = 0 ; j < nparam ; j++ ) {
      param[i][j] = fitparams[j] + lambda[j] * ( (i - 1) == j ? 1 : 0 );
    }
  }

  /* calculate initial values at the simplex nodes */

  for( i = 0 ; i < (nparam + 1) ; i++ ) {
    vals[i] = minfunc(param[i], consts, data, ndata, start, stop);
  }

  /* Begin Nelder-Mead simplex algorithm from Numerical Recipes in C */

  for( j = 0 ; j < nparam ; j++ ) {
    for( sum = 0.0, i = 0 ; i < nparam + 1 ; i++ ) {
      sum += param[i][j];
    }
    psum[j] = sum;
  }


  while( 1 ) {
/*
      determine which point is highest (ihi), next highest (inhi) and
      lowest (ilo) by looping over the points in the simplex
*/
    ilo = 0;

/*  ihi = vals[0] > vals[1] ? (inhi = 1, 0) : (inhi = 0, 1); */
    if(vals[0] > vals[1]) { ihi = 0; inhi = 1; }
    else { ihi = 1; inhi = 0; }

    for( i = 0 ; i < nparam + 1 ; i++) {
      if( vals[i] <= vals[ilo] ) ilo = i;
      if( vals[i] > vals[ihi] ) {
	inhi = ihi;
	ihi = i;
      } else if ( vals[i] > vals[inhi] && i != ihi ) inhi = i;
    }

    /* Are we finished? */

    rtol = 2.0 * fabs(vals[ihi] - vals[ilo]) / 
      (fabs(vals[ihi]) + fabs(vals[ilo]) + TINY);

    if( rtol < TOLERANCE ) {

/* put the best value and best parameters into the first index */

      tmp = vals[0];
      vals[0] = vals[ilo];
      vals[ilo] = tmp;

      for( i = 0 ; i < nparam ; i++ ) {
	tmp = param[0][i];
	param[0][i] = param[ilo][i];
	param[ilo][i] = tmp;
      }

      /* et voila, c'est finis */
      break;
    }

    /* Begin a new iteration */

    /* first, extrapolate by -1 through the face of the simplex across from ihi */

    ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts,
		    data, ndata, start, stop, ihi, -1.0);

    if( ytry <= vals[ilo] ) {

      /* Good result, try additional extrapolation by 2 */

      ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts, 
		      data, ndata, start, stop, ihi, 2.0);

    } else if ( ytry >= vals[inhi] ) {

      /* no good, look for an intermediate lower point by contracting */

      ysave = vals[ihi];
      ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts,
		      data, ndata, start, stop, ihi, 0.5);

      if( ytry >= ysave ) {

	/* Still no good.  Contract around lowest (best) point. */

	for( i = 0 ; i < nparam + 1 ; i++ ) {
	  if( i != ilo ) {
	    for ( j = 0 ; j < nparam ; j++ ) {
	      param[i][j] = psum[j] = 0.5 * (param[i][j] + param[ilo][j]);
	    }
	    vals[i] = minfunc(psum, consts, data, ndata, start, stop);
	  }
	}


	for( j = 0 ; j < nparam ; j++ ) {
	  for( sum = 0.0, i = 0 ; i < nparam + 1 ; i++ ) {
	    sum += param[i][j];
	  }
	  psum[j] = sum;
	}

      }
    }
  }
			   
  for( i = 0 ; i < nparam ; i++ ) {
    fitparams[i] = param[0][i];
  }

  if (ptry!=NULL) {
    free(ptry);
    ptry=NULL;
  }
  free(param[0]);
  free(param);
  free(vals);
  free(psum);
}


double evalfunc(double **param,
		double *vals,
		double *psum,
		double *ptry,
		int nparam,
		double (*minfunc)(double *tryparam, double *consts,
				  struct stat_str *data, int ndata,
				  int start, int stop),
		double *consts,
		void *data,
		int ndata, int start, int stop,
		int ihi,
		double factor) {

  int j;
  double fac1, fac2, ytry;


  fac1 = (1.0 - factor) / nparam;
  fac2 = fac1 - factor;

  for( j = 0 ; j < nparam ; j++ ) {
    ptry[j] = psum[j] * fac1 - param[ihi][j] * fac2;
  }

  ytry = minfunc(ptry, consts, data, ndata, start, stop);

  if( ytry < vals[ihi] ) {
    vals[ihi] = ytry;
    for( j = 0 ; j < nparam ; j++ ) {
      psum[j] += ptry[j] - param[ihi][j];
      param[ihi][j] = ptry[j];
    }
  }

  return ytry;
}

/* end of Nelder-Mead simplex code */

int
proc_hist_n(struct stat_str *sptr, int nstats,
	    struct pstruct *ppst, struct hist_str *histp,
	    int do_trim, struct pstat_str *pu)
{
  int i, j, t_j;
  double s_score, s2_score, ssd, ztrim, mean, var;
  double t_score;
  int max_trim, min_trim;
  int nit, max_hscore;
  char s_string[128];
  char *f_string;

  f_string = histp->stat_info;

  max_hscore = calc_thresh(ppst, sptr, nstats, pu->ngLambda,
			   pu->ngK, pu->ngH, &ztrim);

  s_score = s2_score = 0.0;

  /* calculate mean */
  for (j=0, i=0; i < nstats; i++) {
    if ( sptr[i].score <= max_hscore) {
      s_score += (double)sptr[i].score;
      j++;
    }
    else { sptr[i].n1 = -sptr[i].n1;}
  }

  if (j > 1) {
    pu->r_u.rg.mu = mean = s_score/(double)j;

    /* calculate variance */
    for (i=0; i < nstats; i++) {
      if (sptr[i].n1 < 0) continue;
      ssd = (double)sptr[i].score - mean;
      s2_score += ssd * ssd;
    }
    
    var = pu->r_u.rg.mean_var = s2_score/(double)(j-1);
    pu->r_u.rg.mean_var_sqrt = sqrt(var);
  }
  else {
    pu->r_u.rg.mu = 50.0;
    pu->r_u.rg.mean_var = 10.0;
    pu->r_u.rg.mean_var_sqrt = sqrt(10.0);
  }
  
  if (pu->r_u.rg.mean_var < 0.01) {
    pu->r_u.rg.mean_var = (pu->r_u.rg.mu > 1.0) ? pu->r_u.rg.mu: 1.0;
  }

  /* now remove some scores */

  nit = 5;
  pu->r_u.rg.n_trimmed = 0;
  max_trim = -BIGNUM;
  min_trim = BIGNUM;
  while (nit-- > 0) {
    t_score = 0.0;
    t_j = 0;
    for (i=0; i< nstats; i++) {
      if (sptr[i].n1 < 0) continue;
      ssd = find_zn(sptr[i].score,sptr[i].escore,sptr[i].n1,sptr[i].comp, pu);
      if (ssd > ztrim 
#ifndef NORMAL_DIST
	  || ssd < 20.0
#else
	  || ssd < 0.0
#endif
	  )
	{
	/*      fprintf(stderr,"removing %3d %3d %4.1f\n",
		sptr[i].score, sptr[i].n1,ssd); */

	  ssd = sptr[i].score;
	  if (ssd > max_trim) max_trim = ssd;
	  if (ssd < min_trim) min_trim = ssd;

	  t_score += (double)ssd;
	  t_j++;
	  pu->r_u.rg.n_trimmed++;
	  histp->entries--;
	  sptr[i].n1 = -sptr[i].n1;
	}
    }

    if (j - t_j > 1 ) {
      mean = pu->r_u.rg.mu = (s_score - t_score)/(double)(j - t_j);

    /* calculate variance */
      s2_score = 0.0;
      for (i=0; i < nstats; i++) {
	if ( sptr[i].n1 > 0 && sptr[i].score <= max_hscore) {
	  ssd = (double)sptr[i].score - mean;
	  s2_score += ssd * ssd;
	}
      }