File: codons.c

package info (click to toggle)
codonw 1.4.4-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,060 kB
  • sloc: ansic: 11,227; makefile: 209; sh: 170; perl: 36
file content (1149 lines) | stat: -rwxr-xr-x 49,033 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
/**************************************************************************/
/* CodonW codon usage analysis package                                    */
/* Copyright (C) 2005            John F. Peden                            */
/* This program is free software; you can redistribute                    */
/* it and/or modify it under the terms of the GNU General Public License  */
/* as published by the Free Software Foundation; version 2 of the         */
/* License,                                                               */
/*                                                                        */
/* This program is distributed in the hope that it will be useful, but    */
/* WITHOUT ANY WARRANTY; without even the implied warranty of             */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the           */
/* GNU General Public License for more details.                           */
/* You should have received a copy of the GNU General Public License along*/
/* with this program; if not, write to the Free Software Foundation, Inc.,*/
/* 675 Mass Ave, Cambridge, MA 02139, USA.                                */
/*                                                                        */
/*                                                                        */
/* The author can be contacted by email (jfp#hanson-codonw@yahoo.com Anti-*/
/* Spam please change the # in my email to an _)                          */
/*                                                                        */
/* For the latest version and information see                             */
/* http://codonw.sourceforge.net 					  */
/**************************************************************************/
/*                                                                        */
/* -----------------------        Codons.C       ------------------------ */
/* This file contains main() function and drives CodonW.                  */
/*                                                                        */
/* External subroutines and functions                                     */
/* clearscr           screen clearing Macro defined in CodonW.h           */
/* proc_comm_line     process command line arguments                      */
/* initilize_point    assigns genetic code dependent parameters to structs*/
/* initilize_coa      selects the default codons to exclude from the      */
/*                    Correspondence Analysis                             */
/* main_menu          The interactive menu system                         */
/* clean_up           Re-zeros various internal counters and arrays       */
/* open_file          Open files, checks for existing files               */
/* fileclose          Closes files and returns a NULL pointer or exits    */
/* textbin            Converts codon usage to binary data file            */
/* dot(,X)            prints a period every X times it is called          */
/* PrepAFC            Prepare for the COA                                 */
/* DiagoRC            This routine generates the COA                      */
/* colmout            write the output from COA to file                   */
/* rowout             save as above except records the gene information   */
/* inertialig         analyse row inertia and records the results to file */
/* inertiacol         analyse column inertia and record the results       */
/* suprow             add supplementary genes into COA                    */
/* get_aa             converts a three base codon into a 1 or 3 letter AA */
/* codon_error        Called after all codons read, checks data was OK    */
/* rscu_usage_out     Write out RSCU                                      */
/* codon_usage_out    Write out Codon Usage                               */
/* raau_usage_out     Write out normalised amino acid usage               */
/* dinuc_count        Count the dinucleotide usage                        */
/* dinuc_out          Write out dinucleotide usage                        */
/* aa_usage_out       Write out amino acid usage                          */
/* gc_out             Writes various analyses of base usage               */
/* cutab_out          Write a nice tabulation of the RSCU+CU+AA           */
/* base_sil_us_out    Write out base composition at silent sites          */
/* cai_out            Write out CAI usage                                 */
/* cbi_out            Write out codon bias index                          */
/* fop_out            Write out Frequency of Optimal codons               */
/* enc_out            Write out Effective Number of codons                */
/* hydro_out          Write out Protein hydropathicity                    */
/* aromo_out          Write out Protein aromaticity                       */
/* coa_raw_out        Write out raw codon usage for use by COA analysis   */
/*                                                                        */
/*                                                                        */
/* Internal subroutines to Codon.c                                        */
/* my_exit            Controls exit from CodonW closes any open files     */
/* tidy               reads the input data                                */
/* output             called from tidy to decide what to do with the data */
/* toutput            handles the reformatting and translation of seqs    */
/* output_long        if sequence is very long then process what we know  */
/*                    and write sequence to disk in fragments             */
/* file_close         Closes open files                                   */
/* c_help             Generates help informatio                           */
/* WasHelpCalled      Checks strings to see if help was requested         */
/*                                                                        */
/**************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <string.h>
#include <errno.h>
#include <ctype.h>

#define  ORIG_DEFS    
      /* used to decide whether declarations are external or not          */
      /* Master Header file                                               */
#include "codonW.h"     
#undef   ORIG_DEFS


#if defined(__MWERKS__)
#include <console.h> 
#endif     

/**************************   MAIN   **************************************/
/* The main function processes commandline arguments to decide whether    */
/* CodonW is running in an interactive mode, if so then the menu is called*/
/* CodonW also had the less documented feature of imitating other useful  */
/* codon usage and sequence manipulation program.    If the program is    */
/* called by a recognised name (see proc_comm_line for a list) such as    */
/* rscu then pm->codons is false and it only performs the required tasks  */
/* bypassing the menu system.                                             */
/* Main then calls tidy() to read in the data files, and count codon usage*/
/* depending on the requested output options toutput calls various subrou */
/* tines. If COA has been requested it also calls these subroutines and   */
/* recording useful information to summary.coa.                           */
/**************************************************************************/

