trace.c

hmmer源程序

  for (k = 0; k <= mlen; k++)
    inserts[k] = 0;
  for (idx = 0; idx < nseq; idx++) {
    nins = 0;
    for (tpos = 0; tpos < tr[idx]->tlen; tpos++) {
      switch (tr[idx]->statetype[tpos]) { 
      case STI: nins++; break;
      case STN: if (tr[idx]->statetype[tpos-1] == STN) nins++; break;
      case STC: if (tr[idx]->statetype[tpos-1] == STC) nins++; break;
      case STM:
      case STD:		/* M,D: record max. reset ctr. */
	if (nins > inserts[tr[idx]->nodeidx[tpos]-1])
	  inserts[tr[idx]->nodeidx[tpos]-1] = nins;
	nins = 0;
	break;
      case STB:		/* B; record N-tail max, reset ctr */
	if (nins > inserts[0])
	  inserts[0] = nins;
	nins = 0;
	break;
      case STT:		/* T: record C-tail max */
	if (nins > inserts[mlen])
	  inserts[mlen] = nins;
	break;
      case STS: case STE: break; /* ignore other states */
      case STJ:
	Die("yo! you don't support J in Traces2Alignment(), remember?");
      default:
	Die("Traces2Alignment reports unrecognized statetype %c", 
	    Statetype(tr[idx]->statetype[tpos]));
      }
    }
  }

				/* Insert compression option. */
  if (matchonly) 
    for (k = 0; k <= mlen; k++)
      if (inserts[k] > 1) 
	inserts[k] = 1;

  /***********************************************
   * Construct the alignment
   ***********************************************/
				/* calculate alignment length and matmap */
  matmap= (int *)   MallocOrDie (sizeof(int) * (mlen+1));
  matmap[0] = -1;
  alen = inserts[0];
  for (k = 1; k <= mlen ; k++) {
    matmap[k] = alen;
    alen += inserts[k] + 1;
  }
                                /* allocation for new alignment */
  msa = MSAAlloc(nseq, alen);

  for (idx = 0; idx < nseq; idx++) {
				/* blank an aseq */
    for (apos = 0; apos < alen; apos++)
      msa->aseq[idx][apos] = '.';
    for (k = 1; k <= mlen; k++)
      msa->aseq[idx][matmap[k]] = '-';
    msa->aseq[idx][alen] = '\0';
				/* align the sequence */
    apos = 0;
    for (tpos = 0; tpos < tr[idx]->tlen; tpos++) {
      statetype = tr[idx]->statetype[tpos]; /* just for clarity */
      rpos      = tr[idx]->pos[tpos]; 
      k         = tr[idx]->nodeidx[tpos];

      if (statetype == STM) {
	apos = matmap[k];
	msa->aseq[idx][apos] = Alphabet[(int) dsq[idx][rpos]];
	apos++;
      }
      else if (statetype == STI) {
	if (matchonly) 
	  msa->aseq[idx][apos] = '*'; /* insert compression option */
	else {
	  msa->aseq[idx][apos] = (char) tolower((int) Alphabet[(int) dsq[idx][rpos]]);
	  apos++;
	}
      }
      else if ((statetype == STN || statetype == STC) && rpos > 0) {
	if (matchonly)
	  msa->aseq[idx][apos] = '*'; /* insert compression option */
	else {
	  msa->aseq[idx][apos] = (char) tolower((int) Alphabet[(int) dsq[idx][rpos]]);
	  apos++;
	}
      }
      else if (statetype == STE)
	apos = matmap[mlen]+1;	/* set position for C-term tail */
    }

  /* N-terminal extension is right-justified.
   * Internal inserts are split in half, and C-term is right-justified.
   * C-terminal extension remains left-justified.
   */
    if (! matchonly) {
      rightjustify(msa->aseq[idx], inserts[0]);

      for (k = 1; k < mlen; k++) 
	if (inserts[k] > 1) {
	  for (nins = 0, apos = matmap[k]+1; islower((int) (msa->aseq[idx][apos])); apos++)
	    nins++;
	  nins /= 2;		/* split the insertion in half */
	  rightjustify(msa->aseq[idx]+matmap[k]+1+nins, inserts[k]-nins);
	}
    }

  }
    
  /***********************************************
   * Build the rest of the MSA annotation.
   ***********************************************/
        
  msa->nseq = nseq;
  msa->alen = alen;
  msa->au   = MallocOrDie(sizeof(char) * (strlen(RELEASE)+7));
  sprintf(msa->au, "HMMER %s", RELEASE);
				/* copy sqinfo array and weights */
  for (idx = 0; idx < nseq; idx++)
    {
      msa->sqname[idx] = sre_strdup(sqinfo[idx].name, -1);
      if (sqinfo[idx].flags & SQINFO_ACC)
	MSASetSeqAccession(msa, idx, sqinfo[idx].acc);
      if (sqinfo[idx].flags & SQINFO_DESC)
	MSASetSeqAccession(msa, idx, sqinfo[idx].desc);

      if (sqinfo[idx].flags & SQINFO_SS) {
	if (msa->ss == NULL) msa->ss = MallocOrDie(sizeof(char *) * nseq);
	MakeAlignedString(msa->aseq[idx], alen, 
			  sqinfo[idx].ss, &(msa->ss[idx]));
      }
      if (sqinfo[idx].flags & SQINFO_SA) {
	if (msa->sa == NULL) msa->sa = MallocOrDie(sizeof(char *) * nseq);
	MakeAlignedString(msa->aseq[idx], alen, 
			  sqinfo[idx].sa, &(msa->sa[idx]));
      }
      msa->wgt[idx] = wgt[idx];
    }

  /* #=RF annotation: x for match column, . for insert column
   */
  msa->rf = (char *) MallocOrDie (sizeof(char) * (alen+1));
  for (apos = 0; apos < alen; apos++)
    msa->rf[apos] = '.';
  for (k = 1; k <= mlen; k++)
    msa->rf[matmap[k]] = 'x';
  msa->rf[alen] = '\0';

    /* Currently, we produce no consensus structure. 
     * #=CS, generated from HMM structural annotation, would go here.
     */

  free(inserts);
  free(matmap);
  return msa;
}

