modelmakers.c

hmmer源程序

	      tr[idx]->nodeidx[tpos]   = k;
	      tr[idx]->pos[tpos]       = i;
	      i++;
	      tpos++;
	    }	      
          else if (matassign[apos+1] & ASSIGN_MATCH)
            {                   /* DELETE */
		/* being careful about S/W transitions; no B->D transitions */
	      k++;		/* *always* move on model when ASSIGN_MATCH */
	      if (tr[idx]->statetype[tpos-1] != STB)
		{
		  tr[idx]->statetype[tpos] = STD;
		  tr[idx]->nodeidx[tpos]   = k;
		  tr[idx]->pos[tpos]       = 0;
		  tpos++;
		}
            }
          else if (matassign[apos+1] & EXTERNAL_INSERT_N &&
		   ! isgap(aseq[idx][apos]))
            {                   /* N-TERMINAL TAIL */
              tr[idx]->statetype[tpos] = STN;
              tr[idx]->nodeidx[tpos]   = 0;
              tr[idx]->pos[tpos]       = i;
	      i++;
	      tpos++;
            }
	  else if (matassign[apos+1] & EXTERNAL_INSERT_C &&
		   ! isgap(aseq[idx][apos]))
	    {			/* C-TERMINAL TAIL */
	      tr[idx]->statetype[tpos] = STC;
              tr[idx]->nodeidx[tpos]   = 0;
              tr[idx]->pos[tpos]       = i;
	      i++;
	      tpos++;
	    }
	  else if (! isgap(aseq[idx][apos]))
	    {			/* INSERT */
	      tr[idx]->statetype[tpos] = STI;
              tr[idx]->nodeidx[tpos]   = k;
              tr[idx]->pos[tpos]       = i;
	      i++;
	      tpos++;
	    }

	  if (matassign[apos+1] & LAST_MATCH)
	    {			/* END */
	      /* be careful about S/W transitions; may need to roll
	       * back over some D's because there's no D->E transition
	       */
	      while (tr[idx]->statetype[tpos-1] == STD) 
		tpos--;
	      tr[idx]->statetype[tpos] = STE;
	      tr[idx]->nodeidx[tpos]   = 0;
	      tr[idx]->pos[tpos]       = 0;
	      tpos++;
				/* and then transit E->C;
				   alignments that use J are undefined;
				   C-term tail is emitted on C->C transitions */
	      tr[idx]->statetype[tpos] = STC;
	      tr[idx]->nodeidx[tpos]   = 0;
	      tr[idx]->pos[tpos]       = 0;
	      tpos++;
	    }
        }
                                /* all traces end with T state */
      tr[idx]->statetype[tpos] = STT;
      tr[idx]->nodeidx[tpos]   = 0;
      tr[idx]->pos[tpos]       = 0;
      tr[idx]->tlen = ++tpos;
				/* deal with DI, ID transitions */
				/* k == M here */
      trace_doctor(tr[idx], k, NULL, NULL);

    }    /* end for sequence # idx */

  *ret_tr = tr;
  return;
}

/* Function: trace_doctor()
 * 
 * Purpose:  Plan 7 disallows D->I and I->D "chatter" transitions.
 *           However, these transitions may be implied by many
 *           alignments for hand- or heuristic- built HMMs.
 *           trace_doctor() collapses I->D or D->I into a
 *           single M position in the trace. 
 *           Similarly, B->I and I->E transitions may be implied
 *           by an alignment.
 *           
 *           trace_doctor does not examine any scores when it does
 *           this. In ambiguous situations (D->I->D) the symbol
 *           will be pulled arbitrarily to the left, regardless
 *           of whether that's the best column to put it in or not.
 *           
 * Args:     tr      - trace to doctor
 *           M       - length of model that traces are for 
 *           ret_ndi - number of DI transitions doctored
 *           ret_nid - number of ID transitions doctored
 * 
 * Return:   (void)
 *           tr is modified
 */               