int main(int argc, char *argv[])
{
  FILE       *finput = NULL, *foutput = NULL, *fblkout = NULL;
  FILE       *fcoaout = NULL;
  FILE       *fsummary= NULL;
  int        num_seq  = 0;

  num_sequence     = 0;
  num_seq_int_stop = 0;
  valid_stops      = 0;
  last_aa          = 0;
 
#if defined(__MWERKS__) /* Macintosh code-warrior */
  argc=ccommand(&argv); 
#endif

  pm = &Z_menu;
  pm->totals = FALSE;
  pm->my_err = stderr;


  initilize_point(pm->code, pm->f_type, pm->c_type);
  initilize_coa(pm->code);

  proc_comm_line(&argc, &argv);


  
  /******************** main loop ****************************/

  do {
    if (pm->codonW) {  
                                  /* If the program   chosen is codons    */
      printf(" Welcome to CodonW %.*s for Help type h\n\n",
              (int) strlen(Revision) - 11, Revision +10 ); 
                                  /* Now Run the main menu interface      */      
      if (pm->menu) main_menu(0); 
    } 

          /* if users select human readable output they want nice tables  */
    if (pm->bulk == 'C' && pm->seq_format == 'H')  pm->bulk = 'O';
    if (pm->bulk == 'S' && pm->seq_format == 'H')  pm->bulk = 'O';

     pm->analysis_run = TRUE;       /* codons has started an analysis this*/
                                    /* parameter is checked by my_exit    */

    if (pm->inputfile != NULL)      /* rewind various input files in case */
      rewind(pm->inputfile);        /* this is a second analysis run      */
    if (pm->fopfile != NULL)
      rewind(pm->fopfile);
    if (pm->cbifile != NULL)
      rewind(pm->cbifile);  
    if (pm->caifile != NULL)
      rewind(pm->caifile);

    /* num_sequence                 number of sequences read              */
    /* num_seq_int_stop             number with internal stop codons      */
    /* valid_stops                  No.terminated with a stop codon       */
    /* tot                          total number of codons read           */

    num_sequence = num_seq_int_stop = valid_stops = tot = 0;

    clean_up(ncod, naa);            /*re-zero count of amino and codons   */
    finput = pm->inputfile;
    foutput = pm->outputfile;
    fblkout = pm->tidyoutfile;

    fileclose(&pm->fcoa_out);
    if (pm->coa) 
      if ((pm->fcoa_out = open_file("", "coa_raw", "w", FALSE)) == NULL)
	my_exit(1, "coa_raw");          /*controlled exit from CodonW     */
    fcoaout = pm->fcoa_out;

/*  Tidy                                                                  */
/*  reads input data, returns the number of sequences read in             */
/*  num_sequence is global so I don't really have to assign it here       */    
    num_sequence = tidy(finput, foutput, fblkout, fcoaout);

    fprintf(pm->my_err,"\n\n\t\tNumber of sequences: %i\n",
        num_sequence);

/* num_seq_int_stop  value is calculated in codon_usage_out               */ 
    if (num_seq_int_stop > 0 && pm->warn ) {   
      if (pm->totals && (num_seq_int_stop >= valid_stops ))
	   fprintf(pm->my_err, "\tWARNING\t At least one sequence in your"
                    " input file has\ninternal stop codons (found %i"
                    " internal stops) \tWARNING\n",num_seq_int_stop);
      else
	   fprintf(pm->my_err, "\tWARNING\t %i sequences had internal "
                           "stop codons \tWARNING\n",num_seq_int_stop);
    }
/* don't wait for a pause if no_menu has been set                          */
    if ( pm->codonW && pm->menu ) pause; 

    if ( pm->coa && pm->totals)          /* idiots error catch             */ 
      my_exit(99,"A COA analysis of concatenated sequences is nonsensical\n"
      "I have completed any other requests but not the COA");

/* if COA has been requested then open summary.coa and start the analysis  */
    if (pm->coa) {
     if (fsummary == NULL)
      if ((fsummary = open_file("", "summary.coa", "w", FALSE)) == NULL)
       my_exit(1, "summary.coa");
/* set the number of genes in the analysis to the number read in by tidy   */
     pcoa->rows = num_sequence;
     fileclose(&fcoaout);
/* if COA has been selected then during the reading in phase raw codon usag*/
/* will have been written to the file coa_raw                              */
/* text bin converts this to binary data for the COA analysis program      */
      textbin("coa_raw", "cbrawin");
      printf("Generating correspondence analysis\n"); 
      dot(0,10);
      
         
      fprintf(fsummary, "\t\tSummary of Correspondence Analysis \n\n"
                        "The input file was %s it contained %i genes\n"        
                        "The number of axes generated was %i\n"
                        "A COA was requested of %s%s usage\n\n\n"
                        "Most of the output presented in this file "
                        "has also been written to separate files\n"
                        "genes.coa\tThe position of the genes on the "
                        "first %i axis\n" 
                        "%s.coa\tThe position of the %i %s on the %i "
                        "principle axes\n\n\n",
                        pm->curr_infilename, 
                        pcoa->rows,
                        ((pcoa->rows<pcoa->colm)?pcoa->rows:pcoa->colm)-1,
                        (pm->coa == 'r')        ?"relative synonymous ":"", 
                        (pm->coa == 'a')        ?"amino acid" : "codon",
                        pcoa->axis,
	                    (pm->coa == 'a')        ?"amino" : "codon",
	                    pcoa->colm,
	                    (pm->coa == 'a')        ?"amino acids":"codons", 
	                    pcoa->axis);
/* allocate memory for the rows and columns, scale both, and write out the*/
/* resulting matrix to the file cbrawin                                   */

      PrepAFC("cbrawin");

/* Now do the analysis, calculate the data inertia and all the vectors    */

      DiagoRC(fsummary);

/* colmout records the position of the columns on each of the factors/axes*/

      if (pm->coa == 'a')
          colmout("cbfcco", "amino.coa", paa, fsummary);
      else
          colmout("cbfcco", "codon.coa", paa, fsummary);

/* rowout records the position of the genes on each of the axis           */

      rowout("cbfcli", "genes.coa", "coa_raw", fsummary);

/* pcoa->level == e for exhaustive analysis of inertia                    */  
      if (pcoa->level == 'e') {     

       fprintf(fsummary, "\n\n\nYou requested detailed output from the COA"
           "\n\nThe absolute and relative inertia "
           "of each gene and %s (see also inertia.coa)\n",
           (pm->coa == 'a') ? "amino acids" : "codons");
/* inertialig must preceed inertiacol, records inertia of genes to file   */
/* it opens the raw codon usage file and loads the raw data to memory     */
       inertialig("inertia.coa", "coa_raw" ,fsummary);
/* uses the preloaded raw codon usage, to calculate inertia and other data*/
/* such as contribution of each column to each factor and to the extent   */
/* each column is explained by each factor and what the residual variation*/
/* is                                                                     */
       inertiacol("inertia.coa", fsummary);
      }
      
/* if pcoa->add_row is real string, then it will be the name of the file  */
/* containing additional sequence data, that will be excluded from the COA*/
/* but factored in, using the original COA vectors and then all other     */
/* calculation can proceed as with the original data                      */
      if (strlen(pcoa->add_row)) {
          if ((finput = open_file("", pcoa->add_row, "r", FALSE)) 
              == NULL) my_exit(6, "add_row");
          if ((foutput = tmpfile()) == NULL)
              my_exit(1, "temp file foutput");
          if ((fblkout = tmpfile()) == NULL)
              my_exit(1, "temp file fblkout");

    if ((fcoaout = open_file("", "coa1_raw", "w", FALSE)) == NULL)
      my_exit(1, "coa1_raw");

    clean_up(ncod, naa);
    num_sequence =num_seq_int_stop=valid_stops=tot = 0;
/* load the additional data file and process as normal                    */
/* but don't calculate any indices or write the data to the normal output */
/* files, rather write them to tmp files which will be deleted at end of  */
/* program execution                                                      */
    num_seq = tidy(finput, foutput, fblkout, fcoaout);

/* close the files now we are finished                                    */
    fileclose(&fcoaout);
    fileclose(&foutput);
    fileclose(&fblkout);
    fileclose(&finput);

/* covert to binary, use additional raw data file, note not coa_raw this  */
    textbin("coa1_raw", "cb1raw");
/* now call the routine suprow and add these additional genes, we will    */
/* process this data for inertia and append the gene and col. coordinates */
/* to the original gene.coa and codon.coa (or amino.coa)                  */
    suprow(num_seq, "cbfcvp", "cb1raw", "genes.coa", "coa1_raw", fsummary);
    
/* close these files now that we have finished with them and the COA      */

    fileclose(&foutput);
    fileclose(&fblkout);
    fileclose(&fcoaout);
      }
    }
   printf("\n");
  } while (pm->codonW && pm->menu ); /* OK now we loop back to main_menu  */
/* though only if we are in interactive mode and running as CodonW        */
  my_exit(0,"");                     /* last call to my_exit              */
  return 0;                          /* dummy return to keep pedantic but */
                                     /* brain dead compilers happy        */
}

