postprob.c

hmmer源程序

	  /* match state */
	  mmx[i][k]  = -INFTY;
	  if ((sc = mmx[i-1][k-1]) > mmx[i][k])
	    mmx[i][k] = sc;
	  if ((sc = imx[i-1][k-1]) > mmx[i][k])
	    mmx[i][k] = sc;
	  if ((sc = dmx[i-1][k-1]) > mmx[i][k])
	    mmx[i][k] = sc;
	  if ((sc = xmx[i-1][XMB]) > mmx[i][k])
	    mmx[i][k] = sc;
	  mmx[i][k] = ILogsum(mmx[i][k], posterior->mmx[i][k]);
	  
	  /* delete state */
	  dmx[i][k] = -INFTY;
	  if ((sc = mmx[i][k-1]) > dmx[i][k])
	    dmx[i][k] = sc;
	  if ((sc = dmx[i][k-1]) > dmx[i][k])
	    dmx[i][k] = sc;

	  /* insert state */
	  imx[i][k] = -INFTY;
	  if ((sc = mmx[i-1][k]) > imx[i][k])
	    imx[i][k] = sc;
	  if ((sc = imx[i-1][k]) > imx[i][k])
	    imx[i][k] = sc;
	  imx[i][k] = ILogsum(imx[i][k], posterior->imx[i][k]);
	}
      
      /* Now the special states. Order is important here.
       * remember, C and J emissions are zero score by definition,
       */

      /* N state */
      xmx[i][XMN] = -INFTY;
      if ((sc = ILogsum(xmx[i-1][XMN], posterior->xmx[i][XMN])) > -INFTY)
	xmx[i][XMN] = sc;
      
      /* E state */
      xmx[i][XME] = -INFTY;
      for (k = 1; k <= M; k++)
	if ((sc =  mmx[i][k]) > xmx[i][XME])
	  xmx[i][XME] = sc;
      
      /* J state */
      xmx[i][XMJ] = -INFTY;
      if ((sc = ILogsum(xmx[i-1][XMJ], posterior->xmx[i][XMJ])) > -INFTY)
	xmx[i][XMJ] = sc;
      if ((sc = xmx[i][XME]) > xmx[i][XMJ])      /* no E->J emission */
	xmx[i][XMJ] = sc;
      
      /* B state */
      xmx[i][XMB] = -INFTY;
      if ((sc = xmx[i][XMN]) > -INFTY)
	xmx[i][XMB] = sc;
      if ((sc = xmx[i][XMJ]) > xmx[i][XMB])
	xmx[i][XMB] = sc;
      
      /* C state */
      xmx[i][XMC] = -INFTY;
      if ((sc = ILogsum(xmx[i-1][XMC], posterior->xmx[i][XMC])) > -INFTY)
	xmx[i][XMC] = sc;
      if ((sc = xmx[i][XME]) > xmx[i][XMC])      /* no E->C emission */
	xmx[i][XMC] = sc;
    }

  /* T state (not stored) */
  sc = xmx[L][XMC];
  
  if (ret_tr != NULL) {
    P7OptimalAccuracyTrace(L, M, posterior, mx, &tr);
    *ret_tr = tr;
  }
  
  return Score2Prob(sc,1);	/* the log of the expected accuracy. */
}


/* Function: P7OptimalAccuracyTrace()
 * 
 * Purpose:  Traceback of an optimal accuracy matrix: i.e. retrieval 
 *           of optimum alignment.
 *           
 * Args:     L         - length of sequence
 *           M         - length of HMM
 *           posterior - the posterior matrix
 *           mx        - the matrix to trace back in, (L+1) x M
 *           ret_tr    - RETURN: traceback.
 *           
 * Return:   (void)
 *           ret_tr is allocated here. Free using P7FreeTrace().
 */
void
P7OptimalAccuracyTrace(int L,
		       int M,
		       struct dpmatrix_s *posterior,
		       struct dpmatrix_s *mx,
		       struct p7trace_s **ret_tr)
{
  struct p7trace_s *tr;
  int curralloc;		/* current allocated length of trace */
  int tpos;			/* position in trace */
  int i;			/* position in seq (1..L) */
  int k;			/* position in model (1..M) */
  int **xmx, **mmx, **imx, **dmx;
  int sc;			/* temp var for pre-emission score */

