sfetch_main.c

hmmer源程序

      
				/* get directory name */
      if ((dbdir = getenv(dbenv[dbidx].envname)) == NULL)
	Die("Environment variable %s is not set.\n%s",
	    dbenv[dbidx].envname, usage);
				/* open ssi file */
      ssifile = (char *) MallocOrDie
	((strlen(dbdir) + strlen(dbenv[dbidx].ssiname) + 2) * sizeof(char));
      sprintf(ssifile, "%s/%s", dbdir, dbenv[dbidx].ssiname);
      if ((status = SSIOpen(ssifile, &ssi)) != 0)
	Die("Failed to open SSI index file %s in directory %s\n%s",
	    dbenv[dbidx].ssiname, dbdir, usage);
				/* get seqfile name, file format, and offset */
      if ((status = SSIGetOffsetByName(ssi, getname, &fh, &ssi_offset)) != 0)
	Die("Failed to find key %s in SSI file %s", getname, ssifile);
      if ((status = SSIFileInfo(ssi, fh, &dbfile, &format)) != 0)
	Die("SSI error: %s", SSIErrorString(status));
      free(ssifile);
				/* set up proper seqfile, with path */
      seqfile = (char *) MallocOrDie 
	((strlen(dbdir) + strlen(dbfile) + 2) * sizeof(char));
      sprintf(seqfile, "%s/%s", dbdir, dbfile);
      used_ssi = TRUE;
      SSIClose(ssi);
    }
  else if (! getfirst)	/* Sequence file. SSI index optional. */
    {
      char *dbfile;
      int   fh;

      ssifile = (char *) MallocOrDie ((strlen(seqfile) + 5) * sizeof(char));
      sprintf(ssifile, "%s.ssi", seqfile);
      if ((status = SSIOpen(ssifile, &ssi)) == 0)
	{
	  SQD_DPRINTF1(("Opened SSI index %s...\n", ssifile));
	  if ((status = SSIGetOffsetByName(ssi, getname, &fh, &ssi_offset)) != 0)
	    Die("Failed to find key %s in SSI file %s", getname, ssifile);
	  if ((status = SSIFileInfo(ssi, fh, &dbfile, &format)) != 0)
	    Die("SSI error: %s", SSIErrorString(status));
	  SSIClose(ssi);
	  used_ssi = TRUE;
	}
      free(ssifile);
    }
  
  /***********************************************
   * Open database file
   ***********************************************/

  if ((seqfp = SeqfileOpen(seqfile, format, NULL)) == NULL)
    Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
  if (used_ssi)
    SeqfilePosition(seqfp, &ssi_offset);

  /***********************************************
   * Open output file
   ***********************************************/
  
  /* Determine output format. Default: use same as input. Override: -F option.
   */
  outfmt = seqfp->format;
  if (outformat != NULL)
    {
      outfmt = String2SeqfileFormat(outformat);
      if (outfmt == SQFILE_UNKNOWN)
	Die("Unknown output format %s\n%s", outformat, usage);
      if (IsAlignmentFormat(outfmt))
	Die("Can't output a single sequence in an alignment format (%s)\n", outformat);
    }
				/* open output file for writing;
				   use stdout by default */
  if (outfile == NULL) outfp = stdout;
  else if ((outfp = fopen(outfile, "w")) == NULL)
    Die("cannot open %s for output\n", outfile);


  /***********************************************
   * Main loop
   ***********************************************/

  /* If this is a simple fetch of the complete sequence
   * in native format, and we've been positioned in the file
   * by an SSI index file, we can just read right from the file,
   * partially bypassing the ReadSeq() API, and probably
   * putting our fingers a little too deep into the seqfp object.
   */
  if (getall && used_ssi && outfmt == format && dbname != NULL)
    {
      char *buf   = NULL;
      int  buflen = 0;
      int  endlen;

      if (dbidx == -1) Die("That's weird. No database index available.");
      endlen = strlen(dbenv[dbidx].entryend);
      fputs(seqfp->buf, outfp); /* always do first line */
      /* fputs("\n", outfp); */ /* buf has its /n */
      while (sre_fgets(&buf, &buflen, seqfp->f) != NULL)
	{
	  if (strncmp(buf, dbenv[dbidx].entryend, endlen) == 0)
	    {
	      if (dbenv[dbidx].addend) fputs(buf, outfp);
	      break;
	    }
	  fputs(buf, outfp);
	}
      if (buf != NULL) free(buf);
    }
  else				/* else, the hard way with ReadSeq */
    {
      seq    = NULL;
      frag   = NULL;

      while (ReadSeq(seqfp, format, &seq, &sqinfo))
	{
	  if (used_ssi)		/* GSI file puts us right on our seq. */
	    break;
	  else if (getfirst)	/* Use the first seq in the file. */
	    break;
	  else if (by_accession && 
		   (sqinfo.flags & SQINFO_ACC) &&
		   strcmp(sqinfo.acc, getname) == 0)
	    break;
	  else if (strcmp(sqinfo.name, getname) == 0)
	    break;

	  FreeSequence(seq, &sqinfo);
	  seq = NULL;
	}
      
      if (seq == NULL) 
	Die("failed to extract the subsequence %s\n%s", getname, usage);
      
      if (getall)
	{
	  from = 1;
	  to   = sqinfo.len;
	}
      else if (from == -1) from = 1;
      else if (to   == -1) to   = sqinfo.len;
      
      if (to > sqinfo.len || from > sqinfo.len)
	Warn("Extracting beyond the length of the sequence");
      if (to < 1 || from < 1)
	Warn("Extracting beyond the beginning of the sequence");
      
      /* check for reverse complement */
      if (to != -1 && from > to)
	{
	  int   swapfoo;	/* temp variable for swapping coords */

	  reverse_complement   = TRUE;
	  swapfoo = from; from = to; to = swapfoo;
	}
      if (to > sqinfo.len) to = sqinfo.len;
      if (from < 1)      from = 1;
      
      if ((frag = (char *) calloc (to-from+2, sizeof(char))) == NULL)
	Die("memory error\n");
      
      if (strncpy(frag, seq+from-1, to-from+1) == NULL)
	Die("strncpy() failed\n"); 
      
      if (sqinfo.flags & SQINFO_SS)
	{
	  if ((ss = (char	*) calloc (to-from+2, sizeof(char))) == NULL)
	    Die("memory error\n");
	  if (strncpy(ss, sqinfo.ss+from-1, to-from+1) == NULL)
	    Die("strncpy() failed\n"); 
	  free(sqinfo.ss);
	  sqinfo.ss = ss;
	}
      
      if (reverse_complement)
	{
	  char			*revfrag;		/* temp variable for reverse complement */
	  int   swapfoo;	/* temp variable for swapping coords back */

	  if ((revfrag = calloc ( to-from+2, sizeof(char))) == NULL)
	    Die("memory failure\n"); 
	  revcomp(revfrag, frag);
	  free(frag);
	  frag = revfrag;
	  swapfoo = from; from = to; to = swapfoo;

	  /* reverse complement nullifies secondary structure */
	  if (sqinfo.flags & SQINFO_SS)
	    { free(sqinfo.ss); sqinfo.flags &= ~SQINFO_SS; }
	}
      
      if (! (sqinfo.flags & SQINFO_ID))
	SetSeqinfoString(&sqinfo, sqinfo.name, SQINFO_ID);
      
      if (! (sqinfo.flags & SQINFO_OLEN))
	{ sqinfo.olen = sqinfo.len; sqinfo.flags |= SQINFO_OLEN; }
      
      sqinfo.len = (to > from) ? to-from+1 : from-to+1;
      sqinfo.flags |= SQINFO_LEN;
      
      if (rename != NULL)
	SetSeqinfoString(&sqinfo, rename, SQINFO_NAME);
      
      source_start = (sqinfo.flags & SQINFO_START) ? sqinfo.start : 1;
      source_stop  = (sqinfo.flags & SQINFO_STOP)  ? sqinfo.stop  : sqinfo.len;
      source_orient= (source_stop > source_start)  ? 1 : -1;
      
      sqinfo.start = source_start + (from- 1) * source_orient;
      sqinfo.stop  = source_start + (to  - 1) * source_orient;
      sqinfo.flags |= SQINFO_START | SQINFO_STOP;
      
      WriteSeq(outfp, outfmt, frag, &sqinfo);
      free(frag);
      FreeSequence(seq, &sqinfo);
    }

  if (outfile != NULL)
    printf("Fragment written to file %s\n", outfile);

  SeqfileClose(seqfp);
  fclose(outfp);
  return(0);
}