/**********************  END of MAIN()   **********************************/


/**********************  Subroutines     **********************************/
/* Tidy                                                                   */
/*  reads input data from a sequence file containing fasta like formatted */
/*  sequence discards numbers, but keeps other characters                 */
/*  Each sequence must begin with title line must start with > or ;       */
/*  any following descriptive lines must begin with ; or >.Sequence start */
/*  is the first alphabetic character on the line following the headers   */
/*  There is no limit to sequence length or number of sequences but       */
/*  input lines should be less than 200 char in width                     */
/**************************************************************************/

int tidy(FILE * finput, FILE * foutput, FILE * fblkout, FILE * fcoaout)
{
  char            seq[MAX_GENE + LINE_LENGTH + 1];
  char            in[LINE_LENGTH + 1];
  int             first_line = TRUE, ic = 0;
  int             ii = 0;
  int             i,x;
  long            ic_orig = 0;
/* while still able to read data from the input file keep reading         */
  while ((fgets(in, LINE_LENGTH, finput) != NULL)) {

/* idiot error check to see if the file looks like fasta or PIR format    */ 
    if (!num_sequence && in[0] != ';' && in[0] != '>') {
      fprintf(stderr, "\n Error input file not in a recognised format \n"
                      " you must convert it into FASTA/Pearson format" 
                      " EXITING\n");
      my_exit(99, "input file not in a recognised format:tidy");
    }

    if (in[0] == ';' || in[0] == '>') { /* if true them this is a header   */
      if (first_line) {                 /* if true this is the first header*/

	first_line = FALSE;                 /* will only be reset when reread  */
                                        /* the next sequence               */
	if (num_sequence) {                 /* wait till we have read the first*/
                                        /* before writing to disk          */
/* now if we are concatenating sequence data we need will handle it thus   */
        if (pm->totals) {	    

/* first if translating or reformatting the input file flush the read      */
/* data to the disk                                                        */

         if (strchr("RNT",(int)pm->bulk)!=NULL) output_long(fblkout, seq);
         if (tot) {
   /* if something we have sequence read in, then we need to process this  */
   /* check whether the last codon of the sequence was was a stop          */
            last_aa = codon_usage_tot(seq, tot);
    	    if (pcu->ca[last_aa] == 11) valid_stops++;
	     }
/* rather re-setting everything to zero, we will just blank the array seq  */
         tot = 0;  
	  } else {     
/* else matches if tot; if sequences are not being concatenated we call    */
/* output to decide what to do with all the read data                      */
/* then we blank all the data from memory and start again                  */
        output(seq, foutput, fblkout, fcoaout);
	    clean_up(ncod, naa);      
	  }                 
	}                                    /* matches if(num_sequence)       */

/* If we get here we have read a header line, this then needs to be proc'ed*/
/* first the header is tested to see does it contain spaces the string is  */
/* converted from the first non space character to the title array         */

 for (ii = 1; isspace( (int) in[ii]) && ii < (int) strlen(in); ii++)
     ;
 strncpy(title, in + ii, 99);

/* Titles are cleaned up by removing newline characters and the delimiting */
/* character p,->seperater and also null terminating the title string      */

 for (i = 0; i < (int) strlen(title); i++) {
	  if (title[i] == '\n')	
          title[i] = '\0';          /* chops new line off                  */
	  else if (title[i] == pm->seperator )
	    title[i] = '_';             /* removes the separator if present    */
	  else if (i == (int) (strlen(title) - 1))
	    title[i] = '\0';            /* if we have reached end of title     */
	}

/* if we are reformatting the data, we print a friendly dot just in-case   */
if (strchr("RNT", (int)pm->bulk) ==NULL || pm->totals)
    dot((int) num_sequence, 5);       
/* we have now finished processing our first header line and are reading   */
/* our sequence data                                                       */
num_sequence++; 
      }                            /* matches if first line                */
      continue;                    /* read another line ie. jump to while()*/
    }                              /* if (in[0] == ';' || in[0] == '>')    */
    else{                          /* this must be a line containing seq   */
	first_line = TRUE;             /* so reset the first_line variable     */
      }

/* at this point we have read in the header lines and have been or about to*/
/* process the input data, now we test how much we have read into the array*/
/* seq, tot is equivalent to the last element in the array                 */
/* if tot is greater than or equal to MAX_GENE then the array is quite full*/
/* luckily we made the array seq to be MAX_GENE plus LINE_LENGTH +1        */

    if (tot >= MAX_GENE) {         /* sequence is larger than seq          */
      master_ic += MAX_GENE;       /* now remember how many bases we are   */
      ic_orig = tot;               /* going to write to disk               */
                                   /* and what size the array was to start */

      if (strchr("RNT", (int) pm->bulk) != NULL)
	     output_long(fblkout, seq);/* flush to disk and then continue      */
      else if (pm->bulk == 'D')
       	dinuc_count(seq, tot);     /* then we had better count the dinucs  */

/* Debugging code in-case we are asking for something that we can't handle */
#ifdef DEBUG
      else if (strchr("OCASLDBX", (int) pm->bulk) != NULL) ; /* dummy      */ 
      else if (pm->bulk)
	fprintf(stderr, "ERROR-22 %c pm->bulk undefined\n", pm->bulk);

      if (pm->cai || pm->fop || pm->cbi || pm->enc || pm->gc ||
	  pm->gc3s || pm->sil_base || pm->bulk ||
	  pm->coa);
      else
	fprintf(stderr, "Programming error");
#endif


/* Now count first MAX_GENE bases, luckily MAX_GENE is always a multiple of*/
/* 3, we count the bases and amino acids in codon_usage_tot                */

      last_aa = codon_usage_tot(seq, MAX_GENE);

/* now we move all unprocessed/written/counted bases to the front of seq   */

      for (i = MAX_GENE, x = 0; i < ic_orig; i++, x++)
       seq[x] = seq[i];            /* i is pointing near the end of array  */
      tot = x;                     /* x the front of the array             */
    }                              /* Matches if (tot >= MAX_GENE)         */

    ic = 0;                        /* first base of the input file         */
    while (in[ic] != '\0') {       /* scan input line till we see a Null   */
      if (isalpha((int)in[ic])) ;  /* do nothing if a alpha                */
      else if (pm->bulk == 'R' && in[ic] == '-'); /* do nothing            */
      else if (in[ic] == '*' || in[ic] == '.') ;  /* do nothing            */

      else {
    ic++;                          /* is not one above skip to next letter */
    continue;                     
      }                            /* while( in[ic] != '\0')               */


   in[ic] = (char)toupper((int)in[ic]);/*   converts2capitals              */
   if (strrchr("CG", (int) in[ic]) != NULL)
	GC_TOT++;                          /* is it a G or C                   */
   else if (strrchr("ATU", (int) in[ic]) != NULL)
	AT_TOT++;                          /* is it an A or T                  */
   else if ( in[ic] == '-' )
	GAP_TOT++;                         /* is it a gap character            */
   else
	non_std_char++;                    /* then it isn't a standard base    */

   if (strrchr("ABCDEFGHIKLMNPQRSTVWYZX"
		  ,(int) in[ic]) != NULL)
	AA_TOT++;                          /* it might be an amino acid        */
      if (strrchr("MRWSYKVHDBXN" , (int) in[ic]) != NULL)
	IUBC_TOT++;                        /* it might be a IUBC code          */

  seq[tot] = in[ic];                       /* move base into seq array         */
  seq[tot + 1] = '\0';                     /* make sure array is null term'ed  */


 /* now we test that the first codon is a valid start codon                    */

  if ( tot == 0 && master_ic == 0 ) {
      	
      in[1] = (char)toupper((int)in[1]);  /* Uppercase the first codon         */
      in[2] = (char)toupper((int)in[2]);  

       if ( in[1] == 'T' && (in[0] == 'A' || in[2] == 'G' ))
	  valid_start=TRUE;                /* Yeup it could be a start codon   */
	else
	  valid_start=FALSE;               /* Nope it doesn't seem to be one   */
  }
      ic++;                            /* total No. of sequence bases read */
      tot++;                           /* total currently stored in memory */
    }
  }                                    /* reached end of input file        */

/* Idiot error catch, this file is empty, at least it looks empty to codonW*/

  if ( !num_sequence ) my_exit(99,"The input file was empty");

/* better make sure to write anything left in seq to disk before returning */

  output(seq, foutput, fblkout, fcoaout);
  return (int) num_sequence;
}

