sources/seg/seg.c

   1
   2 /*****************************************************************************/
   3 /***  (seg.c)                                                              ***/
   4 /*** #include precomputed ln(fact(n)) array in "lnfac.h"                   ***/
   5 /*****************************************************************************/
   6
   7 /*--------------------------------------------------------------(includes)---*/
   8
   9 #include "genwin.h"
  10 #include "lnfac.h"
  11
  12 /*---------------------------------------------------------------(defines)---*/
  13
  14 #define LENCORRLIM 120
  15 #define MIN(a,b)        ((a) <= (b) ? (a) : (b))
  16
  17 /*---------------------------------------------------------------(structs)---*/
  18
  19 struct Segment
  20   {
  21    int begin;
  22    int end;
  23    struct Segment *next;
  24   };
  25
  26 /*---------------------------------------------------------------(globals)---*/
  27
  28 int window = 12;
  29 int downset, upset;
  30 double locut = 2.2;
  31 double hicut = 2.5;
  32
  33 int hilenmin = 0;
  34 int overlaps = FALSE;
  35 int hionly = FALSE;
  36 int loonly = FALSE;
  37 int entinfo = TRUE;
  38 int singleseq = FALSE;
  39 int prettyseq = FALSE;
  40 int prettytree = TRUE;
  41 int charline = 60;
  42 int maxtrim = 100;
  43
  44 double getprob(), lnperm(), lnass();
  45
  46 /*------------------------------------------------------------------(main)---*/
  47
  48 main(argc, argv)
  49   int argc;
  50   char *argv[];
  51
  52   {struct Database *db;
  53    struct Sequence *seq;
  54    struct Segment *segs;
  55    int ctime;
  56
  57    genwininit();
  58 /* readlencorr(); */                        /* #include lencorr file */
  59    getparams(argc, argv);
  60    entropy_init(window);
  61
  62    if ((db=opendbase(argv[1]))==NULL)
  63      {
  64       fprintf(stderr, "Error opening file %s\n", argv[1]);
  65       exit(1);
  66      }
  67
  68    for (seq=firstseq(db); seq!=NULL; seq=nextseq(db))
  69      {
  70       segs = (struct Segment *) NULL;
  71       segseq(seq, &segs, 0);
  72       mergesegs(seq, segs);
  73
  74       if (singleseq) singreport(seq, segs);
  75       else if (prettyseq) prettyreport(seq, segs);
  76       else if (prettytree) pretreereport(seq, segs);
  77       else report(seq, segs);
  78
  79       freesegs(segs);
  80       closeseq(seq);
  81      }
  82
  83    closedbase(db);
  84    exit(0);
  85   }
  86
  87 /*-------------------------------------------------------------(getparams)---*/
  88
  89 getparams(argc, argv)
  90   int argc;
  91   char *argv[];
  92
  93   {int i;
  94    int nargc;
  95    char **nargv;
  96    extern char *optarg;
  97    extern int optind;
  98    int c;
  99
 100    if (argc<2)
 101      {
 102       usage();
 103       exit(1);
 104      }
 105
 106    for (i=2; argc>i && argv[i][0]!='-'; i++)
 107      {
 108       if (i==2)
 109         {
 110          window = atoi(argv[2]);
 111         }
 112       else if (i==3)
 113         {
 114          locut = atof(argv[3]);
 115          hicut = locut + 0.3;
 116         }
 117       else if (i==4)
 118         {
 119          hicut = atof(argv[4]);
 120         }
 121      }
 122
 123    if (locut>hicut)
 124      {
 125       fprintf(stderr, "Warning: trigger complexity greater than extension\n");
 126       hicut = locut;
 127      }
 128
 129    downset = (window+1)/2 - 1;
 130    upset = window - downset;
 131
 132    if (i==argc) return;
 133
 134    nargc = argc-i+1;
 135    nargv = argv+(i-1);
 136    while ((c=getopt(nargc, nargv, "m:olhaxpqc:nt:"))!=-1)
 137      {
 138       switch(c)
 139         {
 140          case 'm':
 141             hilenmin = atoi(optarg);
 142             break;
 143          case 'o':
 144             overlaps = TRUE;
 145             hilenmin = 0;
 146             break;
 147          case 'l':
 148             loonly = TRUE;
 149             singleseq = FALSE;
 150             prettyseq = FALSE;
 151             prettytree = FALSE;
 152             break;
 153          case 'h':
 154             hionly = TRUE;
 155             singleseq = FALSE;
 156             prettyseq = FALSE;
 157             prettytree = FALSE;
 158             break;
 159          case 'a':
 160             hionly = FALSE;
 161             loonly = FALSE;
 162             singleseq = FALSE;
 163             prettyseq = FALSE;
 164             prettytree = FALSE;
 165             break;
 166          case 'x':
 167             singleseq = TRUE;
 168             prettyseq = FALSE;
 169             prettytree = FALSE;
 170             hionly = TRUE;
 171             loonly = FALSE;
 172             break;
 173          case 'p':
 174             prettytree = TRUE;
 175             prettyseq = FALSE;
 176             singleseq = FALSE;
 177             hionly = FALSE;
 178             loonly = FALSE;
 179             break;
 180          case 'q':
 181             prettyseq = TRUE;
 182             prettytree = FALSE;
 183             singleseq = FALSE;
 184             hionly = FALSE;
 185             loonly = FALSE;
 186             break;
 187          case 'c':
 188             charline = atoi(optarg);
 189             break;
 190          case 'n':
 191             entinfo = FALSE;
 192             break;
 193          case 't':
 194             maxtrim = atoi(optarg);
 195             break;
 196          case '?':
 197             fprintf(stderr, "Unknown option.\n");
 198             usage();
 199             exit(1);
 200             break;
 201         }
 202      }
 203
 204    return;
 205   }
 206
 207 /*---------------------------------------------------------------(segment)---*/
 208
 209 segseq(seq, segs, offset)
 210   struct Sequence *seq;
 211   struct Segment **segs;
 212   int offset;
 213
 214   {struct Segment *seg, *leftsegs;
 215    struct Sequence *leftseq;
 216    int first, last, lowlim;
 217    int loi, hii, i;
 218    int leftend, rightend, lend, rend;
 219    double *H, *seqent();
 220
 221    H = seqent(seq);
 222    if (H==NULL) return;
 223
 224    first = downset;
 225    last = seq->length - upset;
 226    lowlim = first;
 227
 228    for (i=first; i<=last; i++)
 229      {
 230       if (H[i]<=locut && H[i]!=-1)
 231         {
 232          loi = findlo(i, lowlim, H);
 233          hii = findhi(i, last, H);
 234
 235          leftend = loi - downset;
 236          rightend = hii + upset - 1;
 237
 238          trim(openwin(seq, leftend, rightend-leftend+1), &leftend, &rightend);
 239
 240          if (i+upset-1<leftend)   /* check for trigger window in left trim */
 241            {
 242             lend = loi - downset;
 243             rend = leftend - 1;
 244
 245             leftseq = openwin(seq, lend, rend-lend+1);
 246             leftsegs = (struct Segment *) NULL;
 247             segseq(leftseq, &leftsegs, offset+lend);
 248             if (leftsegs!=NULL)
 249               {
 250                if (*segs==NULL) *segs = leftsegs;
 251                else appendseg(*segs, leftsegs);
 252               }
 253             closewin(leftseq);
 254
 255 /*          trim(openwin(seq, lend, rend-lend+1), &lend, &rend);
 256             seg = (struct Segment *) malloc(sizeof(struct Segment));
 257             seg->begin = lend;
 258             seg->end = rend;
 259             seg->next = (struct Segment *) NULL;
 260             if (segs==NULL) segs = seg;
 261             else appendseg(segs, seg);  */
 262            }
 263
 264          seg = (struct Segment *) malloc(sizeof(struct Segment));
 265          seg->begin = leftend + offset;
 266          seg->end = rightend + offset;
 267          seg->next = (struct Segment *) NULL;
 268
 269          if (*segs==NULL) *segs = seg;
 270          else appendseg(*segs, seg);
 271
 272          i = min(hii, rightend+downset);
 273          lowlim = i + 1;
 274 /*       i = hii;     this ignores the trimmed residues... */
 275         }
 276      }
 277
 278    free(H);
 279    return;
 280   }
 281
 282 /*----------------------------------------------------------------(seqent)---*/
 283
 284 double *seqent(seq)
 285   struct Sequence *seq;
 286
 287   {struct Sequence *win;
 288    double *H;
 289    int i, first, last;
 290
 291    if (window>seq->length)
 292      {
 293       return((double *) NULL);
 294      }
 295
 296    H = (double *) malloc(seq->length*sizeof(double));
 297
 298    for (i=0; i<seq->length; i++)
 299      {
 300       H[i] = -1.;
 301      }
 302
 303    win = openwin(seq, 0, window);
 304    enton(win);
 305
 306    first = downset;
 307    last = seq->length - upset;
 308
 309    for (i=first; i<=last; i++)
 310      {
 311       if (seq->punctuation && hasdash(win))
 312         {H[i] = -1;
 313          shiftwin1(win);
 314          continue;}
 315       H[i] = win->entropy;
 316       shiftwin1(win);
 317      }
 318
 319    closewin(win);
 320    return(H);
 321   }
 322
 323 /*---------------------------------------------------------------(hasdash)---*/
 324
 325 hasdash(win)
 326   struct Sequence *win;
 327 {
 328         register char   *seq, *seqmax;
 329
 330         seq = win->seq;
 331         seqmax = seq + win->length;
 332
 333         while (seq < seqmax) {
 334                 if (*seq++ == '-')
 335                         return TRUE;
 336         }
 337         return FALSE;
 338 }
 339
 340 /*----------------------------------------------------------------(findlo)---*/
 341
 342 findlo(i, limit, H)
 343   int i, limit;
 344   double *H;
 345
 346   {int j;
 347
 348    for (j=i; j>=limit; j--)
 349      {
 350       if (H[j]==-1) break;
 351       if (H[j]>hicut) break;
 352      }
 353
 354    return(j+1);
 355   }
 356
 357 /*----------------------------------------------------------------(findhi)---*/
 358
 359 findhi(i, limit, H)
 360   int i, limit;
 361   double *H;
 362
 363   {int j;
 364
 365    for (j=i; j<=limit; j++)
 366      {
 367       if (H[j]==-1) break;
 368       if (H[j]>hicut) break;
 369      }
 370
 371    return(j-1);
 372   }
 373
 374 /*------------------------------------------------------------------(trim)---*/
 375
 376 trim(seq, leftend, rightend)
 377   struct Sequence *seq;
 378   int *leftend, *rightend;
 379
 380   {struct Sequence *win;
 381    double prob, minprob;
 382    int shift, len, i;
 383    int lend, rend;
 384    int minlen;
 385
 386 /* fprintf(stderr, "%d %d\n", *leftend, *rightend);  */
 387
 388    lend = 0;
 389    rend = seq->length - 1;
 390    minlen = 1;
 391    if ((seq->length-maxtrim)>minlen) minlen = seq->length-maxtrim;
 392
 393    minprob = 1.;
 394    for (len=seq->length; len>minlen; len--)
 395      {
 396       win = openwin(seq, 0, len);
 397       i = 0;
 398
 399       shift = TRUE;
 400       while (shift)
 401         {
 402          prob = getprob(win->state, len);
 403          if (prob<minprob)
 404            {
 405             minprob = prob;
 406             lend = i;
 407             rend = len + i - 1;
 408            }
 409          shift = shiftwin1(win);
 410          i++;
 411         }
 412       closewin(win);
 413      }
 414
 415 /* fprintf(stderr, "%d-%d ", *leftend, *rightend);  */
 416
 417    *leftend = *leftend + lend;
 418    *rightend = *rightend - (seq->length - rend - 1);
 419
 420 /* fprintf(stderr, "%d-%d\n", *leftend, *rightend);  */
 421
 422    closewin(seq);
 423    return;
 424   }
 425
 426 /*---------------------------------------------------------------(getprob)---*/
 427
 428 double getprob(sv, total)
 429   int *sv;
 430   int total;
 431
 432   {double ans, totseq;
 433
 434 #define LN20    2.9957322735539909
 435    totseq = ((double) total) * LN20;
 436
 437    ans = lnass(sv) + lnperm(sv, total) - totseq;
 438
 439    return(ans);
 440   }
 441
 442 /*----------------------------------------------------------------(lnperm)---*/
 443
 444 double lnperm(sv, tot)
 445   int *sv;
 446    int tot;
 447
 448   {double ans;
 449    int i;
 450
 451    ans = lnfac[tot];
 452
 453    for (i=0; sv[i]!=0; i++)
 454      {
 455       ans -= lnfac[sv[i]];
 456      }
 457
 458    return(ans);
 459   }
 460
 461 /*-----------------------------------------------------------------(lnass)---*/
 462
 463 double lnass(sv)
 464         register int    *sv;
 465 {
 466         double  ans;
 467         register int    svi, svim1;
 468         register int    class, total;
 469         register int    i;
 470
 471         ans = lnfac[20];
 472         if (sv[0] == 0)
 473                 return ans;
 474
 475         total = 20;
 476         class = 1;
 477         svim1 = sv[0];
 478         for (i=0;; svim1 = svi) {
 479                 if (++i==20) {
 480                         ans -= lnfac[class];
 481                         break;
 482                       }
 483                 else if ((svi = *++sv) == svim1) {
 484                         class++;
 485                         continue;
 486                 }
 487                 else {
 488                         total -= class;
 489                         ans -= lnfac[class];
 490                         if (svi == 0) {
 491                                 ans -= lnfac[total];
 492                                 break;
 493                         }
 494                         else {
 495                                 class = 1;
 496                                 continue;
 497                         }
 498                 }
 499         }
 500
 501         return ans;
 502 }
 503
 504 /*-------------------------------------------------------------(mergesegs)---*/
 505
 506 mergesegs(seq, segs)
 507   struct Sequence *seq;
 508   struct Segment *segs;
 509
 510   {struct Segment *seg, *nextseg;
 511    int len;
 512
 513    if (overlaps) return;
 514    if (segs==NULL) return;
 515
 516    if (segs->begin<hilenmin) segs->begin = 0;
 517
 518    seg = segs;
 519    nextseg = seg->next;
 520
 521    while (nextseg!=NULL)
 522      {
 523       if (seg->end>=nextseg->begin)               /* overlapping segments */
 524         {
 525          seg->end = nextseg->end;
 526          seg->next = nextseg->next;
 527          free(nextseg);
 528          nextseg = seg->next;
 529          continue;
 530         }
 531       len = nextseg->begin - seg->end - 1;
 532       if (len<hilenmin)                            /* short hient segment */
 533         {
 534          seg->end = nextseg->end;
 535          seg->next = nextseg->next;
 536          free(nextseg);
 537          nextseg = seg->next;
 538          continue;
 539         }
 540       seg = nextseg;
 541       nextseg = seg->next;
 542      }
 543
 544    len = seq->length - seg->end - 1;
 545    if (len<hilenmin) seg->end = seq->length - 1;
 546
 547    return;
 548   }
 549
 550 /*----------------------------------------------------------------(report)---*/
 551
 552 report(seq, segs)
 553   struct Sequence *seq;
 554   struct Segment *segs;
 555
 556   {struct Sequence *subseq;
 557    struct Segment *seg, *nextseg;
 558    static int hi = 1;
 559    static int lo = 0;
 560
 561    if (segs==NULL)
 562      {
 563       enton(seq);
 564       seqout(seq, hi, 1, seq->length);
 565 /*    fputc('\n', stdout);   -- for spacing after each sequence */
 566       return;
 567      }
 568
 569    if (segs->begin>0)
 570      {
 571       subseq = openwin(seq, 0, segs->begin);
 572       enton(subseq);
 573       seqout(subseq, hi, 1, segs->begin);
 574       closewin(subseq);
 575      }
 576
 577    for (seg=segs; seg!=NULL; seg=seg->next)
 578      {
 579       subseq = openwin(seq, seg->begin, seg->end-seg->begin+1);
 580       enton(subseq);
 581       seqout(subseq, lo, seg->begin+1, seg->end+1);
 582       closewin(subseq);
 583
 584       if (seg->next==NULL)
 585         {
 586          break;
 587         }
 588
 589       nextseg = seg->next;
 590
 591       if (nextseg->begin<=seg->end)
 592         {
 593          fprintf(stderr, "Overlapping segments: %s\n", seq->id);
 594          continue;
 595         }
 596
 597       if (nextseg->begin==seg->end+1)
 598         {
 599          continue;
 600         }
 601
 602       subseq = openwin(seq, seg->end+1, nextseg->begin-seg->end-1);
 603       enton(subseq);
 604       seqout(subseq, hi, seg->end+2, nextseg->begin);
 605       closewin(subseq);
 606      }
 607
 608    if (seg->end+1==seq->length)
 609      {
 610 /*    fputc('\n', stdout);   -- for spacing after each sequence */
 611       return;
 612      }
 613
 614    subseq = openwin(seq, seg->end+1, seq->length-seg->end-1);
 615    enton(subseq);
 616    seqout(subseq, hi, seg->end+2, seq->length);
 617    closewin(subseq);
 618
 619 /* fputc('\n', stdout);   -- for spacing after each sequence */
 620    return;
 621   }
 622
 623 /*------------------------------------------------------------(singreport)---*/
 624
 625 singreport(seq, segs)
 626         struct Sequence *seq;
 627         struct Segment  *segs;
 628 {
 629         char    *proseq, *proseqmax;
 630         struct Segment  *seg;
 631         int     begin, end, i, ctr;
 632
 633         proseq = seq->seq;
 634         proseqmax = proseq + seq->length;
 635         upper(proseq, seq->length);
 636
 637         for (seg=segs; seg!=NULL; seg=seg->next) {
 638                 begin = seg->begin;
 639                 end = seg->end;
 640                 memset(proseq + begin, 'x', end - begin +1);
 641         }
 642
 643         fprintf(stdout, "%s\n", seq->header);
 644
 645         for (i=0, ctr=0; proseq < proseqmax; ++i, ++ctr, ++proseq) {
 646                 if (ctr==charline) {
 647                         putc('\n', stdout);
 648                         ctr = 0;
 649                 }
 650                 putc(*proseq, stdout);
 651         }
 652
 653         putc('\n', stdout);
 654         if (putc('\n', stdout) == EOF) {
 655                 fprintf(stderr, "premature EOF on write\n");
 656                 exit(2);
 657         }
 658 }
 659
 660 /*----------------------------------------------------------(prettyreport)---*/
 661
 662 prettyreport(seq, segs)
 663   struct Sequence *seq;
 664   struct Segment *segs;
 665
 666 {
 667         char    *proseq, *proseqmax;
 668         char    format[10];
 669         struct Segment  *seg;
 670         int     begin, end, i, ctr;
 671         int     leftspace;
 672
 673         leftspace = (int) ceil(log10((double)seq->length));
 674         sprintf(format, "%%%dd ", leftspace);
 675
 676         upper(proseq = seq->seq, seq->length);
 677
 678         for (seg=segs; seg!=NULL; seg=seg->next) {
 679                 begin = seg->begin;
 680                 end = seg->end;
 681                 lower(proseq + begin, end - begin + 1);
 682         }
 683
 684         fprintf(stdout, "%s\n", seq->header);
 685
 686         space(leftspace+1);
 687         for (i=0, ctr=1; i<charline; ++i, ++ctr) {
 688                 if (ctr==10) {
 689                         putc('.', stdout);
 690                         ctr = 0;
 691                 }
 692                 else
 693                         putc(' ', stdout);
 694         }
 695         putc('\n', stdout);
 696         fprintf(stdout, format, 1);
 697
 698         proseqmax = proseq + seq->length;
 699         for (i=0, ctr=0; proseq < proseqmax; ++i, ++ctr, ++proseq) {
 700                 if (ctr==charline) {
 701                         putc('\n', stdout);
 702                         fprintf(stdout, format, i+1);
 703                         ctr = 0;
 704                 }
 705                 putc(*proseq, stdout);
 706         }
 707
 708         fprintf(stdout, " %d\n", seq->length);
 709         if (putc('\n', stdout) == EOF) {
 710                 fprintf(stderr, "premature EOF on write\n");
 711                 exit(2);
 712         }
 713 }
 714
 715 /*---------------------------------------------------------(pretreereport)---*/
 716
 717 pretreereport(seq, segs)
 718   struct Sequence *seq;
 719   struct Segment *segs;
 720
 721 {
 722         struct Sequence *win;
 723         struct Segment  *seg;
 724         char    buffer[100], leftfmt[20], rightfmt[20];
 725         char    *curseq;
 726         int     i, left, right, len;
 727         int     current, nextloent;
 728         int     cline;
 729
 730         cline = charline / 2;
 731
 732         fprintf(stdout, "%s\n\n", seq->header);
 733         sprintf(leftfmt, "%%%ds", cline);
 734         sprintf(rightfmt, "%%-%ds", cline);
 735
 736         current = 0;
 737
 738         for (seg=segs; ; seg=seg->next) {
 739                 if (seg==NULL)
 740                         nextloent = seq->length;
 741                 else
 742                         nextloent = seg->begin;
 743
 744                 if (current < nextloent) {
 745                         left = current;
 746                         right = nextloent - 1;
 747                         len = right - left + 1;
 748                         win = openwin(seq, left, len);
 749                         upper(curseq = win->seq, win->length);
 750
 751                         space(cline);
 752                         fprintf(stdout, " %4d-%-4d ", left+1, right+1);
 753
 754                         while (len>0) {
 755                                 if (len<=cline) {
 756                                         fwrite(curseq, len, 1, stdout);
 757                                         putc('\n', stdout);
 758                                         break;
 759                                 }
 760                                 else {
 761                                         fwrite(curseq, cline, 1, stdout);
 762                                         putc('\n', stdout);
 763                                         space(cline+11);
 764                                         curseq += cline;
 765                                         len -= cline;
 766                                 }
 767                         }
 768                         closewin(win);
 769                 }
 770
 771                 if (seg==NULL) break;
 772
 773       left = seg->begin;
 774       right = seg->end;
 775       len = right - left + 1;
 776       win = openwin(seq, left, len);
 777       lower(curseq = win->seq, win->length);
 778
 779                 i = MIN(cline, len);
 780                 if (i < cline)
 781                         fprintf(stdout, "%*s", cline-i, "");
 782                 fwrite(curseq, i, 1, stdout);
 783                 fprintf(stdout, " %4d-%-4d ", left+1, right+1);
 784                 putc('\n', stdout);
 785
 786                 len -= cline;
 787                 if (len>0)
 788                         curseq += cline;
 789
 790                 while (len>0) {
 791                         i = MIN(cline, len);
 792                         if (i < cline)
 793                                 space(cline - i);
 794                         fwrite(curseq, i, 1, stdout);
 795                         putc('\n', stdout);
 796                         len -= i;
 797                         if (len>0)
 798                                 curseq += i;
 799                 }
 800
 801                 closewin(win);
 802                 current = right + 1;
 803         }
 804
 805         putc('\n', stdout);
 806 }
 807
 808 /*-----------------------------------------------------------------(space)---*/
 809
 810 space(len)
 811         register int len;
 812 {
 813         static char     spaces[] =
 814                 {' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',
 815                 ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',
 816                 ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',
 817                 ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',
 818                 ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '
 819                 };
 820         register int i;
 821
 822         while (len > 0) {
 823                 i = MIN(len, sizeof(spaces)/sizeof(spaces[0]));
 824                 fwrite(spaces, i, 1, stdout);
 825                 len -= i;
 826         }
 827 }
 828
 829 /*----------------------------------------------------------------(seqout)---*/
 830
 831 seqout(seq, hilo, begin, end)
 832   struct Sequence *seq;
 833   int hilo;
 834   int begin, end;
 835
 836 #define HDRLEN 60
 837
 838 {
 839         char    *proseq, *proseqmax, *id, *header;
 840    char outbuf[HDRLEN+1];
 841    static int hi = 1;
 842    static int lo = 0;
 843    int i, ctr, iend;
 844
 845    if (hionly && hilo==lo) return;
 846    if (loonly && hilo==hi) return;
 847
 848    proseq = seq->seq;
 849    proseqmax = proseq + seq->length;
 850    id = seq->id;
 851    if (id==NULL) id = seq->parent->id;
 852    header = seq->header;
 853    if (header==NULL) header = seq->parent->header;
 854
 855    iend = findchar(header, ' ');
 856    if (iend!=-1) header = header+iend;
 857
 858    if (entinfo)
 859      {
 860       fprintf(stdout, ">%s(%d-%d)", id, begin, end);
 861 /*    if (iend!=-1 && strlen(header)<=HDRLEN) fprintf(stdout, "%s", header);
 862       else if (iend!=-1) for (i=0; i<HDRLEN; i++) putc(header[i], stdout); */
 863       fprintf(stdout, " complexity=%4.2f (%d/%4.2f/%4.2f)\n",
 864            seq->entropy, window, locut, hicut);
 865      }
 866    else
 867      {
 868       fprintf(stdout, ">%s(%d-%d)", id, begin, end);
 869       if (iend!=-1)   /* fprintf(stdout, "%s\n", header); */
 870         {
 871                  i = MIN(HDRLEN, strlen(header));
 872                  fwrite(header, i, 1, stdout);
 873          putc('\n', stdout);
 874         }
 875       else putc('\n', stdout);
 876      }
 877
 878    if (hilo==lo)
 879      {
 880       lower(proseq, seq->length);
 881      }
 882    else if (hilo==hi && seq->length>=hilenmin)
 883      {
 884       upper(proseq, seq->length);
 885      }
 886    else
 887      {
 888       lower(proseq, seq->length);
 889      }
 890
 891    for (; proseq < proseqmax; proseq+=i) {
 892                 i = MIN(charline, proseqmax - proseq);
 893                 fwrite(proseq, i, 1, stdout);
 894                 putc('\n', stdout);
 895         }
 896
 897         if (putc('\n', stdout) == EOF) {
 898                 fprintf(stderr, "premature EOF on write\n");
 899                 exit(2);
 900         }
 901 }
 902
 903 /*-------------------------------------------------------------(appendseg)---*/
 904
 905 appendseg(segs, seg)
 906   struct Segment *segs, *seg;
 907
 908   {struct Segment *temp;
 909
 910    temp = segs;
 911    while (1)
 912      {
 913       if (temp->next==NULL)
 914         {
 915          temp->next = seg;
 916          break;
 917         }
 918       else
 919         {
 920          temp = temp->next;
 921         }
 922      }
 923
 924    return;
 925   }
 926
 927 /*--------------------------------------------------------------(freesegs)---*/
 928
 929 freesegs(segs)
 930   struct Segment *segs;
 931
 932   {struct Segment *temp;
 933
 934    while (segs!=NULL)
 935      {
 936       temp = segs->next;
 937       free(segs);
 938       segs = temp;
 939      }
 940   }
 941
 942 /*-----------------------------------------------------------------(usage)---*/
 943
 944 usage()
 945
 946   {
 947    fprintf(stderr, "\
 948 Usage: seg <file> <window> <locut> <hicut> <options>\n\
 949          <file>   - filename containing fasta-formatted sequence(s) \n\
 950          <window> - OPTIONAL window size (default 12) \n\
 951          <locut>  - OPTIONAL low (trigger) complexity (default 2.2) \n\
 952          <hicut>  - OPTIONAL high (extension) complexity (default locut + 0.3) \n\
 953          <options> \n\
 954             -x  each input sequence is represented by a single output \n\
 955                 sequence with low-complexity regions replaced by \n\
 956                 strings of 'x' characters \n\
 957             -c <chars> number of sequence characters/line (default 60)\n\
 958             -m <size> minimum length for a high-complexity segment \n\
 959                 (default 0).  Shorter segments are merged with adjacent \n\
 960                 low-complexity segments \n\
 961             -l  show only low-complexity segments (fasta format) \n\
 962             -h  show only high-complexity segments (fasta format) \n\
 963             -a  show all segments (fasta format) \n\
 964             -n  do not add complexity information to the header line \n\
 965             -o  show overlapping low-complexity segments (default merge) \n\
 966             -t <maxtrim> maximum trimming of raw segment (default 100) \n\
 967             -p  prettyprint each segmented sequence (tree format) \n\
 968             -q  prettyprint each segmented sequence (block format) \n");
 969    exit(1);
 970   }
 971