static void
trace_doctor(struct p7trace_s *tr, int mlen, int *ret_ndi, int *ret_nid)
{
  int opos;			/* position in old trace                 */
  int npos;			/* position in new trace (<= opos)       */
  int ndi, nid;			/* number of DI, ID transitions doctored */

				/* overwrite the trace from left to right */
  ndi  = nid  = 0;
  opos = npos = 0;
  while (opos < tr->tlen) {
      /* fix implied D->I transitions; D transforms to M, I pulled in */
    if (tr->statetype[opos] == STD && tr->statetype[opos+1] == STI) {
      tr->statetype[npos] = STM;
      tr->nodeidx[npos]   = tr->nodeidx[opos]; /* D transforms to M      */
      tr->pos[npos]       = tr->pos[opos+1];   /* insert char moves back */
      opos += 2;
      npos += 1;
      ndi++;
    } /* fix implied I->D transitions; D transforms to M, I is pushed in */
    else if (tr->statetype[opos]== STI && tr->statetype[opos+1]== STD) {
      tr->statetype[npos] = STM;
      tr->nodeidx[npos]   = tr->nodeidx[opos+1];/* D transforms to M    */
      tr->pos[npos]       = tr->pos[opos];      /* insert char moves up */
      opos += 2;
      npos += 1;
      nid++; 
    } /* fix implied B->I transitions; pull I back to its M */
    else if (tr->statetype[opos]== STI && tr->statetype[opos-1]== STB) {
      tr->statetype[npos] = STM;
      tr->nodeidx[npos]   = tr->nodeidx[opos]; /* offending I transforms to M */
      tr->pos[npos]       = tr->pos[opos];
      opos++;
      npos++;
    } /* fix implied I->E transitions; push I to next M */
    else if (tr->statetype[opos]== STI && tr->statetype[opos+1]== STE) {
      tr->statetype[npos] = STM;
      tr->nodeidx[npos]   = tr->nodeidx[opos]+1;/* offending I transforms to M */
      tr->pos[npos]       = tr->pos[opos];
      opos++;
      npos++;
    } /* rare: N-N-B-E becomes N-B-M_1-E (swap B,N) */
    else if (tr->statetype[opos]==STB && tr->statetype[opos+1]==STE
	     && tr->statetype[opos-1]==STN && tr->pos[opos-1] > 0) {   
      tr->statetype[npos]   = STM;
      tr->nodeidx[npos]     = 1;
      tr->pos[npos]         = tr->pos[opos-1];
      tr->statetype[npos-1] = STB;
      tr->nodeidx[npos-1]   = 0;
      tr->pos[npos-1]       = 0;
      opos++;
      npos++;
    } /* rare: B-E-C-C-x becomes B-M_M-E-C-x (swap E,C) */
    else if (tr->statetype[opos]==STE && tr->statetype[opos-1]==STB
	     && tr->statetype[opos+1]==STC 
	     && tr->statetype[opos+2]==STC) { 
      tr->statetype[npos]   = STM;
      tr->nodeidx[npos]     = mlen;
      tr->pos[npos]         = tr->pos[opos+2];
      tr->statetype[npos+1] = STE;
      tr->nodeidx[npos+1]   = 0;
      tr->pos[npos+1]       = 0;
      tr->statetype[npos+2] = STC; /* first C must be a nonemitter  */
      tr->nodeidx[npos+2]   = 0;
      tr->pos[npos+2]       = 0;
      opos+=3;
      npos+=3;
    } /* everything else is just copied */
    else {
      tr->statetype[npos] = tr->statetype[opos];
      tr->nodeidx[npos]   = tr->nodeidx[opos];
      tr->pos[npos]       = tr->pos[opos];
      opos++;
      npos++;
    }
  }
  tr->tlen = npos;

  if (ret_ndi != NULL) *ret_ndi = ndi;
  if (ret_nid != NULL) *ret_nid = nid;
  return;
}


/* Function: annotate_model()
 * 
 * Purpose:  Add rf, cs optional annotation to a new model.
 * 
 * Args:     hmm       - new model
 *           matassign - which alignment columns are MAT; [1..alen]
 *           msa       - alignment, including annotation to transfer
 *           
 * Return:   (void)
 */
