hmmio.c

hmmer源程序

   /* specials */
   for (k = 0; k < 4; k++)
     {
       if (! fread((char *) hmm->xt[k], sizeof(float), 2, hmmfp->f))    goto FAILURE;
       if (hmmfp->byteswap) {
	 for (x = 0; x < 2; x++)
	   byteswap((char *)&(hmm->xt[k][x]), sizeof(float));
       }
     }
   
   /* null model */
   if (!fread((char *) &(hmm->p1),sizeof(float), 1, hmmfp->f))        goto FAILURE;
   if (!fread((char *)hmm->null,sizeof(float),Alphabet_size,hmmfp->f))goto FAILURE;

  /* EVD stats */
  if (hmm->flags & PLAN7_STATS) {
    if (! fread((char *) &(hmm->mu),     sizeof(float), 1, hmmfp->f))goto FAILURE;
    if (! fread((char *) &(hmm->lambda), sizeof(float), 1, hmmfp->f))goto FAILURE;

    if (hmmfp->byteswap) {
      byteswap((char *)&(hmm->mu),     sizeof(float));
      byteswap((char *)&(hmm->lambda), sizeof(float));
    }
  }

   /* entry/exit probabilities
    */
   if (! fread((char *)&(hmm->tbd1), sizeof(float), 1, hmmfp->f))       goto FAILURE;
   if (! fread((char *) hmm->begin, sizeof(float), hmm->M+1, hmmfp->f)) goto FAILURE;
   if (! fread((char *) hmm->end,   sizeof(float), hmm->M+1, hmmfp->f)) goto FAILURE;

				/* main model */
   for (k = 1; k <= hmm->M; k++)
     if (! fread((char *) hmm->mat[k], sizeof(float), Alphabet_size, hmmfp->f)) goto FAILURE;
   for (k = 1; k < hmm->M; k++)
     if (! fread((char *) hmm->ins[k], sizeof(float), Alphabet_size, hmmfp->f)) goto FAILURE;
   for (k = 1; k < hmm->M; k++)
     if (! fread((char *) hmm->t[k], sizeof(float), 7, hmmfp->f)) goto FAILURE;

  /* byteswapping
   */
  if (hmmfp->byteswap) {
    for (x = 0; x < Alphabet_size; x++) 
      byteswap((char *) &(hmm->null[x]), sizeof(float));
    byteswap((char *)&(hmm->p1),   sizeof(float));
    byteswap((char *)&(hmm->tbd1), sizeof(float));

    for (k = 1; k <= hmm->M; k++) 
      { 
	for (x = 0; x < Alphabet_size; x++) 
	  byteswap((char *)&(hmm->mat[k][x]), sizeof(float));
	if (k < hmm->M) 
	  for (x = 0; x < Alphabet_size; x++) 
	    byteswap((char *)&(hmm->ins[k][x]), sizeof(float));
	byteswap((char *)&(hmm->begin[k]),  sizeof(float));
	byteswap((char *)&(hmm->end[k]),    sizeof(float));
	if (k < hmm->M)
	  for (x = 0; x < 7; x++) 
	    byteswap((char *)&(hmm->t[k][x]), sizeof(float));
      }
  }

    
  /* set flags and return
   */
  hmm->flags |= PLAN7_HASPROB;	/* probabilities are valid  */
  hmm->flags &= ~PLAN7_HASBITS;	/* scores are not yet valid */
  *ret_hmm = hmm;
  return 1;

FAILURE:
  if (hmm != NULL) FreePlan7(hmm);
  *ret_hmm = NULL;
  return 1;
}





/* Function: read_asc19hmm()
 * Date:     Tue Apr  7 17:11:29 1998 [StL]           
 * 
 * Purpose:  Read ASCII-format tabular (1.9 and later) save files.
 * 
 *           HMMER 1.9 was only used internally at WashU, as far as
 *           I know, so this code shouldn't be terribly important
 *           to anyone.
 */
