core_algorithms.c

hmmer源程序

          endlist[lpos]   = s3;
	}

      if (t2 == STN)
	{			/* if we see STN midpoint, we know the whole N-term is STN */
	  for (; s2 >= s1; s2--) {
	    kassign[s2] = 1;
	    tassign[s2] = STN;
	  }
	}
      if (t2 == STC)
	{			/* if we see STC midpoint, we know whole C-term is STC */
	  for (; s2 <= s3; s2++) {
	    kassign[s2] = hmm->M;
	    tassign[s2] = STC;
	  }
	}
    }

  /*****************************************************************
   * Construct a traceback structure from kassign/tassign by interpolating
   * necessary states. 
   * Trace allocation is as follows. We clearly need L emitting states.
   * We also need nonemitting states as follows:
   * STS,STN,STB,STE,STC,STT = 6
   * STD: count k2-k1-1 in kassign M->M's
   * Also, count N->M's and M->C's (potential wing unfoldings)...
   *   ...and be careful to check wing unfoldings when there aren't
   *      any emitting N or C flanks! (bugfix, 2.1.1b)
   *****************************************************************/ 

  tlen = L + 6;
  for (i = 1; i < L; i++)
    {
      if (tassign[i] == STM && tassign[i+1] == STM)
	tlen += kassign[i+1] - kassign[i] - 1;
      if (tassign[i] == STN && tassign[i+1] == STM)
	tlen += kassign[i+1] - 1;
      if (tassign[i] == STM && tassign[i+1] == STC)
	tlen += hmm->M - kassign[i];
    }
  if (tassign[1] == STM) tlen += kassign[1] - 1; 
  if (tassign[L] == STM) tlen += hmm->M - kassign[L];
  P7AllocTrace(tlen, &tr);

  tr->statetype[0] = STS;
  tr->nodeidx[0]   = 0;
  tr->pos[0]       = 0;
  tr->statetype[1] = STN;
  tr->nodeidx[1]   = 0;
  tr->pos[1]       = 0;
  tpos = 2;
  
  for (i = 1; i <= L; i++)
    {
      switch(tassign[i]) {
      case STM:
				/* check for first match state */
	if (tr->statetype[tpos-1] == STN) {
	  tr->statetype[tpos] = STB;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;
	  tpos++;
				/* check for wing unfolding */
	  if (Prob2Score(hmm->begin[kassign[i]], hmm->p1) + INTSCALE <= hmm->bsc[kassign[i]])
	    for (k = 1; k < kassign[i]; k++) {
	      tr->statetype[tpos] = STD;
	      tr->nodeidx[tpos]   = k;
	      tr->pos[tpos]       = 0;
	      tpos++;
	    }
	}
				/* do the match state itself */
	tr->statetype[tpos] = STM;
	tr->nodeidx[tpos]   = kassign[i];
	tr->pos[tpos]       = i;
	tpos++;
				/* do any deletes necessary 'til next match */
	if (i < L && tassign[i+1] == STM && kassign[i+1] - kassign[i] > 1)
	  for (k = kassign[i] + 1; k < kassign[i+1]; k++)
	    {
	      tr->statetype[tpos] = STD;
	      tr->nodeidx[tpos]   = k;
	      tr->pos[tpos]       = 0;
	      tpos++;
	    }
				/* check for last match state */
	if (i == L || tassign[i+1] == STC) {
				/* check for wing unfolding */
	  if (Prob2Score(hmm->end[kassign[i-1]], 1.) + INTSCALE <=  hmm->esc[kassign[i-1]])
	    for (k = kassign[i]+1; k <= hmm->M; k++)
	      {
		tr->statetype[tpos] = STD;
		tr->nodeidx[tpos]   = k;
		tr->pos[tpos]       = 0;
		tpos++;
	      }
				/* add on the end state */
	  tr->statetype[tpos] = STE;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0; 
	  tpos++;
				/* and a nonemitting C state */
	  tr->statetype[tpos] = STC;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0; 
	  tpos++;
	}
	break;
	
      case STI:
	tr->statetype[tpos] = STI;
	tr->nodeidx[tpos]   = kassign[i];
	tr->pos[tpos]       = i;
	tpos++;
	break;

      case STN:
	tr->statetype[tpos] = STN;
	tr->nodeidx[tpos]   = 0;
	tr->pos[tpos]       = i;
	tpos++;
	break;

      case STC:
	tr->statetype[tpos] = STC;
	tr->nodeidx[tpos]   = 0;
	tr->pos[tpos]       = i; 
	tpos++;
	break;

      default: Die("Bogus state %s", Statetype(tassign[i]));
      }
    }
				/* terminate the trace */
  tr->statetype[tpos] = STT;
  tr->nodeidx[tpos]   = 0;
  tr->pos[tpos]       = 0; 
  tr->tlen = tpos+1;

  *ret_tr = tr;

  free(kassign);
  free(tassign);
  free(startlist);
  free(endlist);
  return ret_sc;
}