  /* Overallocate for the trace.
   * S-N-B- ... - E-C-T  : 6 states + L is minimum trace;
   * add L more as buffer.             
   */
  curralloc = L * 2 + 6; 
  P7AllocTrace(curralloc, &tr);

  xmx = mx->xmx;
  mmx = mx->mmx;
  imx = mx->imx;
  dmx = mx->dmx;

  /* Initialization of trace
   * We do it back to front; ReverseTrace() is called later.
   */
  tr->statetype[0] = STT;
  tr->nodeidx[0]   = 0;
  tr->pos[0]       = 0;
  tr->statetype[1] = STC;
  tr->nodeidx[1]   = 0;
  tr->pos[1]       = 0;
  tpos = 2;
  i    = L;			/* current i (seq pos) we're trying to assign */

  /* Traceback
   */
  while (tr->statetype[tpos-1] != STS) {
    switch (tr->statetype[tpos-1]) {
    case STM:			/* M connects from i-1,k-1, or B */
      sc = mmx[i+1][k+1];
      if (sc == ILogsum(mmx[i][k], posterior->mmx[i+1][k+1]) && i > 0 && k > 0)
	{
	  tr->statetype[tpos] = STM;
	  tr->nodeidx[tpos]   = k--;
	  tr->pos[tpos]       = i--;
	}
      else if (sc == ILogsum(imx[i][k], posterior->mmx[i+1][k+1]) && i > 0 && k > 0)
	{
	  tr->statetype[tpos] = STI;
	  tr->nodeidx[tpos]   = k;
	  tr->pos[tpos]       = i--;
	}
      else if (sc == ILogsum(dmx[i][k], posterior->mmx[i+1][k+1]) && i > 0 && k > 1)
	{
	  tr->statetype[tpos] = STD;
	  tr->nodeidx[tpos]   = k--;
	  tr->pos[tpos]       = 0;
	}
      else if (sc == ILogsum(xmx[i][XMB], posterior->mmx[i+1][k+1]))
	{
	  tr->statetype[tpos] = STB;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;
	}
      else Die("traceback failed");
      break;

    case STD:			/* D connects from M,D */
      if (dmx[i][k+1] == mmx[i][k] && i > 0 && k > 0)
	{
	  tr->statetype[tpos] = STM;
	  tr->nodeidx[tpos]   = k--;
	  tr->pos[tpos]       = i--;
	}
      else if (dmx[i][k+1] == dmx[i][k] && k > 1)
	{
	  tr->statetype[tpos] = STD;
	  tr->nodeidx[tpos]   = k--;
	  tr->pos[tpos]       = 0;
	}
      else Die("traceback failed");
      break;

    case STI:			/* I connects from M,I */
      sc = imx[i+1][k];
      if (sc == ILogsum(mmx[i][k], posterior->imx[i+1][k]) && i > 0 && k > 0)
	{
	  tr->statetype[tpos] = STM;
	  tr->nodeidx[tpos]   = k--;
	  tr->pos[tpos]       = i--;
	}
      else if (sc == ILogsum(imx[i][k], posterior->imx[i+1][k]) && i > 0 && k > 0)
	{
	  tr->statetype[tpos] = STI;
	  tr->nodeidx[tpos]   = k;
	  tr->pos[tpos]       = i--;
	}
      else Die("traceback failed");
      break;

    case STN:			/* N connects from S, N */
      if (i == 0 && xmx[i][XMN] == -INFTY)
	{
	  tr->statetype[tpos] = STS;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;
	}
      else if (i > 0 && xmx[i+1][XMN] == ILogsum(xmx[i][XMN], posterior->xmx[i+1][XMN]) && i > 0)
	{
	  tr->statetype[tpos] = STN;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;    /* note convention adherence:  */
	  tr->pos[tpos-1]     = i--;  /* first N doesn't emit        */
	}
      else Die("traceback failed");
      break;

    case STB:			/* B connects from N, J */
      if (xmx[i][XMB] == xmx[i][XMN])
	{
	  tr->statetype[tpos] = STN;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;
	}
      else if (xmx[i][XMB] == xmx[i][XMJ])
	{
	  tr->statetype[tpos] = STJ;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;
	}
      else Die("traceback failed");
      break;

    case STE:			/* E connects from any M state. k set here */
      for (k = M; k >= 1; k--)
	if (xmx[i][XME] == mmx[i][k] && i > 0)
	  {
	    tr->statetype[tpos] = STM;
	    tr->nodeidx[tpos]   = k--;
	    tr->pos[tpos]       = i--;
	    break;
	  }
      if (k <= 0) Die("traceback failed");
      break;

    case STC:			/* C comes from C, E */
      if (xmx[i][XMC] == ILogsum(xmx[i-1][XMC], posterior->xmx[i][XMC]) && i > 0)
	{
	  tr->statetype[tpos] = STC;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;    /* note convention adherence: */
	  tr->pos[tpos-1]     = i--;  /* first C doesn't emit       */
	}
      else if (xmx[i][XMC] == xmx[i][XME])
	{
	  tr->statetype[tpos] = STE;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0; /* E is a nonemitter */
	}
      else Die("Traceback failed.");
      break;

    case STJ:			/* J connects from E, J */
      if (xmx[i][XMJ] == ILogsum(xmx[i-1][XMJ], posterior->xmx[i][XMJ]) && i > 0)
	{
	  tr->statetype[tpos] = STJ;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0;    /* note convention adherence: */
	  tr->pos[tpos-1]     = i--;  /* first J doesn't emit       */
	}
      else if (xmx[i][XMJ] == xmx[i][XME])
	{
	  tr->statetype[tpos] = STE;
	  tr->nodeidx[tpos]   = 0;
	  tr->pos[tpos]       = 0; /* E is a nonemitter */
	}
      else Die("Traceback failed.");
      break;

    default:
      Die("traceback failed");

    } /* end switch over statetype[tpos-1] */
    
    tpos++;
    if (tpos == curralloc) 
      {				/* grow trace if necessary  */
	curralloc += L;
	P7ReallocTrace(tr, curralloc);
      }

  } /* end traceback, at S state; tpos == tlen now */
  tr->tlen = tpos;
  P7ReverseTrace(tr);
  *ret_tr = tr;

}