static int
read_asc19hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  struct plan7_s *hmm;
  FILE *fp;
  char  buffer[512];
  char *s;
  int   M;			/* length of model  */
  int   k;			/* state number  */
  int   x;			/* symbol number */
  int   atype;			/* Alphabet type */

  hmm = NULL;
  fp  = hmmfp->f;
  if (feof(fp) || fgets(buffer, 512, fp) == NULL) return 0;
  if (strncmp(buffer, "HMMER v1.9", 10) != 0)             goto FAILURE;

  hmm = AllocPlan7Shell();
  				/* read M from first line */
  if ((s = Getword(fp, sqdARG_INT))    == NULL) goto FAILURE;  M = atoi(s);          /* model length */
  if ((s = Getword(fp, sqdARG_INT))    == NULL) goto FAILURE;                        /* ignore alphabet size */
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;  Plan7SetName(hmm, s); /* name */
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE; /* alphabet type */ 
  s2upper(s);           
  if      (strcmp(s, "AMINO") == 0)   atype = hmmAMINO;
  else if (strcmp(s, "NUCLEIC") == 0) atype = hmmNUCLEIC;
  else goto FAILURE;

  if (Alphabet_type == hmmNOTSETYET) SetAlphabet(atype);
  else if (atype != Alphabet_type) 
    Die("Alphabet mismatch error.\nI thought we were working with %s, but tried to read a %s HMM.\n", AlphabetType2String(Alphabet_type), AlphabetType2String(atype));

				/* read alphabet, make sure it's Plan7-compatible... */
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
  if (strncmp(s, Alphabet, Alphabet_size) != 0) goto FAILURE;

				/* whether we have ref, cs info */
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
  if (strcmp(s, "yes") == 0) hmm->flags |= PLAN7_RF;
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
  if (strcmp(s, "yes") == 0) hmm->flags |= PLAN7_CS;

				/* null model. 1.9 has emissions only. invent transitions. */
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
  if (strcmp(s, "null") != 0) goto FAILURE;
  for (x = 0; x < Alphabet_size; x++) {
    if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
    hmm->null[x] = ascii2prob(s, 1.0);
  }
  hmm->p1 = (Alphabet_type == hmmAMINO)? 350./351. : 1000./1001.; 

  /* Done with header; check some stuff before proceeding
   */
  if (feof(hmmfp->f))                goto FAILURE;
  if (M < 1)                         goto FAILURE;
  if (hmm->name == NULL)             goto FAILURE;
  if (Alphabet_type == hmmNOTSETYET) goto FAILURE;

  /* Allocate the model. Set up the probabilities that Plan9
   * doesn't set.
   */
  AllocPlan7Body(hmm, M);
  ZeroPlan7(hmm);
  Plan7LSConfig(hmm);

  /* The zero row has: 4 or 20 unused scores for nonexistent M0 state
   * then: B->M, tbd1, a B->I that Plan7 doesn't have;
   *       three unused D-> transitions; then three I0 transitions that Plan7 doesn't have;
   *       then two unused rf, cs annotations.
   */
  if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE; /* position index ignored */
  for (x = 0; x < Alphabet_size; x++)
    if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE; /* emissions ignored */
  if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
  hmm->begin[1] = ascii2prob(s, 1.0);
  if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
  hmm->tbd1 = ascii2prob(s, 1.0);
				/* renormalize */
  hmm->begin[1] = hmm->begin[1] / (hmm->begin[1] + hmm->tbd1);
  hmm->tbd1     = hmm->tbd1     / (hmm->begin[1] + hmm->tbd1);
				/* skip rest of line, seven integer fields, two char fields */
  for (x = 0; x < 7; x++)
    if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
  if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;

				/* main model: table of emissions, transitions, annotation */
  for (k = 1; k <= hmm->M; k++)
    {
				/* position index ignored */
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
				/* match emissions */
      for (x = 0; x < Alphabet_size; x++) {
	if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
	hmm->mat[k][x] = ascii2prob(s, hmm->null[x]);
      }
				/* nine transitions; two are ignored */
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TMM] = ascii2prob(s, 1.0);
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TMD] = (k == hmm->M) ? 0.0 : ascii2prob(s, 1.0);
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TMI] = ascii2prob(s, 1.0);
      
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TDM] = ascii2prob(s, 1.0);
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TDD] = (k == hmm->M) ? 0.0 : ascii2prob(s, 1.0);
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;/* TDI ignored. */

				/* no insert state at k == M, be careful */
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TIM]  = ascii2prob(s, 1.0);
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE; /* TID ignored. */
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
      if (k < hmm->M) hmm->t[k][TII] = ascii2prob(s, 1.0);
	
				/* annotations */
      if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
      if (hmm->flags & PLAN7_RF) hmm->rf[k] = *s;
      if ((s = Getword(fp, sqdARG_STRING)) == NULL) goto FAILURE;
      if (hmm->flags & PLAN7_CS) hmm->cs[k] = *s;
    }
				/* table of insert emissions; 
                                 * Plan7 has no insert state at 0 or M  */
  for (k = 0; k <= hmm->M; k++)
    {
      if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE; /* position index ignored */
      for (x = 0; x < Alphabet_size; x++) {
	if ((s = Getword(fp, sqdARG_INT)) == NULL) goto FAILURE;
	if (k > 0 && k < hmm->M)
	  hmm->ins[k][x] = ascii2prob(s, hmm->null[x]);
      }
    }

  /* Set flags and return
   */
  hmm->flags |= PLAN7_HASPROB;	/* probabilities are valid */
  hmm->flags &= ~PLAN7_HASBITS;	/* scores are not valid    */
  Plan7Renormalize(hmm);
  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);
  *ret_hmm = hmm;
  return 1;

