home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / science / readseq / readseq.c < prev    next >
C/C++ Source or Header  |  1993-07-17  |  37KB  |  1,141 lines
  1. /* File: readseq.c
  2.  * main() program for ureadseq.c, ureadseq.h
  3.  *
  4.  * Reads and writes nucleic/protein sequence in various
  5.  * formats. Data files may have multiple sequences.
  6.  *
  7.  * Copyright 1990 by d.g.gilbert
  8.  * biology dept., indiana university, bloomington, in 47405
  9.  * e-mail: gilbertd@bio.indiana.edu
  10.  *
  11.  * This program may be freely copied and used by anyone.
  12.  * Developers are encourged to incorporate parts in their
  13.  * programs, rather than devise their own private sequence
  14.  * format.
  15.  *
  16.  * This should compile and run with any ANSI C compiler.
  17.  * Please advise me of any bugs, additions or corrections.
  18.  *
  19.  */
  20.  
  21. const char *title
  22.     = "readSeq (1Feb93), multi-format molbio sequence reader.\n";
  23.  
  24.  /*  History
  25.   27 Feb 90.  1st release to public.
  26.    4 Mar 90.  + Gary Olsen format
  27.               + case change
  28.               * minor corrections to NBRF,EMBL,others
  29.               * output 1 file per sequence for gcg, unknown
  30.               * define -DNOSTR for c-libraries w/o strstr
  31.               - readseq.p, pascal version, becomes out-of-date
  32.   24 May 90.  + Phylip 3.2 output format (no input)
  33.   20 Jul 90.  + Phylip 3.3 output (no input yet)
  34.               + interactive output re-direction
  35.               + verbose progress info
  36.               * interactive help output
  37.               * dropped line no.s on NBRF output
  38.               * patched in HyperGCG XCMD corrections,
  39.                 - except for seq. documentation handling
  40.               * dropped the IG special nuc codes, as IG has
  41.                 adopted the standard IUB codes (now if only
  42.                 everyone would adopt a standard format !)
  43.   11 Oct 90.  * corrected bug in reading/writing of EMBL format
  44.  
  45.   17 Oct 91.  * corrected bug in reading Olsen format
  46.                 (serious-deletion)
  47.   10 Nov 91.  * corrected bug in reading some GCG format files
  48.                 (serious-last line duplicated)
  49.               + add format name parsing (-fgb, -ffasta, ...)
  50.               + Phylip v3.4 output format (== v3.2, sequential)
  51.               + add checksum output to all forms that have document
  52.               + skip mail headers in seq file
  53.               + add pipe for standard input == seq file (with -p)
  54.               * fold in parts of MacApp Seq object
  55.               * strengthen format detection
  56.               * clarify program structure
  57.               * remove fixed sequence size limit (now dynamic, sizeof memory)
  58.               * check and fold in accumulated bug reports:
  59.               *   Now ANSI-C fopen(..,"w") & check open failure
  60.               *   Define -DFIXTOUPPER for nonANSI C libraries that mess
  61.                   up toupper/tolower
  62.               = No command-line changes; callers of readseq main() should be okay
  63.               - ureadseq.h functions have changed; client programs need to note.
  64.               + added Unix and VMS Make scripts, including validation tests
  65.  
  66.    4 May 92.  + added 32 bit CRC checksum as alternative to GCG 6.5bit checksum
  67.                 (-DBIGCHECKSUM)
  68.     Aug 92    = fixed Olsen format input to handle files w/ more sequences,
  69.                 not to mess up when more than one seq has same identifier,
  70.                 and to convert number masks to symbols.
  71.               = IG format fix to understand ^L
  72.  
  73.   25-30 Dec 92
  74.               * revised command-line & interactive interface.  Suggested form is now
  75.                   readseq infile -format=genbank -output=outfile -item=1,3,4 ...
  76.                 but remains compatible with prior commandlines:
  77.                   readseq infile -f2 -ooutfile -i3 ...
  78.               + added GCG MSF multi sequence file format
  79.               + added PIR/CODATA format
  80.               + added NCBI ASN.1 sequence file format
  81.               + added Pretty, multi sequence pretty output (only)
  82.               + added PAUP multi seq format
  83.               + added degap option
  84.               + added Gary Williams (GWW, G.Williams@CRC.AC.UK) reverse-complement option.
  85.               + added support for reading Phylip formats (interleave & sequential)
  86.               * string fixes, dropped need for compiler flags NOSTR, FIXTOUPPER, NEEDSTRCASECMP
  87.               * changed 32bit checksum to default, -DSMALLCHECKSUM for GCG version
  88.  
  89.    1Feb93
  90.               = revert GenBank output to a fixed left number width which 
  91.                other software depends on.
  92.           = fix for MSF input to handle symbols in names
  93.           = fix bug for possible memory overrun when truncating seqs for
  94.         Phylip or Paup formats (thanks Anthony Persechini)
  95.  
  96.  */
  97.  
  98.  
  99.  
  100. /*
  101.    Readseq has been tested with:
  102.       Macintosh MPW C
  103.       GNU gcc
  104.       SGI cc
  105.       VAX-VMS cc
  106.    Any ANSI C compiler should be able to handle this.
  107.    Old-style C compilers barf all over the source.
  108.  
  109.  
  110. How do I build the readseq program if I have an Ansi C compiler?
  111. #--------------------
  112. # Unix ANSI C
  113. # Use the supplied Makefile this way:
  114. %  make CC=name-of-c-compiler
  115. # OR do this...
  116. % gcc readseq.c ureadseq.c -o readseq
  117.  
  118. #--------------------
  119. $!VAX-VMS cc
  120. $! Use the supplied Make.Com this way:
  121. $  @make
  122. $! OR, do this:
  123. $ cc readseq, ureadseq
  124. $ link readseq, ureadseq, sys$library:vaxcrtl/lib
  125. $ readseq :== $ MyDisk:[myacct]readseq
  126.  
  127. #--------------------
  128. # Macintosh Simple Input/Output Window application
  129. # requires MPW-C and SIOW library (from APDA)
  130. # also uses files macinit.c, macinit.r, readseqSIOW.make
  131. #
  132. Buildprogram readseqSIOW
  133.  
  134. #--------------------
  135. #MPW-C v3 tool
  136. C  ureadseq.c
  137. C  readseq.c
  138. link -w -o readseq -t MPST -c 'MPS ' À
  139.    readseq.c.o Ureadseq.c.o À
  140.     "{Libraries}"Interface.o À
  141.     "{Libraries}"ToolLibs.o À
  142.     "{Libraries}"Runtime.o À
  143.     "{CLibraries}"StdClib.o
  144. readseq -i1 ig.seq
  145.  
  146. # MPW-C with NCBI tools
  147.  
  148. set NCBI "{Boot}@molbio:ncbi:"; EXPORT NCBI
  149. set NCBILIB1  "{NCBI}"lib:libncbi.o; export NCBILIB1
  150. set NCBILIB2  "{NCBI}"lib:libncbiobj.o; export NCBILIB2
  151. set NCBILIB3  "{NCBI}"lib:libncbicdr.o; export NCBILIB3
  152. set NCBILIB4  "{NCBI}"lib:libvibrant.o; export NCBILIB4
  153.  
  154. C  ureadseq.c
  155. C  -d NCBI -i "{NCBI}"include: ureadasn.c
  156. C  -d NCBI -i "{NCBI}"include: readseq.c
  157. link -w -o readseq -t MPST -c 'MPS ' À
  158.    ureadseq.c.o ureadasn.c.o readseq.c.o  À
  159.     {NCBILIB4} {NCBILIB2} {NCBILIB1} À
  160.     "{Libraries}"Interface.o À
  161.     "{Libraries}"ToolLibs.o À
  162.     "{Libraries}"Runtime.o À
  163.     "{CLibraries}"CSANELib.o À
  164.     "{CLibraries}"Math.o À
  165.     "{CLibraries}"StdClib.o
  166.  
  167. ===========================================================*/
  168.  
  169.  
  170.  
  171. #include <stdio.h>
  172. #include <stdlib.h>    /* for diverse functions */
  173. #include <string.h>
  174. #include <ctype.h>
  175.  
  176. #include "ureadseq.h"
  177.  
  178. #pragma segment readseq
  179.  
  180.  
  181.  
  182. static char inputfilestore[256], *inputfile = inputfilestore;
  183.  
  184. const char *formats[kMaxFormat+1] = {
  185.     " 1. IG/Stanford",
  186.     " 2. GenBank/GB",
  187.     " 3. NBRF",
  188.     " 4. EMBL",
  189.     " 5. GCG",
  190.     " 6. DNAStrider",
  191.     " 7. Fitch",
  192.     " 8. Pearson/Fasta",
  193.     " 9. Zuker (in-only)",
  194.     "10. Olsen (in-only)",
  195.     "11. Phylip3.2",
  196.     "12. Phylip",
  197.     "13. Plain/Raw",
  198.     "14. PIR/CODATA",
  199.     "15. MSF",
  200.     "16. ASN.1",
  201.     "17. PAUP/NEXUS",
  202.     "18. Pretty (out-only)",
  203.     "" };
  204.  
  205. #define kFormCount  30
  206. #define kMaxFormName 15
  207.  
  208. const  struct formatTable {
  209.   char  *name;
  210.   short num;
  211.   } formname[] = {
  212.     {"ig",  kIG},
  213.     {"stanford", kIG},
  214.     {"genbank", kGenBank},
  215.     {"gb", kGenBank},
  216.     {"nbrf", kNBRF},
  217.     {"embl", kEMBL},
  218.     {"gcg", kGCG},
  219.     {"uwgcg", kGCG},
  220.     {"dnastrider", kStrider},
  221.     {"strider", kStrider},
  222.     {"fitch", kFitch},
  223.     {"pearson", kPearson},
  224.     {"fasta", kPearson},
  225.     {"zuker", kZuker},
  226.     {"olsen", kOlsen},
  227.     {"phylip", kPhylip},
  228.     {"phylip3.2", kPhylip2},
  229.     {"phylip3.3", kPhylip3},
  230.     {"phylip3.4", kPhylip4},
  231.     {"phylip-interleaved", kPhylip4},
  232.     {"phylip-sequential", kPhylip2},
  233.     {"plain", kPlain},
  234.     {"raw", kPlain},
  235.     {"pir", kPIR},
  236.     {"codata", kPIR},
  237.     {"asn.1", kASN1},
  238.     {"msf", kMSF},
  239.     {"paup", kPAUP},
  240.     {"nexus", kPAUP},
  241.     {"pretty", kPretty},
  242.   };
  243.  
  244. const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
  245.  
  246. /* GWW table for getting the complement of a nucleotide (IUB codes) */
  247. /*                     ! "#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[ \]^_`abcdefghijklmnopqrstuvwxyz{|}~ */
  248. const char compl[] = " !\"#$%&'()*+,-./0123456789:;<=>?@TVGHNNCDNNMNKNNYRYSAABWNRN[\\]^_`tvghnncdnnmnknnyrysaabwnrn{|}~";
  249.  
  250.  
  251.  
  252. char *formatstr( short format)
  253. {
  254.   if (format < 1 || format > kMaxFormat) {
  255.     switch (format) {
  256.       case kASNseqentry :
  257.       case kASNseqset   : return formats[kASN1-1];
  258.       case kPhylipInterleave:
  259.       case kPhylipSequential: return formats[kPhylip-1];
  260.       default: return "(unknown)";
  261.       }
  262.     }
  263.   else return formats[format-1];
  264. }
  265.  
  266. int parseformat( char *name)
  267. {
  268. #define kDupmatch  -2
  269.   int   namelen, maxlen, i, match, matchat;
  270.   char  lname[kMaxFormName+1];
  271.  
  272.   skipwhitespace(name);
  273.   namelen = strlen(name);
  274.   if (namelen == 0)
  275.     return kNoformat;
  276.   else if (isdigit(*name)) {
  277.     i = atol( name);
  278.     if (i < kMinFormat | i > kMaxFormat) return kNoformat;
  279.     else return i;
  280.     }
  281.  
  282.   /* else match character name */
  283.   maxlen = min( kMaxFormName, namelen);
  284.   for (i=0; i<maxlen; i++) lname[i] = to_lower(name[i]);
  285.   lname[maxlen]=0;
  286.   matchat = kNoformat;
  287.  
  288.   for (i=0; i<kFormCount; i++) {
  289.     match = strncmp( lname, formname[i].name, maxlen);
  290.     if (match == 0) {
  291.       if (strlen(formname[i].name) == namelen) return (formname[i].num);
  292.       else if (matchat == kNoformat) matchat = i;
  293.       else matchat = kDupmatch; /* 2 or more partial matches */
  294.       }
  295.     }
  296.   if (matchat == kNoformat || matchat == kDupmatch)
  297.     return kNoformat;
  298.   else
  299.     return formname[matchat].num;
  300. }
  301.  
  302.  
  303.  
  304. static void dumpSeqList(char *list, short format)
  305. {
  306.   long i, l, listlen;
  307.   char s[256];
  308.  
  309.   listlen = strlen(list);
  310.   printf("Sequences in %s  (format is %s)\n", inputfile, formatstr(format));
  311.   for (i=0, l=0; i < listlen; i++) {
  312.     if (list[i] == (char)NEWLINE) {
  313.       s[l] = '\0'; l = 0;
  314.       puts(s);
  315.       }
  316.     else if (l < 255)
  317.       s[l++] = list[i];
  318.     }
  319.   putchar('\n');
  320. }
  321.  
  322.  
  323.  
  324. void usage()
  325. {
  326.   short   i, midi;
  327.  
  328.   fprintf(stderr,title);
  329.   fprintf(stderr,
  330.   "usage: readseq [-options] in.seq > out.seq\n");
  331.   fprintf(stderr," options\n");
  332. /* ? add -d[igits] to allow digits in sequence data, &/or option to specify seq charset !? */
  333.   fprintf(stderr, "    -a[ll]         select All sequences\n");
  334.   fprintf(stderr, "    -c[aselower]   change to lower case\n");
  335.   fprintf(stderr, "    -C[ASEUPPER]   change to UPPER CASE\n");
  336.   fprintf(stderr, "    -degap[=-]     remove gap symbols\n");
  337.   fprintf(stderr, "    -i[tem=2,3,4]  select Item number(s) from several\n");
  338.   fprintf(stderr, "    -l[ist]        List sequences only\n");
  339.   fprintf(stderr, "    -o[utput=]out.seq  redirect Output\n");
  340.   fprintf(stderr, "    -p[ipe]        Pipe (command line, <stdin, >stdout)\n");
  341.   fprintf(stderr, "    -r[everse]     change to Reverse-complement\n");
  342.   fprintf(stderr, "    -v[erbose]     Verbose progress\n");
  343.   fprintf(stderr, "    -f[ormat=]#    Format number for output,  or\n");
  344.   fprintf(stderr, "    -f[ormat=]Name Format name for output:\n");
  345.   midi = (kMaxFormat+1) / 2;
  346.   for (i = kMinFormat-1; i < midi; i++)
  347.    fprintf( stderr, "        %-20s      %-20s\n",
  348.     formats[i], formats[midi+i]);
  349.  
  350.   /* new output format options, esp. for pretty format: */
  351.   fprintf(stderr, "     \n");
  352.   fprintf(stderr, "   Pretty format options: \n");
  353.   fprintf(stderr, "    -wid[th]=#            sequence line width\n");
  354.   fprintf(stderr, "    -tab=#                left indent\n");
  355.   fprintf(stderr, "    -col[space]=#         column space within sequence line on output\n");
  356.   fprintf(stderr, "    -gap[count]           count gap chars in sequence numbers\n");
  357.   fprintf(stderr, "    -nameleft, -nameright[=#]   name on left/right side [=max width]\n");
  358.   fprintf(stderr, "    -nametop              name at top/bottom\n");
  359.   fprintf(stderr, "    -numleft, -numright   seq index on left/right side\n");
  360.   fprintf(stderr, "    -numtop, -numbot      index on top/bottom\n");
  361.   fprintf(stderr, "    -match[=.]            use match base for 2..n species\n");
  362.   fprintf(stderr, "    -inter[line=#]        blank line(s) between sequence blocks\n");
  363.  
  364.   /******  not ready yet
  365.   fprintf(stderr, "    -code=none,rtf,postscript,ps   code syntax\n");
  366.   fprintf(stderr, "    -namefont=, -numfont=, -seqfont=font   font choice\n");
  367.   fprintf(stderr, "       font suggestions include times,courier,helvetica\n");
  368.   fprintf(stderr, "    -namefontsize=, -numfontsize=, -seqfontsize=#\n");
  369.   fprintf(stderr, "       fontsize suggestions include 9,10,12,14\n");
  370.   fprintf(stderr, "    -namefontstyle=, -numfontstyle=, -seqfontstyle= style  fontstyle for names\n");
  371.   fprintf(stderr, "       fontstyle options are plain,italic,bold,bold-italic\n");
  372.   ******/
  373. }
  374.  
  375. void erralert(short err)
  376. {
  377.   switch (err) {
  378.     case 0  :
  379.       break;
  380.     case eFileNotFound: fprintf(stderr, "File not found: %s\n", inputfile);
  381.       break;
  382.     case eFileCreate: fprintf(stderr, "Can't open output file.\n");
  383.       break;
  384.     case eASNerr: fprintf(stderr, "Error in ASN.1 sequence routines.\n");
  385.       break;
  386.     case eNoData: fprintf(stderr, "No data in file.\n");
  387.       break;
  388.     case eItemNotFound: fprintf(stderr, "Specified item not in file.\n");
  389.       break;
  390.     case eUnequalSize:  fprintf(stderr,
  391.       "This format requires equal length sequences.\nSequence truncated or padded to fit.\n");
  392.       break;
  393.     case eUnknownFormat: fprintf(stderr, "Error: this format is unknown to me.\n");
  394.       break;
  395.     case eOneFormat: fprintf(stderr,
  396.       "Warning: This format permits only 1 sequence per file.\n");
  397.       break;
  398.     case eMemFull: fprintf(stderr, "Out of storage memory. Sequence truncated.\n");
  399.       break;
  400.     default: fprintf(stderr, "readSeq error = %d\n", err);
  401.       break;
  402.     }
  403. } /* erralert */
  404.  
  405.  
  406. int chooseFormat( boolean quietly)
  407. {
  408.   char  sform[128];
  409.   int   midi, i, outform;
  410.  
  411.     if (quietly)
  412.       return kPearson;  /* default */
  413.     else {
  414.       midi = (kMaxFormat+1) / 2;
  415.       for (i = kMinFormat-1; i < midi; i++)
  416.         fprintf( stderr, "        %-20s      %-20s\n",
  417.                         formats[i], formats[midi+i]);
  418.       fprintf(stderr,"\nChoose an output format (name or #): \n");
  419.       gets(sform);
  420.       outform = parseformat(sform);
  421.       if (outform == kNoformat) outform = kPearson;
  422.       return outform;
  423.       }
  424. }
  425.  
  426.  
  427.  
  428. /* read paramater(s) */
  429.  
  430. boolean checkopt( boolean casesense, char *sopt, const char *smatch, short minword)
  431. {
  432.   long  lenopt, lenmatch;
  433.   boolean result;
  434.   short minmaxw;
  435.  
  436.   lenopt = strlen(sopt);
  437.   lenmatch= strlen(smatch);
  438.   minmaxw= max(minword, min(lenopt, lenmatch));
  439.  
  440.   if (casesense)
  441.     result= (!strncmp( sopt, smatch, minmaxw));
  442.   else
  443.     result= (!Strncasecmp( sopt, smatch, minmaxw ));
  444.   /* if (result) { */
  445.     /* fprintf(stderr,"true checkopt(opt=%s,match=%s,param=%s)\n", sopt, smatch, *sparam); */
  446.   /*  } */
  447.   return result;
  448. }
  449.  
  450.  
  451. #define   kMaxwhichlist  50
  452.  
  453. /* global for readopt(), main() */
  454. boolean   chooseall = false, quietly = false, gotinputfile = false,
  455.           listonly = false, closeout = false, verbose = false,
  456.           manyout = false, dolower = false, doupper = false, doreverse= false,
  457.           askout  = true, dopipe= false, interleaved = false;
  458. short     nfile = 0, iwhichlist=0, nwhichlist = 0;
  459. short     whichlist[kMaxwhichlist+1];
  460. long      whichSeq = 0, outform = kNoformat;
  461. char      onamestore[128], *oname = onamestore;
  462. FILE      *foo = NULL;
  463.  
  464. void resetGlobals()
  465. /* need this when used from SIOW, as these globals are not reinited automatically
  466. between calls to local main() */
  467. {
  468.   chooseall = false; quietly = false; gotinputfile = false;
  469.   listonly = false; closeout = false; verbose = false;
  470.   manyout = false; dolower = false; doupper = false; doreverse= false;
  471.   askout  = true; dopipe= false; interleaved = false;
  472.   nfile = 0; iwhichlist=0; nwhichlist = 0;
  473.   whichSeq = 0; outform = kNoformat;
  474.   oname = onamestore;
  475.   foo = NULL;
  476.  
  477.   gPrettyInit(gPretty);
  478. }
  479.  
  480.  
  481. #define kOptOkay  1
  482. #define kOptNone  0
  483.  
  484. int readopt( char *sopt)
  485. {
  486.   char    sparamstore[256], *sparam= sparamstore;
  487.   short   n, slen= strlen(sopt);
  488.  
  489.   /* fprintf(stderr,"readopt( %s) == ", sopt); */
  490.  
  491.   if (*sopt == '?') {
  492.     usage();
  493.     return kOptNone;   /*? eOptionBad or kOptNone */
  494.     }
  495.  
  496.   else if (*sopt == '-') {
  497.  
  498.     char *cp= strchr(sopt,'=');
  499.     *sparam= '\0';
  500.     if (cp) {
  501.       strcpy(sparam, cp+1);
  502.       *cp= 0;
  503.       }
  504.  
  505.     if (checkopt( false, sopt, "-help", 2)) {
  506.       usage();
  507.       return kOptNone;
  508.       }
  509.  
  510.     if (checkopt( false, sopt, "-all", 2)) {
  511.       whichSeq= 1; chooseall= true;
  512.       return kOptOkay;
  513.       }
  514.  
  515.     if (checkopt( false, sopt, "-colspace", 4)) { /* test before -c[ase] */
  516.       n= atoi( sparam);
  517.       gPretty.spacer = n;
  518.       return kOptOkay;
  519.       }
  520.  
  521.     if (checkopt( true, sopt, "-caselower", 2)) {
  522.       dolower= true;
  523.       return kOptOkay;
  524.       }
  525.     if (checkopt( true, sopt, "-CASEUPPER", 2)) {
  526.       doupper= true;
  527.       return kOptOkay;
  528.       }
  529.  
  530.     if (checkopt( false, sopt, "-pipe", 2)) {
  531.       dopipe= true; askout= false;
  532.       return kOptOkay;
  533.       }
  534.  
  535.     if (checkopt( false, sopt, "-list", 2)) {
  536.       listonly = true; askout = false;
  537.       return kOptOkay;
  538.       }
  539.  
  540.     if (checkopt( false, sopt, "-reverse", 2)) {
  541.       doreverse = true;
  542.       return kOptOkay;
  543.       }
  544.  
  545.     if (checkopt( false, sopt, "-verbose", 2)) {
  546.       verbose = true;
  547.       return kOptOkay;
  548.       }
  549.  
  550.     if (checkopt( false, sopt, "-match", 5)) {
  551.       gPretty.domatch= true;
  552.       if (*sparam >= ' ') gPretty.matchchar= *sparam;
  553.       return kOptOkay;
  554.       }
  555.     if (checkopt( false, sopt, "-degap", 4)) {
  556.       gPretty.degap= true;
  557.       if (*sparam >= ' ') gPretty.gapchar= *sparam;
  558.       return kOptOkay;
  559.       }
  560.  
  561.     if (checkopt( false, sopt, "-interline", 4)) {
  562.       gPretty.interline= atoi( sparam);
  563.       return kOptOkay;
  564.       }
  565.  
  566.     if (checkopt( false, sopt, "-item", 2)) {
  567.       char  *cp = sparam;
  568.       nwhichlist= 0;
  569.       whichlist[0]= 0;
  570.       if (*cp == 0) cp= sopt+2; /* compatible w/ old way */
  571.       do {
  572.         while (*cp!=0 && !isdigit(*cp)) cp++;
  573.         if (*cp!=0) {
  574.           n = atoi( cp);
  575.           whichlist[nwhichlist++]= n;
  576.           while (*cp!=0 && isdigit(*cp)) cp++;
  577.           }
  578.       } while (*cp!=0 && n>0 && nwhichlist<kMaxwhichlist);
  579.       whichlist[nwhichlist++]= 0; /* 0 == stopsign for loop */
  580.       whichSeq= max(1,whichlist[0]); iwhichlist= 1;
  581.       return kOptOkay;
  582.       }
  583.  
  584.     if (checkopt( false, sopt, "-format", 5)) {/* -format=phylip, -f2, -form=phylip */
  585.       if (*sparam==0) { for (sparam= sopt+2; isalpha(*sparam); sparam++) ; }
  586.       outform = parseformat( sparam);
  587.       return kOptOkay;
  588.       }
  589.     if (checkopt( false, sopt, "-f", 2)) { /* compatible w/ -fphylip prior version */
  590.       if (*sparam==0) sparam= sopt+2;
  591.       outform = parseformat( sparam);
  592.       return kOptOkay;
  593.       }
  594.  
  595.     if (checkopt( false, sopt, "-output", 3)) {/* -output=myseq */
  596.       if (*sparam==0) { for (sparam= sopt+3; isalpha(*sparam); sparam++) ; }
  597.       strcpy( oname, sparam);
  598.       foo = fopen( oname, "w");
  599.       if (!foo) { erralert(eFileCreate); return eFileCreate; }
  600.       closeout = true;
  601.       askout = false;
  602.       return kOptOkay;
  603.       }
  604.     if (checkopt( false, sopt, "-o", 2)) {  /* compatible w/ -omyseq prior version */
  605.       if (*sparam==0) sparam= sopt+2;
  606.       strcpy( oname, sparam);
  607.       foo = fopen( oname, "w");
  608.       if (!foo) { erralert(eFileCreate); return eFileCreate; }
  609.       closeout = true;
  610.       askout = false;
  611.       return kOptOkay;
  612.       }
  613.  
  614.     if (checkopt( false, sopt, "-width", 2)) {
  615.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  616.       n= atoi( sparam);
  617.       if (n>0) gPretty.seqwidth = n;
  618.       return kOptOkay;
  619.       }
  620.  
  621.     if (checkopt( false, sopt, "-tab", 4)) {
  622.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  623.       n= atoi( sparam);
  624.       gPretty.tab = n;
  625.       return kOptOkay;
  626.       }
  627.  
  628.     if (checkopt( false, sopt, "-gapcount", 4)) {
  629.       gPretty.baseonlynum = false;
  630.       /* if (*sparam >= ' ') gPretty.gapchar= *sparam; */
  631.       return kOptOkay;
  632.       }
  633.     if (checkopt( false, sopt, "-nointerleave", 8)) {
  634.       gPretty.noleaves = true;
  635.       return kOptOkay;
  636.       }
  637.  
  638.     if (checkopt( false, sopt, "-nameleft", 7)) {
  639.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  640.       n= atoi( sparam);
  641.       if (n>0 && n<50) gPretty.namewidth =  n;
  642.       gPretty.nameleft= true;
  643.       return kOptOkay;
  644.       }
  645.     if (checkopt( false, sopt, "-nameright", 7)) {
  646.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  647.       n= atoi( sparam);
  648.       if (n>0 && n<50) gPretty.namewidth =  n;
  649.       gPretty.nameright= true;
  650.       return kOptOkay;
  651.       }
  652.     if (checkopt( false, sopt, "-nametop", 6)) {
  653.       gPretty.nametop= true;
  654.       return kOptOkay;
  655.       }
  656.  
  657.     if (checkopt( false, sopt, "-numleft", 6)) {
  658.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  659.       n= atoi( sparam);
  660.       if (n>0 && n<50) gPretty.numwidth =  n;
  661.       gPretty.numleft= true;
  662.       return kOptOkay;
  663.       }
  664.     if (checkopt( false, sopt, "-numright", 6)) {
  665.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  666.       n= atoi( sparam);
  667.       if (n>0 && n<50) gPretty.numwidth =  n;
  668.       gPretty.numright= true;
  669.       return kOptOkay;
  670.       }
  671.  
  672.     if (checkopt( false, sopt, "-numtop", 6)) {
  673.       gPretty.numtop= true;
  674.       return kOptOkay;
  675.       }
  676.     if (checkopt( false, sopt, "-numbottom", 6)) {
  677.       gPretty.numbot= true;
  678.       return kOptOkay;
  679.       }
  680.  
  681.     else {
  682.       usage();
  683.       return eOptionBad;
  684.       }
  685.     }
  686.  
  687.   else {
  688.     strcpy( inputfile, sopt);
  689.     gotinputfile = (*inputfile != 0);
  690.     nfile++;
  691.     return kOptOkay;
  692.     }
  693.  
  694.  /* return kOptNone; -- never here */
  695. }
  696.  
  697.  
  698.  
  699.  
  700. /* this program suffers some as it tries to be a quiet translator pipe
  701.    _and_ a noisy user interactor
  702. */
  703.  
  704. /* return is best for SIOW, okay for others */
  705. #ifdef SIOW
  706. #define Exit(a)   return(a)
  707. siow_main( int argc, char *argv[])
  708.  
  709. #else
  710. #define Exit(a)   exit(a)
  711.  
  712. main( int argc, char *argv[])
  713. #endif
  714. {
  715. boolean   closein = false;
  716. short     ifile, nseq, atseq, format, err = 0, seqtype = kDNA,
  717.           nlines, seqout = 0, phylvers = 2;
  718. long      i, skiplines, seqlen, seqlen0;
  719. unsigned long  checksum= 0, checkall= 0;
  720. char      *seq, *cp, *firstseq = NULL, *seqlist, *progname, tempname[256];
  721. char      seqid[256], *seqidptr = seqid;
  722. char      stempstore[256], *stemp = stempstore;
  723. FILE      *ftmp, *fin, *fout;
  724. long      outindexmax= 0, noutindex= 0, *outindex = NULL;
  725.  
  726. #define exit_main(err) {        \
  727.   if (closeout) fclose(fout);   \
  728.   if (closein) fclose(fin);   \
  729.   if (*tempname!=0) remove(tempname);\
  730.   Exit(err); }
  731.  
  732. #define indexout()  if (interleaved) {\
  733.   if (noutindex>=outindexmax) {\
  734.     outindexmax= noutindex + 20;\
  735.     outindex= (long*) realloc(outindex, sizeof(long)*outindexmax);\
  736.     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }\
  737.     }\
  738.   outindex[noutindex++]= ftell(fout);\
  739.   }
  740.  
  741.  
  742.   resetGlobals();
  743.   foo = stdout;
  744.   progname = argv[0];
  745.   *oname = 0;
  746.   *tempname = 0;
  747.   /* initialize gPretty ?? -- done in header */
  748.  
  749.   for (i=1; i < argc; i++) {
  750.     err= readopt( argv[i]);
  751.     if (err <= 0) exit_main(err);
  752.     }
  753.  
  754.                             /* pipe input from stdin !? */
  755.   if (dopipe && !gotinputfile) {
  756.     int c;
  757.     tmpnam(tempname);
  758.     inputfile = tempname;
  759.     ftmp = fopen( inputfile, "w");
  760.     if (!ftmp) { erralert(eFileCreate); exit_main(eFileCreate); }
  761.     while ((c = getc(stdin)) != EOF) fputc(c, ftmp);
  762.     fclose(ftmp);
  763.     gotinputfile= true;
  764.     }
  765.  
  766.   quietly = (dopipe || (gotinputfile && (listonly || whichSeq != 0)));
  767.  
  768.   if (verbose || (!quietly && !gotinputfile)) fprintf( stderr, title);
  769.   ifile = 1;
  770.  
  771.                             /* UI: Choose output */
  772.   if (askout && !closeout && !quietly) {
  773.     askout = false;
  774.     fprintf(stderr,"\nName of output file (?=help, defaults to display): \n");
  775.     gets(oname= onamestore);
  776.     skipwhitespace(oname);
  777.     if (*oname == '?') { usage(); exit_main(0); }
  778.     else if (*oname != 0) {
  779.       closeout = true;
  780.       foo = fopen( oname, "w");
  781.       if (!foo) { erralert(eFileCreate); exit_main(eFileCreate); }
  782.       }
  783.     }
  784.  
  785.   fout = foo;
  786.   if (outform == kNoformat) outform = chooseFormat(quietly);
  787.  
  788.                           /* set up formats ... */
  789.   switch (outform) {
  790.     case kPhylip2:
  791.       interleaved= false;
  792.       phylvers = 2;
  793.       outform = kPhylip;
  794.       break;
  795.  
  796.     case kPhylip4:
  797.       interleaved= true;
  798.       phylvers = 4;
  799.       outform = kPhylip;
  800.       break;
  801.  
  802.     case kMSF:
  803.     case kPAUP:
  804.       interleaved= true;
  805.       break;
  806.  
  807.     case kPretty:
  808.       gPretty.isactive= true;
  809.       interleaved= true;
  810.       break;
  811.  
  812.     }
  813.  
  814.   if (gPretty.isactive && gPretty.noleaves) interleaved= false;
  815.   if (interleaved) {
  816.     fout = ftmp = tmpfile();
  817.     outindexmax= 30; noutindex= 0;
  818.     outindex = (long*) malloc(outindexmax*sizeof(long));
  819.     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }
  820.     }
  821.  
  822.                         /* big loop over all input files */
  823.   do {
  824.                         /* select next input file */
  825.     gotinputfile = (*tempname != 0);
  826.     while ((ifile < argc) && (!gotinputfile)) {
  827.       if (*argv[ifile] != '-') {
  828.         strcpy( inputfile, argv[ifile]);
  829.         gotinputfile = (*inputfile != 0);
  830.         --nfile;
  831.         }
  832.       ifile++;
  833.       }
  834.  
  835.     while (!gotinputfile) {
  836.       fprintf(stderr,"\nName an input sequence or -option: \n");
  837.       inputfile= inputfilestore;
  838.  
  839.       gets(stemp= stempstore);
  840.       if (*stemp==0) goto fini;  /* !! need this to finish work during interactive use */
  841.       stemp= strtok(stempstore, " \n\r\t");
  842.       while (stemp) {
  843.         err= readopt( stemp); /* will read inputfile if it exists */
  844.         if (err<0) exit_main(err);
  845.         stemp= strtok( NULL, " \n\r\t");
  846.         }
  847.       }
  848.               /* thanks to AJB@UK.AC.DARESBURY.DLVH for this PHYLIP3 fix: */
  849.               /* head for end (interleave if needed) */
  850.     if (*inputfile == 0) break;
  851.  
  852.     format = seqFileFormat( inputfile, &skiplines, &err);
  853.  
  854.     if (err == 0)  {
  855. #ifdef NCBI
  856.       if (format == kASNseqentry || format == kASNseqset)
  857.         seqlist = listASNSeqs( inputfile, skiplines, format, &nseq, &err);
  858.       else
  859. #endif
  860.         seqlist = listSeqs( inputfile, skiplines, format, &nseq, &err);
  861.       }
  862.  
  863.     if (err != 0)
  864.       erralert(err);
  865.  
  866.     else if (listonly) {
  867.       dumpSeqList(seqlist,format);
  868.       free( seqlist);
  869.       }
  870.  
  871.     else {
  872.                                 /* choose whichSeq if needed */
  873.       if (nseq == 1 || chooseall || (quietly && whichSeq == 0)) {
  874.         chooseall= true;
  875.         whichSeq = 1;
  876.         quietly = true; /* no loop */
  877.         }
  878.       else if (whichSeq > nseq && quietly) {
  879.         erralert(eItemNotFound);
  880.         err= eItemNotFound;
  881.         }
  882.       else if (whichSeq > nseq || !quietly) {
  883.         dumpSeqList(seqlist, format);
  884.         fprintf(stderr,"\nChoose a sequence (# or All): \n");
  885.         gets(stemp= stempstore);
  886.         skipwhitespace(stemp);
  887.         if (to_lower(*stemp) == 'a') {
  888.           chooseall= true;
  889.           whichSeq = 1;
  890.           quietly = true; /* !? this means we don't ask for another file 
  891.                             as well as no more whichSeqs... */
  892.           }
  893.         else if (isdigit(*stemp)) whichSeq= atol(stemp);
  894.         else whichSeq= 1; /* default */
  895.         }
  896.       free( seqlist);
  897.  
  898.       if (false /*chooseall*/) {  /* this isn't debugged yet...*/
  899.         fin = fopen(inputfile, "r");
  900.         closein= true;
  901.         }
  902.  
  903.       while (whichSeq > 0 && whichSeq <= nseq) {
  904.                                 /* need to open multiple output files ? */
  905.         manyout = ((chooseall || nwhichlist>1) && nseq > 1
  906.                   && (outform == kPlain || outform == kGCG));
  907.         if (manyout) {
  908.           if ( whichSeq == 1 ) erralert(eOneFormat);
  909.           else if (closeout) {
  910.             sprintf( stemp,"%s_%d", oname, whichSeq);
  911.             freopen( stemp, "w", fout);
  912.             fprintf( stderr,"Writing sequence %d to file %s\n", whichSeq, stemp);
  913.             }
  914.           }
  915.  
  916.         if (closein) {
  917.           /* !! this fails... skips most seqs... */
  918.           /* !! in sequential read, must count seqs already read from whichSeq ... */
  919.           /* need major revision of ureadseq before we can do this */
  920.           atseq= whichSeq-1;
  921.           seqidptr= seqid;
  922.           seq = readSeqFp( whichSeq, fin, skiplines, format,
  923.                           &seqlen, &atseq, &err, seqidptr);
  924.           skiplines= 0;
  925.           }
  926.         else {
  927.           atseq= 0;
  928.           seqidptr= seqid;
  929. #ifdef NCBI
  930.           if (format == kASNseqentry || format == kASNseqset) {
  931.             seqidptr= NULL;
  932.             seq = readASNSeq( whichSeq, inputfile, skiplines, format,
  933.                      &seqlen, &atseq, &err, &seqidptr);
  934.             }
  935.           else
  936. #endif
  937.           seq = readSeq( whichSeq, inputfile, skiplines, format,
  938.                           &seqlen, &atseq, &err, seqidptr);
  939.           }
  940.  
  941.  
  942.         if (gPretty.degap) {
  943.           char *newseq;
  944.           long newlen;
  945.           newseq= compressSeq( gPretty.gapchar, seq, seqlen, &newlen);
  946.           if (newseq) {
  947.             free(seq); seq= newseq; seqlen= newlen;
  948.             }
  949.           }
  950.  
  951.         if (outform == kMSF) checksum= GCGchecksum(seq, seqlen, &checkall);
  952.         else if (verbose) checksum= seqchecksum(seq, seqlen, &checkall);
  953.         if (verbose)
  954.           fprintf( stderr, "Sequence %d, length= %d, checksum= %X, format= %s, id= %s\n",
  955.                 whichSeq, seqlen, checksum, formatstr(format), seqidptr);
  956.  
  957.         if (err != 0) erralert(err);
  958.         else {
  959.                                   /* format fixes that writeseq doesn't do */
  960.           switch (outform) {
  961.             case kPIR:
  962.               if (seqout == 0) fprintf( foo,"\\\\\\\n");
  963.               break;
  964.             case kASN1:
  965.               if (seqout == 0) fprintf( foo, kASN1headline);
  966.               break;
  967.  
  968.             case kPhylip:
  969.               if (seqout == 0) {
  970.                 if (!interleaved) {  /*  bug, nseq is for 1st infile only */
  971.                   if (chooseall) i= nseq; else i=1;
  972.                   if (phylvers >= 4) fprintf(foo," %d %d\n", i, seqlen);
  973.                   else fprintf(foo," %d %d YF\n", i, seqlen);
  974.                   }
  975.                 seqlen0 = seqlen;
  976.                 }
  977.               else if (seqlen != seqlen0) {
  978.                 erralert(eUnequalSize);
  979.                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
  980.                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
  981.                 seqlen = seqlen0;
  982.                 seq[seqlen] = 0; 
  983.                 }
  984.               break;
  985.  
  986.             case kPAUP:
  987.               if (seqout == 0) {
  988.                 seqtype= getseqtype(seq, seqlen);
  989.                 seqlen0 = seqlen;
  990.                 }
  991.               else if (seqlen != seqlen0) {
  992.                 erralert(eUnequalSize);
  993.                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0); 
  994.                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
  995.                 seqlen = seqlen0;
  996.                 seq[seqlen] = 0; 
  997.                 }
  998.               break;
  999.  
  1000.             }
  1001.  
  1002.           if (doupper)
  1003.             for (i = 0; i<seqlen; i++) seq[i] = to_upper(seq[i]);
  1004.           else if (dolower)
  1005.             for (i = 0; i<seqlen; i++) seq[i] = to_lower(seq[i]);
  1006.  
  1007.           if (doreverse) {
  1008.             long  j, k;
  1009.             char  ctemp;
  1010.             for (j=0, k=seqlen-1; j <= k; j++, k--) {
  1011.               ctemp = compl[seq[j] - ' '];
  1012.               seq[j] = compl[seq[k] - ' '];
  1013.               seq[k] = ctemp;
  1014.               }
  1015.             }
  1016.  
  1017.           if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq != NULL) {
  1018.             for (i=0; i<seqlen; i++)
  1019.               if (seq[i]==firstseq[i]) seq[i]= gPretty.matchchar;
  1020.             }
  1021.  
  1022.  
  1023.           if (gPretty.isactive && gPretty.numtop && seqout == 0) {
  1024.             gPretty.numline = 1;
  1025.             indexout();
  1026.             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1027.             gPretty.numline = 2;
  1028.             indexout();
  1029.             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1030.             gPretty.numline = 0;
  1031.             }
  1032.  
  1033.           indexout();
  1034.           nlines = writeSeq( fout, seq, seqlen, outform, seqidptr);
  1035.           seqout++;
  1036.           }
  1037.  
  1038.         if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq == NULL) {
  1039.           firstseq= seq;
  1040.           seq = NULL;
  1041.           }
  1042.         else if (seq!=NULL) { free(seq); seq = NULL; }
  1043.  
  1044. #ifdef NCBI
  1045.        if ( (format == kASNseqentry || format == kASNseqset)
  1046.           && seqidptr && seqidptr!= seqid)
  1047.             free(seqidptr);
  1048. #endif
  1049.         if (chooseall) whichSeq++;
  1050.         else if (iwhichlist<nwhichlist) whichSeq= whichlist[iwhichlist++];
  1051.         else whichSeq= 0;
  1052.         }
  1053.       if (closein) { fclose(fin); closein= false; }
  1054.       }
  1055.     whichSeq  = 0;
  1056.   } while (nfile > 0 || !quietly);
  1057.  
  1058.  
  1059. fini:
  1060.   if (firstseq) { free(firstseq); firstseq= NULL; }
  1061.   if (err || listonly) exit_main(err);
  1062.  
  1063.   if (gPretty.isactive && gPretty.numbot) {
  1064.     gPretty.numline = 2;
  1065.     indexout();
  1066.     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1067.     gPretty.numline = 1;
  1068.     indexout();
  1069.     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1070.     gPretty.numline = 0;
  1071.     }
  1072.  
  1073.   if (outform == kMSF) {
  1074.     if (*oname) cp= oname; else cp= inputfile;
  1075.     fprintf(foo,"\n %s  MSF: %d  Type: N  January 01, 1776  12:00  Check: %d ..\n\n",
  1076.                   cp, seqlen, checkall);
  1077.     }
  1078.  
  1079.   if (outform == kPAUP) {
  1080.     fprintf(foo,"#NEXUS\n");
  1081.     if (*oname) cp= oname; else cp= inputfile;
  1082.     fprintf(foo,"[%s -- data title]\n\n", cp);
  1083.     /* ! now have header lines for each sequence... put them before "begin data;... */
  1084.     }
  1085.  
  1086.   if (outform==kPhylip && interleaved) {
  1087.     if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
  1088.     else fprintf(foo," %d %d YF\n", seqout, seqlen);
  1089.     }
  1090.  
  1091.   if (interleaved) {
  1092.     /* interleave species lines in true output */
  1093.     /* nlines is # lines / sequence */
  1094.     short iline, j, leaf, iseq;
  1095.     char  *s = stempstore;
  1096.  
  1097.     indexout();  noutindex--; /* mark eof */
  1098.  
  1099.     for (leaf=0; leaf<nlines; leaf++) {
  1100.       if (outform == kMSF && leaf == 1) {
  1101.         fputs("//\n\n", foo);
  1102.         }
  1103.       if (outform == kPAUP && leaf==1) {
  1104.         switch (seqtype) {
  1105.           case kDNA     : cp= "dna"; break;
  1106.           case kRNA     : cp= "rna"; break;
  1107.           case kNucleic : cp= "dna"; break;
  1108.           case kAmino   : cp= "protein"; break;
  1109.           case kOtherSeq: cp= "dna"; break;
  1110.           }
  1111.         fprintf(foo,"\nbegin data;\n");
  1112.         fprintf(foo," dimensions ntax=%d nchar=%d;\n", seqout, seqlen);
  1113.         fprintf(foo," format datatype=%s interleave missing=%c", cp, gPretty.gapchar);
  1114.         if (gPretty.domatch) fprintf(foo," matchchar=%c", gPretty.matchchar);
  1115.         fprintf(foo,";\n  matrix\n");
  1116.         }
  1117.  
  1118.       for (iseq=0; iseq<noutindex; iseq++) {
  1119.         fseek(ftmp, outindex[iseq], 0);
  1120.         for (iline=0; iline<=leaf; iline++)
  1121.           if (!fgets(s, 256, ftmp)) *s= 0;
  1122.         if (ftell(ftmp) <= outindex[iseq+1])
  1123.           fputs( s, foo);
  1124.         }
  1125.  
  1126.       for (j=0; j<gPretty.interline; j++)
  1127.         fputs( "\n", foo);  /* some want spacer line */
  1128.       }
  1129.     fclose(ftmp); /* tmp disappears */
  1130.     fout= foo;
  1131.     }
  1132.  
  1133.   if (outform == kASN1)  fprintf( foo, "} }\n");
  1134.   if (outform == kPAUP)  fprintf( foo,";\n  end;\n");
  1135.  
  1136.   if (outindex != NULL) free(outindex);
  1137.   exit_main(0);
  1138. }
  1139.  
  1140.  
  1141.