trace.c

hmmer源程序

  int   i;

  buffer[ALILENGTH] = '\0';
  endi = ali->sqfrom - 1;
  for (pos = 0; pos < ali->len; pos += ALILENGTH)
    {
				/* coords of target seq for this line */
      starti = endi + 1;
      for (i = pos; ali->aseq[i] != '\0' && i < pos + ALILENGTH; i++)
	if (!isgap(ali->aseq[i])) endi++;
	    
      if (ali->csline != NULL) {
	strncpy(buffer, ali->csline+pos, ALILENGTH);
	fprintf(fp, "  %16s %s\n", "CS", buffer);
      }
      if (ali->rfline != NULL) {
	strncpy(buffer, ali->rfline+pos, ALILENGTH);
	fprintf(fp, "  %16s %s\n", "RF", buffer);
      }
      if (ali->model  != NULL) {
	strncpy(buffer, ali->model+pos, ALILENGTH);
	fprintf(fp, "  %16s %s\n", " ", buffer);
      }
      if (ali->mline  != NULL) {
	strncpy(buffer, ali->mline+pos, ALILENGTH);
	fprintf(fp, "  %16s %s\n", " ", buffer);
      }
      if (ali->aseq   != NULL) { 
	strncpy(buffer, ali->aseq+pos, ALILENGTH);
	if (endi >= starti)
	  fprintf(fp, "  %10.10s %5d %s %-5d\n\n", ali->target, starti, buffer, endi);
	else
	  fprintf(fp, "  %10.10s %5s %s %-5s\n\n", ali->target, "-", buffer, "-"); 
      }
    }

  /* Cleanup and return
   */
  fflush(fp);
  return;
} 



/* Function: TraceDecompose()
 * Date:     Sat Aug 30 11:18:40 1997 (Denver CO)
 * 
 * Purpose:  Decompose a long multi-hit trace into zero or more
 *           traces without N,C,J transitions: for consistent
 *           scoring and statistical evaluation of single domain
 *           hits.
 *           
 * Args:     otr    - original trace structure
 *           ret_tr - RETURN: array of simpler traces        
 *           ret_ntr- RETURN: number of traces.
 *           
 * Return:   (void)
 *           ret_tr alloc'ed here; free individuals with FreeTrace().
 */            
void
TraceDecompose(struct p7trace_s *otr, struct p7trace_s ***ret_tr, int *ret_ntr)
{
  struct p7trace_s **tr;        /* array of new traces          */
  int ntr;			/* number of traces             */
  int i,j;			/* position counters in traces  */
  int idx;			/* index over ntr subtraces     */

  /* First pass: count begin states to get ntr.
   */
  for (ntr = 0, i = 0; i < otr->tlen; i++)
    if (otr->statetype[i] == STB) ntr++;

  /* Allocations.
   */
  if (ntr == 0) {
    *ret_ntr = 0;
    *ret_tr  = NULL;
    return;
  }
  tr = (struct p7trace_s **) MallocOrDie (sizeof(struct p7trace_s *) * ntr);

  for (idx = 0, i = 0; i < otr->tlen; i++) /* i = position in old trace */
    if (otr->statetype[i] == STB)
      {
	for (j = i+1; j < otr->tlen; j++) /* j = tmp; get length of subtrace */
	  if (otr->statetype[j] == STE) break;
			/* trace = S-N-(B..E)-C-T : len + 4 : j-i+1 + 4*/
	P7AllocTrace(j-i+5, &(tr[idx]));
	tr[idx]->tlen = j-i+5;

	tr[idx]->statetype[0] = STS;
	tr[idx]->nodeidx[0]   = 0;
	tr[idx]->pos[0]       = 0;
	tr[idx]->statetype[1] = STN;
	tr[idx]->nodeidx[1]   = 0;
	tr[idx]->pos[1]       = 0;
	j = 2;			/* now j = position in new subtrace */
	while (1)		/* copy subtrace */
	  {
	    tr[idx]->statetype[j] = otr->statetype[i];
	    tr[idx]->nodeidx[j]   = otr->nodeidx[i];
	    tr[idx]->pos[j]       = otr->pos[i];
	    if (otr->statetype[i] == STE) break;
	    i++; j++;
	  }
	j++;
	tr[idx]->statetype[j] = STC;
	tr[idx]->nodeidx[j]   = 0;
	tr[idx]->pos[j]       = 0;
	j++;
	tr[idx]->statetype[j] = STT;
	tr[idx]->nodeidx[j]   = 0;
	tr[idx]->pos[j]       = 0;
	idx++;
      }

  *ret_tr  = tr;
  *ret_ntr = ntr;
  return;
}