static void
annotate_model(struct plan7_s *hmm, int *matassign, MSA *msa)
{                      
  int   apos;			/* position in matassign, 1.alen  */
  int   k;			/* position in model, 1.M         */
  char *pri;			/* X-PRM, X-PRI, X-PRT annotation */

  /* Transfer reference coord annotation from alignment,
   * if available
   */
  if (msa->rf != NULL) {
    hmm->rf[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos] & ASSIGN_MATCH) /* ainfo is off by one from HMM */
	hmm->rf[k++] = (msa->rf[apos-1] == ' ') ? '.' : msa->rf[apos-1];
    hmm->rf[k] = '\0';
    hmm->flags |= PLAN7_RF;
  }

  /* Transfer consensus structure annotation from alignment, 
   * if available
   */
  if (msa->ss_cons != NULL) {
    hmm->cs[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos] & ASSIGN_MATCH)
	hmm->cs[k++] = (msa->ss_cons[apos-1] == ' ') ? '.' : msa->ss_cons[apos-1];
    hmm->cs[k] = '\0';
    hmm->flags |= PLAN7_CS;
  }

  /* Transfer surface accessibility annotation from alignment,
   * if available
   */
  if (msa->sa_cons != NULL) {
    hmm->ca[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos] & ASSIGN_MATCH)
	hmm->ca[k++] = (msa->sa_cons[apos-1] == ' ') ? '.' : msa->sa_cons[apos-1];
    hmm->ca[k] = '\0';
    hmm->flags |= PLAN7_CA;
  }

  /* Store the alignment map
   */
  for (apos = k = 1; apos <= msa->alen; apos++)
    if (matassign[apos] & ASSIGN_MATCH)
      hmm->map[k++] = apos;
  hmm->flags |= PLAN7_MAP;

  /* Translate and transfer X-PRM annotation. 
   * 0-9,[a-zA-Z] are legal; translate as 0-9,10-35 into hmm->mpri.
   * Any other char is translated as -1, and this will be interpreted
   * as a flag that means "unknown", e.g. use the normal mixture Dirichlet
   * procedure for this column.
   */
  if ((pri = MSAGetGC(msa, "X-PRM")) != NULL)
    {
      hmm->mpri = MallocOrDie(sizeof(int) * (hmm->M+1));
      for (apos = k = 1; apos <= msa->alen; apos++)
	if (matassign[apos] & ASSIGN_MATCH)
	  {
	    if      (isdigit((int) pri[apos-1])) hmm->mpri[k] = pri[apos-1] - '0';
	    else if (islower((int) pri[apos-1])) hmm->mpri[k] = pri[apos-1] - 'a' + 10;
	    else if (isupper((int) pri[apos-1])) hmm->mpri[k] = pri[apos-1] - 'A' + 10;
	    else hmm->mpri[k] = -1;
	    k++;
	  }
    }
  /* And again for X-PRI annotation on insert priors:
   */
  if ((pri = MSAGetGC(msa, "X-PRI")) != NULL)
    {
      hmm->ipri = MallocOrDie(sizeof(int) * (hmm->M+1));
      for (apos = k = 1; apos <= msa->alen; apos++)
	if (matassign[apos] & ASSIGN_MATCH)
	  {
	    if      (isdigit((int) pri[apos-1])) hmm->ipri[k] = pri[apos-1] - '0';
	    else if (islower((int) pri[apos-1])) hmm->ipri[k] = pri[apos-1] - 'a' + 10;
	    else if (isupper((int) pri[apos-1])) hmm->ipri[k] = pri[apos-1] - 'A' + 10;
	    else hmm->ipri[k] = -1;
	    k++;
	  }
    }
  /* And one last time for X-PRT annotation on transition priors:
   */
  if ((pri = MSAGetGC(msa, "X-PRT")) != NULL)
    {
      hmm->tpri = MallocOrDie(sizeof(int) * (hmm->M+1));
      for (apos = k = 1; apos <= msa->alen; apos++)
	if (matassign[apos] & ASSIGN_MATCH)
	  {
	    if      (isdigit((int) pri[apos-1])) hmm->tpri[k] = pri[apos-1] - '0';
	    else if (islower((int) pri[apos-1])) hmm->tpri[k] = pri[apos-1] - 'a' + 10;
	    else if (isupper((int) pri[apos-1])) hmm->tpri[k] = pri[apos-1] - 'A' + 10;
	    else hmm->tpri[k] = -1;
	    k++;
	  }
    }

}

static void
print_matassign(int *matassign, int alen)
{
  int apos;

  for (apos = 0; apos <= alen; apos++) {
    printf("%3d %c %c %c\n", 
	   apos,
	   (matassign[apos] & ASSIGN_MATCH) ? 'x':' ',
	   (matassign[apos] & FIRST_MATCH || matassign[apos] & LAST_MATCH) ? '<' : ' ',
	   (matassign[apos] & EXTERNAL_INSERT_N ||
	    matassign[apos] & EXTERNAL_INSERT_C) ? '|':' ');
  }
}