forester/archive/RIO/others/hmmer/squid/translate_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 /* translate_main.c
  12  *
  13  * translate - create a file of all possible protein ORFs, given
  14  *             an input nucleic acid sequence
  15  *
  16  *
  17  * Not currently compliant w/ HMMER API.
  18  *
  19  * 1.02 Thu Apr 20 16:12:41 1995
  20  *     + incorporated into squid
  21  *     +  -a, -s options added
  22  *
  23  * CVS $Id: translate_main.c,v 1.1.1.1 2005/03/22 08:34:27 cmzmasek Exp $
  24  */
  25
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29 #include "squid.h"
  30 #include "version.h"
  31
  32 #ifdef NEED_GETOPTH
  33 #include <getopt.h>
  34 #endif
  35
  36 #define OPTIONS "ahl:o:qs:"
  37
  38 static char usage[] = "\
  39 Usage: translate [-options] <seqfile>\n\
  40    Translate a nucleic acid sequence into protein ORFs.\n\
  41    Available options are:\n\
  42    -a            : translate in full, with stops; no individual ORFs\n\
  43    -h            : help; show brief usage and version info\n\
  44    -l <minlen>   : report only ORFs greater than minlen (default 20)\n\
  45    -o <outfile>  : save results in output file\n\
  46    -q            : quiet; silence banner, for piping or redirection\n\
  47    -s <stopchar> : with -a, set stop character to <stopchar>\n";
  48
  49 int
  50 main(int argc, char **argv)
  51 {
  52   char        *seqfile;         /* name of seq file to read             */
  53   SQFILE      *seqfp;           /* ptr to opened seq file               */
  54   int          format;          /* format of sequence file              */
  55   char        *seq;             /* ptr to current sequence              */
  56   SQINFO       sqinfo;          /* sequence information                 */
  57   char        *revseq;          /* reverse complement of seq            */
  58   int          start, end;      /* coords of ORF in current seq         */
  59   int          orfnumber;       /* counter for ORFs in current seq      */
  60   char        *aaseq[6];        /* full translations in all 6 frames    */
  61   char        *orf;             /* ptr to translated ORF sequence       */
  62   char        *sptr;            /* ptr into orf                         */
  63   int          len;             /* length of an ORF                     */
  64   int          frame;           /* counter for frames (3..5 are reverse)*/
  65
  66   int          minimum_len;     /* minimum length of ORFs to print out  */
  67   char        *outfile;         /* file to save output in               */
  68   FILE        *ofp;             /* where to direct output               */
  69   char         stopchar;        /* what to use as a stop character      */
  70   int          keepstops;       /* TRUE to do six big ORFs              */
  71   int          quiet;           /* TRUE to silence banner               */
  72
  73   int          optchar;         /* option character */
  74   extern char *optarg;          /* for getopt() */
  75   extern int   optind;          /* for getopt() */
  76
  77   /***********************************************
  78    * Parse the command line
  79    ***********************************************/
  80
  81   format      = SQFILE_UNKNOWN; /* autodetect by default */
  82   minimum_len = 20;
  83   outfile     = NULL;
  84   stopchar    = '*';
  85   keepstops   = FALSE;
  86   quiet       = FALSE;
  87
  88   while ((optchar = getopt(argc, argv, OPTIONS)) != -1)
  89     switch (optchar) {
  90
  91     case 'a': keepstops   = TRUE;         break;
  92     case 'l': minimum_len = atoi(optarg); break;
  93     case 'o': outfile     = optarg;       break;
  94     case 'q': quiet       = TRUE;         break;
  95     case 's': stopchar    = *optarg;      break;
  96
  97     case 'h':
  98       printf("translate %s, %s\n%s\n", RELEASE, RELEASEDATE, usage);
  99       exit(EXIT_SUCCESS);
 100     default:
 101       Die("%s\n", usage);
 102     }
 103
 104   if (argc - optind != 1)
 105     Die("Incorrect number of command line arguments\n%s\n", usage);
 106
 107   seqfile = argv[optind];
 108
 109   /***********************************************
 110    * Open sequence file and output file
 111    ***********************************************/
 112
 113   seqfp = SeqfileOpen(seqfile, format, NULL);
 114   if (seqfp == NULL)
 115     Die("Failed to open sequence file %s\n%s\n",
 116         seqfile, usage);
 117
 118   if (outfile != NULL)
 119     {
 120       if ((ofp = fopen(outfile, "w")) == NULL)
 121         Die("Failed to open output file %s\n", outfile);
 122     }
 123   else
 124     ofp = stdout;
 125
 126
 127   /***********************************************
 128    * Main routine
 129    ***********************************************/
 130
 131   if (! quiet) printf("translate %s, %s\n", RELEASE, RELEASEDATE);
 132
 133   while (ReadSeq(seqfp, seqfp->format, &seq, &sqinfo))
 134     {
 135       s2upper(seq);
 136       revseq = (char *) malloc (sqinfo.len + 1);
 137       revcomp(revseq, seq);
 138       orfnumber = 1;
 139
 140                                 /* Translate seq in all six frames */
 141       aaseq[0] = Translate(seq, stdcode1);
 142       aaseq[1] = Translate(seq + 1, stdcode1);
 143       aaseq[2] = Translate(seq + 2, stdcode1);
 144       aaseq[3] = Translate(revseq, stdcode1);
 145       aaseq[4] = Translate(revseq + 1, stdcode1);
 146       aaseq[5] = Translate(revseq + 2, stdcode1);
 147
 148
 149
 150       if (keepstops)
 151         {                       /* full translation including stops */
 152           for (frame = 0; frame < 6; frame++)
 153             {
 154               fprintf(ofp, "> %s:%d", sqinfo.name, frame);
 155               for (sptr = aaseq[frame]; *sptr; sptr++)
 156                 {
 157                   if (*sptr == '*') *sptr = stopchar;
 158                   if (! ((sptr - aaseq[frame]) % 50)) putc('\n', ofp);
 159                   putc((int) *sptr, ofp);
 160                 }
 161               putc('\n', ofp);
 162             }
 163         }
 164       else
 165         {                       /* Print all decent ORF's in FASTA format */
 166           for (frame = 0; frame < 6; frame++)
 167             {
 168                                 /* initialize strtok on the first ORF;
 169                                    termination codons are '*' symbols */
 170               orf = strtok(aaseq[frame], "*");
 171               while (orf != NULL)
 172                 {
 173                   len = strlen(orf);
 174                   if (len > minimum_len)
 175                     {
 176                                 /* calculate coords */
 177                       start = (orf - aaseq[frame]) * 3 + 1;
 178                       if (frame < 3) start += frame; /* frame corrections */
 179                       else       start -= frame-3;
 180
 181                       if (frame < 3)
 182                         end = start + len * 3;
 183                       else
 184                         {
 185                           start = -1 * (start - sqinfo.len - 1);
 186                           end = start - len * 3;
 187                         }
 188
 189                       fprintf(ofp, "> %s.%d    length %d, nt %d..%d",
 190                               sqinfo.name,
 191                               orfnumber,
 192                               len,
 193                               start,
 194                               end);
 195
 196                       for (sptr = orf; *sptr; sptr++)
 197                         {
 198                           if (! ((sptr - orf) % 50))
 199                             putc('\n', ofp);
 200                           putc((int) *sptr, ofp);
 201                         }
 202                       putc('\n', ofp);
 203
 204                       orfnumber++;
 205                     }
 206
 207                                 /* pick off next orf */
 208                   orf = strtok(NULL, "*");
 209                 }
 210             }
 211         }
 212
 213       for (frame = 0; frame < 6; frame++)
 214         free(aaseq[frame]);
 215       FreeSequence(seq, &sqinfo);
 216       free(revseq);
 217     }
 218
 219   SeqfileClose(seqfp);
 220
 221   /**************************************************
 222    * Successful return to invocation environment
 223    **************************************************/
 224   return 0;
 225 }
 226