/* Function: TraceDomainNumber()
 * 
 * Purpose:  Count how many times we traverse the
 *           model in a single Plan7 trace -- equivalent
 *           to counting the number of domains.
 *           
 *           (A weakness is that we might discard some of
 *           those domains because they have low scores
 *           below E or T threshold.)
 */
int
TraceDomainNumber(struct p7trace_s *tr)
{
  int i;
  int ndom = 0;
  
  for (i = 0; i < tr->tlen; i++)
    if (tr->statetype[i] == STB) ndom++;
  return ndom;
}


/* Function: TraceSimpleBounds()
 * 
 * Purpose:  For a trace that contains only a single
 *           traverse of the model (i.e. something that's
 *           come from TraceDecompose(), or a global
 *           alignment), determine the bounds of
 *           the match on both the sequence [1..L] and the
 *           model [1..M].
 *           
 * Args:     tr   - trace to look at
 *           i1   - RETURN: start point in sequence [1..L]
 *           i2   - RETURN: end point in sequence [1..L]
 *           k1   - RETURN: start point in model [1..M]
 *           k2   - RETURN: end point in model [1..M]
 */
void
TraceSimpleBounds(struct p7trace_s *tr, int *ret_i1, int *ret_i2, 
		  int *ret_k1,  int *ret_k2)
{
  int i1, i2, k1, k2, tpos;

  i1 = k1 = i2 = k2 = -1;

				/* Look forwards to find start of match */
  for (tpos = 0; tpos < tr->tlen; tpos++)
    {
      if (k1 == -1 && (tr->statetype[tpos] == STM || tr->statetype[tpos] == STD))
	k1 = tr->nodeidx[tpos];
      if (tr->statetype[tpos] == STM)
	{
	  i1 = tr->pos[tpos];
	  break;
	}
    }
  if (tpos == tr->tlen || i1 == -1 || k1 == -1)
    Die("sanity check failed: didn't find a match state in trace");

				/* Look backwards to find end of match */
  for (tpos = tr->tlen-1; tpos >= 0; tpos--)
    {
      if (k2 == -1 && (tr->statetype[tpos] == STM || tr->statetype[tpos] == STD))
	k2 = tr->nodeidx[tpos];
      if (tr->statetype[tpos] == STM)
	{
	  i2 = tr->pos[tpos];
	  break;
	}
    }
  if (tpos == tr->tlen || i2 == -1 || k2 == -1)
    Die("sanity check failed: didn't find a match state in trace");

  *ret_k1 = k1;
  *ret_i1 = i1;
  *ret_k2 = k2;
  *ret_i2 = i2;
}


/* Function: MasterTraceFromMap()
 * Date:     SRE, Tue Jul  7 18:51:11 1998 [St. Louis]
 *
 * Purpose:  Convert an alignment map (e.g. hmm->map) to
 *           a master trace. Used for mapping an alignment
 *           onto an HMM. Generally precedes a call to
 *           ImposeMasterTrace(). Compare P7ViterbiAlignAlignment(),
 *           which aligns an alignment to the model using a
 *           Viterbi algorithm to get a master trace. 
 *           MasterTraceFromMap() only works if the alignment
 *           is exactly the one used to train the model.
 *
 * Args:     map  - the map (usually hmm->map is passed) 1..M
 *           M    - length of map (model; usually hmm->M passed)
 *           alen - length of alignment that map refers to
 *
 * Returns:  ptr to master trace
 *           Caller must free: P7FreeTrace().
 */
struct p7trace_s *
MasterTraceFromMap(int *map, int M, int alen)
{
  struct p7trace_s *tr;         /* RETURN: master trace */ 
  int tpos;			/* position in trace */
  int apos;			/* position in alignment, 1..alen */
  int k;			/* position in model */

  /* Allocate for the trace.
   * S-N-B- ... - E-C-T  : 6 states + alen is maximum trace,
   * because each of alen columns is an N*, M*, I*, or C* metastate.
   * No D* metastates possible.
   */
  P7AllocTrace(alen+6, &tr);

  /* Initialize the trace
   */
  tpos = 0;
  TraceSet(tr, tpos, STS, 0, 0); tpos++;
  TraceSet(tr, tpos, STN, 0, 0); tpos++;

  /* Leading N's
   */
  for (apos = 1; apos < map[1]; apos++) {
    TraceSet(tr, tpos, STN, 0, apos); tpos++;
  } /* now apos == map[1] */
  TraceSet(tr, tpos, STB, 0, 0); tpos++;

  for (k = 1; k < M; k++)
    {
      TraceSet(tr, tpos, STM, k, apos); tpos++;
      apos++;

      for (; apos < map[k+1]; apos++) {
	TraceSet(tr, tpos, STI, k, apos); tpos++;
      }
    } /* now apos == map[M] and k == M*/
      
  TraceSet(tr, tpos, STM, M, apos); tpos++;
  apos++;

  /* Trailing C's
   */
  TraceSet(tr, tpos, STE, 0, 0); tpos++;
  TraceSet(tr, tpos, STC, 0, 0); tpos++;
  for (; apos <= alen; apos++) {
    TraceSet(tr, tpos, STC, 0, apos); tpos++;
  }

  /* Terminate and return
   */
  TraceSet(tr, tpos, STT, 0, 0); tpos++;
  tr->tlen = tpos;
  return tr;
}



