source: tags/arb_5.2/READSEQ/readseq.c

Last change on this file was 5390, checked in by westram, 16 years ago
  • TAB-Ex
  • 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, "File not found: %s\n", inputfile);
385      break;
386    case eFileCreate: fprintf(stderr, "Can't open output file.\n");
387      break;
388    case eASNerr: fprintf(stderr, "Error in ASN.1 sequence routines.\n");
389      break;
390    case eNoData: fprintf(stderr, "No data in file.\n");
391      break;
392    case eItemNotFound: fprintf(stderr, "Specified item not in file.\n");
393      break;
394    case eUnequalSize:  fprintf(stderr,
395      "This format requires equal length sequences.\nSequence truncated or padded to fit.\n");
396      break;
397    case eUnknownFormat: fprintf(stderr, "Error: this format is unknown to me.\n");
398      break;
399    case eOneFormat: fprintf(stderr,
400      "Warning: This format permits only 1 sequence per file.\n");
401      break;
402    case eMemFull: fprintf(stderr, "Out of storage memory. Sequence truncated.\n");
403      break;
404    default: fprintf(stderr, "readSeq error = %d\n", err);
405      break;
406    }
407} /* erralert */
408
409
410int chooseFormat( boolean quietly)
411{
412  char  sform[128];
413  int   midi, i, outform;
414
415    if (quietly)
416      return kPearson;  /* default */
417    else {
418      midi = (kMaxFormat+1) / 2;
419      for (i = kMinFormat-1; i < midi; i++)
420        fprintf( stderr, "        %-20s      %-20s\n",
421                        formats[i], formats[midi+i]);
422      fprintf(stderr,"\nChoose an output format (name or #): \n");
423      gets(sform);
424      outform = parseformat(sform);
425      if (outform == kNoformat) outform = kPearson;
426      return outform;
427      }
428}
429
430
431
432/* read paramater(s) */
433
434boolean checkopt( boolean casesense, char *sopt, const char *smatch, short minword)
435{
436  long  lenopt, lenmatch;
437  boolean result;
438  short minmaxw;
439
440  lenopt = strlen(sopt);
441  lenmatch= strlen(smatch);
442  minmaxw= max(minword, min(lenopt, lenmatch));
443
444  if (casesense)
445    result= (!strncmp( sopt, smatch, minmaxw));
446  else
447    result= (!Strncasecmp( sopt, smatch, minmaxw ));
448  /* if (result) { */
449    /* fprintf(stderr,"true checkopt(opt=%s,match=%s,param=%s)\n", sopt, smatch, *sparam); */
450  /*  } */
451  return result;
452}
453
454
455#define   kMaxwhichlist  50
456
457/* global for readopt(), main() */
458boolean   chooseall = false, quietly = false, gotinputfile = false,
459          listonly = false, closeout = false, verbose = false,
460          manyout = false, dolower = false, doupper = false, doreverse= false,
461          askout  = true, dopipe= false, interleaved = false;
462short     nfile = 0, iwhichlist=0, nwhichlist = 0;
463short     whichlist[kMaxwhichlist+1];
464long      whichSeq = 0, outform = kNoformat;
465char      onamestore[128], *oname = onamestore;
466FILE      *foo = NULL;
467
468void resetGlobals()
469/* need this when used from SIOW, as these globals are not reinited automatically
470between calls to local main() */
471{
472  chooseall = false; quietly = false; gotinputfile = false;
473  listonly = false; closeout = false; verbose = false;
474  manyout = false; dolower = false; doupper = false; doreverse= false;
475  askout  = true; dopipe= false; interleaved = false;
476  nfile = 0; iwhichlist=0; nwhichlist = 0;
477  whichSeq = 0; outform = kNoformat;
478  oname = onamestore;
479  foo = NULL;
480
481  gPrettyInit(gPretty);
482}
483
484
485#define kOptOkay  1
486#define kOptNone  0
487
488int readopt( char *sopt)
489{
490  char    sparamstore[256], *sparam= sparamstore;
491  short   n, slen= strlen(sopt);
492
493  /* fprintf(stderr,"readopt( %s) == ", sopt); */
494
495  if (*sopt == '?') {
496    usage();
497    return kOptNone;   /*? eOptionBad or kOptNone */
498    }
499
500  else if (*sopt == '-') {
501
502    char *cp= strchr(sopt,'=');
503    *sparam= '\0';
504    if (cp) {
505      strcpy(sparam, cp+1);
506      *cp= 0;
507      }
508
509    if (checkopt( false, sopt, "-help", 2)) {
510      usage();
511      return kOptNone;
512      }
513
514    if (checkopt( false, sopt, "-all", 2)) {
515      whichSeq= 1; chooseall= true;
516      return kOptOkay;
517      }
518
519    if (checkopt( false, sopt, "-colspace", 4)) { /* test before -c[ase] */
520      n= atoi( sparam);
521      gPretty.spacer = n;
522      return kOptOkay;
523      }
524
525    if (checkopt( true, sopt, "-caselower", 2)) {
526      dolower= true;
527      return kOptOkay;
528      }
529    if (checkopt( true, sopt, "-CASEUPPER", 2)) {
530      doupper= true;
531      return kOptOkay;
532      }
533
534    if (checkopt( false, sopt, "-pipe", 2)) {
535      dopipe= true; askout= false;
536      return kOptOkay;
537      }
538
539    if (checkopt( false, sopt, "-list", 2)) {
540      listonly = true; askout = false;
541      return kOptOkay;
542      }
543
544    if (checkopt( false, sopt, "-reverse", 2)) {
545      doreverse = true;
546      return kOptOkay;
547      }
548
549    if (checkopt( false, sopt, "-verbose", 2)) {
550      verbose = true;
551      return kOptOkay;
552      }
553
554    if (checkopt( false, sopt, "-match", 5)) {
555      gPretty.domatch= true;
556      if (*sparam >= ' ') gPretty.matchchar= *sparam;
557      return kOptOkay;
558      }
559    if (checkopt( false, sopt, "-degap", 4)) {
560      gPretty.degap= true;
561      if (*sparam >= ' ') gPretty.gapchar= *sparam;
562      return kOptOkay;
563      }
564
565    if (checkopt( false, sopt, "-interline", 4)) {
566      gPretty.interline= atoi( sparam);
567      return kOptOkay;
568      }
569
570    if (checkopt( false, sopt, "-item", 2)) {
571      char  *cp = sparam;
572      nwhichlist= 0;
573      whichlist[0]= 0;
574      if (*cp == 0) cp= sopt+2; /* compatible w/ old way */
575      do {
576        while (*cp!=0 && !rs_isdigit(*cp)) cp++;
577        if (*cp!=0) {
578          n = atoi( cp);
579          whichlist[nwhichlist++]= n;
580          while (*cp!=0 && rs_isdigit(*cp)) cp++;
581          }
582      } while (*cp!=0 && n>0 && nwhichlist<kMaxwhichlist);
583      whichlist[nwhichlist++]= 0; /* 0 == stopsign for loop */
584      whichSeq= max(1,whichlist[0]); iwhichlist= 1;
585      return kOptOkay;
586      }
587
588    if (checkopt( false, sopt, "-format", 5)) {/* -format=phylip, -f2, -form=phylip */
589      if (*sparam==0) { for (sparam= sopt+2; isalpha(*sparam); sparam++) ; }
590      outform = parseformat( sparam);
591      return kOptOkay;
592      }
593    if (checkopt( false, sopt, "-f", 2)) { /* compatible w/ -fphylip prior version */
594      if (*sparam==0) sparam= sopt+2;
595      outform = parseformat( sparam);
596      return kOptOkay;
597      }
598
599    if (checkopt( false, sopt, "-output", 3)) {/* -output=myseq */
600      if (*sparam==0) { for (sparam= sopt+3; isalpha(*sparam); sparam++) ; }
601      strcpy( oname, sparam);
602      foo = fopen( oname, "w");
603      if (!foo) { erralert(eFileCreate); return eFileCreate; }
604      closeout = true;
605      askout = false;
606      return kOptOkay;
607      }
608    if (checkopt( false, sopt, "-o", 2)) {  /* compatible w/ -omyseq prior version */
609      if (*sparam==0) sparam= sopt+2;
610      strcpy( oname, sparam);
611      foo = fopen( oname, "w");
612      if (!foo) { erralert(eFileCreate); return eFileCreate; }
613      closeout = true;
614      askout = false;
615      return kOptOkay;
616      }
617
618    if (checkopt( false, sopt, "-width", 2)) {
619      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
620      n= atoi( sparam);
621      if (n>0) gPretty.seqwidth = n;
622      return kOptOkay;
623      }
624
625    if (checkopt( false, sopt, "-tab", 4)) {
626      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
627      n= atoi( sparam);
628      gPretty.tab = n;
629      return kOptOkay;
630      }
631
632    if (checkopt( false, sopt, "-gapcount", 4)) {
633      gPretty.baseonlynum = false;
634      /* if (*sparam >= ' ') gPretty.gapchar= *sparam; */
635      return kOptOkay;
636      }
637    if (checkopt( false, sopt, "-nointerleave", 8)) {
638      gPretty.noleaves = true;
639      return kOptOkay;
640      }
641
642    if (checkopt( false, sopt, "-nameleft", 7)) {
643      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
644      n= atoi( sparam);
645      if (n>0 && n<50) gPretty.namewidth =  n;
646      gPretty.nameleft= true;
647      return kOptOkay;
648      }
649    if (checkopt( false, sopt, "-nameright", 7)) {
650      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
651      n= atoi( sparam);
652      if (n>0 && n<50) gPretty.namewidth =  n;
653      gPretty.nameright= true;
654      return kOptOkay;
655      }
656    if (checkopt( false, sopt, "-nametop", 6)) {
657      gPretty.nametop= true;
658      return kOptOkay;
659      }
660
661    if (checkopt( false, sopt, "-numleft", 6)) {
662      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
663      n= atoi( sparam);
664      if (n>0 && n<50) gPretty.numwidth =  n;
665      gPretty.numleft= true;
666      return kOptOkay;
667      }
668    if (checkopt( false, sopt, "-numright", 6)) {
669      if (*sparam==0) { for (sparam= sopt+2; !rs_isdigit(*sparam) && *sparam!=0; sparam++) ; }
670      n= atoi( sparam);
671      if (n>0 && n<50) gPretty.numwidth =  n;
672      gPretty.numright= true;
673      return kOptOkay;
674      }
675
676    if (checkopt( false, sopt, "-numtop", 6)) {
677      gPretty.numtop= true;
678      return kOptOkay;
679      }
680    if (checkopt( false, sopt, "-numbottom", 6)) {
681      gPretty.numbot= true;
682      return kOptOkay;
683      }
684
685    else {
686      usage();
687      return eOptionBad;
688      }
689    }
690
691  else {
692    strcpy( inputfile, sopt);
693    gotinputfile = (*inputfile != 0);
694    nfile++;
695    return kOptOkay;
696    }
697
698 /* return kOptNone; -- never here */
699}
700
701
702
703
704/* this program suffers some as it tries to be a quiet translator pipe
705   _and_ a noisy user interactor
706*/
707
708/* return is best for SIOW, okay for others */
709#ifdef SIOW
710#define Exit(a)   return(a)
711siow_main( int argc, char *argv[])
712
713#else
714#define Exit(a)   exit(a)
715
716main( int argc, char *argv[])
717#endif
718{
719boolean   closein = false;
720short     ifile, nseq, atseq, format, err = 0, seqtype = kDNA,
721          nlines, seqout = 0, phylvers = 2;
722long      i, skiplines, seqlen, seqlen0;
723unsigned long  checksum= 0, checkall= 0;
724char      *seq, *cp, *firstseq = NULL, *seqlist, *progname, tempname[256];
725char      seqid[256], *seqidptr = seqid;
726char      stempstore[256], *stemp = stempstore;
727FILE      *ftmp, *fin, *fout;
728long      outindexmax= 0, noutindex= 0, *outindex = NULL;
729
730#define exit_main(err) {        \
731  if (closeout) fclose(fout);   \
732  if (closein) fclose(fin);   \
733  if (*tempname!=0) remove(tempname);\
734  Exit(err); }
735
736#define indexout()  if (interleaved) {\
737  if (noutindex>=outindexmax) {\
738    outindexmax= noutindex + 20;\
739    outindex= (long*) realloc(outindex, sizeof(long)*outindexmax);\
740    if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }\
741    }\
742  outindex[noutindex++]= ftell(fout);\
743  }
744
745
746  resetGlobals();
747  foo = stdout;
748  progname = argv[0];
749  *oname = 0;
750  *tempname = 0;
751  /* initialize gPretty ?? -- done in header */
752
753  for (i=1; i < argc; i++) {
754    err= readopt( argv[i]);
755    if (err <= 0) exit_main(err);
756    }
757
758                            /* pipe input from stdin !? */
759  if (dopipe && !gotinputfile) {
760    int c;
761    tmpnam(tempname);
762    inputfile = tempname;
763    ftmp = fopen( inputfile, "w");
764    if (!ftmp) { erralert(eFileCreate); exit_main(eFileCreate); }
765    while ((c = getc(stdin)) != EOF) fputc(c, ftmp);
766    fclose(ftmp);
767    gotinputfile= true;
768    }
769
770  quietly = (dopipe || (gotinputfile && (listonly || whichSeq != 0)));
771
772  if (verbose || (!quietly && !gotinputfile)) fprintf( stderr, title);
773  ifile = 1;
774
775                            /* UI: Choose output */
776  if (askout && !closeout && !quietly) {
777    askout = false;
778    fprintf(stderr,"\nName of output file (?=help, defaults to display): \n");
779    gets(oname= onamestore);
780    skipwhitespace(oname);
781    if (*oname == '?') { usage(); exit_main(0); }
782    else if (*oname != 0) {
783      closeout = true;
784      foo = fopen( oname, "w");
785      if (!foo) { erralert(eFileCreate); exit_main(eFileCreate); }
786      }
787    }
788
789  fout = foo;
790  if (outform == kNoformat) outform = chooseFormat(quietly);
791
792                          /* set up formats ... */
793  switch (outform) {
794    case kPhylip2:
795      interleaved= false;
796      phylvers = 2;
797      outform = kPhylip;
798      break;
799
800    case kPhylip4:
801      interleaved= true;
802      phylvers = 4;
803      outform = kPhylip;
804      break;
805
806    case kMSF:
807    case kPAUP:
808      interleaved= true;
809      break;
810
811    case kPretty:
812      gPretty.isactive= true;
813      interleaved= true;
814      break;
815
816    }
817
818  if (gPretty.isactive && gPretty.noleaves) interleaved= false;
819  if (interleaved) {
820    fout = ftmp = tmpfile();
821    outindexmax= 30; noutindex= 0;
822    outindex = (long*) malloc(outindexmax*sizeof(long));
823    if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }
824    }
825
826                        /* big loop over all input files */
827  do {
828                        /* select next input file */
829    gotinputfile = (*tempname != 0);
830    while ((ifile < argc) && (!gotinputfile)) {
831      if (*argv[ifile] != '-') {
832        strcpy( inputfile, argv[ifile]);
833        gotinputfile = (*inputfile != 0);
834        --nfile;
835        }
836      ifile++;
837      }
838
839    while (!gotinputfile) {
840      fprintf(stderr,"\nName an input sequence or -option: \n");
841      inputfile= inputfilestore;
842
843      gets(stemp= stempstore);
844      if (*stemp==0) goto fini;  /* !! need this to finish work during interactive use */
845      stemp= strtok(stempstore, " \n\r\t");
846      while (stemp) {
847        err= readopt( stemp); /* will read inputfile if it exists */
848        if (err<0) exit_main(err);
849        stemp= strtok( NULL, " \n\r\t");
850        }
851      }
852              /* thanks to AJB@UK.AC.DARESBURY.DLVH for this PHYLIP3 fix: */
853              /* head for end (interleave if needed) */
854    if (*inputfile == 0) break;
855
856    format = seqFileFormat( inputfile, &skiplines, &err);
857
858    if (err == 0)  {
859#ifdef NCBI
860      if (format == kASNseqentry || format == kASNseqset)
861        seqlist = listASNSeqs( inputfile, skiplines, format, &nseq, &err);
862      else
863#endif
864        seqlist = listSeqs( inputfile, skiplines, format, &nseq, &err);
865      }
866
867    if (err != 0)
868      erralert(err);
869
870    else if (listonly) {
871      dumpSeqList(seqlist,format);
872      free( seqlist);
873      }
874
875    else {
876                                /* choose whichSeq if needed */
877      if (nseq == 1 || chooseall || (quietly && whichSeq == 0)) {
878        chooseall= true;
879        whichSeq = 1;
880        quietly = true; /* no loop */
881        }
882      else if (whichSeq > nseq && quietly) {
883        erralert(eItemNotFound);
884        err= eItemNotFound;
885        }
886      else if (whichSeq > nseq || !quietly) {
887        dumpSeqList(seqlist, format);
888        fprintf(stderr,"\nChoose a sequence (# or All): \n");
889        gets(stemp= stempstore);
890        skipwhitespace(stemp);
891        if (to_lower(*stemp) == 'a') {
892          chooseall= true;
893          whichSeq = 1;
894          quietly = true; /* !? this means we don't ask for another file
895                            as well as no more whichSeqs... */
896          }
897        else if (rs_isdigit(*stemp)) whichSeq= atol(stemp);
898        else whichSeq= 1; /* default */
899        }
900      free( seqlist);
901
902      if (false /*chooseall*/) {  /* this isn't debugged yet...*/
903        fin = fopen(inputfile, "r");
904        closein= true;
905        }
906
907      while (whichSeq > 0 && whichSeq <= nseq) {
908                                /* need to open multiple output files ? */
909        manyout = ((chooseall || nwhichlist>1) && nseq > 1
910                  && (outform == kPlain || outform == kGCG));
911        if (manyout) {
912          if ( whichSeq == 1 ) erralert(eOneFormat);
913          else if (closeout) {
914            sprintf( stemp,"%s_%ld", oname, whichSeq);
915            freopen( stemp, "w", fout);
916            fprintf( stderr,"Writing sequence %ld to file %s\n", whichSeq, stemp);
917            }
918          }
919
920        if (closein) {
921          /* !! this fails... skips most seqs... */
922          /* !! in sequential read, must count seqs already read from whichSeq ... */
923          /* need major revision of ureadseq before we can do this */
924          atseq= whichSeq-1;
925          seqidptr= seqid;
926          seq = readSeqFp( whichSeq, fin, skiplines, format,
927                          &seqlen, &atseq, &err, seqidptr);
928          skiplines= 0;
929          }
930        else {
931          atseq= 0;
932          seqidptr= seqid;
933#ifdef NCBI
934          if (format == kASNseqentry || format == kASNseqset) {
935            seqidptr= NULL;
936            seq = readASNSeq( whichSeq, inputfile, skiplines, format,
937                     &seqlen, &atseq, &err, &seqidptr);
938            }
939          else
940#endif
941          seq = readSeq( whichSeq, inputfile, skiplines, format,
942                          &seqlen, &atseq, &err, seqidptr);
943          }
944
945
946        if (gPretty.degap) {
947          char *newseq;
948          long newlen;
949          newseq= compressSeq( gPretty.gapchar, seq, seqlen, &newlen);
950          if (newseq) {
951            free(seq); seq= newseq; seqlen= newlen;
952            }
953          }
954
955        if (outform == kMSF) checksum= GCGchecksum(seq, seqlen, &checkall);
956        else if (verbose) checksum= seqchecksum(seq, seqlen, &checkall);
957        if (verbose)
958          fprintf( stderr, "Sequence %ld, length= %ld, checksum= %lX, format= %s, id= %s\n",
959                whichSeq, seqlen, checksum, formatstr(format), seqidptr);
960
961        if (err != 0) erralert(err);
962        else {
963                                  /* format fixes that writeseq doesn't do */
964          switch (outform) {
965            case kPIR:
966              if (seqout == 0) fprintf( foo,"\\\\\\\n");
967              break;
968            case kASN1:
969              if (seqout == 0) fprintf( foo, kASN1headline);
970              break;
971
972            case kPhylip:
973              if (seqout == 0) {
974                if (!interleaved) {  /*  bug, nseq is for 1st infile only */
975                  if (chooseall) i= nseq; else i=1;
976                  if (phylvers >= 4) fprintf(foo," %ld %ld\n", i, seqlen);
977                  else fprintf(foo," %ld %ld YF\n", i, seqlen);
978                  }
979                seqlen0 = seqlen;
980                }
981              else if (seqlen != seqlen0) {
982                erralert(eUnequalSize);
983                if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
984                for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
985                seqlen = seqlen0;
986                seq[seqlen] = 0;
987                }
988              break;
989
990            case kPAUP:
991              if (seqout == 0) {
992                seqtype= getseqtype(seq, seqlen);
993                seqlen0 = seqlen;
994                }
995              else if (seqlen != seqlen0) {
996                erralert(eUnequalSize);
997                if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
998                for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
999                seqlen = seqlen0;
1000                seq[seqlen] = 0;
1001                }
1002              break;
1003
1004            }
1005
1006          if (doupper)
1007            for (i = 0; i<seqlen; i++) seq[i] = to_upper(seq[i]);
1008          else if (dolower)
1009            for (i = 0; i<seqlen; i++) seq[i] = to_lower(seq[i]);
1010
1011          if (outform==kPhylip){
1012            for (i = 0; i<seqlen; i++) if (seq[i] == '.') seq[i] = '?';
1013          }
1014
1015          if (doreverse) {
1016            long  j, k;
1017            char  ctemp;
1018            for (j=0, k=seqlen-1; j <= k; j++, k--) {
1019              ctemp = compl[seq[j] - ' '];
1020              seq[j] = compl[seq[k] - ' '];
1021              seq[k] = ctemp;
1022              }
1023            }
1024
1025          if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq != NULL) {
1026            for (i=0; i<seqlen; i++){
1027                if (seq[i] == gPretty.matchchar) seq[i] = 'o';
1028              if (seq[i]==firstseq[i]) seq[i]= gPretty.matchchar;
1029            }
1030          }
1031
1032
1033          if (gPretty.isactive && gPretty.numtop && seqout == 0) {
1034            gPretty.numline = 1;
1035            indexout();
1036            (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1037            gPretty.numline = 2;
1038            indexout();
1039            (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1040            gPretty.numline = 0;
1041            }
1042
1043          indexout();
1044          nlines = writeSeq( fout, seq, seqlen, outform, seqidptr);
1045          seqout++;
1046          }
1047
1048        if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq == NULL) {
1049          firstseq= seq;
1050          seq = NULL;
1051          }
1052        else if (seq!=NULL) { free(seq); seq = NULL; }
1053
1054#ifdef NCBI
1055       if ( (format == kASNseqentry || format == kASNseqset)
1056          && seqidptr && seqidptr!= seqid)
1057            free(seqidptr);
1058#endif
1059        if (chooseall) whichSeq++;
1060        else if (iwhichlist<nwhichlist) whichSeq= whichlist[iwhichlist++];
1061        else whichSeq= 0;
1062        }
1063      if (closein) { fclose(fin); closein= false; }
1064      }
1065    whichSeq  = 0;
1066  } while (nfile > 0 || !quietly);
1067
1068
1069fini:
1070  if (firstseq) { free(firstseq); firstseq= NULL; }
1071  if (err || listonly) exit_main(err);
1072
1073  if (gPretty.isactive && gPretty.numbot) {
1074    gPretty.numline = 2;
1075    indexout();
1076    (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1077    gPretty.numline = 1;
1078    indexout();
1079    (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
1080    gPretty.numline = 0;
1081    }
1082
1083  if (outform == kMSF) {
1084    if (*oname) cp= oname; else cp= inputfile;
1085    fprintf(foo,"\n %s  MSF: %ld  Type: N  January 01, 1776  12:00  Check: %lu ..\n\n",
1086                  cp, seqlen, checkall);
1087    }
1088
1089  if (outform == kPAUP) {
1090    fprintf(foo,"#NEXUS\n");
1091    if (*oname) cp= oname; else cp= inputfile;
1092    fprintf(foo,"[%s -- data title]\n\n", cp);
1093    /* ! now have header lines for each sequence... put them before "begin data;... */
1094    }
1095
1096  if (outform==kPhylip && interleaved) {
1097    if (phylvers >= 4) fprintf(foo," %d %ld\n", seqout, seqlen);
1098    else fprintf(foo," %d %ld YF\n", seqout, seqlen);
1099    }
1100
1101  if (interleaved) {
1102    /* interleave species lines in true output */
1103    /* nlines is # lines / sequence */
1104    short iline, j, leaf, iseq;
1105    char  *s = stempstore;
1106
1107    indexout();  noutindex--; /* mark eof */
1108
1109    for (leaf=0; leaf<nlines; leaf++) {
1110      if (outform == kMSF && leaf == 1) {
1111        fputs("//\n\n", foo);
1112        }
1113      if (outform == kPAUP && leaf==1) {
1114        switch (seqtype) {
1115          case kDNA     : cp= "dna"; break;
1116          case kRNA     : cp= "rna"; break;
1117          case kNucleic : cp= "dna"; break;
1118          case kAmino   : cp= "protein"; break;
1119          case kOtherSeq: cp= "dna"; break;
1120          }
1121        fprintf(foo,"\nbegin data;\n");
1122        fprintf(foo," dimensions ntax=%d nchar=%ld;\n", seqout, seqlen);
1123        /* fix by Ralf Westram (ARB): '-' means 'gap' ('.' means 'missing') */
1124        fprintf(foo," format datatype=%s interleave=yes missing=. gap=%c", cp, gPretty.gapchar);
1125        if (gPretty.domatch) fprintf(foo," matchchar=%c", gPretty.matchchar);
1126        fprintf(foo,";\n  matrix\n");
1127        }
1128
1129      for (iseq=0; iseq<noutindex; iseq++) {
1130        fseek(ftmp, outindex[iseq], 0);
1131        for (iline=0; iline<=leaf; iline++)
1132          if (!fgets(s, 256, ftmp)) *s= 0;
1133        if (ftell(ftmp) <= outindex[iseq+1])
1134          fputs( s, foo);
1135        }
1136
1137      for (j=0; j<gPretty.interline; j++)
1138        fputs( "\n", foo);  /* some want spacer line */
1139      }
1140    fclose(ftmp); /* tmp disappears */
1141    fout= foo;
1142    }
1143
1144  if (outform == kASN1)  fprintf( foo, "} }\n");
1145  if (outform == kPAUP)  fprintf( foo,";\n  end;\n");
1146
1147  if (outindex != NULL) free(outindex);
1148  exit_main(0);
1149}
1150
1151
Note: See TracBrowser for help on using the repository browser.