stockholm.c

这是一个基于HMM 模型的生物多序列比对算法的linux实现版本。hmmer

       * These are stored, for example, as:
       *     msa->gs[0][0] = "PDB; 1xxx;\nPDB; 2yyy;"
       * and must be decomposed.
       */
      for (j = 0; j < msa->nseq; j++)
	if (msa->gs[i][j] != NULL)
	  {
	    s = msa->gs[i][j];
	    while ((tok = sre_strtok(&s, "\n", NULL)) != NULL)
	      fprintf(fp, "#=GS %*.*s %5s %s\n", namewidth, namewidth,
		      msa->sqname[j], msa->gs_tag[i], tok);
	  }
      fprintf(fp, "\n");
    }

  /* Alignment section:
   * contains aligned sequence, #=GR annotation, and #=GC annotation
   */
  for (currpos = 0; currpos < msa->alen; currpos += cpl)
    {
      if (currpos > 0) fprintf(fp, "\n");
      for (i = 0; i < msa->nseq; i++)
	{
	  strncpy(buf, msa->aseq[i] + currpos, cpl);
	  buf[cpl] = '\0';	      
	  fprintf(fp, "%-*.*s  %s\n", namewidth+typewidth+markupwidth, namewidth+typewidth+markupwidth, 
		  msa->sqname[i], buf);

	  if (msa->ss != NULL && msa->ss[i] != NULL) {
	    strncpy(buf, msa->ss[i] + currpos, cpl);
	    buf[cpl] = '\0';	 
	    fprintf(fp, "#=GR %-*.*s SS     %s\n", namewidth, namewidth, msa->sqname[i], buf);
	  }
	  if (msa->sa != NULL && msa->sa[i] != NULL) {
	    strncpy(buf, msa->sa[i] + currpos, cpl);
	    buf[cpl] = '\0';
	    fprintf(fp, "#=GR %-*.*s SA     %s\n", namewidth, namewidth, msa->sqname[i], buf);
	  }
	  for (j = 0; j < msa->ngr; j++)
	    if (msa->gr[j][i] != NULL) {
	      strncpy(buf, msa->gr[j][i] + currpos, cpl);
	      buf[cpl] = '\0';
	      fprintf(fp, "#=GR %-*.*s %5s  %s\n", 
		      namewidth, namewidth, msa->sqname[i], msa->gr_tag[j], buf);
	    }
	}
      if (msa->ss_cons != NULL) {
	strncpy(buf, msa->ss_cons + currpos, cpl);
	buf[cpl] = '\0';
	fprintf(fp, "#=GC %-*.*s %s\n", namewidth+typewidth, namewidth+typewidth, "SS_cons", buf);
      }

      if (msa->sa_cons != NULL) {
	strncpy(buf, msa->sa_cons + currpos, cpl);
	buf[cpl] = '\0';
	fprintf(fp, "#=GC %-*.*s %s\n", namewidth+typewidth, namewidth+typewidth, "SA_cons", buf);
      }

      if (msa->rf != NULL) {
	strncpy(buf, msa->rf + currpos, cpl);
	buf[cpl] = '\0';
	fprintf(fp, "#=GC %-*.*s %s\n", namewidth+typewidth, namewidth+typewidth, "RF", buf);
      }
      for (j = 0; j < msa->ngc; j++) {
	strncpy(buf, msa->gc[j] + currpos, cpl);
	buf[cpl] = '\0';
	fprintf(fp, "#=GC %-*.*s %s\n", namewidth+typewidth, namewidth+typewidth, 
		msa->gc_tag[j], buf);
      }
    }
  fprintf(fp, "//\n");
  free(buf);
}





/* Format of a GF line:
 *    #=GF <featurename> <text>
 */
static int
parse_gf(MSA *msa, char *buf)
{
  char *gf;
  char *featurename;
  char *text;
  char *s;

  s = buf;
  if ((gf          = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((featurename = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((text        = sre_strtok(&s, "\n",       NULL)) == NULL) return 0;
  while (*text && (*text == ' ' || *text == '\t')) text++;

  if      (strcmp(featurename, "ID") == 0) 
    msa->name                 = sre_strdup(text, -1);
  else if (strcmp(featurename, "AC") == 0) 
    msa->acc                  = sre_strdup(text, -1);
  else if (strcmp(featurename, "DE") == 0) 
    msa->desc                 = sre_strdup(text, -1);
  else if (strcmp(featurename, "AU") == 0) 
    msa->au                   = sre_strdup(text, -1);
  else if (strcmp(featurename, "GA") == 0) 
    {				/* Pfam has GA1, GA2. Rfam just has GA1. */
      s = text;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
      msa->cutoff[MSA_CUTOFF_GA1]        = atof(text);
      msa->cutoff_is_set[MSA_CUTOFF_GA1] = TRUE;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) != NULL) {
	msa->cutoff[MSA_CUTOFF_GA2]        = atof(text);
	msa->cutoff_is_set[MSA_CUTOFF_GA2] = TRUE;
      }
    }
  else if (strcmp(featurename, "NC") == 0) 
    {
      s = text;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
      msa->cutoff[MSA_CUTOFF_NC1]        = atof(text);
      msa->cutoff_is_set[MSA_CUTOFF_NC1] = TRUE;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) != NULL) {
	msa->cutoff[MSA_CUTOFF_NC2]        = atof(text);
	msa->cutoff_is_set[MSA_CUTOFF_NC2] = TRUE;
      }
    }
  else if (strcmp(featurename, "TC") == 0) 
    {
      s = text;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
      msa->cutoff[MSA_CUTOFF_TC1]        = atof(text);
      msa->cutoff_is_set[MSA_CUTOFF_TC1] = TRUE;
      if ((text = sre_strtok(&s, WHITESPACE, NULL)) != NULL) {
	msa->cutoff[MSA_CUTOFF_TC2]        = atof(text);
	msa->cutoff_is_set[MSA_CUTOFF_TC2] = TRUE;
      }
    }
  else 
    MSAAddGF(msa, featurename, text);

  return 1;
}