/* Function: TransitionScoreLookup()
 * 
 * Purpose:  Convenience function used in PrintTrace() and TraceScore();
 *           given state types and node indices for a transition,
 *           return the integer score for that transition. 
 */
int
TransitionScoreLookup(struct plan7_s *hmm, char st1, int k1, 
		      char st2, int k2)
{
  switch (st1) {
  case STS: return 0;	/* S never pays */
  case STN:
    switch (st2) {
    case STB: return hmm->xsc[XTN][MOVE]; 
    case STN: return hmm->xsc[XTN][LOOP]; 
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STB:
    switch (st2) {
    case STM: return hmm->bsc[k2]; 
    case STD: return Prob2Score(hmm->tbd1, 1.);
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STM:
    switch (st2) {
    case STM: return hmm->tsc[k1][TMM];
    case STI: return hmm->tsc[k1][TMI];
    case STD: return hmm->tsc[k1][TMD];
    case STE: return hmm->esc[k1];
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STI:
    switch (st2) {
    case STM: return hmm->tsc[k1][TIM];
    case STI: return hmm->tsc[k1][TII];
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STD:
    switch (st2) {
    case STM: return hmm->tsc[k1][TDM]; 
    case STD: return hmm->tsc[k1][TDD];
    case STE: return 0;	/* D_m->E has probability 1.0 by definition in Plan7 */
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STE:
    switch (st2) {
    case STC: return hmm->xsc[XTE][MOVE]; 
    case STJ: return hmm->xsc[XTE][LOOP]; 
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STJ:
    switch (st2) {
    case STB: return hmm->xsc[XTJ][MOVE]; 
    case STJ: return hmm->xsc[XTJ][LOOP]; 
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STC:
    switch (st2) {
    case STT: return hmm->xsc[XTC][MOVE]; 
    case STC: return hmm->xsc[XTC][LOOP]; 
    default:      Die("illegal %s->%s transition", Statetype(st1), Statetype(st2));
    }
    break;
  case STT:   return 0;	/* T makes no transitions */
  default:        Die("illegal state %s in traceback", Statetype(st1));
  }
  /*NOTREACHED*/
  return 0;
}


/* Function: CreateFancyAli()
 * Date:     SRE, Mon Oct 27 06:49:44 1997 [Sanger Centre UK]
 * 
 * Purpose:  Output of an HMM/sequence alignment, using a
 *           traceback structure. Deliberately similar to 
 *           the output of BLAST, to make it easier for
 *           people to adapt their Perl parsers (or what have
 *           you) from BLAST to HMMER.
 *           
 * Args:     tr  - traceback structure that gives the alignment      
 *           hmm - the model 
 *           dsq - the sequence (digitized form)       
 *           name- name of the sequence  
 *                 
 * Return:   allocated, filled fancy alignment structure.
 */
struct fancyali_s *
CreateFancyAli(struct p7trace_s *tr, struct plan7_s *hmm,
	       char *dsq, char *name)
{
  struct fancyali_s *ali;       /* alignment to create                */
  int   tpos;			/* position in trace and alignment    */
  int   bestsym;		/* index of best symbol at this pos   */
  float mthresh;		/* above this P(x), display uppercase */

  /* Allocate and initialize the five lines of display
   */
  ali         = AllocFancyAli();
  ali->rfline = NULL;
  ali->csline = NULL;
  ali->model  = MallocOrDie (sizeof(char) * (tr->tlen+1));
  ali->mline  = MallocOrDie (sizeof(char) * (tr->tlen+1));
  ali->aseq   = MallocOrDie (sizeof(char) * (tr->tlen+1));

  memset(ali->model,  ' ', tr->tlen);
  memset(ali->mline,  ' ', tr->tlen);
  memset(ali->aseq,   ' ', tr->tlen);

  if (hmm->flags & PLAN7_RF)
    {
      ali->rfline = (char *) MallocOrDie (sizeof(char) * (tr->tlen+1));
      memset(ali->rfline, ' ', tr->tlen);
    }
  if (hmm->flags & PLAN7_CS)
    {
      ali->csline = (char *) MallocOrDie (sizeof(char) * (tr->tlen+1));
      memset(ali->csline, ' ', tr->tlen);  
    }

  ali->query  = Strdup(hmm->name);
  ali->target = Strdup(name);

  if (Alphabet_type == hmmAMINO) mthresh = 0.5;
  else                           mthresh = 0.9;

  /* Find first, last seq position
   * HMM start/end positions currently not recorded, because there
   * might be multiple HMM hits per sequence.
   */
  for (tpos = 0; tpos < tr->tlen; tpos++) 
    if (tr->pos[tpos] > 0) {
	ali->sqfrom = tr->pos[tpos];
	break;
    }
  for (tpos = tr->tlen-1; tpos >= 0; tpos--)
    if (tr->pos[tpos] > 0) {
      ali->sqto = tr->pos[tpos];
      break;
    }

  /* Fill in the five lines of display
   */
  for (tpos = 0; tpos < tr->tlen; tpos++) {
    switch (tr->statetype[tpos]) {
    case STS: 
    case STT:
      ali->model[tpos] = '*';
      break;

    case STN:
    case STJ:
    case STC:
      ali->model[tpos] = '-';
      if (tr->pos[tpos] > 0) { 
	ali->aseq[tpos] = tolower(Alphabet[(int) dsq[tr->pos[tpos]]]);
      }
      break;

    case STB: 
      ali->model[tpos] = '>';
      break;

    case STE:
      ali->model[tpos] = '<';
      break;

    case STM:
      if (hmm->flags & PLAN7_RF) ali->rfline[tpos] = hmm->rf[tr->nodeidx[tpos]];
      if (hmm->flags & PLAN7_CS) ali->csline[tpos] = hmm->cs[tr->nodeidx[tpos]];
      bestsym = FMax(hmm->mat[tr->nodeidx[tpos]], Alphabet_size);
      ali->model[tpos] = Alphabet[bestsym];
      if (hmm->mat[tr->nodeidx[tpos]][bestsym] < mthresh)
	ali->model[tpos] = tolower(ali->model[tpos]);
      if (dsq[tr->pos[tpos]] == bestsym)
	{
	  ali->mline[tpos] = Alphabet[(int) dsq[tr->pos[tpos]]];
	  if (hmm->mat[tr->nodeidx[tpos]][bestsym] < mthresh)
	    ali->mline[tpos] = tolower(ali->mline[tpos]);
	}
      else if (hmm->msc[(int) dsq[tr->pos[tpos]]] [tr->nodeidx[tpos]] > 0)
	ali->mline[tpos] = '+';
      ali->aseq[tpos]  = Alphabet[(int) dsq[tr->pos[tpos]]];
      break;

    case STD:
      if (hmm->flags & PLAN7_RF) ali->rfline[tpos] = hmm->rf[tr->nodeidx[tpos]];
      if (hmm->flags & PLAN7_CS) ali->csline[tpos] = hmm->cs[tr->nodeidx[tpos]];
      bestsym = FMax(hmm->mat[tr->nodeidx[tpos]], Alphabet_size);
      ali->model[tpos] = Alphabet[bestsym];
      if (hmm->mat[tr->nodeidx[tpos]][bestsym] < mthresh)
	ali->model[tpos] = tolower(ali->model[tpos]);
      ali->aseq[tpos]  = '-';
      break;

    case STI:
      ali->model[tpos] = '.';
      if (hmm->isc[(int) dsq[tr->pos[tpos]]] [tr->nodeidx[tpos]] > 0)
	ali->mline[tpos] = '+';
      ali->aseq[tpos]  = (char) tolower((int) Alphabet[(int) dsq[tr->pos[tpos]]]);
      break;

    default:
      Die("bogus statetype");
    } /* end switch over statetypes */
  }  /* end loop over tpos */

  ali->len          = tpos;
  if (hmm->flags & PLAN7_RF) ali->rfline[tpos] = '\0';
  if (hmm->flags & PLAN7_CS) ali->csline[tpos] = '\0';
  ali->model[tpos]  = '\0';
  ali->mline[tpos]  = '\0';
  ali->aseq[tpos]   = '\0';
  return ali;
} 


/* Function: PrintFancyAli()
 * Date:     SRE, Mon Oct 27 06:56:42 1997 [Sanger Centre UK]
 * 
 * Purpose:  Print an HMM/sequence alignment from a fancyali_s 
 *           structure. Line length controlled by ALILENGTH in
 *           config.h (set to 50).
 *           
 * Args:     fp  - where to print it (stdout or open FILE)
 *           ali - alignment to print 
 *                 
 * Return:   (void)                
 */
void
PrintFancyAli(FILE *fp, struct fancyali_s *ali)
{
  char  buffer[ALILENGTH+1];	/* output line buffer                 */
  int   starti, endi;
  int   pos;