/* Function: ImposeMasterTrace()
 * Date:     SRE, Sun Jul  5 14:27:16 1998 [St. Louis]
 *
 * Purpose:  Goes with P7ViterbiAlignAlignment(), which gives us
 *           a "master trace" for a whole alignment. Now, given
 *           the alignment and the master trace, construct individual
 *           tracebacks for each sequence. Later we'll hand these
 *           (and presumably other traces) to P7Traces2Alignment().
 *           
 *           It is possible to generate individual traces that
 *           are not consistent with Plan7 (e.g. D->I and I->D 
 *           transitions may be present). P7Traces2Alignment()
 *           can handle such traces; other functions may not.
 *           See modelmaker.c:trace_doctor() if this is a problem.
 * 
 *           Akin to modelmaker.c:fake_tracebacks().
 *
 * Args:     aseq  - aligned seqs
 *           nseq  - number of aligned seqs 
 *           mtr   - master traceback
 *           ret_tr- RETURN: array of individual tracebacks, one for each aseq
 *
 * Returns:  (void)
 */
void
ImposeMasterTrace(char **aseq, int nseq, struct p7trace_s *mtr, struct p7trace_s ***ret_tr)
{
  struct p7trace_s **tr;
  int  idx;			/* counter over sequences */
  int  i;                       /* position in raw sequence (1..L) */
  int  tpos;			/* position in traceback           */
  int  mpos;			/* position in master trace        */

  tr = (struct p7trace_s **) MallocOrDie (sizeof(struct p7trace_s *) * nseq);
  
  for (idx = 0; idx < nseq; idx++)
    {
      P7AllocTrace(mtr->tlen, &tr[idx]); /* we're guaranteed that individuals len < master len */
      
      tpos = 0;
      i    = 1;
      for (mpos = 0; mpos < mtr->tlen; mpos++)
	{
	  switch (mtr->statetype[mpos]) 
	    {
	    case STS:		/* straight copies w/ no emission: S, B, D, E, T*/
	    case STB:
	    case STD:
	    case STE:
	    case STT:
	      TraceSet(tr[idx], tpos, mtr->statetype[mpos], mtr->nodeidx[mpos], 0);
	      tpos++;
	      break;

	    case STM:		/* M* implies M or D */
	      if (isgap(aseq[idx][mtr->pos[mpos]-1])) 
		TraceSet(tr[idx], tpos, STD, mtr->nodeidx[mpos], 0);
	      else {
		TraceSet(tr[idx], tpos, STM, mtr->nodeidx[mpos], i);
		i++;
	      }
	      tpos++;
	      break;

	    case STI:		/* I* implies I or nothing */
	      if (!isgap(aseq[idx][mtr->pos[mpos]-1])) {
		TraceSet(tr[idx], tpos, STI, mtr->nodeidx[mpos], i);
		i++;
		tpos++;
	      }
	      break;

	    case STJ:		/* N,J,C: first N* -> N. After that, N* -> N or nothing. */
	    case STN:
	    case STC:
	      if (mtr->pos[mpos] == 0) { 
		TraceSet(tr[idx], tpos, mtr->statetype[mpos], 0, 0);
		tpos++;
	      } else if (!isgap(aseq[idx][mtr->pos[mpos]-1])) {
		TraceSet(tr[idx], tpos, mtr->statetype[mpos], 0, i);
		i++;
		tpos++; 
	      }
	      break;

	    case STBOGUS:
	      Die("never happens. Trust me.");
	    }
	}
      tr[idx]->tlen = tpos;
    }	  
  *ret_tr = tr;
}


/* Function: rightjustify()
 * 
 * Purpose:  Given a gap-containing string of length n,
 *           pull all the non-gap characters as far as
 *           possible to the right, leaving gaps on the
 *           left side. Used to rearrange the positions
 *           of insertions in HMMER alignments.
 */
static void
rightjustify(char *s, int n)
{
  int npos;
  int opos;

  npos = n-1;
  opos = n-1;
  while (opos >= 0) {
    if (isgap(s[opos])) opos--;
    else                s[npos--]=s[opos--];  
  }
  while (npos >= 0) 
    s[npos--] = '.';
}