sources/readseq/readseq.c

   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 <string.h>
 173 #include <ctype.h>
 174
 175 #include "ureadseq.h"
 176
 177 #pragma segment readseq
 178
 179
 180
 181 static char inputfilestore[256], *inputfile = inputfilestore;
 182
 183 const char *formats[kMaxFormat+1] = {
 184     " 1. IG/Stanford",
 185     " 2. GenBank/GB",
 186     " 3. NBRF",
 187     " 4. EMBL",
 188     " 5. GCG",
 189     " 6. DNAStrider",
 190     " 7. Fitch",
 191     " 8. Pearson/Fasta",
 192     " 9. Zuker (in-only)",
 193     "10. Olsen (in-only)",
 194     "11. Phylip3.2",
 195     "12. Phylip",
 196     "13. Plain/Raw",
 197     "14. PIR/CODATA",
 198     "15. MSF",
 199     "16. ASN.1",
 200     "17. PAUP/NEXUS",
 201     "18. Pretty (out-only)",
 202     "" };
 203
 204 #define kFormCount  30
 205 #define kMaxFormName 15
 206
 207 const  struct formatTable {
 208   char  *name;
 209   short num;
 210   } formname[] = {
 211     {"ig",  kIG},
 212     {"stanford", kIG},
 213     {"genbank", kGenBank},
 214     {"gb", kGenBank},
 215     {"nbrf", kNBRF},
 216     {"embl", kEMBL},
 217     {"gcg", kGCG},
 218     {"uwgcg", kGCG},
 219     {"dnastrider", kStrider},
 220     {"strider", kStrider},
 221     {"fitch", kFitch},
 222     {"pearson", kPearson},
 223     {"fasta", kPearson},
 224     {"zuker", kZuker},
 225     {"olsen", kOlsen},
 226     {"phylip", kPhylip},
 227     {"phylip3.2", kPhylip2},
 228     {"phylip3.3", kPhylip3},
 229     {"phylip3.4", kPhylip4},
 230     {"phylip-interleaved", kPhylip4},
 231     {"phylip-sequential", kPhylip2},
 232     {"plain", kPlain},
 233     {"raw", kPlain},
 234     {"pir", kPIR},
 235     {"codata", kPIR},
 236     {"asn.1", kASN1},
 237     {"msf", kMSF},
 238     {"paup", kPAUP},
 239     {"nexus", kPAUP},
 240     {"pretty", kPretty},
 241   };
 242
 243 const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
 244
 245 /* GWW table for getting the complement of a nucleotide (IUB codes) */
 246 /*                     ! "#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[ \]^_`abcdefghijklmnopqrstuvwxyz{|}~ */
 247 const char compl[] = " !\"#$%&'()*+,-./0123456789:;<=>?@TVGHNNCDNNMNKNNYRYSAABWNRN[\\]^_`tvghnncdnnmnknnyrysaabwnrn{|}~";
 248
 249
 250
 251 char *formatstr( short format)
 252 {
 253   if (format < 1 || format > kMaxFormat) {
 254     switch (format) {
 255       case kASNseqentry :
 256       case kASNseqset   : return formats[kASN1-1];
 257       case kPhylipInterleave:
 258       case kPhylipSequential: return formats[kPhylip-1];
 259       default: return "(unknown)";
 260       }
 261     }
 262   else return formats[format-1];
 263 }
 264
 265 int parseformat( char *name)
 266 {
 267 #define kDupmatch  -2
 268   int   namelen, maxlen, i, match, matchat;
 269   char  lname[kMaxFormName+1];
 270
 271   skipwhitespace(name);
 272   namelen = strlen(name);
 273   if (namelen == 0)
 274     return kNoformat;
 275   else if (isdigit(*name)) {
 276     i = atol( name);
 277     if (i < kMinFormat | i > kMaxFormat) return kNoformat;
 278     else return i;
 279     }
 280
 281   /* else match character name */
 282   maxlen = min( kMaxFormName, namelen);
 283   for (i=0; i<maxlen; i++) lname[i] = to_lower(name[i]);
 284   lname[maxlen]=0;
 285   matchat = kNoformat;
 286
 287   for (i=0; i<kFormCount; i++) {
 288     match = strncmp( lname, formname[i].name, maxlen);
 289     if (match == 0) {
 290       if (strlen(formname[i].name) == namelen) return (formname[i].num);
 291       else if (matchat == kNoformat) matchat = i;
 292       else matchat = kDupmatch; /* 2 or more partial matches */
 293       }
 294     }
 295   if (matchat == kNoformat || matchat == kDupmatch)
 296     return kNoformat;
 297   else
 298     return formname[matchat].num;
 299 }
 300
 301
 302
 303 static void dumpSeqList(char *list, short format)
 304 {
 305   long i, l, listlen;
 306   char s[256];
 307
 308   listlen = strlen(list);
 309   printf("Sequences in %s  (format is %s)\n", inputfile, formatstr(format));
 310   for (i=0, l=0; i < listlen; i++) {
 311     if (list[i] == (char)NEWLINE) {
 312       s[l] = '\0'; l = 0;
 313       puts(s);
 314       }
 315     else if (l < 255)
 316       s[l++] = list[i];
 317     }
 318   putchar('\n');
 319 }
 320
 321
 322
 323 void usage()
 324 {
 325   short   i, midi;
 326
 327   fprintf(stderr,title);
 328   fprintf(stderr,
 329   "usage: readseq [-options] in.seq > out.seq\n");
 330   fprintf(stderr," options\n");
 331 /* ? add -d[igits] to allow digits in sequence data, &/or option to specify seq charset !? */
 332   fprintf(stderr, "    -a[ll]         select All sequences\n");
 333   fprintf(stderr, "    -c[aselower]   change to lower case\n");
 334   fprintf(stderr, "    -C[ASEUPPER]   change to UPPER CASE\n");
 335   fprintf(stderr, "    -degap[=-]     remove gap symbols\n");
 336   fprintf(stderr, "    -i[tem=2,3,4]  select Item number(s) from several\n");
 337   fprintf(stderr, "    -l[ist]        List sequences only\n");
 338   fprintf(stderr, "    -o[utput=]out.seq  redirect Output\n");
 339   fprintf(stderr, "    -p[ipe]        Pipe (command line, <stdin, >stdout)\n");
 340   fprintf(stderr, "    -r[everse]     change to Reverse-complement\n");
 341   fprintf(stderr, "    -v[erbose]     Verbose progress\n");
 342   fprintf(stderr, "    -f[ormat=]#    Format number for output,  or\n");
 343   fprintf(stderr, "    -f[ormat=]Name Format name for output:\n");
 344   midi = (kMaxFormat+1) / 2;
 345   for (i = kMinFormat-1; i < midi; i++)
 346    fprintf( stderr, "        %-20s      %-20s\n",
 347     formats[i], formats[midi+i]);
 348
 349   /* new output format options, esp. for pretty format: */
 350   fprintf(stderr, "     \n");
 351   fprintf(stderr, "   Pretty format options: \n");
 352   fprintf(stderr, "    -wid[th]=#            sequence line width\n");
 353   fprintf(stderr, "    -tab=#                left indent\n");
 354   fprintf(stderr, "    -col[space]=#         column space within sequence line on output\n");
 355   fprintf(stderr, "    -gap[count]           count gap chars in sequence numbers\n");
 356   fprintf(stderr, "    -nameleft, -nameright[=#]   name on left/right side [=max width]\n");
 357   fprintf(stderr, "    -nametop              name at top/bottom\n");
 358   fprintf(stderr, "    -numleft, -numright   seq index on left/right side\n");
 359   fprintf(stderr, "    -numtop, -numbot      index on top/bottom\n");
 360   fprintf(stderr, "    -match[=.]            use match base for 2..n species\n");
 361   fprintf(stderr, "    -inter[line=#]        blank line(s) between sequence blocks\n");
 362
 363   /******  not ready yet
 364   fprintf(stderr, "    -code=none,rtf,postscript,ps   code syntax\n");
 365   fprintf(stderr, "    -namefont=, -numfont=, -seqfont=font   font choice\n");
 366   fprintf(stderr, "       font suggestions include times,courier,helvetica\n");
 367   fprintf(stderr, "    -namefontsize=, -numfontsize=, -seqfontsize=#\n");
 368   fprintf(stderr, "       fontsize suggestions include 9,10,12,14\n");
 369   fprintf(stderr, "    -namefontstyle=, -numfontstyle=, -seqfontstyle= style  fontstyle for names\n");
 370   fprintf(stderr, "       fontstyle options are plain,italic,bold,bold-italic\n");
 371   ******/
 372 }
 373
 374 void erralert(short err)
 375 {
 376   switch (err) {
 377     case 0  :
 378       break;
 379     case eFileNotFound: fprintf(stderr, "File not found: %s\n", inputfile);
 380       break;
 381     case eFileCreate: fprintf(stderr, "Can't open output file.\n");
 382       break;
 383     case eASNerr: fprintf(stderr, "Error in ASN.1 sequence routines.\n");
 384       break;
 385     case eNoData: fprintf(stderr, "No data in file.\n");
 386       break;
 387     case eItemNotFound: fprintf(stderr, "Specified item not in file.\n");
 388       break;
 389     case eUnequalSize:  fprintf(stderr,
 390       "This format requires equal length sequences.\nSequence truncated or padded to fit.\n");
 391       break;
 392     case eUnknownFormat: fprintf(stderr, "Error: this format is unknown to me.\n");
 393       break;
 394     case eOneFormat: fprintf(stderr,
 395       "Warning: This format permits only 1 sequence per file.\n");
 396       break;
 397     case eMemFull: fprintf(stderr, "Out of storage memory. Sequence truncated.\n");
 398       break;
 399     default: fprintf(stderr, "readSeq error = %d\n", err);
 400       break;
 401     }
 402 } /* erralert */
 403
 404
 405 int chooseFormat( boolean quietly)
 406 {
 407   char  sform[128];
 408   int   midi, i, outform;
 409
 410     if (quietly)
 411       return kPearson;  /* default */
 412     else {
 413       midi = (kMaxFormat+1) / 2;
 414       for (i = kMinFormat-1; i < midi; i++)
 415         fprintf( stderr, "        %-20s      %-20s\n",
 416                         formats[i], formats[midi+i]);
 417       fprintf(stderr,"\nChoose an output format (name or #): \n");
 418       gets(sform);
 419       outform = parseformat(sform);
 420       if (outform == kNoformat) outform = kPearson;
 421       return outform;
 422       }
 423 }
 424
 425
 426
 427 /* read paramater(s) */
 428
 429 boolean checkopt( boolean casesense, char *sopt, const char *smatch, short minword)
 430 {
 431   long  lenopt, lenmatch;
 432   boolean result;
 433   short minmaxw;
 434
 435   lenopt = strlen(sopt);
 436   lenmatch= strlen(smatch);
 437   minmaxw= max(minword, min(lenopt, lenmatch));
 438
 439   if (casesense)
 440     result= (!strncmp( sopt, smatch, minmaxw));
 441   else
 442     result= (!Strncasecmp( sopt, smatch, minmaxw ));
 443   /* if (result) { */
 444     /* fprintf(stderr,"true checkopt(opt=%s,match=%s,param=%s)\n", sopt, smatch, *sparam); */
 445   /*  } */
 446   return result;
 447 }
 448
 449
 450 #define   kMaxwhichlist  50
 451
 452 /* global for readopt(), main() */
 453 boolean   chooseall = false, quietly = false, gotinputfile = false,
 454           listonly = false, closeout = false, verbose = false,
 455           manyout = false, dolower = false, doupper = false, doreverse= false,
 456           askout  = true, dopipe= false, interleaved = false;
 457 short     nfile = 0, iwhichlist=0, nwhichlist = 0;
 458 short     whichlist[kMaxwhichlist+1];
 459 long      whichSeq = 0, outform = kNoformat;
 460 char      onamestore[128], *oname = onamestore;
 461 FILE      *foo = NULL;
 462
 463 void resetGlobals()
 464 /* need this when used from SIOW, as these globals are not reinited automatically
 465 between calls to local main() */
 466 {
 467   chooseall = false; quietly = false; gotinputfile = false;
 468   listonly = false; closeout = false; verbose = false;
 469   manyout = false; dolower = false; doupper = false; doreverse= false;
 470   askout  = true; dopipe= false; interleaved = false;
 471   nfile = 0; iwhichlist=0; nwhichlist = 0;
 472   whichSeq = 0; outform = kNoformat;
 473   oname = onamestore;
 474   foo = NULL;
 475
 476   gPrettyInit(gPretty);
 477 }
 478
 479
 480 #define kOptOkay  1
 481 #define kOptNone  0
 482
 483 int readopt( char *sopt)
 484 {
 485   char    sparamstore[256], *sparam= sparamstore;
 486   short   n, slen= strlen(sopt);
 487
 488   /* fprintf(stderr,"readopt( %s) == ", sopt); */
 489
 490   if (*sopt == '?') {
 491     usage();
 492     return kOptNone;   /*? eOptionBad or kOptNone */
 493     }
 494
 495   else if (*sopt == '-') {
 496
 497     char *cp= strchr(sopt,'=');
 498     *sparam= '\0';
 499     if (cp) {
 500       strcpy(sparam, cp+1);
 501       *cp= 0;
 502       }
 503
 504     if (checkopt( false, sopt, "-help", 2)) {
 505       usage();
 506       return kOptNone;
 507       }
 508
 509     if (checkopt( false, sopt, "-all", 2)) {
 510       whichSeq= 1; chooseall= true;
 511       return kOptOkay;
 512       }
 513
 514     if (checkopt( false, sopt, "-colspace", 4)) { /* test before -c[ase] */
 515       n= atoi( sparam);
 516       gPretty.spacer = n;
 517       return kOptOkay;
 518       }
 519
 520     if (checkopt( true, sopt, "-caselower", 2)) {
 521       dolower= true;
 522       return kOptOkay;
 523       }
 524     if (checkopt( true, sopt, "-CASEUPPER", 2)) {
 525       doupper= true;
 526       return kOptOkay;
 527       }
 528
 529     if (checkopt( false, sopt, "-pipe", 2)) {
 530       dopipe= true; askout= false;
 531       return kOptOkay;
 532       }
 533
 534     if (checkopt( false, sopt, "-list", 2)) {
 535       listonly = true; askout = false;
 536       return kOptOkay;
 537       }
 538
 539     if (checkopt( false, sopt, "-reverse", 2)) {
 540       doreverse = true;
 541       return kOptOkay;
 542       }
 543
 544     if (checkopt( false, sopt, "-verbose", 2)) {
 545       verbose = true;
 546       return kOptOkay;
 547       }
 548
 549     if (checkopt( false, sopt, "-match", 5)) {
 550       gPretty.domatch= true;
 551       if (*sparam >= ' ') gPretty.matchchar= *sparam;
 552       return kOptOkay;
 553       }
 554     if (checkopt( false, sopt, "-degap", 4)) {
 555       gPretty.degap= true;
 556       if (*sparam >= ' ') gPretty.gapchar= *sparam;
 557       return kOptOkay;
 558       }
 559
 560     if (checkopt( false, sopt, "-interline", 4)) {
 561       gPretty.interline= atoi( sparam);
 562       return kOptOkay;
 563       }
 564
 565     if (checkopt( false, sopt, "-item", 2)) {
 566       char  *cp = sparam;
 567       nwhichlist= 0;
 568       whichlist[0]= 0;
 569       if (*cp == 0) cp= sopt+2; /* compatible w/ old way */
 570       do {
 571         while (*cp!=0 && !isdigit(*cp)) cp++;
 572         if (*cp!=0) {
 573           n = atoi( cp);
 574           whichlist[nwhichlist++]= n;
 575           while (*cp!=0 && isdigit(*cp)) cp++;
 576           }
 577       } while (*cp!=0 && n>0 && nwhichlist<kMaxwhichlist);
 578       whichlist[nwhichlist++]= 0; /* 0 == stopsign for loop */
 579       whichSeq= max(1,whichlist[0]); iwhichlist= 1;
 580       return kOptOkay;
 581       }
 582
 583     if (checkopt( false, sopt, "-format", 5)) {/* -format=phylip, -f2, -form=phylip */
 584       if (*sparam==0) { for (sparam= sopt+2; isalpha(*sparam); sparam++) ; }
 585       outform = parseformat( sparam);
 586       return kOptOkay;
 587       }
 588     if (checkopt( false, sopt, "-f", 2)) { /* compatible w/ -fphylip prior version */
 589       if (*sparam==0) sparam= sopt+2;
 590       outform = parseformat( sparam);
 591       return kOptOkay;
 592       }
 593
 594     if (checkopt( false, sopt, "-output", 3)) {/* -output=myseq */
 595       if (*sparam==0) { for (sparam= sopt+3; isalpha(*sparam); sparam++) ; }
 596       strcpy( oname, sparam);
 597       foo = fopen( oname, "w");
 598       if (!foo) { erralert(eFileCreate); return eFileCreate; }
 599       closeout = true;
 600       askout = false;
 601       return kOptOkay;
 602       }
 603     if (checkopt( false, sopt, "-o", 2)) {  /* compatible w/ -omyseq prior version */
 604       if (*sparam==0) sparam= sopt+2;
 605       strcpy( oname, sparam);
 606       foo = fopen( oname, "w");
 607       if (!foo) { erralert(eFileCreate); return eFileCreate; }
 608       closeout = true;
 609       askout = false;
 610       return kOptOkay;
 611       }
 612
 613     if (checkopt( false, sopt, "-width", 2)) {
 614       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 615       n= atoi( sparam);
 616       if (n>0) gPretty.seqwidth = n;
 617       return kOptOkay;
 618       }
 619
 620     if (checkopt( false, sopt, "-tab", 4)) {
 621       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 622       n= atoi( sparam);
 623       gPretty.tab = n;
 624       return kOptOkay;
 625       }
 626
 627     if (checkopt( false, sopt, "-gapcount", 4)) {
 628       gPretty.baseonlynum = false;
 629       /* if (*sparam >= ' ') gPretty.gapchar= *sparam; */
 630       return kOptOkay;
 631       }
 632     if (checkopt( false, sopt, "-nointerleave", 8)) {
 633       gPretty.noleaves = true;
 634       return kOptOkay;
 635       }
 636
 637     if (checkopt( false, sopt, "-nameleft", 7)) {
 638       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 639       n= atoi( sparam);
 640       if (n>0 && n<50) gPretty.namewidth =  n;
 641       gPretty.nameleft= true;
 642       return kOptOkay;
 643       }
 644     if (checkopt( false, sopt, "-nameright", 7)) {
 645       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 646       n= atoi( sparam);
 647       if (n>0 && n<50) gPretty.namewidth =  n;
 648       gPretty.nameright= true;
 649       return kOptOkay;
 650       }
 651     if (checkopt( false, sopt, "-nametop", 6)) {
 652       gPretty.nametop= true;
 653       return kOptOkay;
 654       }
 655
 656     if (checkopt( false, sopt, "-numleft", 6)) {
 657       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 658       n= atoi( sparam);
 659       if (n>0 && n<50) gPretty.numwidth =  n;
 660       gPretty.numleft= true;
 661       return kOptOkay;
 662       }
 663     if (checkopt( false, sopt, "-numright", 6)) {
 664       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
 665       n= atoi( sparam);
 666       if (n>0 && n<50) gPretty.numwidth =  n;
 667       gPretty.numright= true;
 668       return kOptOkay;
 669       }
 670
 671     if (checkopt( false, sopt, "-numtop", 6)) {
 672       gPretty.numtop= true;
 673       return kOptOkay;
 674       }
 675     if (checkopt( false, sopt, "-numbottom", 6)) {
 676       gPretty.numbot= true;
 677       return kOptOkay;
 678       }
 679
 680     else {
 681       usage();
 682       return eOptionBad;
 683       }
 684     }
 685
 686   else {
 687     strcpy( inputfile, sopt);
 688     gotinputfile = (*inputfile != 0);
 689     nfile++;
 690     return kOptOkay;
 691     }
 692
 693  /* return kOptNone; -- never here */
 694 }
 695
 696
 697
 698
 699 /* this program suffers some as it tries to be a quiet translator pipe
 700    _and_ a noisy user interactor
 701 */
 702
 703 /* return is best for SIOW, okay for others */
 704 #ifdef SIOW
 705 #define Exit(a)   return(a)
 706 siow_main( int argc, char *argv[])
 707
 708 #else
 709 #define Exit(a)   exit(a)
 710
 711 main( int argc, char *argv[])
 712 #endif
 713 {
 714 boolean   closein = false;
 715 short     ifile, nseq, atseq, format, err = 0, seqtype = kDNA,
 716           nlines, seqout = 0, phylvers = 2;
 717 long      i, skiplines, seqlen, seqlen0;
 718 unsigned long  checksum= 0, checkall= 0;
 719 char      *seq, *cp, *firstseq = NULL, *seqlist, *progname, tempname[256];
 720 char      seqid[256], *seqidptr = seqid;
 721 char      stempstore[256], *stemp = stempstore;
 722 FILE      *ftmp, *fin, *fout;
 723 long      outindexmax= 0, noutindex= 0, *outindex = NULL;
 724
 725 #define exit_main(err) {        \
 726   if (closeout) fclose(fout);   \
 727   if (closein) fclose(fin);   \
 728   if (*tempname!=0) remove(tempname);\
 729   Exit(err); }
 730
 731 #define indexout()  if (interleaved) {\
 732   if (noutindex>=outindexmax) {\
 733     outindexmax= noutindex + 20;\
 734     outindex= (long*) realloc(outindex, sizeof(long)*outindexmax);\
 735     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }\
 736     }\
 737   outindex[noutindex++]= ftell(fout);\
 738   }
 739
 740
 741   resetGlobals();
 742   foo = stdout;
 743   progname = argv[0];
 744   *oname = 0;
 745   *tempname = 0;
 746   /* initialize gPretty ?? -- done in header */
 747
 748   for (i=1; i < argc; i++) {
 749     err= readopt( argv[i]);
 750     if (err <= 0) exit_main(err);
 751     }
 752
 753                             /* pipe input from stdin !? */
 754   if (dopipe && !gotinputfile) {
 755     int c;
 756     tmpnam(tempname);
 757     inputfile = tempname;
 758     ftmp = fopen( inputfile, "w");
 759     if (!ftmp) { erralert(eFileCreate); exit_main(eFileCreate); }
 760     while ((c = getc(stdin)) != EOF) fputc(c, ftmp);
 761     fclose(ftmp);
 762     gotinputfile= true;
 763     }
 764
 765   quietly = (dopipe || (gotinputfile && (listonly || whichSeq != 0)));
 766
 767   if (verbose || (!quietly && !gotinputfile)) fprintf( stderr, title);
 768   ifile = 1;
 769
 770                             /* UI: Choose output */
 771   if (askout && !closeout && !quietly) {
 772     askout = false;
 773     fprintf(stderr,"\nName of output file (?=help, defaults to display): \n");
 774     gets(oname= onamestore);
 775     skipwhitespace(oname);
 776     if (*oname == '?') { usage(); exit_main(0); }
 777     else if (*oname != 0) {
 778       closeout = true;
 779       foo = fopen( oname, "w");
 780       if (!foo) { erralert(eFileCreate); exit_main(eFileCreate); }
 781       }
 782     }
 783
 784   fout = foo;
 785   if (outform == kNoformat) outform = chooseFormat(quietly);
 786
 787                           /* set up formats ... */
 788   switch (outform) {
 789     case kPhylip2:
 790       interleaved= false;
 791       phylvers = 2;
 792       outform = kPhylip;
 793       break;
 794
 795     case kPhylip4:
 796       interleaved= true;
 797       phylvers = 4;
 798       outform = kPhylip;
 799       break;
 800
 801     case kMSF:
 802     case kPAUP:
 803       interleaved= true;
 804       break;
 805
 806     case kPretty:
 807       gPretty.isactive= true;
 808       interleaved= true;
 809       break;
 810
 811     }
 812
 813   if (gPretty.isactive && gPretty.noleaves) interleaved= false;
 814   if (interleaved) {
 815     fout = ftmp = tmpfile();
 816     outindexmax= 30; noutindex= 0;
 817     outindex = (long*) malloc(outindexmax*sizeof(long));
 818     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }
 819     }
 820
 821                         /* big loop over all input files */
 822   do {
 823                         /* select next input file */
 824     gotinputfile = (*tempname != 0);
 825     while ((ifile < argc) && (!gotinputfile)) {
 826       if (*argv[ifile] != '-') {
 827         strcpy( inputfile, argv[ifile]);
 828         gotinputfile = (*inputfile != 0);
 829         --nfile;
 830         }
 831       ifile++;
 832       }
 833
 834     while (!gotinputfile) {
 835       fprintf(stderr,"\nName an input sequence or -option: \n");
 836       inputfile= inputfilestore;
 837
 838       gets(stemp= stempstore);
 839       if (*stemp==0) goto fini;  /* !! need this to finish work during interactive use */
 840       stemp= strtok(stempstore, " \n\r\t");
 841       while (stemp) {
 842         err= readopt( stemp); /* will read inputfile if it exists */
 843         if (err<0) exit_main(err);
 844         stemp= strtok( NULL, " \n\r\t");
 845         }
 846       }
 847               /* thanks to AJB@UK.AC.DARESBURY.DLVH for this PHYLIP3 fix: */
 848               /* head for end (interleave if needed) */
 849     if (*inputfile == 0) break;
 850
 851     format = seqFileFormat( inputfile, &skiplines, &err);
 852
 853     if (err == 0)  {
 854 #ifdef NCBI
 855       if (format == kASNseqentry || format == kASNseqset)
 856         seqlist = listASNSeqs( inputfile, skiplines, format, &nseq, &err);
 857       else
 858 #endif
 859         seqlist = listSeqs( inputfile, skiplines, format, &nseq, &err);
 860       }
 861
 862     if (err != 0)
 863       erralert(err);
 864
 865     else if (listonly) {
 866       dumpSeqList(seqlist,format);
 867       free( seqlist);
 868       }
 869
 870     else {
 871                                 /* choose whichSeq if needed */
 872       if (nseq == 1 || chooseall || (quietly && whichSeq == 0)) {
 873         chooseall= true;
 874         whichSeq = 1;
 875         quietly = true; /* no loop */
 876         }
 877       else if (whichSeq > nseq && quietly) {
 878         erralert(eItemNotFound);
 879         err= eItemNotFound;
 880         }
 881       else if (whichSeq > nseq || !quietly) {
 882         dumpSeqList(seqlist, format);
 883         fprintf(stderr,"\nChoose a sequence (# or All): \n");
 884         gets(stemp= stempstore);
 885         skipwhitespace(stemp);
 886         if (to_lower(*stemp) == 'a') {
 887           chooseall= true;
 888           whichSeq = 1;
 889           quietly = true; /* !? this means we don't ask for another file
 890                             as well as no more whichSeqs... */
 891           }
 892         else if (isdigit(*stemp)) whichSeq= atol(stemp);
 893         else whichSeq= 1; /* default */
 894         }
 895       free( seqlist);
 896
 897       if (false /*chooseall*/) {  /* this isn't debugged yet...*/
 898         fin = fopen(inputfile, "r");
 899         closein= true;
 900         }
 901
 902       while (whichSeq > 0 && whichSeq <= nseq) {
 903                                 /* need to open multiple output files ? */
 904         manyout = ((chooseall || nwhichlist>1) && nseq > 1
 905                   && (outform == kPlain || outform == kGCG));
 906         if (manyout) {
 907           if ( whichSeq == 1 ) erralert(eOneFormat);
 908           else if (closeout) {
 909             sprintf( stemp,"%s_%d", oname, whichSeq);
 910             freopen( stemp, "w", fout);
 911             fprintf( stderr,"Writing sequence %d to file %s\n", whichSeq, stemp);
 912             }
 913           }
 914
 915         if (closein) {
 916           /* !! this fails... skips most seqs... */
 917           /* !! in sequential read, must count seqs already read from whichSeq ... */
 918           /* need major revision of ureadseq before we can do this */
 919           atseq= whichSeq-1;
 920           seqidptr= seqid;
 921           seq = readSeqFp( whichSeq, fin, skiplines, format,
 922                           &seqlen, &atseq, &err, seqidptr);
 923           skiplines= 0;
 924           }
 925         else {
 926           atseq= 0;
 927           seqidptr= seqid;
 928 #ifdef NCBI
 929           if (format == kASNseqentry || format == kASNseqset) {
 930             seqidptr= NULL;
 931             seq = readASNSeq( whichSeq, inputfile, skiplines, format,
 932                      &seqlen, &atseq, &err, &seqidptr);
 933             }
 934           else
 935 #endif
 936           seq = readSeq( whichSeq, inputfile, skiplines, format,
 937                           &seqlen, &atseq, &err, seqidptr);
 938           }
 939
 940
 941         if (gPretty.degap) {
 942           char *newseq;
 943           long newlen;
 944           newseq= compressSeq( gPretty.gapchar, seq, seqlen, &newlen);
 945           if (newseq) {
 946             free(seq); seq= newseq; seqlen= newlen;
 947             }
 948           }
 949
 950         if (outform == kMSF) checksum= GCGchecksum(seq, seqlen, &checkall);
 951         else if (verbose) checksum= seqchecksum(seq, seqlen, &checkall);
 952         if (verbose)
 953           fprintf( stderr, "Sequence %d, length= %d, checksum= %X, format= %s, id= %s\n",
 954                 whichSeq, seqlen, checksum, formatstr(format), seqidptr);
 955
 956         if (err != 0) erralert(err);
 957         else {
 958                                   /* format fixes that writeseq doesn't do */
 959           switch (outform) {
 960             case kPIR:
 961               if (seqout == 0) fprintf( foo,"\\\\\\\n");
 962               break;
 963             case kASN1:
 964               if (seqout == 0) fprintf( foo, kASN1headline);
 965               break;
 966
 967             case kPhylip:
 968               if (seqout == 0) {
 969                 if (!interleaved) {  /*  bug, nseq is for 1st infile only */
 970                   if (chooseall) i= nseq; else i=1;
 971                   if (phylvers >= 4) fprintf(foo," %d %d\n", i, seqlen);
 972                   else fprintf(foo," %d %d YF\n", i, seqlen);
 973                   }
 974                 seqlen0 = seqlen;
 975                 }
 976               else if (seqlen != seqlen0) {
 977                 erralert(eUnequalSize);
 978                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
 979                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
 980                 seqlen = seqlen0;
 981                 seq[seqlen] = 0;
 982                 }
 983               break;
 984
 985             case kPAUP:
 986               if (seqout == 0) {
 987                 seqtype= getseqtype(seq, seqlen);
 988                 seqlen0 = seqlen;
 989                 }
 990               else if (seqlen != seqlen0) {
 991                 erralert(eUnequalSize);
 992                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
 993                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
 994                 seqlen = seqlen0;
 995                 seq[seqlen] = 0;
 996                 }
 997               break;
 998
 999             }
1000
1001           if (doupper)
1002             for (i = 0; i<seqlen; i++) seq[i] = to_upper(seq[i]);
1003           else if (dolower)
1004             for (i = 0; i<seqlen; i++) seq[i] = to_lower(seq[i]);
1005
1006           if (doreverse) {
1007             long  j, k;
1008             char  ctemp;
1009             for (j=0, k=seqlen-1; j <= k; j++, k--) {
1010               ctemp = compl[seq[j] - ' '];
1011               seq[j] = compl[seq[k] - ' '];
1012               seq[k] = ctemp;
1013               }
1014             }
1015
1016           if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq != NULL) {
1017             for (i=0; i<seqlen; i++)
1018               if (seq[i]==firstseq[i]) seq[i]= gPretty.matchchar;
1019             }
1020
1021
1022           if (gPretty.isactive && gPretty.numtop && seqout == 0) {
1023             gPretty.numline = 1;
1024             indexout();
1025             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1026             gPretty.numline = 2;
1027             indexout();
1028             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1029             gPretty.numline = 0;
1030             }
1031
1032           indexout();
1033           nlines = writeSeq( fout, seq, seqlen, outform, seqidptr);
1034           seqout++;
1035           }
1036
1037         if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq == NULL) {
1038           firstseq= seq;
1039           seq = NULL;
1040           }
1041         else if (seq!=NULL) { free(seq); seq = NULL; }
1042
1043 #ifdef NCBI
1044        if ( (format == kASNseqentry || format == kASNseqset)
1045           && seqidptr && seqidptr!= seqid)
1046             free(seqidptr);
1047 #endif
1048         if (chooseall) whichSeq++;
1049         else if (iwhichlist<nwhichlist) whichSeq= whichlist[iwhichlist++];
1050         else whichSeq= 0;
1051         }
1052       if (closein) { fclose(fin); closein= false; }
1053       }
1054     whichSeq  = 0;
1055   } while (nfile > 0 || !quietly);
1056
1057
1058 fini:
1059   if (firstseq) { free(firstseq); firstseq= NULL; }
1060   if (err || listonly) exit_main(err);
1061
1062   if (gPretty.isactive && gPretty.numbot) {
1063     gPretty.numline = 2;
1064     indexout();
1065     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1066     gPretty.numline = 1;
1067     indexout();
1068     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1069     gPretty.numline = 0;
1070     }
1071
1072   if (outform == kMSF) {
1073     if (*oname) cp= oname; else cp= inputfile;
1074     fprintf(foo,"\n %s  MSF: %d  Type: N  January 01, 1776  12:00  Check: %d ..\n\n",
1075                   cp, seqlen, checkall);
1076     }
1077
1078   if (outform == kPAUP) {
1079     fprintf(foo,"#NEXUS\n");
1080     if (*oname) cp= oname; else cp= inputfile;
1081     fprintf(foo,"[%s -- data title]\n\n", cp);
1082     /* ! now have header lines for each sequence... put them before "begin data;... */
1083     }
1084
1085   if (outform==kPhylip && interleaved) {
1086     if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
1087     else fprintf(foo," %d %d YF\n", seqout, seqlen);
1088     }
1089
1090   if (interleaved) {
1091     /* interleave species lines in true output */
1092     /* nlines is # lines / sequence */
1093     short iline, j, leaf, iseq;
1094     char  *s = stempstore;
1095
1096     indexout();  noutindex--; /* mark eof */
1097
1098     for (leaf=0; leaf<nlines; leaf++) {
1099       if (outform == kMSF && leaf == 1) {
1100         fputs("//\n\n", foo);
1101         }
1102       if (outform == kPAUP && leaf==1) {
1103         switch (seqtype) {
1104           case kDNA     : cp= "dna"; break;
1105           case kRNA     : cp= "rna"; break;
1106           case kNucleic : cp= "dna"; break;
1107           case kAmino   : cp= "protein"; break;
1108           case kOtherSeq: cp= "dna"; break;
1109           }
1110         fprintf(foo,"\nbegin data;\n");
1111         fprintf(foo," dimensions ntax=%d nchar=%d;\n", seqout, seqlen);
1112         fprintf(foo," format datatype=%s interleave missing=%c", cp, gPretty.gapchar);
1113         if (gPretty.domatch) fprintf(foo," matchchar=%c", gPretty.matchchar);
1114         fprintf(foo,";\n  matrix\n");
1115         }
1116
1117       for (iseq=0; iseq<noutindex; iseq++) {
1118         fseek(ftmp, outindex[iseq], 0);
1119         for (iline=0; iline<=leaf; iline++)
1120           if (!fgets(s, 256, ftmp)) *s= 0;
1121         if (ftell(ftmp) <= outindex[iseq+1])
1122           fputs( s, foo);
1123         }
1124
1125       for (j=0; j<gPretty.interline; j++)
1126         fputs( "\n", foo);  /* some want spacer line */
1127       }
1128     fclose(ftmp); /* tmp disappears */
1129     fout= foo;
1130     }
1131
1132   if (outform == kASN1)  fprintf( foo, "} }\n");
1133   if (outform == kPAUP)  fprintf( foo,";\n  end;\n");
1134
1135   if (outindex != NULL) free(outindex);
1136   exit_main(0);
1137 }
1138
1139