forester/archive/RIO/others/hmmer/squid/revcomp_main.c

   1 /*****************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  *****************************************************************/
  10
  11 /* main for revcomp
  12  *
  13  * revcomp - generate reverse complement of sequences
  14  * SRE, Thu Aug  5 17:36:57 1993
  15  * RCS $Id: revcomp_main.c,v 1.1.1.1 2005/03/22 08:34:31 cmzmasek Exp $
  16  */
  17
  18 #include <stdio.h>
  19 #include <string.h>
  20 #include "squid.h"
  21 #include "version.h"
  22
  23 #define OPTIONS "h"
  24
  25 char usage[]  = "Usage: revcomp [-options] <seqfile>\n\
  26   Reverse complement a nucleic acid sequence.\n\
  27   Available options:\n\
  28   -h    : help; print version and usage info\n";
  29
  30 int
  31 main(int argc, char **argv)
  32 {
  33   char  *seqfile;               /* name of sequence file */
  34   SQFILE *dbfp;                 /* open sequence file */
  35   int    fmt;                   /* format of seqfile  */
  36   char  *seq;                   /* sequence */
  37   SQINFO sqinfo;                /* additional sequence info */
  38   char  *rev;                   /* reverse complement */
  39   int    swap;
  40
  41   int    optchar;               /* option character, command line */
  42   extern int   optind;
  43
  44   /***********************************************
  45    * Parse command line
  46    ***********************************************/
  47
  48   fmt = SQFILE_UNKNOWN;
  49
  50   while ((optchar = getopt(argc, argv, OPTIONS)) != -1)
  51     switch (optchar) {
  52     case 'h':
  53       printf("revcomp %s, %s\n%s\n", RELEASE, RELEASEDATE, usage);
  54       exit(EXIT_SUCCESS);
  55     default:
  56       Die("%s\n", usage);
  57     }
  58
  59   if (argc - optind != 1) Die("%s\n", usage);
  60   seqfile = argv[optind];
  61
  62   if ((dbfp = SeqfileOpen(seqfile, fmt, NULL)) == NULL)
  63     Die("Failed to open sequence file %s for reading", seqfile);
  64
  65   while (ReadSeq(dbfp, dbfp->format, &seq, &sqinfo))
  66     {
  67       if ((rev = (char *) malloc ((sqinfo.len + 1) * sizeof(char))) == NULL)
  68         Die("malloc failed");
  69
  70       revcomp(rev, seq);
  71       if (sqinfo.flags & (SQINFO_START | SQINFO_STOP))
  72         {
  73           swap         = sqinfo.start;
  74           sqinfo.start = sqinfo.stop;
  75           sqinfo.stop  = swap;
  76         }
  77         /* secondary structure of reverse strand is nonsense
  78          */
  79       if (sqinfo.flags & SQINFO_SS)
  80         {
  81           sqinfo.flags = sqinfo.flags & ~SQINFO_SS;
  82           free(sqinfo.ss);
  83         }
  84
  85       WriteSeq(stdout, SQFILE_FASTA, rev, &sqinfo);
  86
  87       free(rev);
  88       FreeSequence(seq, &sqinfo);
  89     }
  90
  91   SeqfileClose(dbfp);
  92   return 0;
  93 }