/* Function: Plan7ESTViterbi()
 * 
 * Purpose:  Frameshift-tolerant alignment of protein model to cDNA EST.
 *           
 * 
 */
float
Plan7ESTViterbi(char *dsq, int L, struct plan7_s *hmm, struct dpmatrix_s **ret_mx)
{
  struct dpmatrix_s *mx;
  int **xmx;
  int **mmx;
  int **imx;
  int **dmx;
  int   i,k;
  int   sc;
  int   codon;
  
  /* Allocate a DP matrix with 0..L rows, 0..M+1 columns.
   */ 
  mx = AllocPlan7Matrix(L+1, hmm->M, &xmx, &mmx, &imx, &dmx);

  /* Initialization of the zero row (DNA sequence of length 0)
   * Note that xmx[i][stN] = 0 by definition for all i,
   *    and xmx[i][stT] = xmx[i][stC], so neither stN nor stT need
   *    to be calculated in DP matrices.
   */
  xmx[0][XMN] = 0;		                     /* S->N, p=1            */
  xmx[0][XMB] = hmm->xsc[XTN][MOVE];                 /* S->N->B, no N-tail   */
  xmx[0][XME] = xmx[0][XMC] = xmx[0][XMJ] = -INFTY;  /* need seq to get here */
  for (k = 0; k <= hmm->M; k++)
    mmx[0][k] = imx[0][k] = dmx[0][k] = -INFTY;      /* need seq to get here */

  /* Initialization of the first row (DNA sequence of length 1);
   * only N state can make this nucleotide.
   */
  xmx[1][XMN] = xmx[0][XMN] + hmm->xsc[XTN][LOOP];
  xmx[1][XMB] = xmx[1][XMN] + hmm->xsc[XTN][MOVE]; 
  xmx[0][XME] = xmx[0][XMC] = xmx[0][XMJ] = -INFTY;  /* need 2 nt to get here */
  for (k = 0; k <= hmm->M; k++)
    mmx[0][k] = imx[0][k] = dmx[0][k] = -INFTY;      /* need 2 nt to get into model */

  /* Recursion. Done as a pull.
   * Note some slightly wasteful boundary conditions:  
   *    tsc[0] = -INFTY for all eight transitions (no node 0)
   *    D_M and I_M are wastefully calculated (they don't exist)
   */
  for (i = 2; i <= L; i++) {
    mmx[i][0] = imx[i][0] = dmx[i][0] = -INFTY;

				/* crude calculation of lookup value for codon */
    if (i > 2) {
      if (dsq[i-2] < 4 && dsq[i-1] < 4 && dsq[i] < 4)
	codon = dsq[i-2] * 16 + dsq[i-1] * 4 + dsq[i];
      else
	codon = 64;		/* ambiguous codon; punt */
    }

    for (k = 1; k <= hmm->M; k++) {
				/* match state */
      if (i > 2) {
	mmx[i][k]  = mmx[i-3][k-1] + hmm->tsc[k-1][TMM];
	if ((sc = imx[i-3][k-1] + hmm->tsc[k-1][TIM]) > mmx[i][k])
	  mmx[i][k] = sc;
	if ((sc = xmx[i-3][XMB] + hmm->bsc[k]) > mmx[i][k])
	  mmx[i][k] = sc;
	if ((sc = dmx[i-3][k-1] + hmm->tsc[k-1][TDM]) > mmx[i][k])
	  mmx[i][k] = sc;
	mmx[i][k] += hmm->dnam[codon][k];
      }
				/* -1 frameshifts into match state */
      if ((sc = mmx[i-2][k-1] + hmm->tsc[k-1][TMM] + hmm->dna2) > mmx[i][k])
	mmx[i][k] = sc;
      if ((sc = imx[i-2][k-1] + hmm->tsc[k-1][TIM] + hmm->dna2) > mmx[i][k])
	mmx[i][k] = sc;
      if ((sc = xmx[i-2][XMB] + hmm->bsc[k] + hmm->dna2) > mmx[i][k])
	mmx[i][k] = sc;
      if ((sc = dmx[i-2][k-1] + hmm->tsc[k-1][TDM] + hmm->dna2) > mmx[i][k])
	mmx[i][k] = sc;
      
				/* +1 frameshifts into match state */
      if (i > 3) {
	if ((sc = mmx[i-4][k-1] + hmm->tsc[k-1][TMM] + hmm->dna4) > mmx[i][k])
	  mmx[i][k] = sc;
	if ((sc = imx[i-4][k-1] + hmm->tsc[k-1][TIM] + hmm->dna4) > mmx[i][k])
	  mmx[i][k] = sc;
	if ((sc = xmx[i-4][XMB] + hmm->bsc[k] + hmm->dna4) > mmx[i][k])
	  mmx[i][k] = sc;
	if ((sc = dmx[i-4][k-1] + hmm->tsc[k-1][TDM] + hmm->dna4) > mmx[i][k])
	  mmx[i][k] = sc;
      }
      				/* delete state */
      dmx[i][k]  = mmx[i][k-1] + hmm->tsc[k-1][TMD];
      if ((sc = dmx[i][k-1] + hmm->tsc[k-1][TDD]) > dmx[i][k])
	dmx[i][k] = sc;

				/* insert state */
      if (i > 2) {
	imx[i][k] = mmx[i-3][k] + hmm->tsc[k][TMI];
	if ((sc = imx[i-3][k] + hmm->tsc[k][TII]) > imx[i][k])
	  imx[i][k] = sc;
	imx[i][k] += hmm->dnai[codon][k];
      }

				/* -1 frameshifts into insert state */
      if ((sc = mmx[i-2][k] + hmm->tsc[k][TMI] + hmm->dna2) > imx[i][k])
	imx[i][k] = sc;
      if ((sc = imx[i-2][k] + hmm->tsc[k][TII] + hmm->dna2) > imx[i][k])
	imx[i][k] = sc;

				/* +1 frameshifts into insert state */
      if (i > 4) {
	if ((sc = mmx[i-4][k] + hmm->tsc[k][TMI] + hmm->dna4) > imx[i][k])
	  imx[i][k] = sc;
	if ((sc = imx[i-4][k] + hmm->tsc[k][TII] + hmm->dna4) > imx[i][k])
	  imx[i][k] = sc;
      }
    }
    /* Now the special states. Order is important here.
     * remember, C and J emissions are zero score by definition,
     */
				/* N state: +1 nucleotide */
    xmx[i][XMN] = xmx[i-1][XMN] + hmm->xsc[XTN][LOOP];
                                /* E state: collect from M's, and last D  */
    xmx[i][XME] = dmx[i][hmm->M];    /* transition prob from last D = 1.0 */
    for (k = 1; k <= hmm->M; k++)
      if ((sc =  mmx[i][k] + hmm->esc[k]) > xmx[i][XME])
        xmx[i][XME] = sc;
                                /* J state: +1 nucleotide */
    xmx[i][XMJ] = xmx[i-1][XMJ] + hmm->xsc[XTJ][LOOP];
    if ((sc = xmx[i][XME]   + hmm->xsc[XTE][LOOP]) > xmx[i][XMJ])
      xmx[i][XMJ] = sc;
                                /* B state: collect from N,J */
    xmx[i][XMB] = xmx[i][XMN] + hmm->xsc[XTN][MOVE];
    if ((sc = xmx[i][XMJ] + hmm->xsc[XTJ][MOVE]) > xmx[i][XMB])
      xmx[i][XMB] = sc;
				/* C state: +1 nucleotide */
    xmx[i][XMC] = xmx[i-1][XMC] + hmm->xsc[XTC][LOOP];
    if ((sc = xmx[i][XME] + hmm->xsc[XTE][MOVE]) > xmx[i][XMC])
      xmx[i][XMC] = sc;
  }

  sc = xmx[L][XMC] + hmm->xsc[XTC][MOVE];

  if (ret_mx != NULL) *ret_mx = mx;
  else                FreePlan7Matrix(mx);

  return Scorify(sc);            /* the total Viterbi score. */
}