/************************  TOUTPUT       **********************************/
/* toutput                                                                */
/*                                                                        */
/* This subroutine is very similar to output_long, basically it reformats */
/* or translates sequences less than MAX_GENE in length as a single read  */
/* It writes in reader format "ACG ATT ATC" i.e writes the sequence in    */
/* codons. Because it works with output_long it needs to know whether     */
/* the sequence being written to disk is a fragment or a complete gene    */
/**************************************************************************/
int toutput(FILE * fblkout, char *seq) {
  long int        ic = 0;
  int             space = 3;
  char            codon[4];
  int i,x;
  
  if (long_seq == FALSE) {         /* then this must be a complete genes  */
    switch (pm->bulk) {
    case 'T':                      /* tidy or fasta formatted header      */      
      fprintf(fblkout, ">%-20.20s%6li\n",
          title, (long int) tot + master_ic);
      break;
    case 'R':                      /* reader header .. don't ask          */ 
      fprintf(fblkout, ">%6li %-70.70s\n",
          (long int) tot + master_ic, title);
      break;
    case 'N':                      /* Conceptually translated DNA header  */
      fprintf(fblkout, ">%-20.20s%6li\n", 
           title, (long int) ((tot + master_ic) / 3));
      break;
    default:                       /* whoops                              */
      printf("\nProgramming error type A2 check code \n");
      my_exit(99, "toutput");
      break;
    }
  } else {                        
      
/* then long_seq must be true, this means we are about to finish writing a*/
/* sequence that has already been written in MAX_GENE chunks to disk)     */
/* when we wrote the original header line, we didn't know the size of the */
/* sequence, but now we do so we are going to update that bit of info     */
/* luckily remembered to record where the header line is in the file      */
/* its at fl_pos_start                                                    */

    fl_pos_curr = ftell(fblkout);   /* record where we are at present     */
    fseek(fblkout, fl_pos_start, 0);/* find the header line for this seq  */
    switch (pm->bulk) {
    case 'T':                       /* Now update the info                */
      fprintf(fblkout, ">%-20.20s%6li",     
          title, (long int) tot + master_ic);   
      break;
    case 'R':                  
      fprintf(fblkout, ">%6li %-70.70s",    
          (long int) tot + master_ic, title);   
      break;
    case 'N':
      fprintf(fblkout, ">%-20.20s%6li", title, 
              (long int) ((tot + master_ic) / 3));
      break;
    default:
      printf("\nProgramming error type A3 check code \n");
      my_exit(99, "output");
    }
    fseek(fblkout, fl_pos_curr, 0);/* now we move back to where we were   */
  }


  while (ic < tot) {               /* keep writing till the array is empty*/
    switch (pm->bulk) {           
    case 'T':
      fprintf(fblkout, "%c", seq[ic++]);
      reg++;
      break;
    case 'R':
      if (space == 3) {            /* Its reader format so print a space   */                                
    fprintf(fblkout, " ");         /* every third base                     */
    space = 0;                      
      } else {                     /* not the 3rd base yet so just print   */
    fprintf(fblkout, "%c", seq[ic++]);
    space++;
    reg++;
      }
      break;
    case 'N':
      for (i = (int) ic, x = 0; i < (int) ic + 3 && i < tot; i++, x++)  
          codon[x] = *(seq + i);   /* get the next three bases if there    */
      codon[x] = '\0';             /* null terminate the codon array       */
      ic += 3;                     /* remember that we have read 3 bases   */
      /* use the function get_aa to return the amino acid for the codon    */
      /* 1 = is for the one letter code of the codon                       */
      fprintf(fblkout, "%c", *get_aa(1, codon));   
      reg++;
      break;
    }
    if (!(reg % 61)) {             /* every 60 bases print a new line char */
      reg = 1;                  
      fprintf(fblkout, "\n");
    }
  }

  if (reg != 1) {                  /* reached the end of sequence so we    */
    fprintf(fblkout, "\n");        /* print a \n char unless we just did   */
    reg = 1;                       /* reset  number of bases printed       */
  }

/* Now that we have finished writing this sequence to disk lets have a     */
/* closer look at it, and do a few diagnostics about the bases used        */

  if (AT_TOT + GC_TOT > AA_TOT*0.5) {/* Assume its DNA then                */
    fprintf(pm->my_err, "%3li>\t%6li %-40.40s\tDNA\tGC%"                   
        " =%5.3f\n"                 /* with G+C content and length of gene */
        ,num_sequence
        ,(long int) tot + master_ic, title
        ,(float) GC_TOT / (GC_TOT + AT_TOT));

    if (non_std_char - IUBC_TOT && pm->warn )  /* any non IUBC characters */
      fprintf(pm->my_err, "\t\t WARNING %d non IUBC standard characters "
          "in sequence %i\n"
          ,non_std_char - IUBC_TOT
          ,num_sequence);
  } else {                         /* if not DNA then it must be a protein */
    fprintf(pm->my_err, "\t%3i>\t%6li %-40.40s\tPROTEIN\n"
        ,num_sequence
        ,(long int) tot + master_ic
        ,title);
    if ( (tot+master_ic)-AA_TOT && pm->warn)  /* non IUBC AA chars        */ 
      fprintf(pm->my_err, "\t\t WARNING %d non "                               
          "standard AA characters "
          "in sequence %i\n"
          ,non_std_char
          ,num_sequence);
  }
  return 1;                        /* return to calling function           */
}