/* Format of a GS line:
 *    #=GS <seqname> <featurename> <text>
 */
static int
parse_gs(MSA *msa, char *buf)
{
  char *gs;
  char *seqname;
  char *featurename;
  char *text; 
  int   seqidx;
  char *s;

  s = buf;
  if ((gs          = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((seqname     = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((featurename = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((text        = sre_strtok(&s, "\n",       NULL)) == NULL) return 0;
  while (*text && (*text == ' ' || *text == '\t')) text++;
  
  /* GS usually follows another GS; guess lastidx+1
   */
  seqidx = MSAGetSeqidx(msa, seqname, msa->lastidx+1);
  msa->lastidx = seqidx;

  if (strcmp(featurename, "WT") == 0)
    {
      msa->wgt[seqidx]          = atof(text);
      msa->flags |= MSA_SET_WGT;
    }

  else if (strcmp(featurename, "AC") == 0)
    MSASetSeqAccession(msa, seqidx, text);

  else if (strcmp(featurename, "DE") == 0)
    MSASetSeqDescription(msa, seqidx, text);

  else				
    MSAAddGS(msa, featurename, seqidx, text);

  return 1;
}

/* Format of a GC line:
 *    #=GC <featurename> <text>
 */
static int 
parse_gc(MSA *msa, char *buf)
{
  char *gc;
  char *featurename;
  char *text; 
  char *s;
  int   len;

  s = buf;
  if ((gc          = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((featurename = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((text        = sre_strtok(&s, WHITESPACE, &len)) == NULL) return 0;
  
  if (strcmp(featurename, "SS_cons") == 0)
    sre_strcat(&(msa->ss_cons), -1, text, len);
  else if (strcmp(featurename, "SA_cons") == 0)
    sre_strcat(&(msa->sa_cons), -1, text, len);
  else if (strcmp(featurename, "RF") == 0)
    sre_strcat(&(msa->rf), -1, text, len);
  else
    MSAAppendGC(msa, featurename, text);

  return 1;
}

/* Format of a GR line:
 *    #=GR <seqname> <featurename> <text>
 */
static int
parse_gr(MSA *msa, char *buf)
{
  char *gr;
  char *seqname;
  char *featurename;
  char *text;
  int   seqidx;
  int   len;
  int   j;
  char *s;

  s = buf;
  if ((gr          = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((seqname     = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((featurename = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((text        = sre_strtok(&s, WHITESPACE, &len)) == NULL) return 0;

  /* GR usually follows sequence it refers to; guess msa->lastidx */
  seqidx = MSAGetSeqidx(msa, seqname, msa->lastidx);
  msa->lastidx = seqidx;

  if (strcmp(featurename, "SS") == 0) 
    {
      if (msa->ss == NULL)
	{
	  msa->ss    = MallocOrDie(sizeof(char *) * msa->nseqalloc);
	  msa->sslen = MallocOrDie(sizeof(int)    * msa->nseqalloc);
	  for (j = 0; j < msa->nseqalloc; j++)
	    {
	      msa->ss[j]    = NULL;
	      msa->sslen[j] = 0;
	    }
	}
      msa->sslen[seqidx] = sre_strcat(&(msa->ss[seqidx]), msa->sslen[seqidx], text, len);
    }
  else if (strcmp(featurename, "SA") == 0)
    {
      if (msa->sa == NULL)
	{
	  msa->sa    = MallocOrDie(sizeof(char *) * msa->nseqalloc);
	  msa->salen = MallocOrDie(sizeof(int)    * msa->nseqalloc);
	  for (j = 0; j < msa->nseqalloc; j++) 
	    {
	      msa->sa[j]    = NULL;
	      msa->salen[j] = 0;
	    }
	}
      msa->salen[seqidx] = sre_strcat(&(msa->sa[seqidx]), msa->salen[seqidx], text, len);
    }
  else 
    MSAAppendGR(msa, featurename, seqidx, text);

  return 1;
}


/* comments are simply stored verbatim, not parsed
 */
static int
parse_comment(MSA *msa, char *buf)
{
  char *s;
  char *comment;

  s = buf + 1;			               /* skip leading '#' */
  if (*s == '\n') { *s = '\0'; comment = s; }  /* deal with blank comment */
  else if ((comment = sre_strtok(&s, "\n", NULL)) == NULL) return 0;
  
  MSAAddComment(msa, comment);
  return 1;
}

static int
parse_sequence(MSA *msa, char *buf)
{
  char *s;
  char *seqname;
  char *text;
  int   seqidx;
  int   len;

  s = buf;
  if ((seqname     = sre_strtok(&s, WHITESPACE, NULL)) == NULL) return 0;
  if ((text        = sre_strtok(&s, WHITESPACE, &len)) == NULL) return 0; 
  
  /* seq usually follows another seq; guess msa->lastidx +1 */
  seqidx = MSAGetSeqidx(msa, seqname, msa->lastidx+1);
  msa->lastidx = seqidx;

  msa->sqlen[seqidx] = sre_strcat(&(msa->aseq[seqidx]), msa->sqlen[seqidx], text, len);
  return 1;
}