/* Function: get_wee_midpt()
 * Date:     SRE, Wed Mar  4 08:27:11 1998 [St. Louis]
 *
 * Purpose:  The heart of the divide and conquer algorithm
 *           for P7WeeViterbi(). This function is called
 *           recursively to find successive optimal midpoints
 *           in the alignment matrix. See P7WeeViterbi() for
 *           further comments on the assumptions of this algorithm.
 *           
 * Args:     hmm   - the model, set up for integer scores
 *           dsq   - the sequence, digitized
 *           L     - length of the sequence
 *           k1    - model node to start with, 1..M
 *           t1    - state type to start with, STM | STI | STN | STC; STS to start
 *           s1    - sequence position to start with, 1..L; 1 to start
 *           k3    - model node to end with, 1..M
 *           t3    - state type to end with, STM | STI | STN | STC; STT to start
 *           s3    - sequence position to end with, 1..L; L to start
 *          ret_k2 - RETURN: optimal midpoint, node position in model
 *          ret_t2 - RETURN: optimal midpoint, state type 
 *          ret_s2 - RETURN: optimal midpoint, sequence position
 *
 * Returns: score of optimal alignment, in bits. 
 */
static float
get_wee_midpt(struct plan7_s *hmm, char *dsq, int L, 
	      int k1, char t1, int s1,
	      int k3, char t3, int s3,
	      int *ret_k2, char *ret_t2, int *ret_s2)
{
  struct dpmatrix_s *fwd;
  struct dpmatrix_s *bck;
  int        **xmx;             /* convenience ptr into special states */
  int        **mmx;             /* convenience ptr into match states   */
  int        **imx;             /* convenience ptr into insert states  */
  int        **dmx;             /* convenience ptr into delete states  */
  int          k2;
  char         t2;
  int          s2;
  int          cur, prv, nxt;	/* current, previous, next row index (0 or 1)*/
  int          i,k;		/* indices for seq, model */
  int          sc;		/* integer score */
  int          max;		/* maximum integer score */
  int          start;		/* s1 to start at (need, for STS special case) */

 
  /* Choose our midpoint.
   * Special cases: s1, s3 adjacent and t1 == STS: s2 = s1
   *                s1, s3 adjacent and t3 == STT: s2 = s3
   *                (where we must replace STS, STT eventually)
   */
  s2 = s1 + (s3-s1) / 2;
  if (s3-s1 == 1 && t1 == STS) s2 = s1;
  if (s3-s1 == 1 && t3 == STT) s2 = s3;

  /* STS is a special case. STS aligns to row zero by convention,
   * but we'll be passed s1=1, t1=STS. We have to init on row
   * zero then start DP on row 1.
   */
  start = (t1 == STS) ? 0 : s1;

  /* Allocate our forward two rows.
   * Initialize row zero.
   */
  fwd = AllocPlan7Matrix(2, hmm->M, &xmx, &mmx, &imx, &dmx);
  cur = start%2;
  xmx[cur][XMN] = xmx[cur][XMB] = -INFTY;
  xmx[cur][XME] = xmx[cur][XMC] = -INFTY;  
  for (k = k1; k <= k3; k++)
    mmx[cur][k] = imx[cur][k] = dmx[cur][k] = -INFTY;      

  /* Where to put our zero for our start point...
   * (only possible to start on an emitting state; J disallowed)
   */
  switch (t1) {
  case STM: mmx[cur][k1]  = 0; break;
  case STI: imx[cur][k1]  = 0; break;
  case STN: xmx[cur][XMN] = 0; break;
  case STC: xmx[cur][XMC] = 0; break;
  case STS: xmx[cur][XMN] = 0; break;
  default:  Die("you can't init get_wee_midpt with a %s\n", Statetype(t1));
  }

  /* Still initializing.
   * Deal with pulling horizontal matrix moves in initial row.
   * These are any transitions to nonemitters:
   *    STM-> E, D
   *    STI-> none
   *    STN-> B
   *    STC-> (T, but we never observe this in the forward pass of a d&c)
   *    STE-> C