/************************* output_long   **********************************/
/* called to write a block of a sequence that has exceeded the MAX_GENE   */
/* limit. If this is the first time it has been called for this sequence  */
/* (ie. long_seq is false) it write a dummy header line which is updated  */
/* by toutput when the last fragment of the sequence is written to disk   */
/**************************************************************************/

int output_long(FILE * fblkout, char *seq)
{
  long int        ic = 0;
  char            space = 3;
  char            codon[4];
  int i,x;
  
  if (long_seq == FALSE) {         
/* First call to output_long for seq. So record where the header line is  */
/* and then write the dummy header line.                                  */

      fl_pos_start = ftell(fblkout);
    if (pm->bulk == 'R')
      fprintf(fblkout, ">%6s %-72.72s\n", "      ", title); 
    else                   
      fprintf(fblkout, ">%-20.20s%9s\n", title, "    ");    
    long_seq = TRUE;               
  }
/* see toutput for explanation of the switch statement                    */                       
  while (ic < MAX_GENE && ic < tot) {     
    switch (pm->bulk) {
    case 'T':
      fprintf(fblkout, "%c", seq[ic++]);
      reg++;
      break;
    case 'R':
      if (space == 3) {
    fprintf(fblkout, " ");
    space = 0;
      } else {
    fprintf(fblkout, "%c", seq[ic++]);
    space++;
    reg++;
      }
      break;
    case 'N':
      for (i = (int) ic, x = 0; i < (int) ic + 3 && i < tot; i++, x++)
	codon[x] = *(seq + i);
      codon[x] = '\0';
      fprintf(fblkout, "%c", *get_aa(1, codon));
      ic += 3;
      reg++;
      break;
    default:
      printf("\nProgramming error type A1 check code \n");
      my_exit(99, "output_long");
    }
    if (!(reg % 61)) {
      reg = 1;
      fprintf(fblkout, "\n");
    }
  }
  return 1;                        /* return to tidy                      */
}