FAILURE:
  if (hmm  != NULL) FreePlan7(hmm);
  *ret_hmm = NULL;
  return 1;
}

static int  
read_bin19hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  unsigned int     magic;
  struct plan7_s  *hmm;     /* plan7 HMM */ 
  struct plan9_s  *p9hmm;   /* old style 1.x HMM */
  
  /* Read the magic number; if we don't see it, then we
   * must be out of data in the file.
   */
  if (feof(hmmfp->f)) return 0;
  if (! fread((char *) &magic, sizeof(unsigned int), 1, hmmfp->f)) return 0;

  p9hmm = read_plan9_binhmm(hmmfp->f, HMMER1_9B, hmmfp->byteswap);
  if (p9hmm == NULL) { *ret_hmm = NULL; return 1; }

  Plan9toPlan7(p9hmm, &hmm);

  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);

  P9FreeHMM(p9hmm);
 *ret_hmm = hmm;
  return 1;
}
static int  
read_asc17hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  struct plan7_s  *hmm;     /* plan7 HMM */ 
  struct plan9_s  *p9hmm;   /* old style 1.x HMM */
  char   buffer[512];

  /* Read the magic header; if we don't see it, then
   * we must be out of data in the file.
   */
  if (feof(hmmfp->f) || fgets(buffer, 512, hmmfp->f) == NULL) return 0;

  p9hmm = read_plan9_aschmm(hmmfp->f, HMMER1_7F);
  if (p9hmm == NULL) { *ret_hmm = NULL; return 1; }

  Plan9toPlan7(p9hmm, &hmm);

  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);

  P9FreeHMM(p9hmm);
  Plan7Renormalize(hmm);
 *ret_hmm = hmm;
  return 1;
}

static int  
read_bin17hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  unsigned int     magic;
  struct plan7_s  *hmm;     /* plan7 HMM */ 
  struct plan9_s  *p9hmm;   /* old style 1.x HMM */
  
  /* Read the magic number; if we don't see it, then we
   * must be out of data in the file.
   */
  if (feof(hmmfp->f)) return 0;
  if (! fread((char *) &magic, sizeof(unsigned int), 1, hmmfp->f)) return 0;

  p9hmm = read_plan9_binhmm(hmmfp->f, HMMER1_7B, hmmfp->byteswap);
  if (p9hmm == NULL) { *ret_hmm = NULL; return 1; }

  Plan9toPlan7(p9hmm, &hmm);

  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);

  P9FreeHMM(p9hmm);
 *ret_hmm = hmm;
  return 1;
}

