sources/readseq/ureadseq.c

   1 /* File: ureadseq.c
   2  *
   3  * Reads and writes nucleic/protein sequence in various
   4  * formats. Data files may have multiple sequences.
   5  *
   6  * Copyright 1990 by d.g.gilbert
   7  * biology dept., indiana university, bloomington, in 47405
   8  * e-mail: gilbertd@bio.indiana.edu
   9  *
  10  * This program may be freely copied and used by anyone.
  11  * Developers are encourged to incorporate parts in their
  12  * programs, rather than devise their own private sequence
  13  * format.
  14  *
  15  * This should compile and run with any ANSI C compiler.
  16  *
  17  * Alexander Sherstnev, 18/11/2013
  18  * renamed getline to locgetline
  19  *
  20  */
  21
  22
  23 #include <stdio.h>
  24 #include <ctype.h>
  25 #include <string.h>
  26
  27 #define UREADSEQ_G
  28 #include "ureadseq.h"
  29
  30 #pragma segment ureadseq
  31
  32
  33 int Strcasecmp(const char *a, const char *b)  /* from Nlm_StrICmp */
  34 {
  35   int diff, done;
  36   if (a == b)  return 0;
  37   done = 0;
  38   while (! done) {
  39     diff = to_upper(*a) - to_upper(*b);
  40     if (diff) return diff;
  41     if (*a == '\0') done = 1;
  42     else { a++; b++; }
  43     }
  44   return 0;
  45 }
  46
  47 int Strncasecmp(const char *a, const char *b, long maxn) /* from Nlm_StrNICmp */
  48 {
  49   int diff, done;
  50   if (a == b)  return 0;
  51   done = 0;
  52   while (! done) {
  53     diff = to_upper(*a) - to_upper(*b);
  54     if (diff) return diff;
  55     if (*a == '\0') done = 1;
  56     else {
  57       a++; b++; maxn--;
  58       if (! maxn) done = 1;
  59       }
  60     }
  61   return 0;
  62 }
  63
  64
  65
  66
  67
  68 #ifndef Local
  69 # define Local      static    /* local functions */
  70 #endif
  71
  72 #define kStartLength  500
  73
  74 const char *aminos      = "ABCDEFGHIKLMNPQRSTVWXYZ*";
  75 const char *primenuc    = "ACGTU";
  76 const char *protonly    = "EFIPQZ";
  77
  78 const char kNocountsymbols[5]  = "_.-?";
  79 const char stdsymbols[6]  = "_.-*?";
  80 const char allsymbols[32] = "_.-*?<>{}[]()!@#$%^&=+;:'/|`~\"\\";
  81 static const char *seqsymbols   = allsymbols;
  82
  83 const char nummask[11]   = "0123456789";
  84 const char nonummask[11] = "~!@#$%^&*(";
  85
  86 /*
  87     use general form of isseqchar -- all chars + symbols.
  88     no formats except nbrf (?) use symbols in data area as
  89     anything other than sequence chars.
  90 */
  91
  92
  93
  94                           /* Local variables for readSeq: */
  95 struct ReadSeqVars {
  96   short choice, err, nseq;
  97   long  seqlen, maxseq, seqlencount;
  98   short topnseq;
  99   long  topseqlen;
 100   const char *fname;
 101   char *seq, *seqid, matchchar;
 102   boolean allDone, done, filestart, addit;
 103   FILE  *f;
 104   long  linestart;
 105   char  s[256], *sp;
 106
 107   int (*isseqchar)();
 108   /* int  (*isseqchar)(int c);  << sgi cc hates (int c) */
 109 };
 110
 111
 112
 113 int isSeqChar(int c)
 114 {
 115   return (isalpha(c) || strchr(seqsymbols,c));
 116 }
 117
 118 int isSeqNumChar(int c)
 119 {
 120   return (isalnum(c) || strchr(seqsymbols,c));
 121 }
 122
 123
 124 int isAnyChar(int c)
 125 {
 126   return isascii(c); /* wrap in case isascii is macro */
 127 }
 128
 129 Local void readline(FILE *f, char *s, long *linestart)
 130 {
 131   char  *cp;
 132
 133   *linestart= ftell(f);
 134   if (NULL == fgets(s, 256, f))
 135     *s = 0;
 136   else {
 137     cp = strchr(s, '\n');
 138     if (cp != NULL) *cp = 0;
 139     }
 140 }
 141
 142 Local void locgetline(struct ReadSeqVars *V)
 143 {
 144   readline(V->f, V->s, &V->linestart);
 145 }
 146
 147 Local void ungetline(struct ReadSeqVars *V)
 148 {
 149   fseek(V->f, V->linestart, 0);
 150 }
 151
 152
 153 Local void addseq(char *s, struct ReadSeqVars *V)
 154 {
 155   char  *ptr;
 156
 157   if (V->addit) while (*s != 0) {
 158     if ((V->isseqchar)(*s)) {
 159       if (V->seqlen >= V->maxseq) {
 160         V->maxseq += kStartLength;
 161         ptr = (char*) realloc(V->seq, V->maxseq+1);
 162         if (ptr==NULL) {
 163           V->err = eMemFull;
 164           return;
 165           }
 166         else V->seq = ptr;
 167         }
 168       V->seq[(V->seqlen)++] = *s;
 169       }
 170     s++;
 171     }
 172 }
 173
 174 Local void countseq(char *s, struct ReadSeqVars *V)
 175  /* this must count all valid seq chars, for some formats (paup-sequential) even
 176     if we are skipping seq... */
 177 {
 178   while (*s != 0) {
 179     if ((V->isseqchar)(*s)) {
 180       (V->seqlencount)++;
 181       }
 182     s++;
 183     }
 184 }
 185
 186
 187 Local void addinfo(char *s, struct ReadSeqVars *V)
 188 {
 189   char s2[256], *si;
 190   boolean saveadd;
 191
 192   si = s2;
 193   while (*s == ' ') s++;
 194   sprintf(si, " %d)  %s\n", V->nseq, s);
 195
 196   saveadd = V->addit;
 197   V->addit = true;
 198   V->isseqchar = isAnyChar;
 199   addseq( si, V);
 200   V->addit = saveadd;
 201   V->isseqchar = isSeqChar;
 202 }
 203
 204
 205
 206
 207 Local void readLoop(short margin, boolean addfirst,
 208             boolean (*endTest)(boolean *addend, boolean *ungetend, struct ReadSeqVars *V),
 209             struct ReadSeqVars *V)
 210 {
 211   boolean addend = false;
 212   boolean ungetend = false;
 213
 214   V->nseq++;
 215   if (V->choice == kListSequences) V->addit = false;
 216   else V->addit = (V->nseq == V->choice);
 217   if (V->addit) V->seqlen = 0;
 218
 219   if (addfirst) addseq(V->s, V);
 220   do {
 221     locgetline(V);
 222     V->done = feof(V->f);
 223     V->done |= (*endTest)( &addend, &ungetend, V);
 224     if (V->addit && (addend || !V->done) && (strlen(V->s) > margin)) {
 225       addseq( (V->s)+margin, V);
 226     }
 227   } while (!V->done);
 228
 229   if (V->choice == kListSequences) addinfo(V->seqid, V);
 230   else {
 231     V->allDone = (V->nseq >= V->choice);
 232     if (V->allDone && ungetend) ungetline(V);
 233     }
 234 }
 235
 236
 237
 238 Local boolean endIG( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 239 {
 240   *addend = true; /* 1 or 2 occur in line w/ bases */
 241   *ungetend= false;
 242   return((strchr(V->s,'1')!=NULL) || (strchr(V->s,'2')!=NULL));
 243 }
 244
 245 Local void readIG(struct ReadSeqVars *V)
 246 {
 247 /* 18Aug92: new IG format -- ^L between sequences in place of ";" */
 248   char  *si;
 249
 250   while (!V->allDone) {
 251     do {
 252       locgetline(V);
 253       for (si= V->s; *si != 0 && *si < ' '; si++) *si= ' '; /* drop controls */
 254       if (*si == 0) *V->s= 0; /* chop line to empty */
 255     } while (! (feof(V->f) || ((*V->s != 0) && (*V->s != ';') ) ));
 256     if (feof(V->f))
 257       V->allDone = true;
 258     else {
 259       strcpy(V->seqid, V->s);
 260       readLoop(0, false, endIG, V);
 261       }
 262   }
 263 }
 264
 265
 266
 267 Local boolean endStrider( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 268 {
 269   *addend = false;
 270   *ungetend= false;
 271   return (strstr( V->s, "//") != NULL);
 272 }
 273
 274 Local void readStrider(struct ReadSeqVars *V)
 275 { /* ? only 1 seq/file ? */
 276
 277   while (!V->allDone) {
 278     locgetline(V);
 279     if (strstr(V->s,"; DNA sequence  ") == V->s)
 280       strcpy(V->seqid, (V->s)+16);
 281     else
 282       strcpy(V->seqid, (V->s)+1);
 283     while ((!feof(V->f)) && (*V->s == ';')) {
 284       locgetline(V);
 285       }
 286     if (feof(V->f)) V->allDone = true;
 287     else readLoop(0, true, endStrider, V);
 288   }
 289 }
 290
 291
 292 Local boolean endPIR( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 293 {
 294   *addend = false;
 295   *ungetend= (strstr(V->s,"ENTRY") == V->s);
 296   return ((strstr(V->s,"///") != NULL) || *ungetend);
 297 }
 298
 299 Local void readPIR(struct ReadSeqVars *V)
 300 { /*PIR -- many seqs/file */
 301
 302   while (!V->allDone) {
 303     while (! (feof(V->f) || strstr(V->s,"ENTRY")  || strstr(V->s,"SEQUENCE")) )
 304       locgetline(V);
 305     strcpy(V->seqid, (V->s)+16);
 306     while (! (feof(V->f) || strstr(V->s,"SEQUENCE") == V->s))
 307       locgetline(V);
 308     readLoop(0, false, endPIR, V);
 309
 310     if (!V->allDone) {
 311      while (! (feof(V->f) || ((*V->s != 0)
 312        && (strstr( V->s,"ENTRY") == V->s))))
 313         locgetline(V);
 314       }
 315     if (feof(V->f)) V->allDone = true;
 316   }
 317 }
 318
 319
 320 Local boolean endGB( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 321 {
 322   *addend = false;
 323   *ungetend= (strstr(V->s,"LOCUS") == V->s);
 324   return ((strstr(V->s,"//") != NULL) || *ungetend);
 325 }
 326
 327 Local void readGenBank(struct ReadSeqVars *V)
 328 { /*GenBank -- many seqs/file */
 329
 330   while (!V->allDone) {
 331     strcpy(V->seqid, (V->s)+12);
 332     while (! (feof(V->f) || strstr(V->s,"ORIGIN") == V->s))
 333       locgetline(V);
 334     readLoop(0, false, endGB, V);
 335
 336     if (!V->allDone) {
 337      while (! (feof(V->f) || ((*V->s != 0)
 338        && (strstr( V->s,"LOCUS") == V->s))))
 339         locgetline(V);
 340       }
 341     if (feof(V->f)) V->allDone = true;
 342   }
 343 }
 344
 345
 346 Local boolean endNBRF( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 347 {
 348   char  *a;
 349
 350   if ((a = strchr(V->s, '*')) != NULL) { /* end of 1st seq */
 351     /* "*" can be valid base symbol, drop it here */
 352     *a = 0;
 353     *addend = true;
 354     *ungetend= false;
 355     return(true);
 356     }
 357   else if (*V->s == '>') { /* start of next seq */
 358     *addend = false;
 359     *ungetend= true;
 360     return(true);
 361     }
 362   else
 363     return(false);
 364 }
 365
 366 Local void readNBRF(struct ReadSeqVars *V)
 367 {
 368   while (!V->allDone) {
 369     strcpy(V->seqid, (V->s)+4);
 370     locgetline(V);   /*skip title-junk line*/
 371     readLoop(0, false, endNBRF, V);
 372     if (!V->allDone) {
 373      while (!(feof(V->f) || (*V->s != 0 && *V->s == '>')))
 374         locgetline(V);
 375       }
 376     if (feof(V->f)) V->allDone = true;
 377   }
 378 }
 379
 380
 381
 382 Local boolean endPearson( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 383 {
 384   *addend = false;
 385   *ungetend= true;
 386   return(*V->s == '>');
 387 }
 388
 389 Local void readPearson(struct ReadSeqVars *V)
 390 {
 391   while (!V->allDone) {
 392     strcpy(V->seqid, (V->s)+1);
 393     readLoop(0, false, endPearson, V);
 394     if (!V->allDone) {
 395      while (!(feof(V->f) || ((*V->s != 0) && (*V->s == '>'))))
 396         locgetline(V);
 397       }
 398     if (feof(V->f)) V->allDone = true;
 399   }
 400 }
 401
 402
 403
 404 Local boolean endEMBL( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 405 {
 406   *addend = false;
 407   *ungetend= (strstr(V->s,"ID   ") == V->s);
 408   return ((strstr(V->s,"//") != NULL) || *ungetend);
 409 }
 410
 411 Local void readEMBL(struct ReadSeqVars *V)
 412 {
 413   while (!V->allDone) {
 414     strcpy(V->seqid, (V->s)+5);
 415     do {
 416       locgetline(V);
 417     } while (!(feof(V->f) | (strstr(V->s,"SQ   ") == V->s)));
 418
 419     readLoop(0, false, endEMBL, V);
 420     if (!V->allDone) {
 421       while (!(feof(V->f) |
 422          ((*V->s != '\0') & (strstr(V->s,"ID   ") == V->s))))
 423       locgetline(V);
 424     }
 425     if (feof(V->f)) V->allDone = true;
 426   }
 427 }
 428
 429
 430
 431 Local boolean endZuker( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 432 {
 433   *addend = false;
 434   *ungetend= true;
 435   return( *V->s == '(' );
 436 }
 437
 438 Local void readZuker(struct ReadSeqVars *V)
 439 {
 440   /*! 1st string is Zuker's Fortran format */
 441
 442   while (!V->allDone) {
 443     locgetline(V);  /*s == "seqLen seqid string..."*/
 444     strcpy(V->seqid, (V->s)+6);
 445     readLoop(0, false, endZuker, V);
 446     if (!V->allDone) {
 447       while (!(feof(V->f) |
 448         ((*V->s != '\0') & (*V->s == '('))))
 449           locgetline(V);
 450       }
 451     if (feof(V->f)) V->allDone = true;
 452   }
 453 }
 454
 455
 456
 457 Local boolean endFitch( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 458 {
 459   /* this is a somewhat shaky end,
 460     1st char of line is non-blank for seq. title
 461   */
 462   *addend = false;
 463   *ungetend= true;
 464   return( *V->s != ' ' );
 465 }
 466
 467 Local void readFitch(struct ReadSeqVars *V)
 468 {
 469   boolean first;
 470
 471   first = true;
 472   while (!V->allDone) {
 473     if (!first) strcpy(V->seqid, V->s);
 474     readLoop(0, first, endFitch, V);
 475     if (feof(V->f)) V->allDone = true;
 476     first = false;
 477     }
 478 }
 479
 480
 481 Local void readPlain(struct ReadSeqVars *V)
 482 {
 483   V->nseq++;
 484   V->addit = (V->choice > 0);
 485   if (V->addit) V->seqlen = 0;
 486   addseq(V->seqid, V);   /*from above..*/
 487   if (V->fname!=NULL) sprintf(V->seqid, "%s  [Unknown form]", V->fname);
 488   else sprintf(V->seqid, "  [Unknown form]");
 489   do {
 490     addseq(V->s, V);
 491     V->done = feof(V->f);
 492     locgetline(V);
 493   } while (!V->done);
 494   if (V->choice == kListSequences) addinfo(V->seqid, V);
 495   V->allDone = true;
 496 }
 497
 498
 499 Local void readUWGCG(struct ReadSeqVars *V)
 500 {
 501 /*
 502 10nov91: Reading GCG files casued duplication of last line when
 503          EOF followed that line !!!
 504     fix: locgetline now sets *V->s = 0
 505 */
 506   char  *si;
 507
 508   V->nseq++;
 509   V->addit = (V->choice > 0);
 510   if (V->addit) V->seqlen = 0;
 511   strcpy(V->seqid, V->s);
 512   /*writeseq: "    %s  Length: %d  (today)  Check: %d  ..\n" */
 513   /*drop above or ".." from id*/
 514   if (si = strstr(V->seqid,"  Length: ")) *si = 0;
 515   else if (si = strstr(V->seqid,"..")) *si = 0;
 516   do {
 517     V->done = feof(V->f);
 518     locgetline(V);
 519     if (!V->done) addseq((V->s), V);
 520   } while (!V->done);
 521   if (V->choice == kListSequences) addinfo(V->seqid, V);
 522   V->allDone = true;
 523 }
 524
 525
 526 Local void readOlsen(struct ReadSeqVars *V)
 527 { /* G. Olsen /print output from multiple sequence editor */
 528
 529   char    *si, *sj, *sk, *sm, sid[40], snum[20];
 530   boolean indata = false;
 531   int snumlen;
 532
 533   V->addit = (V->choice > 0);
 534   if (V->addit) V->seqlen = 0;
 535   rewind(V->f); V->nseq= 0;
 536   do {
 537     locgetline(V);
 538     V->done = feof(V->f);
 539
 540     if (V->done && !(*V->s)) break;
 541     else if (indata) {
 542       if ( (si= strstr(V->s, sid))
 543         /* && (strstr(V->s, snum) == si - snumlen - 1) ) { */
 544         && (sm= strstr(V->s, snum)) && (sm < si - snumlen) ) {
 545
 546         /* Spaces are valid alignment data !! */
 547 /* 17Oct91: Error, the left margin is 21 not 22! */
 548 /* dropped some nucs up to now -- my example file was right shifted ! */
 549 /* variable right id margin, drop id-2 spaces at end */
 550 /*
 551   VMS CC COMPILER (VAXC031) mess up:
 552   -- Index of 21 is chopping 1st nuc on VMS systems Only!
 553   Byte-for-byte same ame rnasep.olsen sequence file !
 554 */
 555
 556         /* si = (V->s)+21; < was this before VMS CC wasted my time */
 557         si += 10;  /* use strstr index plus offset to outfox VMS CC bug */
 558
 559         if (sk = strstr(si, sid)) *(sk-2) = 0;
 560         for (sk = si; *sk != 0; sk++) {
 561            if (*sk == ' ') *sk = '.';
 562            /* 18aug92: !! some olsen masks are NUMBERS !! which addseq eats */
 563            else if (isdigit(*sk)) *sk= nonummask[*sk - '0'];
 564            }
 565
 566         addseq(si, V);
 567         }
 568       }
 569
 570     else if (sk = strstr(V->s, "): ")) {  /* seq info header line */
 571   /* 18aug92: correct for diff seqs w/ same name -- use number, e.g. */
 572   /*   3 (Agr.tume):  agrobacterium.prna  18-JUN-1987 16:12 */
 573   /* 328 (Agr.tume):  agrobacterium.prna XYZ  19-DEC-1992   */
 574       (V->nseq)++;
 575       si = 1 + strchr(V->s,'(');
 576       *sk = ' ';
 577       if (V->choice == kListSequences) addinfo( si, V);
 578       else if (V->nseq == V->choice) {
 579         strcpy(V->seqid, si);
 580         sj = strchr(V->seqid, ':');
 581         while (*(--sj) == ' ') ;
 582         while (--sj != V->seqid) { if (*sj == ' ') *sj = '_'; }
 583
 584         *sk = 0;
 585         while (*(--sk) == ' ') *sk = 0;
 586         strcpy(sid, si);
 587
 588         si= V->s;
 589         while ((*si <= ' ') && (*si != 0)) si++;
 590         snumlen=0;
 591         while (si[snumlen] > ' ' && snumlen<20)
 592          { snum[snumlen]= si[snumlen]; snumlen++; }
 593         snum[snumlen]= 0;
 594         }
 595
 596       }
 597
 598     else if (strstr(V->s,"identity:   Data:")) {
 599       indata = true;
 600       if (V->choice == kListSequences) V->done = true;
 601       }
 602
 603   } while (!V->done);
 604
 605   V->allDone = true;
 606 } /*readOlsen*/
 607
 608
 609 Local void readMSF(struct ReadSeqVars *V)
 610 { /* gcg's MSF, mult. sequence format, interleaved ! */
 611
 612   char    *si, *sj, sid[128];
 613   boolean indata = false;
 614   int     atseq= 0, iline= 0;
 615
 616   V->addit = (V->choice > 0);
 617   if (V->addit) V->seqlen = 0;
 618   rewind(V->f); V->nseq= 0;
 619   do {
 620     locgetline(V);
 621     V->done = feof(V->f);
 622
 623     if (V->done && !(*V->s)) break;
 624     else if (indata) {
 625       /*somename  ...gpvedai .......t.. aaigr..vad tvgtgptnse aipaltaaet */
 626       /*       E  gvenae.kgv tentna.tad fvaqpvylpe .nqt...... kv.affynrs */
 627
 628       si= V->s;
 629       skipwhitespace(si);
 630       /* for (sj= si; isalnum(*sj); sj++) ; bug -- cdelwiche uses "-", "_" and others in names*/
 631       for (sj= si; *sj > ' '; sj++) ;
 632       *sj= 0;
 633       if ( *si ) {
 634         if ( (0==strcmp(si, sid)) ) {
 635           addseq(sj+1, V);
 636           }
 637         iline++;
 638         }
 639       }
 640
 641     else if (NULL != (si = strstr(V->s, "Name: "))) {  /* seq info header line */
 642       /* Name: somename      Len:   100  Check: 7009  Weight:  1.00 */
 643
 644       (V->nseq)++;
 645       si += 6;
 646       if (V->choice == kListSequences) addinfo( si, V);
 647       else if (V->nseq == V->choice) {
 648         strcpy(V->seqid, si);
 649         si = V->seqid;
 650         skipwhitespace(si);
 651         /* for (sj= si; isalnum(*sj); sj++) ; -- bug */
 652         for (sj= si; *sj > ' '; sj++) ;
 653         *sj= 0;
 654         strcpy(sid, si);
 655         }
 656       }
 657
 658     else if ( strstr(V->s,"//") /*== V->s*/ )  {
 659       indata = true;
 660       iline= 0;
 661       if (V->choice == kListSequences) V->done = true;
 662       }
 663
 664   } while (!V->done);
 665
 666
 667   V->allDone = true;
 668 } /*readMSF*/
 669
 670
 671
 672 Local void readPAUPinterleaved(struct ReadSeqVars *V)
 673 { /* PAUP mult. sequence format, interleaved or sequential! */
 674
 675   char    *si, *sj, *send, sid[40], sid1[40], saveseq[255];
 676   boolean first = true, indata = false, domatch;
 677   int     atseq= 0, iline= 0, ifmc, saveseqlen=0;
 678
 679 #define fixmatchchar(s) { \
 680   for (ifmc=0; ifmc<saveseqlen; ifmc++) \
 681     if (s[ifmc] == V->matchchar) s[ifmc]= saveseq[ifmc]; }
 682
 683   V->addit = (V->choice > 0);
 684   V->seqlencount = 0;
 685   if (V->addit) V->seqlen = 0;
 686   /* rewind(V->f); V->nseq= 0;  << do in caller !*/
 687   indata= true; /* call here after we find "matrix" */
 688   domatch= (V->matchchar > 0);
 689
 690   do {
 691     locgetline(V);
 692     V->done = feof(V->f);
 693
 694     if (V->done && !(*V->s)) break;
 695     else if (indata) {
 696       /* [         1                    1                    1         ]*/
 697       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
 698       /* chimp     ................a.t. .c.................a ..........*/
 699       /* !! need to correct for V->matchchar */
 700       si= V->s;
 701       skipwhitespace(si);
 702       if (strchr(si,';')) indata= false;
 703
 704       if (isalnum(*si))  {
 705         /* valid data line starts w/ a left-justified seq name in columns [0..8] */
 706         if (first) {
 707           (V->nseq)++;
 708           if (V->nseq >= V->topnseq) first= false;
 709           for (sj = si; isalnum(*sj); sj++) ;
 710           send= sj;
 711           skipwhitespace(sj);
 712           if (V->choice == kListSequences) {
 713             *send= 0;
 714             addinfo( si, V);
 715             }
 716           else if (V->nseq == V->choice) {
 717             if (domatch) {
 718               if (V->nseq == 1) { strcpy( saveseq, sj); saveseqlen= strlen(saveseq); }
 719               else fixmatchchar( sj);
 720               }
 721             addseq(sj, V);
 722             *send= 0;
 723             strcpy(V->seqid, si);
 724             strcpy(sid, si);
 725             if (V->nseq == 1) strcpy(sid1, sid);
 726             }
 727           }
 728
 729         else if ( (strstr(si, sid) == si) ){
 730           while (isalnum(*si)) si++;
 731           skipwhitespace(si);
 732           if (domatch) {
 733             if (V->nseq == 1) { strcpy( saveseq, si); saveseqlen= strlen(saveseq); }
 734             else fixmatchchar( si);
 735             }
 736           addseq(si, V);
 737           }
 738
 739         else if (domatch && (strstr(si, sid1) == si)) {
 740           strcpy( saveseq, si);
 741           saveseqlen= strlen(saveseq);
 742           }
 743
 744         iline++;
 745         }
 746       }
 747
 748     else if ( strstr(V->s,"matrix") )  {
 749       indata = true;
 750       iline= 0;
 751       if (V->choice == kListSequences) V->done = true;
 752       }
 753
 754   } while (!V->done);
 755
 756   V->allDone = true;
 757 } /*readPAUPinterleaved*/
 758
 759
 760
 761 Local void readPAUPsequential(struct ReadSeqVars *V)
 762 { /* PAUP mult. sequence format, interleaved or sequential! */
 763   char    *si, *sj;
 764   boolean atname = true, indata = false;
 765
 766   V->addit = (V->choice > 0);
 767   if (V->addit) V->seqlen = 0;
 768   V->seqlencount = 0;
 769   /* rewind(V->f); V->nseq= 0;  << do in caller !*/
 770   indata= true; /* call here after we find "matrix" */
 771   do {
 772     locgetline(V);
 773     V->done = feof(V->f);
 774
 775     if (V->done && !(*V->s)) break;
 776     else if (indata) {
 777       /* [         1                    1                    1         ]*/
 778       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
 779       /*           aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
 780       /* chimp     ................a.t. .c.................a ..........*/
 781       /*           ................a.t. .c.................a ..........*/
 782
 783       si= V->s;
 784       skipwhitespace(si);
 785       if (strchr(si,';')) indata= false;
 786       if (isalnum(*si))  {
 787         /* valid data line starts w/ a left-justified seq name in columns [0..8] */
 788         if (atname) {
 789           (V->nseq)++;
 790           V->seqlencount = 0;
 791           atname= false;
 792           sj= si+1;
 793           while (isalnum(*sj)) sj++;
 794           if (V->choice == kListSequences) {
 795             /* !! we must count bases to know when topseqlen is reached ! */
 796             countseq(sj, V);
 797             if (V->seqlencount >= V->topseqlen) atname= true;
 798             *sj= 0;
 799             addinfo( si, V);
 800             }
 801           else if (V->nseq == V->choice) {
 802             addseq(sj, V);
 803             V->seqlencount= V->seqlen;
 804             if (V->seqlencount >= V->topseqlen) atname= true;
 805             *sj= 0;
 806             strcpy(V->seqid, si);
 807             }
 808           else {
 809             countseq(sj, V);
 810             if (V->seqlencount >= V->topseqlen) atname= true;
 811             }
 812           }
 813
 814         else if (V->nseq == V->choice) {
 815           addseq(V->s, V);
 816           V->seqlencount= V->seqlen;
 817           if (V->seqlencount >= V->topseqlen) atname= true;
 818           }
 819         else {
 820           countseq(V->s, V);
 821           if (V->seqlencount >= V->topseqlen) atname= true;
 822           }
 823         }
 824       }
 825
 826     else if ( strstr(V->s,"matrix") )  {
 827       indata = true;
 828       atname= true;
 829       if (V->choice == kListSequences) V->done = true;
 830       }
 831
 832   } while (!V->done);
 833
 834   V->allDone = true;
 835 } /*readPAUPsequential*/
 836
 837
 838 Local void readPhylipInterleaved(struct ReadSeqVars *V)
 839 {
 840   char    *si, *sj;
 841   boolean first = true;
 842   int     iline= 0;
 843
 844   V->addit = (V->choice > 0);
 845   if (V->addit) V->seqlen = 0;
 846   V->seqlencount = 0;
 847   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); << topnseq == 0 !!! bad scan !! */
 848   si= V->s;
 849   skipwhitespace(si);
 850   V->topnseq= atoi(si);
 851   while (isdigit(*si)) si++;
 852   skipwhitespace(si);
 853   V->topseqlen= atol(si);
 854   /* fprintf(stderr,"Phylip-ileaf: topnseq=%d  topseqlen=%d\n",V->topnseq, V->topseqlen); */
 855
 856   do {
 857     locgetline(V);
 858     V->done = feof(V->f);
 859
 860     if (V->done && !(*V->s)) break;
 861     si= V->s;
 862     skipwhitespace(si);
 863     if (*si != 0) {
 864
 865       if (first) {  /* collect seq names + seq, as fprintf(outf,"%-10s  ",seqname); */
 866         (V->nseq)++;
 867         if (V->nseq >= V->topnseq) first= false;
 868         sj= V->s+10;  /* past name, start of data */
 869         if (V->choice == kListSequences) {
 870           *sj= 0;
 871           addinfo( si, V);
 872           }
 873         else if (V->nseq == V->choice) {
 874           addseq(sj, V);
 875           *sj= 0;
 876           strcpy(V->seqid, si);
 877           }
 878         }
 879       else if ( iline % V->nseq == V->choice -1 ) {
 880         addseq(si, V);
 881         }
 882       iline++;
 883     }
 884   } while (!V->done);
 885
 886   V->allDone = true;
 887 } /*readPhylipInterleaved*/
 888
 889
 890
 891 Local boolean endPhylipSequential( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
 892 {
 893   *addend = false;
 894   *ungetend= false;
 895   countseq( V->s, V);
 896   return V->seqlencount >= V->topseqlen;
 897 }
 898
 899 Local void readPhylipSequential(struct ReadSeqVars *V)
 900 {
 901   short  i;
 902   char  *si;
 903   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); < ? bad sscan ? */
 904   si= V->s;
 905   skipwhitespace(si);
 906   V->topnseq= atoi(si);
 907   while (isdigit(*si)) si++;
 908   skipwhitespace(si);
 909   V->topseqlen= atol(si);
 910   locgetline(V);
 911   while (!V->allDone) {
 912     V->seqlencount= 0;
 913     strncpy(V->seqid, (V->s), 10);
 914     V->seqid[10]= 0;
 915     for (i=0; i<10 && V->s[i]; i++) V->s[i]= ' ';
 916     readLoop(0, true, endPhylipSequential, V);
 917     if (feof(V->f)) V->allDone = true;
 918     }
 919 }
 920
 921
 922
 923
 924 Local void readSeqMain(
 925       struct ReadSeqVars *V,
 926       const long  skiplines_,
 927       const short format_)
 928 {
 929 #define tolowerstr(s) { long Itlwr, Ntlwr= strlen(s); \
 930   for (Itlwr=0; Itlwr<Ntlwr; Itlwr++) s[Itlwr]= to_lower(s[Itlwr]); }
 931
 932   boolean gotuw;
 933   long l;
 934
 935   V->linestart= 0;
 936   V->matchchar= 0;
 937   if (V->f == NULL)
 938     V->err = eFileNotFound;
 939   else {
 940
 941     for (l = skiplines_; l > 0; l--) locgetline( V);
 942
 943     do {
 944       locgetline( V);
 945       for (l= strlen(V->s); (l > 0) && (V->s[l] == ' '); l--) ;
 946     } while ((l == 0) && !feof(V->f));
 947
 948     if (feof(V->f)) V->err = eNoData;
 949
 950     else switch (format_) {
 951       case kPlain : readPlain(V); break;
 952       case kIG    : readIG(V); break;
 953       case kStrider: readStrider(V); break;
 954       case kGenBank: readGenBank(V); break;
 955       case kPIR   : readPIR(V); break;
 956       case kNBRF  : readNBRF(V); break;
 957       case kPearson: readPearson(V); break;
 958       case kEMBL  : readEMBL(V); break;
 959       case kZuker : readZuker(V); break;
 960       case kOlsen : readOlsen(V); break;
 961       case kMSF   : readMSF(V); break;
 962
 963       case kPAUP    : {
 964         boolean done= false;
 965         boolean interleaved= false;
 966         char  *cp;
 967         /* rewind(V->f); V->nseq= 0; ?? assume it is at top ?? skiplines ... */
 968         while (!done) {
 969           locgetline( V);
 970           tolowerstr( V->s);
 971           if (strstr( V->s, "matrix")) done= true;
 972           if (strstr( V->s, "interleav")) interleaved= true;
 973           if (NULL != (cp=strstr( V->s, "ntax=")) )  V->topnseq= atoi(cp+5);
 974           if (NULL != (cp=strstr( V->s, "nchar=")) )  V->topseqlen= atoi(cp+6);
 975           if (NULL != (cp=strstr( V->s, "matchchar=")) )  {
 976             cp += 10;
 977             if (*cp=='\'') cp++;
 978             else if (*cp=='"') cp++;
 979             V->matchchar= *cp;
 980             }
 981           }
 982         if (interleaved) readPAUPinterleaved(V);
 983         else readPAUPsequential(V);
 984         }
 985         break;
 986
 987       /* kPhylip: ! can't determine in middle of file which type it is...*/
 988       /* test for interleave or sequential and use Phylip4(ileave) or Phylip2 */
 989       case kPhylip2:
 990         readPhylipSequential(V);
 991         break;
 992       case kPhylip4: /* == kPhylip3 */
 993         readPhylipInterleaved(V);
 994         break;
 995
 996       default:
 997         V->err = eUnknownFormat;
 998         break;
 999
1000       case kFitch :
1001         strcpy(V->seqid, V->s); locgetline(V);
1002         readFitch(V);
1003         break;
1004
1005       case kGCG:
1006         do {
1007           gotuw = (strstr(V->s,"..") != NULL);
1008           if (gotuw) readUWGCG(V);
1009           locgetline(V);
1010         } while (!(feof(V->f) || V->allDone));
1011         break;
1012       }
1013     }
1014
1015   V->filestart= false;
1016   V->seq[V->seqlen] = 0; /* stick a string terminator on it */
1017 }
1018
1019
1020 char *readSeqFp(
1021       const short whichEntry_,  /* index to sequence in file */
1022       FILE  *fp_,   /* pointer to open seq file */
1023       const long  skiplines_,
1024       const short format_,      /* sequence file format */
1025       long  *seqlen_,     /* return seq size */
1026       short *nseq_,       /* number of seqs in file, for listSeqs() */
1027       short *error_,      /* return error */
1028       char  *seqid_)      /* return seq name/info */
1029 {
1030   struct ReadSeqVars V;
1031
1032   if (format_ < kMinFormat || format_ > kMaxFormat) {
1033     *error_ = eUnknownFormat;
1034     *seqlen_ = 0;
1035     return NULL;
1036     }
1037
1038   V.choice = whichEntry_;
1039   V.fname  = NULL;  /* don't know */
1040   V.seq    = (char*) calloc(1, kStartLength+1);
1041   V.maxseq = kStartLength;
1042   V.seqlen = 0;
1043   V.seqid  = seqid_;
1044
1045   V.f = fp_;
1046   V.filestart= (ftell( fp_) == 0);
1047   /* !! in sequential read, must remove current seq position from choice/whichEntry_ counter !! ... */
1048   if (V.filestart)  V.nseq = 0;
1049   else V.nseq= *nseq_;  /* track where we are in file...*/
1050
1051   *V.seqid = '\0';
1052   V.err = 0;
1053   V.nseq = 0;
1054   V.isseqchar = isSeqChar;
1055   if (V.choice == kListSequences) ; /* leave as is */
1056   else if (V.choice <= 0) V.choice = 1; /* default ?? */
1057   V.addit = (V.choice > 0);
1058   V.allDone = false;
1059
1060   readSeqMain(&V, skiplines_, format_);
1061
1062   *error_ = V.err;
1063   *seqlen_ = V.seqlen;
1064   *nseq_ = V.nseq;
1065   return V.seq;
1066 }
1067
1068 char *readSeq(
1069       const short whichEntry_,  /* index to sequence in file */
1070       const char  *filename_,   /* file name */
1071       const long  skiplines_,
1072       const short format_,      /* sequence file format */
1073       long  *seqlen_,     /* return seq size */
1074       short *nseq_,       /* number of seqs in file, for listSeqs() */
1075       short *error_,      /* return error */
1076       char  *seqid_)      /* return seq name/info */
1077 {
1078   struct ReadSeqVars V;
1079
1080   if (format_ < kMinFormat || format_ > kMaxFormat) {
1081     *error_ = eUnknownFormat;
1082     *seqlen_ = 0;
1083     return NULL;
1084     }
1085
1086   V.choice = whichEntry_;
1087   V.fname  = filename_;  /* don't need to copy string, just ptr to it */
1088   V.seq    = (char*) calloc(1, kStartLength+1);
1089   V.maxseq = kStartLength;
1090   V.seqlen = 0;
1091   V.seqid  = seqid_;
1092
1093   V.f = NULL;
1094   *V.seqid = '\0';
1095   V.err = 0;
1096   V.nseq = 0;
1097   V.isseqchar = isSeqChar;
1098   if (V.choice == kListSequences) ; /* leave as is */
1099   else if (V.choice <= 0) V.choice = 1; /* default ?? */
1100   V.addit = (V.choice > 0);
1101   V.allDone = false;
1102
1103   V.f = fopen(V.fname, "r");
1104   V.filestart= true;
1105
1106   readSeqMain(&V, skiplines_, format_);
1107
1108   if (V.f != NULL) fclose(V.f);
1109   *error_ = V.err;
1110   *seqlen_ = V.seqlen;
1111   *nseq_ = V.nseq;
1112   return V.seq;
1113 }
1114
1115
1116
1117
1118
1119 char *listSeqs(
1120       const char  *filename_,   /* file name */
1121       const long skiplines_,
1122       const short format_,      /* sequence file format */
1123       short *nseq_,       /* number of seqs in file, for listSeqs() */
1124       short *error_)      /* return error */
1125 {
1126   char  seqid[256];
1127   long  seqlen;
1128
1129   return readSeq( kListSequences, filename_, skiplines_, format_,
1130                   &seqlen, nseq_, error_, seqid);
1131 }
1132
1133
1134
1135
1136 short seqFileFormat(    /* return sequence format number, see ureadseq.h */
1137     const char *filename,
1138     long  *skiplines,   /* return #lines to skip any junk like mail header */
1139     short *error)       /* return any error value or 0 */
1140 {
1141   FILE      *fseq;
1142   short      format;
1143
1144   fseq  = fopen(filename, "r");
1145   format= seqFileFormatFp( fseq, skiplines, error);
1146   if (fseq!=NULL) fclose(fseq);
1147   return format;
1148 }
1149
1150 short seqFileFormatFp(
1151     FILE *fseq,
1152     long  *skiplines,   /* return #lines to skip any junk like mail header */
1153     short *error)       /* return any error value or 0 */
1154 {
1155   boolean   foundDNA= false, foundIG= false, foundStrider= false,
1156             foundGB= false, foundPIR= false, foundEMBL= false, foundNBRF= false,
1157             foundPearson= false, foundFitch= false, foundPhylip= false, foundZuker= false,
1158             gotolsen= false, gotpaup = false, gotasn1 = false, gotuw= false, gotMSF= false,
1159             isfitch= false,  isphylip= false, done= false;
1160   short     format= kUnknown;
1161   int       nlines= 0, k, splen= 0, otherlines= 0, aminolines= 0, dnalines= 0;
1162   char      sp[256];
1163   long      linestart=0;
1164   int     maxlines2check=500;
1165
1166 #define ReadOneLine(sp)   \
1167   { done |= (feof(fseq)); \
1168     readline( fseq, sp, &linestart);  \
1169     if (!done) { splen = strlen(sp); ++nlines; } }
1170
1171   *skiplines = 0;
1172   *error = 0;
1173   if (fseq == NULL) { *error = eFileNotFound;  return kNoformat; }
1174
1175   while ( !done ) {
1176     ReadOneLine(sp);
1177
1178     /* check for mailer head & skip past if found */
1179     if (nlines < 4 && !done) {
1180       if ((strstr(sp,"From ") == sp) || (strstr(sp,"Received:") == sp)) {
1181         do {
1182           /* skip all lines until find one blank line */
1183           ReadOneLine(sp);
1184           if (!done) for (k=0; (k<splen) && (sp[k]==' '); k++) ;
1185           } while ((!done) && (k < splen));
1186         *skiplines = nlines; /* !? do we want #lines or #bytes ?? */
1187         }
1188       }
1189
1190     if (sp==NULL || *sp==0)
1191       ; /* nada */
1192
1193     /* high probability identities: */
1194
1195     else if ( strstr(sp,"MSF:") && strstr(sp,"Type:") && strstr(sp,"Check:") )
1196       gotMSF= true;
1197
1198     else if ((strstr(sp,"..") != NULL) && (strstr(sp,"Check:") != NULL))
1199       gotuw= true;
1200
1201     else if (strstr(sp,"identity:   Data:") != NULL)
1202       gotolsen= true;
1203
1204     else if ( strstr(sp,"::=") &&
1205       (strstr(sp,"Bioseq") ||       /* Bioseq or Bioseq-set */
1206        strstr(sp,"Seq-entry") ||
1207        strstr(sp,"Seq-submit") ) )  /* can we read submit format? */
1208           gotasn1= true;
1209
1210     else if ( strstr(sp,"#NEXUS") == sp )
1211       gotpaup= true;
1212
1213     /* uncertain identities: */
1214
1215     else if (*sp ==';') {
1216       if (strstr(sp,"Strider") !=NULL) foundStrider= true;
1217       else foundIG= true;
1218       }
1219
1220     else if (strstr(sp,"LOCUS") == sp)
1221       foundGB= true;
1222     else if (strstr(sp,"ORIGIN") == sp)
1223       foundGB= true;
1224
1225     else if (strstr(sp,"ENTRY   ") == sp)  /* ? also (strcmp(sp,"\\\\\\")==0) */
1226       foundPIR= true;
1227     else if (strstr(sp,"SEQUENCE") == sp)
1228       foundPIR= true;
1229
1230     else if (*sp == '>') {
1231       if (sp[3] == ';') foundNBRF= true;
1232       else foundPearson= true;
1233       }
1234
1235     else if (strstr(sp,"ID   ") == sp)
1236       foundEMBL= true;
1237     else if (strstr(sp,"SQ   ") == sp)
1238       foundEMBL= true;
1239
1240     else if (*sp == '(')
1241       foundZuker= true;
1242
1243     else {
1244       if (nlines - *skiplines == 1) {
1245         int ispp= 0, ilen= 0;
1246         sscanf( sp, "%d%d", &ispp, &ilen);
1247         if (ispp > 0 && ilen > 0) isphylip= true;
1248         }
1249       else if (isphylip && nlines - *skiplines == 2) {
1250         int  tseq;
1251         tseq= getseqtype(sp+10, strlen(sp+10));
1252         if ( isalpha(*sp)     /* 1st letter in 2nd line must be of a name */
1253          && (tseq != kOtherSeq))  /* sequence section must be okay */
1254             foundPhylip= true;
1255         }
1256
1257       for (k=0, isfitch= true; isfitch & (k < splen); k++) {
1258         if (k % 4 == 0) isfitch &= (sp[k] == ' ');
1259         else isfitch &= (sp[k] != ' ');
1260         }
1261       if (isfitch & (splen > 20)) foundFitch= true;
1262
1263       /* kRNA && kDNA are fairly certain...*/
1264       switch (getseqtype( sp, splen)) {
1265         case kOtherSeq: otherlines++; break;
1266         case kAmino   : if (splen>20) aminolines++; break;
1267         case kDNA     :
1268         case kRNA     : if (splen>20) dnalines++; break;
1269         case kNucleic : break; /* not much info ? */
1270         }
1271
1272       }
1273
1274                     /* pretty certain */
1275     if (gotolsen) {
1276       format= kOlsen;
1277       done= true;
1278       }
1279     else if (gotMSF) {
1280       format= kMSF;
1281       done= true;
1282       }
1283     else if (gotasn1) {
1284       /* !! we need to look further and return  kASNseqentry | kASNseqset */
1285       /*
1286         seqentry key is Seq-entry ::=
1287         seqset key is Bioseq-set ::=
1288         ?? can't read these yet w/ ncbi tools ??
1289           Seq-submit ::=
1290           Bioseq ::=  << fails both bioseq-seq and seq-entry parsers !
1291       */
1292       if (strstr(sp,"Bioseq-set")) format= kASNseqset;
1293       else if (strstr(sp,"Seq-entry")) format= kASNseqentry;
1294       else format= kASN1;  /* other form, we can't yet read... */
1295       done= true;
1296       }
1297     else if (gotpaup) {
1298       format= kPAUP;
1299       done= true;
1300       }
1301
1302     else if (gotuw) {
1303       if (foundIG) format= kIG;  /* a TOIG file from GCG for certain */
1304       else format= kGCG;
1305       done= true;
1306       }
1307
1308     else if ((dnalines > 1) || done || (nlines > maxlines2check)) {
1309           /* decide on most likely format */
1310           /* multichar idents: */
1311       if (foundStrider) format= kStrider;
1312       else if (foundGB) format= kGenBank;
1313       else if (foundPIR) format= kPIR;
1314       else if (foundEMBL) format= kEMBL;
1315       else if (foundNBRF) format= kNBRF;
1316           /* single char idents: */
1317       else if (foundIG) format= kIG;
1318       else if (foundPearson) format= kPearson;
1319       else if (foundZuker) format= kZuker;
1320           /* digit ident: */
1321       else if (foundPhylip) format= kPhylip;
1322           /* spacing ident: */
1323       else if (foundFitch) format= kFitch;
1324           /* no format chars: */
1325       else if (otherlines > 0) format= kUnknown;
1326       else if (dnalines > 1) format= kPlain;
1327       else if (aminolines > 1) format= kPlain;
1328       else format= kUnknown;
1329
1330       done= true;
1331       }
1332
1333     /* need this for possible long header in olsen format */
1334      else if (strstr(sp,"): ") != NULL)
1335        maxlines2check++;
1336     }
1337
1338   if (format == kPhylip) {
1339     /* check for interleaved or sequential -- really messy */
1340     int tname, tseq;
1341     long i, j, nspp= 0, nlen= 0, ilen, leaf= 0, seq= 0;
1342     char  *ps;
1343
1344     rewind(fseq);
1345     for (i=0; i < *skiplines; i++) ReadOneLine(sp);
1346     nlines= 0;
1347     ReadOneLine(sp);
1348     sscanf( sp, "%d%d", &nspp, &nlen);
1349     ReadOneLine(sp); /* 1st seq line */
1350     for (ps= sp+10, ilen=0; *ps!=0; ps++) if (isprint(*ps)) ilen++;
1351
1352     for (i= 1; i<nspp; i++) {
1353       ReadOneLine(sp);
1354
1355       tseq= getseqtype(sp+10, strlen(sp+10));
1356       tname= getseqtype(sp, 10);
1357       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++);
1358       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
1359
1360       /* find probable interleaf or sequential ... */
1361       if (j>=9) seq += 10; /* pretty certain not ileaf */
1362       else {
1363         if (tseq != tname) leaf++; else seq++;
1364         if (tname == kDNA || tname == kRNA) seq++; else leaf++;
1365         }
1366
1367       if (ilen <= nlen && j<9) {
1368         if (tname == kOtherSeq) leaf += 10;
1369         else if (tname == kAmino || tname == kDNA || tname == kRNA) seq++; else leaf++;
1370         }
1371       else if (ilen > nlen) {
1372         ilen= 0;
1373         }
1374       }
1375     for ( nspp *= 2 ; i<nspp; i++) {  /* this should be only bases if interleaf */
1376       ReadOneLine(sp);
1377
1378       tseq= getseqtype(sp+10, strlen(sp+10));
1379       tname= getseqtype(sp, 10);
1380       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
1381       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++);
1382       if (j<9) {
1383         if (tname == kOtherSeq) seq += 10;
1384         if (tseq != tname) seq++; else leaf++;
1385         if (tname == kDNA || tname == kRNA) leaf++; else seq++;
1386         }
1387       if (ilen > nlen) {
1388         if (j>9) leaf += 10; /* must be a name here for sequent */
1389         else if (tname == kOtherSeq) seq += 10;
1390         ilen= 0;
1391         }
1392       }
1393
1394     if (leaf > seq) format= kPhylip4;
1395     else format= kPhylip2;
1396     }
1397
1398   return(format);
1399 #undef  ReadOneLine
1400 } /* SeqFileFormat */
1401
1402
1403
1404
1405 unsigned long GCGchecksum( const char *seq, const long seqlen, unsigned long *checktotal)
1406 /* GCGchecksum */
1407 {
1408   register long  i, check = 0, count = 0;
1409
1410   for (i = 0; i < seqlen; i++) {
1411     count++;
1412     check += count * to_upper(seq[i]);
1413     if (count == 57) count = 0;
1414     }
1415   check %= 10000;
1416   *checktotal += check;
1417   *checktotal %= 10000;
1418   return check;
1419 }
1420
1421 /* Table of CRC-32's of all single byte values (made by makecrc.c of ZIP source) */
1422 const unsigned long crctab[] = {
1423   0x00000000L, 0x77073096L, 0xee0e612cL, 0x990951baL, 0x076dc419L,
1424   0x706af48fL, 0xe963a535L, 0x9e6495a3L, 0x0edb8832L, 0x79dcb8a4L,
1425   0xe0d5e91eL, 0x97d2d988L, 0x09b64c2bL, 0x7eb17cbdL, 0xe7b82d07L,
1426   0x90bf1d91L, 0x1db71064L, 0x6ab020f2L, 0xf3b97148L, 0x84be41deL,
1427   0x1adad47dL, 0x6ddde4ebL, 0xf4d4b551L, 0x83d385c7L, 0x136c9856L,
1428   0x646ba8c0L, 0xfd62f97aL, 0x8a65c9ecL, 0x14015c4fL, 0x63066cd9L,
1429   0xfa0f3d63L, 0x8d080df5L, 0x3b6e20c8L, 0x4c69105eL, 0xd56041e4L,
1430   0xa2677172L, 0x3c03e4d1L, 0x4b04d447L, 0xd20d85fdL, 0xa50ab56bL,
1431   0x35b5a8faL, 0x42b2986cL, 0xdbbbc9d6L, 0xacbcf940L, 0x32d86ce3L,
1432   0x45df5c75L, 0xdcd60dcfL, 0xabd13d59L, 0x26d930acL, 0x51de003aL,
1433   0xc8d75180L, 0xbfd06116L, 0x21b4f4b5L, 0x56b3c423L, 0xcfba9599L,
1434   0xb8bda50fL, 0x2802b89eL, 0x5f058808L, 0xc60cd9b2L, 0xb10be924L,
1435   0x2f6f7c87L, 0x58684c11L, 0xc1611dabL, 0xb6662d3dL, 0x76dc4190L,
1436   0x01db7106L, 0x98d220bcL, 0xefd5102aL, 0x71b18589L, 0x06b6b51fL,
1437   0x9fbfe4a5L, 0xe8b8d433L, 0x7807c9a2L, 0x0f00f934L, 0x9609a88eL,
1438   0xe10e9818L, 0x7f6a0dbbL, 0x086d3d2dL, 0x91646c97L, 0xe6635c01L,
1439   0x6b6b51f4L, 0x1c6c6162L, 0x856530d8L, 0xf262004eL, 0x6c0695edL,
1440   0x1b01a57bL, 0x8208f4c1L, 0xf50fc457L, 0x65b0d9c6L, 0x12b7e950L,
1441   0x8bbeb8eaL, 0xfcb9887cL, 0x62dd1ddfL, 0x15da2d49L, 0x8cd37cf3L,
1442   0xfbd44c65L, 0x4db26158L, 0x3ab551ceL, 0xa3bc0074L, 0xd4bb30e2L,
1443   0x4adfa541L, 0x3dd895d7L, 0xa4d1c46dL, 0xd3d6f4fbL, 0x4369e96aL,
1444   0x346ed9fcL, 0xad678846L, 0xda60b8d0L, 0x44042d73L, 0x33031de5L,
1445   0xaa0a4c5fL, 0xdd0d7cc9L, 0x5005713cL, 0x270241aaL, 0xbe0b1010L,
1446   0xc90c2086L, 0x5768b525L, 0x206f85b3L, 0xb966d409L, 0xce61e49fL,
1447   0x5edef90eL, 0x29d9c998L, 0xb0d09822L, 0xc7d7a8b4L, 0x59b33d17L,
1448   0x2eb40d81L, 0xb7bd5c3bL, 0xc0ba6cadL, 0xedb88320L, 0x9abfb3b6L,
1449   0x03b6e20cL, 0x74b1d29aL, 0xead54739L, 0x9dd277afL, 0x04db2615L,
1450   0x73dc1683L, 0xe3630b12L, 0x94643b84L, 0x0d6d6a3eL, 0x7a6a5aa8L,
1451   0xe40ecf0bL, 0x9309ff9dL, 0x0a00ae27L, 0x7d079eb1L, 0xf00f9344L,
1452   0x8708a3d2L, 0x1e01f268L, 0x6906c2feL, 0xf762575dL, 0x806567cbL,
1453   0x196c3671L, 0x6e6b06e7L, 0xfed41b76L, 0x89d32be0L, 0x10da7a5aL,
1454   0x67dd4accL, 0xf9b9df6fL, 0x8ebeeff9L, 0x17b7be43L, 0x60b08ed5L,
1455   0xd6d6a3e8L, 0xa1d1937eL, 0x38d8c2c4L, 0x4fdff252L, 0xd1bb67f1L,
1456   0xa6bc5767L, 0x3fb506ddL, 0x48b2364bL, 0xd80d2bdaL, 0xaf0a1b4cL,
1457   0x36034af6L, 0x41047a60L, 0xdf60efc3L, 0xa867df55L, 0x316e8eefL,
1458   0x4669be79L, 0xcb61b38cL, 0xbc66831aL, 0x256fd2a0L, 0x5268e236L,
1459   0xcc0c7795L, 0xbb0b4703L, 0x220216b9L, 0x5505262fL, 0xc5ba3bbeL,
1460   0xb2bd0b28L, 0x2bb45a92L, 0x5cb36a04L, 0xc2d7ffa7L, 0xb5d0cf31L,
1461   0x2cd99e8bL, 0x5bdeae1dL, 0x9b64c2b0L, 0xec63f226L, 0x756aa39cL,
1462   0x026d930aL, 0x9c0906a9L, 0xeb0e363fL, 0x72076785L, 0x05005713L,
1463   0x95bf4a82L, 0xe2b87a14L, 0x7bb12baeL, 0x0cb61b38L, 0x92d28e9bL,
1464   0xe5d5be0dL, 0x7cdcefb7L, 0x0bdbdf21L, 0x86d3d2d4L, 0xf1d4e242L,
1465   0x68ddb3f8L, 0x1fda836eL, 0x81be16cdL, 0xf6b9265bL, 0x6fb077e1L,
1466   0x18b74777L, 0x88085ae6L, 0xff0f6a70L, 0x66063bcaL, 0x11010b5cL,
1467   0x8f659effL, 0xf862ae69L, 0x616bffd3L, 0x166ccf45L, 0xa00ae278L,
1468   0xd70dd2eeL, 0x4e048354L, 0x3903b3c2L, 0xa7672661L, 0xd06016f7L,
1469   0x4969474dL, 0x3e6e77dbL, 0xaed16a4aL, 0xd9d65adcL, 0x40df0b66L,
1470   0x37d83bf0L, 0xa9bcae53L, 0xdebb9ec5L, 0x47b2cf7fL, 0x30b5ffe9L,
1471   0xbdbdf21cL, 0xcabac28aL, 0x53b39330L, 0x24b4a3a6L, 0xbad03605L,
1472   0xcdd70693L, 0x54de5729L, 0x23d967bfL, 0xb3667a2eL, 0xc4614ab8L,
1473   0x5d681b02L, 0x2a6f2b94L, 0xb40bbe37L, 0xc30c8ea1L, 0x5a05df1bL,
1474   0x2d02ef8dL
1475 };
1476
1477 unsigned long CRC32checksum(const char *seq, const long seqlen, unsigned long *checktotal)
1478 /*CRC32checksum: modified from CRC-32 algorithm found in ZIP compression source */
1479 {
1480   register unsigned long c = 0xffffffffL;
1481   register long n = seqlen;
1482
1483   while (n--) {
1484     c = crctab[((int)c ^ (to_upper(*seq))) & 0xff] ^ (c >> 8);
1485     seq++; /* fixed aug'98 finally */
1486     }
1487   c= c ^ 0xffffffffL;
1488   *checktotal += c;
1489   return c;
1490 }
1491
1492
1493
1494
1495 short getseqtype( const char *seq, const long seqlen)
1496 { /* return sequence kind: kDNA, kRNA, kProtein, kOtherSeq, ??? */
1497   char  c;
1498   short i, maxtest;
1499   short na = 0, aa = 0, po = 0, nt = 0, nu = 0, ns = 0, no = 0;
1500
1501   maxtest = min(300, seqlen);
1502   for (i = 0; i < maxtest; i++) {
1503     c = to_upper(seq[i]);
1504     if (strchr(protonly, c)) po++;
1505     else if (strchr(primenuc,c)) {
1506       na++;
1507       if (c == 'T') nt++;
1508       else if (c == 'U') nu++;
1509       }
1510     else if (strchr(aminos,c)) aa++;
1511     else if (strchr(seqsymbols,c)) ns++;
1512     else if (isalpha(c)) no++;
1513     }
1514
1515   if ((no > 0) || (po+aa+na == 0)) return kOtherSeq;
1516   /* ?? test for probability of kOtherSeq ?, e.g.,
1517   else if (po+aa+na / maxtest < 0.70) return kOtherSeq;
1518   */
1519   else if (po > 0) return kAmino;
1520   else if (aa == 0) {
1521     if (nu > nt) return kRNA;
1522     else return kDNA;
1523     }
1524   else if (na > aa) return kNucleic;
1525   else return kAmino;
1526 } /* getseqtype */
1527
1528
1529 char* compressSeq( const char gapc, const char *seq, const long seqlen, long *newlen)
1530 {
1531   register char *a, *b;
1532   register long i;
1533   char  *newseq;
1534
1535   *newlen= 0;
1536   if (!seq) return NULL;
1537   newseq = (char*) malloc(seqlen+1);
1538   if (!newseq) return NULL;
1539   for (a= (char*)seq, b=newseq, i=0; *a!=0; a++)
1540     if (*a != gapc) {
1541       *b++= *a;
1542       i++;
1543       }
1544   *b= '\0';
1545   newseq = (char*) realloc(newseq, i+1);
1546   *newlen= i;
1547   return newseq;
1548 }
1549
1550
1551
1552 /***
1553 char *rtfhead = "{\\rtf1\\defformat\\mac\\deff2 \
1554 {\\fonttbl\
1555   {\\f1\\fmodern Courier;}{\\f2\\fmodern Monaco;}\
1556   {\\f3\\fswiss Helvetica;}{\\f4\\fswiss Geneva;}\
1557   {\\f5\\froman Times;}{\\f6\\froman Palatino;}\
1558   {\\f7\\froman New Century Schlbk;}{\\f8\\ftech Symbol;}}\
1559 {\\stylesheet\
1560   {\\s1 \\f5\\fs20 \\sbasedon0\\snext1 name;}\
1561   {\\s2 \\f3\\fs20 \\sbasedon0\\snext2 num;}\
1562   {\\s3 \\f1\\f21 \\sbasedon0\\snext3 seq;}}";
1563
1564 char *rtftail = "}";
1565 ****/
1566
1567 short writeSeq(FILE *outf, const char *seq, const long seqlen,
1568                 const short outform, const char *seqid)
1569 /* dump sequence to standard output */
1570 {
1571   const short kSpaceAll = -9;
1572 #define kMaxseqwidth  250
1573
1574   boolean baseonlynum= false; /* nocountsymbols -- only count true bases, not "-" */
1575   short  numline = 0; /* only true if we are writing seq number line (for interleave) */
1576   boolean numright = false, numleft = false;
1577   boolean nameright = false, nameleft = false;
1578   short   namewidth = 8, numwidth = 8;
1579   short   spacer = 0, width  = 50, tab = 0;
1580   /* new parameters: width, spacer, those above... */
1581
1582   short linesout = 0, seqtype = kNucleic;
1583   long  i, j, l, l1, ibase;
1584   char  idword[31], endstr[10];
1585   char  seqnamestore[128], *seqname = seqnamestore;
1586   char  s[kMaxseqwidth], *cp;
1587   char  nameform[10], numform[10], nocountsymbols[10];
1588   unsigned long checksum = 0, checktotal = 0;
1589
1590   gPretty.atseq++;
1591   skipwhitespace(seqid);
1592   l = min(128, strlen(seqid));
1593   strncpy( seqnamestore, seqid, l);
1594   seqname[l] = 0;
1595
1596   sscanf( seqname, "%30s", idword);
1597   sprintf(numform, "%d", seqlen);
1598   numwidth= strlen(numform)+1;
1599   nameform[0]= '\0';
1600
1601   if (strstr(seqname,"checksum") != NULL) {
1602     cp = strstr(seqname,"bases");
1603     if (cp!=NULL) {
1604       for ( ; (cp!=seqname) && (*cp!=','); cp--) ;
1605       if (cp!=seqname) *cp=0;
1606       }
1607     }
1608
1609   strcpy( endstr,"");
1610   l1 = 0;
1611
1612   if (outform == kGCG || outform == kMSF)
1613     checksum = GCGchecksum(seq, seqlen, &checktotal);
1614   else
1615     checksum = seqchecksum(seq, seqlen, &checktotal);
1616
1617   switch (outform) {
1618
1619     case kPlain:
1620     case kUnknown:    /* no header, just sequence */
1621       strcpy(endstr,"\n"); /* end w/ extra blank line */
1622       break;
1623
1624     case kOlsen:  /* Olsen seq. editor takes plain nucs OR Genbank  */
1625     case kGenBank:
1626       fprintf(outf,"LOCUS       %s       %d bp\n", idword, seqlen);
1627       fprintf(outf,"DEFINITION  %s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1628    /* fprintf(outf,"ACCESSION   %s\n", accnum); */
1629       fprintf(outf,"ORIGIN      \n");
1630       spacer = 11;
1631       numleft = true;
1632       numwidth = 8;  /* dgg. 1Feb93, patch for GDE fail to read short numwidth */
1633       strcpy(endstr, "\n//");
1634       linesout += 4;
1635       break;
1636
1637     case kPIR:
1638       /* somewhat like genbank... \\\*/
1639       /* fprintf(outf,"\\\\\\\n"); << only at top of file, not each entry... */
1640       fprintf(outf,"ENTRY           %s \n", idword);
1641       fprintf(outf,"TITLE           %s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1642    /* fprintf(outf,"ACCESSION       %s\n", accnum); */
1643       fprintf(outf,"SEQUENCE        \n");
1644       numwidth = 7;
1645       width= 30;
1646       spacer = kSpaceAll;
1647       numleft = true;
1648       strcpy(endstr, "\n///");
1649       /* run a top number line for PIR */
1650       for (j=0; j<numwidth; j++) fputc(' ',outf);
1651       for (j= 5; j<=width; j += 5) fprintf(outf,"%10d",j);
1652       fputc('\n',outf);
1653       linesout += 5;
1654       break;
1655
1656     case kNBRF:
1657       if (getseqtype(seq, seqlen) == kAmino)
1658         fprintf(outf,">P1;%s\n", idword);
1659       else
1660         fprintf(outf,">DL;%s\n", idword);
1661       fprintf(outf,"%s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1662       spacer = 11;
1663       strcpy(endstr,"*\n");
1664       linesout += 3;
1665       break;
1666
1667     case kEMBL:
1668       fprintf(outf,"ID   %s\n", idword);
1669   /*  fprintf(outf,"AC   %s\n", accnum); */
1670       fprintf(outf,"DE   %s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1671       fprintf(outf,"SQ             %d BP\n", seqlen);
1672       strcpy(endstr, "\n//"); /* 11Oct90: bug fix*/
1673       tab = 4;     /** added 31jan91 */
1674       spacer = 11; /** added 31jan91 */
1675       width = 60;
1676       linesout += 4;
1677       break;
1678
1679     case kGCG:
1680       fprintf(outf,"%s\n", seqname);
1681    /* fprintf(outf,"ACCESSION   %s\n", accnum); */
1682       fprintf(outf,"    %s  Length: %d  (today)  Check: %d  ..\n", idword, seqlen, checksum);
1683       spacer = 11;
1684       numleft = true;
1685       strcpy(endstr, "\n");  /* this is insurance to help prevent misreads at eof */
1686       linesout += 3;
1687       break;
1688
1689     case kStrider: /* ?? map ?*/
1690       fprintf(outf,"; ### from DNA Strider ;-)\n");
1691       fprintf(outf,"; DNA sequence  %s, %d bases, %X checksum.\n;\n", seqname, seqlen, checksum);
1692       strcpy(endstr, "\n//");
1693       linesout += 3;
1694       break;
1695
1696     case kFitch:
1697       fprintf(outf,"%s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1698       spacer = 4;
1699       width = 60;
1700       linesout += 1;
1701       break;
1702
1703     case kPhylip2:
1704     case kPhylip4:
1705       /* this is version 3.2/3.4 -- simplest way to write
1706         version 3.3 is to write as version 3.2, then
1707         re-read file and interleave the species lines */
1708       if (strlen(idword)>10) idword[10] = 0;
1709       fprintf(outf,"%-10s  ",idword);
1710       l1  = -1;
1711       tab = 12;
1712       spacer = 11;
1713       break;
1714
1715     case kASN1:
1716       seqtype= getseqtype(seq, seqlen);
1717       switch (seqtype) {
1718         case kDNA     : cp= "dna"; break;
1719         case kRNA     : cp= "rna"; break;
1720         case kNucleic : cp= "na"; break;
1721         case kAmino   : cp= "aa"; break;
1722         case kOtherSeq: cp= "not-set"; break;
1723         }
1724       fprintf(outf,"  seq {\n");
1725       fprintf(outf,"    id { local id %d },\n", gPretty.atseq);
1726       fprintf(outf,"    descr { title \"%s\" },\n", seqid);
1727       fprintf(outf,"    inst {\n");
1728       fprintf(outf,"      repr raw, mol %s, length %d, topology linear,\n", cp, seqlen);
1729       fprintf(outf,"      seq-data\n");
1730       if (seqtype == kAmino)
1731         fprintf(outf,"        iupacaa \"");
1732       else
1733         fprintf(outf,"        iupacna \"");
1734       l1  = 17;
1735       spacer = 0;
1736       width  = 78;
1737       tab  = 0;
1738       strcpy(endstr,"\"\n      } } ,");
1739       linesout += 7;
1740       break;
1741
1742     case kPAUP:
1743       nameleft= true;
1744       namewidth = 9;
1745       spacer = 21;
1746       width  = 100;
1747       tab  = 0; /* 1; */
1748       /* strcpy(endstr,";\nend;"); << this is end of all seqs.. */
1749       /* do a header comment line for paup */
1750       fprintf(outf,"[Name: %-16s  Len:%6d  Check: %8X]\n", idword, seqlen, checksum);
1751       linesout += 1;
1752       break;
1753
1754     case kPretty:
1755       numline= gPretty.numline;
1756       baseonlynum= gPretty.baseonlynum;
1757       namewidth = gPretty.namewidth;
1758       numright = gPretty.numright;
1759       numleft = gPretty.numleft;
1760       nameright = gPretty.nameright;
1761       nameleft = gPretty.nameleft;
1762       spacer = gPretty.spacer + 1;
1763       width  = gPretty.seqwidth;
1764       tab  = gPretty.tab;
1765       /* also add rtf formatting w/ font, size, style */
1766       if (gPretty.nametop) {
1767         fprintf(outf,"Name: %-16s  Len:%6d  Check: %8X\n", idword, seqlen, checksum);
1768         linesout++;
1769         }
1770       break;
1771
1772     case kMSF:
1773       fprintf(outf," Name: %-16s Len:%6d  Check: %5d  Weight:  1.00\n",
1774                     idword, seqlen, checksum);
1775       linesout++;
1776       nameleft= true;
1777       namewidth= 15; /* need MAX namewidth here... */
1778       sprintf(nameform, "%%+%ds ",namewidth);
1779       spacer = 11;
1780       width  = 50;
1781       tab = 0; /* 1; */
1782       break;
1783
1784     case kIG:
1785       fprintf(outf,";%s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1786       fprintf(outf,"%s\n", idword);
1787       strcpy(endstr,"1"); /* == linear dna */
1788       linesout += 2;
1789       break;
1790
1791     default :
1792     case kZuker: /* don't attempt Zuker's ftn format */
1793     case kPearson:
1794       fprintf(outf,">%s, %d bases, %X checksum.\n", seqname, seqlen, checksum);
1795       linesout += 1;
1796       break;
1797     }
1798
1799   if (*nameform==0) sprintf(nameform, "%%%d.%ds ",namewidth,namewidth);
1800   if (numline) sprintf(numform, "%%%ds ",numwidth);
1801   else sprintf(numform, "%%%dd ",numwidth);
1802   strcpy( nocountsymbols, kNocountsymbols);
1803   if (baseonlynum) {
1804     if (strchr(nocountsymbols,gPretty.gapchar)==NULL) {
1805       strcat(nocountsymbols," ");
1806       nocountsymbols[strlen(nocountsymbols)-1]= gPretty.gapchar;
1807       }
1808     if (gPretty.domatch && (cp=strchr(nocountsymbols,gPretty.matchchar))!=NULL) {
1809       *cp= ' ';
1810       }
1811     }
1812
1813   if (numline) {
1814    *idword= 0;
1815    }
1816
1817   width = min(width,kMaxseqwidth);
1818   for (i=0, l=0, ibase = 1; i < seqlen; ) {
1819
1820     if (l1 < 0) l1 = 0;
1821     else if (l1 == 0) {
1822       if (nameleft) fprintf(outf, nameform, idword);
1823       if (numleft) { if (numline) fprintf(outf, numform, "");
1824                     else fprintf(outf, numform, ibase);}
1825       for (j=0; j<tab; j++) fputc(' ',outf);
1826       }
1827
1828     l1++;                 /* don't count spaces for width*/
1829     if (numline) {
1830       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1)) {
1831         if (numline==1) fputc(' ',outf);
1832         s[l++] = ' ';
1833         }
1834       if (l1 % 10 == 1 || l1 == width) {
1835         if (numline==1) fprintf(outf,"%-9d ",i+1);
1836         s[l++]= '|'; /* == put a number here */
1837         }
1838       else s[l++]= ' ';
1839       i++;
1840       }
1841
1842     else {
1843       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
1844         s[l++] = ' ';
1845       if (!baseonlynum) ibase++;
1846       else if (0==strchr(nocountsymbols,seq[i])) ibase++;
1847       s[l++] = seq[i++];
1848       }
1849
1850     if (l1 == width || i == seqlen) {
1851       if (outform==kPretty) for ( ; l1<width; l1++) {
1852         if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
1853           s[l++] = ' ';
1854         s[l++]=' '; /* pad w/ blanks */
1855         }
1856       s[l] = '\0';
1857       l = 0; l1 = 0;
1858
1859       if (numline) {
1860         if (numline==2) fprintf(outf,"%s",s); /* finish numberline ! and | */
1861         }
1862       else {
1863         if (i == seqlen) fprintf(outf,"%s%s",s,endstr);
1864         else fprintf(outf,"%s",s);
1865         if (numright || nameright) fputc(' ',outf);
1866         if (numright)  fprintf(outf,numform, ibase-1);
1867         if (nameright) fprintf(outf, nameform,idword);
1868         }
1869       fputc('\n',outf);
1870       linesout++;
1871       }
1872     }
1873   return linesout;
1874 }  /*writeSeq*/
1875
1876
1877
1878 /* End file: ureadseq.c */