Logo Search packages:      
Sourcecode: hmmer version File versions  Download package

debug.c

/*****************************************************************
 * HMMER - Biological sequence analysis with profile HMMs
 * Copyright (C) 1992-1998 Washington University School of Medicine
 * All Rights Reserved
 * 
 *     This source code is distributed under the terms of the
 *     GNU General Public License. See the files COPYING and LICENSE
 *     for details.
 *****************************************************************/

/* debug.c
 * Thu Nov 21 09:58:05 1996
 * 
 * Printing out or naming various useful things from HMMER
 * innards.
 * 
 * RCS $Id: debug.c,v 1.9 1998/10/09 17:40:25 eddy Exp $
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <ctype.h>

#include "structs.h"
#include "config.h"
#include "funcs.h" 
#include "squid.h"

#ifdef MEMDEBUG
#include "dbmalloc.h"
#endif



/* Function: Statetype()
 * 
 * Purpose:  Returns the state type in text.
 * Example:  Statetype(S) = "S"
 */
char *
Statetype(char st)
{
  switch (st) {
  case STS: return "S";
  case STN: return "N";
  case STB: return "B";
  case STM: return "M";
  case STD: return "D";
  case STI: return "I";
  case STE: return "E";
  case STJ: return "J";
  case STC: return "C";
  case STT: return "T";
  default: return "BOGUS";
  }
}

/* Function: P7PrintTrace()
 * 
 * Purpose:  Print out a traceback structure.
 *           If hmm is non-NULL, also print transition and emission scores.
 *           
 * Args:     fp  - stderr or stdout, often
 *           tr  - trace structure to print
 *           hmm - NULL or hmm containing scores to print
 *           dsq - NULL or digitized sequence trace refers to.                
 */
void
P7PrintTrace(FILE *fp, struct p7trace_s *tr, struct plan7_s *hmm, char *dsq)
{
  int tpos;             /* counter for trace position */
  int sym;
  int sc; 

  if (hmm == NULL) {
    fprintf(fp, "st  node   rpos  - traceback len %d\n", tr->tlen);
    fprintf(fp, "--  ---- ------\n");
    for (tpos = 0; tpos < tr->tlen; tpos++) {
      fprintf(fp, "%1s  %4d %6d\n", 
            Statetype(tr->statetype[tpos]),
            tr->nodeidx[tpos],
            tr->pos[tpos]);
    } 
  } else {
    if (!(hmm->flags & PLAN7_HASBITS))
      Die("oi, you can't print scores from that hmm, it's not ready.");

    sc = 0;
    fprintf(fp, "st  node   rpos  transit emission - traceback len %d\n", tr->tlen);
    fprintf(fp, "--  ---- ------  ------- --------\n");
    for (tpos = 0; tpos < tr->tlen; tpos++) {
      if (dsq != NULL) sym = (int) dsq[tr->pos[tpos]];

      fprintf(fp, "%1s  %4d %6d  %7d", 
            Statetype(tr->statetype[tpos]),
            tr->nodeidx[tpos],
            tr->pos[tpos],
            (tpos < tr->tlen-1) ? 
            TransitionScoreLookup(hmm, tr->statetype[tpos], tr->nodeidx[tpos],
                            tr->statetype[tpos+1], tr->nodeidx[tpos+1]) : 0);

      if (tpos < tr->tlen-1)
      sc += TransitionScoreLookup(hmm, tr->statetype[tpos], tr->nodeidx[tpos],
                            tr->statetype[tpos+1], tr->nodeidx[tpos+1]);

      if (dsq != NULL) {
      if (tr->statetype[tpos] == STM)  
        {
          fprintf(fp, " %8d %c", hmm->msc[sym][tr->nodeidx[tpos]], 
                Alphabet[sym]);
          sc += hmm->msc[sym][tr->nodeidx[tpos]];
        }
      else if (tr->statetype[tpos] == STI) 
        {
          fprintf(fp, " %8d %c", hmm->isc[sym][tr->nodeidx[tpos]], 
                tolower(Alphabet[sym]));
          sc += hmm->isc[sym][tr->nodeidx[tpos]];
        }
      else if ((tr->statetype[tpos] == STN && tr->statetype[tpos-1] == STN) ||
             (tr->statetype[tpos] == STC && tr->statetype[tpos-1] == STC) ||
             (tr->statetype[tpos] == STJ && tr->statetype[tpos-1] == STJ))
        {
          fprintf(fp, " %8d %c", 0, tolower(Alphabet[sym]));
        }
      } else {
      fprintf(fp, " %8s %c", "-", '-');
      }


      fputs("\n", fp);
    }
    fprintf(fp, "                 ------- --------\n");
    fprintf(fp, "           total: %6d\n\n", sc);
  }
}

/* Function: P7PrintPrior()
 * 
 * Purpose:  Print out a Plan 7 prior structure.
 */
void
P7PrintPrior(FILE *fp, struct p7prior_s *pri)
{
  int q, x;             /* counters for mixture component, element */
  
  if      (pri->strategy == PRI_DCHLET) fputs("Dirichlet\n", fp);
  else if (pri->strategy == PRI_PAM)    fputs("PAM\n", fp);
  else Die("No such strategy.");
  
  if      (Alphabet_type == hmmAMINO)   fputs("Amino\n", fp);
  else if (Alphabet_type == hmmNUCLEIC) fputs("Nucleic\n", fp);

  /* Transitions
   */
  fprintf(fp, "\n%d\n", pri->tnum);
  for (q = 0; q < pri->tnum; q++)
    {
      fprintf(fp, "%.4f\n", pri->tq[q]);
      for (x = 0; x < 7; x++)
      fprintf(fp, "%.4f ", pri->t[q][x]);
      fputs("\n", fp);
    }

  /* Match emissions
   */
  fprintf(fp, "\n%d\n", pri->mnum);
  for (q = 0; q < pri->mnum; q++)
    {
      fprintf(fp, "%.4f\n", pri->mq[q]);
      for (x = 0; x < Alphabet_size; x++)
      fprintf(fp, "%.4f ", pri->m[q][x]);
      fputs("\n", fp);
    }

  /* Insert emissions
   */
  fprintf(fp, "\n%d\n", pri->inum);
  for (q = 0; q < pri->inum; q++)
    {
      fprintf(fp, "%.4f\n", pri->iq[q]);
      for (x = 0; x < Alphabet_size; x++)
      fprintf(fp, "%.4f ", pri->i[q][x]);
      fputs("\n", fp);
    }
}