static int  
read_asc11hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  Die("1.1 ASCII HMMs unsupported");
  return 1;
}
static int  
read_bin11hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  unsigned int     magic;
  struct plan7_s  *hmm;     /* plan7 HMM */ 
  struct plan9_s  *p9hmm;   /* old style 1.x HMM */
  
  /* Read the magic number; if we don't see it, then we
   * must be out of data in the file.
   */
  if (feof(hmmfp->f)) return 0;
  if (! fread((char *) &magic, sizeof(unsigned int), 1, hmmfp->f)) return 0;

  p9hmm = read_plan9_binhmm(hmmfp->f, HMMER1_1B, hmmfp->byteswap);
  if (p9hmm == NULL) { *ret_hmm = NULL; return 1; }

  Plan9toPlan7(p9hmm, &hmm);

  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);

  P9FreeHMM(p9hmm);
 *ret_hmm = hmm;
  return 1;
}

static int  
read_asc10hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  Die("1.0 ASCII HMMs unsupported");
  return 1;
}

static int  
read_bin10hmm(HMMFILE *hmmfp, struct plan7_s **ret_hmm)
{
  unsigned int     magic;
  struct plan7_s  *hmm;     /* plan7 HMM */ 
  struct plan9_s  *p9hmm;   /* old style 1.x HMM */
  
  /* Read the magic number; if we don't see it, then we
   * must be out of data in the file.
   */
  if (feof(hmmfp->f)) return 0;
  if (! fread((char *) &magic, sizeof(unsigned int), 1, hmmfp->f)) return 0;

  p9hmm = read_plan9_binhmm(hmmfp->f, HMMER1_0B, hmmfp->byteswap);
  if (p9hmm == NULL) { *ret_hmm = NULL; return 1; }

  Plan9toPlan7(p9hmm, &hmm);

  hmm->comlog = Strdup("[converted from an old Plan9 HMM]");
  Plan7SetCtime(hmm);

  P9FreeHMM(p9hmm);
 *ret_hmm = hmm;
  return 1; 
}

/*****************************************************************
 * Some miscellaneous utility functions
 *****************************************************************/

/* Function: prob2ascii()
 * 
 * Purpose:  Format a probability for output to an ASCII save
 *           file. Returns a ptr to a static internal buffer.
 *              
 */
static char *
prob2ascii(float p, float null)
{
  static char buffer[8];

  if (p == 0.0) return "*";
  sprintf(buffer, "%6d", Prob2Score(p, null));
  return buffer;
}


/* Function: ascii2prob()
 * 
 * Purpose:  Convert a saved string back to a probability.
 */
static float
ascii2prob(char *s, float null)
{
  return (*s == '*') ? 0. : Score2Prob(atoi(s), null);
}

/* Function: byteswap()
 * 
 * Purpose:  Swap between big-endian and little-endian.
 *           For example:
 *               int foo = 0x12345678;
 *               byteswap((char *) &foo, sizeof(int));
 *               printf("%x\n", foo)
 *           gives 78563412.
 *           
 *           I don't fully understand byte-swapping issues.
 *           However, I have tested this on chars through floats,
 *           on various machines:
 *               SGI IRIX 4.0.5, SunOS 4.1.3, DEC Alpha OSF/1, Alliant
 *               
 *           Note: this is only a partial solution to the problem of
 *           binary file portability. 32 bit integers are assumed by HMMER,
 *           for instance. This should be true for all UNIX, VAX, and WinNT
 *           platforms, I believe.     
 *
 * Date: Sun Feb 12 10:26:22 1995              
 */
static void
byteswap(char *swap, int nbytes)
{