/*************************  output       **********************************/
/* Called from after subroutine tidy has read the sequence into memory    */
/* or  more accurately counted the codon and amino acid usage. This sub-  */
/* routine, via a switch checks which parameters and indices have been    */
/* requested and write these to file, it handles all output except for COA*/
/**************************************************************************/


void output(char *seq, FILE * foutput, FILE * fblkout, FILE * fcoaout)
{
  char sp;
 
  /* set the column delimiter to something shorter than pm->seperator     */
  sp = (char) (pm->seq_format=='H')? (char) '\t': (char) pm->seperator;

  if (tot) {                       /* still data in array seq..           */
    last_aa = codon_usage_tot(seq, tot);
    if (pcu->ca[last_aa] == 11)
      valid_stops++;                /* check the last codon was a stop    */
  }
  
  /* codon_error, if 4th parameter is 1, then checks for valid start and  */
  /* internal stop codon, if 4th parmater is 2, checks that the last codon*/
  /* is a stop or was partial, and for non-translatable codons            */
  codon_error(last_aa, valid_stops, title, (char) 1); 
  codon_error(last_aa, valid_stops, title, (char) 2);

  /* if we are concatenating sequences then change the title to avger_of  */ 
  if(pm->totals)                   
    (pm->seq_format=='M')?
      strcpy(title, "Average_of_genes"):
      strcpy(title, "Average of genes");
  
                
  if (strchr("RNT", (int) pm->bulk) != NULL) {
    /* better write the remaing sequence in seq to disk                   */
    toutput(fblkout, seq);                              
  } else if (strchr("OCASDLDBX", (int) pm->bulk) != NULL) {

/* These subroutines are self explanatory (see the top of this file)      */
/* are called such that only one can be called for each sequence read     */
/* all these calls are written to the bulk output file                    */

    switch ((int) pm->bulk) {
    case 'S':
      rscu_usage_out(fblkout, ncod, naa);
      break;
    case 'C':
      codon_usage_out(fblkout, ncod, last_aa, valid_stops, title);
      break;
    case 'L':
      raau_usage_out(fblkout, naa);
      break;
    case 'D':
      dinuc_count(seq, tot);
      dinuc_out(fblkout, title);
      break;
    case 'A':
      aa_usage_out(fblkout, naa);
      break;
    case 'B':
      gc_out(foutput, fblkout, 1);
      break;
    case 'O':
      cutab_out(fblkout, ncod, naa);
      break;
    case 'X':
            /* X is no bulk output written to file */
      break;
    default:
      fprintf(stderr, "ERROR-23 %s bulk undefined\n",  pm->prog);
      my_exit(99, "output");
      break;
    }
  } else if (pm->bulk) {            /* just a programming error catch     */
    fprintf(stderr, "ERROR-24 %s -prog undefined\n", pm->prog);
    my_exit(99, "output");
  }
  
  
  /* if an index has been requested then this is true                     */                    
  if (pm->sil_base || pm->cai || pm->fop   || pm->enc  || pm->gc3s ||
            pm->gc || pm->cbi || pm->L_sym || pm->L_aa || pm->coa  || 
            pm->hyd|| pm->aro) {
      /* if this is the first sequence then write a header line           */

    if (num_sequence == 1 || pm->totals) {

      fprintf(foutput, (pm->seq_format == 'H')?
	      "%-25.25s%c":"%-.25s%c"
	      ,"title",sp);
      if (pm->sil_base)
	fprintf(foutput, "%s%c%s%c%s%c%s%c", "T3s",sp,"C3s",sp,"A3s",sp,
"G3s",sp);
      if (pm->cai)
	fprintf(foutput, "%s%c", "CAI",sp);
      if (pm->cbi)
	fprintf(foutput, "%s%c", "CBI",sp);
      if (pm->fop)
	fprintf(foutput, "%s%c", "Fop",sp);
      if (pm->enc)
	fprintf(foutput, "%s%c", "Nc",sp);
      if (pm->gc3s)
	fprintf(foutput, "%s%c", "GC3s" ,sp);
      if (pm->gc)
	fprintf(foutput, "%s%c", "GC"   ,sp);
      if (pm->L_sym)
	fprintf(foutput, "%s%c", "L_sym",sp);
      if (pm->L_aa)
	fprintf(foutput, "%s%c", "L_aa" ,sp);
      if (pm->hyd)
	fprintf(foutput, "%s%c", "Gravy",sp);
      if (pm->aro)
	fprintf(foutput, "%s%c", "Aromo",sp);

      fprintf(foutput, "\n");
    }
    
    /* if output format is human readable print the fixed width sequence  */
    /* name, else print only the name of the sequence                     */
    fprintf(foutput, (pm->seq_format == 'H')?
	      "%-25.25s%c":"%-.25s%c"
	      ,title,sp);

    /*Need to use if statements as we allow more than one index to be calc*/
    /* per sequence read in                                               */
    if (pm->sil_base)
      base_sil_us_out(foutput, ncod, naa);
    if (pm->cai)
      cai_out(foutput, ncod);
    if (pm->cbi)
      cbi_out(foutput, ncod, naa);  
    if (pm->fop)
      fop_out(foutput, ncod);
    if (pm->enc)
      enc_out(foutput, ncod, naa);
    if (pm->gc3s)
      gc_out(foutput, fblkout, 3);
    if (pm->gc)
      gc_out(foutput, fblkout, 2);
    if (pm->L_sym)
      gc_out(foutput, fblkout, 4);
    if (pm->L_aa)
      gc_out(foutput, fblkout, 5);
    if (pm->hyd)
      hydro_out(foutput, naa);
    if (pm->aro)
      aromo_out(foutput, naa);
    if (pm->coa)
      coa_raw_out(fcoaout, ncod, naa, title);

    fprintf(foutput, "\n");

  }
  return;
}

