source: branches/stable/READSEQ/readseq.c

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