/* Function: TraceVerify()
 * Date:     SRE, Mon Feb  2 07:48:52 1998 [St. Louis]
 *
 * Purpose:  Check a traceback structure for internal consistency.
 *           Used in Shiva testsuite, for example.
 *
 * Args:     tr  - traceback to verify
 *           M   - length of HMM
 *           N   - length of sequence      
 *
 * Returns:  1 if OK. 0 if not.
 */
int
TraceVerify(struct p7trace_s *tr, int M, int N)
{
  int tpos;             /* position in trace                  */
  int k;                /* current position in HMM nodes 1..M */
  int i;                /* current position in seq 1..N       */
  int nn, nc, nj;       /* number of STN's, STC's, STJ's seen */
  int nm;               /* number of STM's seen */
  
  /* Basic checks on ends.
   */
  if (tr->statetype[0] != STS)          return 0;
  if (tr->statetype[1] != STN)          return 0;
  if (tr->statetype[tr->tlen-2] != STC) return 0;
  if (tr->statetype[tr->tlen-1] != STT) return 0;
  if (tr->pos[1] != 0)                  return 0;

  /* Check for consistency throughout trace
   */
  k = i = nn = nc = nj = nm = 0;
  for (tpos = 0; tpos < tr->tlen; tpos++)
    {
      switch (tr->statetype[tpos]) {
      case STS:
      if (tr->nodeidx[tpos] != 0) return 0;
      if (tr->pos[tpos]     != 0) return 0;
      if (k != 0)                 return 0;
      if (i != 0)                 return 0;
      if (tpos != 0)              return 0;
      break;

      case STN:               /* first N doesn't emit. */
      if (tr->nodeidx[tpos] != 0) return 0;
      if (k != 0)                 return 0;
      if (nn > 0)
        {
          if (tr->pos[tpos] != i+1) return 0;
          i++;
        }
      else 
        {
          if (tr->pos[tpos] != 0) return 0;
          if (i != 0)             return 0;
        }
      nn++;
      break;

      case STB:
      if (tr->nodeidx[tpos] != 0) return 0;
      if (tr->pos[tpos]     != 0) return 0;
      nm = 0;
      break;

      case STM:               /* can enter anywhere on first M */
      if (tr->pos[tpos] != i+1) return 0;
      if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M) return 0;
      i++;
      if (nm == 0)  k = tr->nodeidx[tpos];
      else {
        if (tr->nodeidx[tpos] != k+1) return 0;
        k++;
      }
      nm++;
      break;

      case STI:
      if (tr->pos[tpos] != i+1)   return 0;
      if (tr->nodeidx[tpos] != k) return 0;
      if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M-1) return 0;
      if (k >= M)                 return 0;
      i++;
      break;

      case STD:
      if (tr->pos[tpos] != 0)       return 0;
      if (tr->nodeidx[tpos] != k+1) return 0;
      if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M) return 0;
      k++;
      break;

      case STE:
      if (tr->nodeidx[tpos] != 0) return 0;
      if (tr->pos[tpos]     != 0) return 0;
      nj = 0;
      break;

      case STJ:
      if (tr->nodeidx[tpos] != 0) return 0;
      if (nj > 0)
        {
          if (tr->pos[tpos] != i+1) return 0;
          i++;
        }
      else if (tr->pos[tpos] != 0) return 0;
      nj++;
      break;

      case STC:
      if (tr->nodeidx[tpos] != 0) return 0;
      if (nc > 0)
        {
          if (tr->pos[tpos] != i+1) return 0;
          i++;
        }
      else if (tr->pos[tpos] != 0)  return 0;
      nc++;
      break;

      case STT:
      if (tpos != tr->tlen - 1)   return 0;
      if (tr->nodeidx[tpos] != 0) return 0;
      if (tr->pos[tpos]     != 0) return 0;
      if (i != N)                 return 0;
      break;

      case STBOGUS:
      default:
      return 0;
      }     /* end switch over statetypes */
    } /* end loop over trace positions */

  return 1;
}


/* Function: TraceCompare()
 * Date:     SRE, Wed Mar  4 17:26:49 1998 [St. Louis]
 *
 * Purpose:  Compare two tracebacks; return 1 if they're
 *           identical, else 0. Written for Shiva testsuite.
 *
 * Args:     t1 - first trace
 *           t2 - second trace     
 *
 * Returns:  1 if identical; 0 elsewise
 */
int
TraceCompare(struct p7trace_s *t1, struct p7trace_s *t2)
{
  int tpos;

  if (t1->tlen != t2->tlen) return 0;

  for (tpos = 0; tpos < t1->tlen; tpos++)
    {
      if (t1->statetype[tpos] != t2->statetype[tpos]) return 0;
      if (t1->nodeidx[tpos]   != t2->nodeidx[tpos])   return 0;
      if (t1->pos[tpos]       != t2->pos[tpos])       return 0;
    }
  return 1;
}


Generated by  Doxygen 1.6.0   Back to index