/************************* my_exit       **********************************/
/* Called to clean up open files and generate an intelligent exit message */
/* Also warns if no analysis has been run, the user did not select R from */
/* the main menu. If COA was selected then it reminds the user to look    */
/* at the file summary.coa, and deletes any stray binary files            */
/**************************************************************************/

int my_exit(int error_num, char *message)
{

  fileclose(&pm->inputfile);
  
  /* if we are masuquarading as another program we assign both outputfile */
  /* and tidyout the same filehandle (we don't want to close this twice   */
  if ( pm->outputfile == pm->tidyoutfile ){
   fileclose(&pm->outputfile);
  }else{
   fileclose(&pm->outputfile);
   fileclose(&pm->tidyoutfile);
  }
  
  fileclose(&pm->cuout);
  fileclose(&pm->fopfile);
  fileclose(&pm->cbifile);
  fileclose(&pm->caifile);
  fileclose(&pm->logfile);
  fileclose(&pm->fcoa_in);
  fileclose(&pm->fcoa_out);

  if (pm->inputfile = fopen("cbrawin", "r")) {
    fclose(pm->inputfile);
    deletefile("cbrawin");
  }
  if (pm->inputfile = fopen("cbfcco", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcco");
  }
  if (pm->inputfile = fopen("cbfcli", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcli");
  }
  if (pm->inputfile = fopen("cbfcpc", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcpc");
  }
  if (pm->inputfile = fopen("cbfcpl", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcpl");
  }
  if (pm->inputfile = fopen("cbfcta", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcta");
  }
  if (pm->inputfile = fopen("cbfcvp", "r")) {
    fclose(pm->inputfile);
    deletefile("cbfcvp");
  }
  if (pm->inputfile = fopen("cb1rawin", "r")) {
    fclose(pm->inputfile);
    deletefile("cb1rawin");
  }
  if (error_num == 2 || error_num == 0 ) {
    if (pm->analysis_run) {
      fprintf(stderr, "Files used:\n");
      if (strlen(pm->curr_infilename))
	fprintf(pm->my_err, " Input file was\t %s \n", 
            	pm->curr_infilename);

      if (strlen(pm->curr_outfilename)){
	fprintf(pm->my_err, " Output file was\t %s %s",
		pm->curr_outfilename,
		(pm->codonW) ? " (codon usage indices, e.g. gc3s)\n":"\n");
      }

      if (strlen(pm->curr_tidyoutname)){
	fprintf(pm->my_err, " Output file was\t %s %s",
		pm->curr_tidyoutname,
		(pm->codonW) ? " (bulk output e.g. raw codon usage)\n":"\n");	
      }

      if (pm->coa)
	fprintf(pm->my_err, " For more information about the COrrespondence "
		"Analysis see summary.coa\n");
    } else if ( pm->codonW )          
      fprintf(stderr, " \n\n WARNING You are exiting before codonW has generated any results\n"
	      "  Select 'r' from the main menu to run\n");
  }

  if ( pm->codonW )  printf("\n CodonW has finished\n");

  switch ((int) error_num) {

  case 0:
    /* silent exit */
    exit(0);
    break;
  case 1:
    printf("failed to open file for output <%s>\n", message);
    exit(1);
    break;
  case 2:
    printf("user requested exit <%s>\n", message);
    exit(0);
    break;
  case 3:
    printf("failed to allocate memory <%s>\n", message);
    exit(1);
    break;
  case 4:
    printf("Write to disk failed ! <%s>\n", message);
    exit(1);
    break;
  case 5:
    printf("Read from disk failed! <%s>\n", message );
    exit(1);
    break;
  case 6:
    printf("failed to open file for reading <%s>\n", message);
    exit(1);
    break;
  case 7: 
    printf("failed to close file <%s>\n", message);
    exit(1);
  case 99:
    printf(" Controlled exit <%s>\n",message);
    exit(0);
    break;
  default:
    printf("for unknown reason\n");
    exit(1);
    break;
  }
  return 0;
}

/************************** file_close   **********************************/
/* Fileclose function checks whether the filepointer is open, if so it    */
/* attempts to close the open file handle and assigns a null pointer      */
/* to that  handle                                                        */
/**************************************************************************/

int  fileclose(FILE ** file_pointer)
{
   if (*file_pointer != NULL ) {
     if (fclose(*file_pointer) == EOF ) {
       fprintf(stderr,"Failed to close file %i \n",errno);
       perror ("Unexpected condition in fileclose");
       exit(7);
     }    
     *file_pointer = NULL;           /* make sure file_pointer is null*/
   }
  return 1;
}

/************************** Chelp    **************************************/
/* Chelp scans opens the help file and returns text associated with that  */
/* help keyword. Help keywords are surrounded by hashs, starting in the   */
/* first column of the ASCII help file and are terminated by //           */
/**************************************************************************/

int chelp ( char *help_keyword )
{
 char helplib [MAX_FILENAME_LEN]="";
 char *p=NULL, inhelp=FALSE;
 char QueryString[120];            /* limit for help phrase is 120 chars  */
 char HelpMessage[121];
 int  line_counter=2;              /* assume 2 blank lines to start with  */
 FILE *hfp=NULL;
/* Inital steps is to locate help file                                    */
/* First check if CODONW_H has been set as an environment variable        */
/* If not then assume that the help file is in the current directory      */       
   
 p=getenv( "CODONW_H" );
 if ( p != NULL ) 
     strcpy ( helplib , p );   
 else {
     strcpy ( helplib , "/usr/share/codonw/codonW.hlp");
  }

 hfp=open_file("",helplib, "r", FALSE);

/* if we can't open the help file then explain what we where trying to do  */

 if ( hfp == NULL ) {
           fprintf ( stderr , 
                  "Could not open help file codonw.hlp\n" 
                   "Expected to find this file in %s\n" 
                   "This can be overridden by setting the" 
                   "environmental variable\n"
                   "CODONW_H to the help file location\n", 
                   helplib);
                   pause;           /* make sure they Ack. the error mesg  */
                   return 0;        /* abort                               */
     }   
/* Now that we have opened the help file, assemble the help keyword string */

strcpy (QueryString , "#");
strcat (QueryString , help_keyword );
strcat (QueryString , "#");
fprintf(stderr,"\n\n");             

/* now scan the help file looking for this keyword                         */ 

while ( fgets ( HelpMessage, 120, hfp ) )   {

 if ( strstr (HelpMessage,QueryString) != NULL ) 
     inhelp=TRUE;                                          /* we found it  */

 else if  ( inhelp && strstr ( HelpMessage  , "//") )  {   /* found the end*/
           fileclose(&hfp );                            
	       if ( line_counter )pause;
           return 1;                                     
     }
 /* if inhelp is true we have found the help keyword but not reached EOF   */
 else if  ( inhelp ) {       
     if ( strchr(HelpMessage,'\n') )
        fprintf ( stderr, "%s",HelpMessage );    
                                           /*stderr,it must be interactive */
     else
        fprintf ( stderr, "%s\n",HelpMessage );
                                           /*make sure there are line feeds*/ 


  /* count how many lines I have printed to the terminal and compare it    */
  /* with the length of the terminal screen as defined by pm->term_length  */

  if (line_counter++ >= pm->term_length-3 && line_counter ) {
    line_counter=0;
    pause;
    fprintf(stderr, "%s",HelpMessage);
   }
 }   
}   

/* Error catches for problems with help file                              */
if ( HelpMessage == NULL && inhelp == FALSE ){
   fprintf ( stderr ," Error in help file, %s not found ", QueryString);
   pause;
  }
else {
    fprintf (stderr , "Premature end of help file ...  \n");
    pause;
 }   
return 0;                           /* failed for some reason             */
}


/******************** WasHelpCalled     ***********************************/
/* Checks the string input to see if the user asked for help              */
/**************************************************************************/

char WasHelpCalled ( char * input ) {
    char ans = FALSE;

    if ( strlen ( input) == 1 && (char)toupper((int)input[0]) == 'H') 
        ans = TRUE;
    else if ( !strcmp ( input , "help") ) 
        ans = TRUE;
    else if ( !strcmp ( input , "HELP") ) 
        ans = TRUE;

    return ans;
}