/* Function: PostalCode()
 * Date:     SRE, Sun Nov  7 15:31:35 1999 [Cold Spring Harbor]
 *
 * Purpose:  Given a traceback and one of Ian's posterior
 *           probability matrices, calculate a string that
 *           represents the confidence values on each 
 *           residue in the sequence.
 *           
 *           The code string is 0..L-1  (L = len of target seq),
 *           so it's in the coordinate system of the sequence string;
 *           off by one from dsq; and convertible to the coordinate
 *           system of aseq using MakeAlignedString().
 *           
 *           Values are 0-9,*  
 *           for example, 9 means with >=90% posterior probabiility,
 *           residue i is aligned to the state k that it
 *           is assigned to in the given trace.
 *
 * Args:     L  - length of seq
 *           mx - posterior prob matrix: see P7EmitterPosterior()
 *           tr - a traceback to get a Postal code string for.   
 *
 * Returns:  char * array of codes, 0..L-1
 *           Caller is responsible for free'ing it.
 */
static char
score2postcode(int sc)
{
  char i;
  i = (char) (Score2Prob(sc, 1.) * 10.);
  return ((i > 9) ? '*' : '0'+i);
}
char *
PostalCode(int L, struct dpmatrix_s *mx, struct p7trace_s *tr)
{
  int tpos;
  int i;
  int k;
  char *postcode;

  postcode = MallocOrDie((L+1) * sizeof(char)); 
  for (tpos = 0; tpos < tr->tlen; tpos++)
    {
      i = tr->pos[tpos];
      k = tr->nodeidx[tpos];
      if (i == 0) continue;

      switch (tr->statetype[tpos]) {
      case STM: postcode[i-1] = score2postcode(mx->mmx[i][k]);   break;
      case STI: postcode[i-1] = score2postcode(mx->imx[i][k]);   break;
      case STN:	postcode[i-1] = score2postcode(mx->xmx[i][XMN]); break;
      case STC: postcode[i-1] = score2postcode(mx->xmx[i][XMC]); break;
      case STJ:	postcode[i-1] = score2postcode(mx->xmx[i][XMJ]); break;
      }
    }
  postcode[L] = '\0';

  return postcode;
}