source: branches/help/GDE/CLUSTAL/amenu.c

Last change on this file was 19518, checked in by westram, 3 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.0 KB
Line 
1/* Menus and command line interface for CLUSTAL V  */
2
3#include <stdio.h>
4#include <string.h>
5#include <ctype.h>
6#include <stdlib.h>
7#include "clustalv.h"
8
9#include <gets_noOverflow.h>
10
11/*
12 *      Prototypes
13 */
14
15extern void     getstr(const char *,char *, int);
16extern double   getreal(char *,double,double);
17extern int              getint(const char *,int,int,int);
18extern void             do_system(void);
19extern void             make_pamo(int);
20extern int          readseqs(int);
21extern void             get_path(char *,char *);
22extern void             show_pair(void);
23extern void             upgma(int);
24extern void             myers(int);
25extern void     phylogenetic_tree(void);
26extern void         bootstrap_tree(void);
27extern void             error(const char *,...);
28extern int              SeqGCGCheckSum(char *, int);
29
30
31void init_amenu(void);
32void parse_params(void);
33void main_menu(void);
34FILE * open_output_file(const char *, const char *, char *, const char *);
35#if UNIX
36FILE * open_path(const char *);
37#endif
38
39static Boolean check_param(const char *, const char *, char *, int *);
40static void get_help(int);                      /* Help procedure */
41static void pair_menu(void);
42static void multi_menu(void);
43static void multiple_align_menu(void);          /* multiple alignments menu */
44static void profile_align_menu(void);           /* profile       "      "   */
45static void phylogenetic_tree_menu(void);       /* NJ trees/distances menu  */
46static void read_matrix(void);
47static void seq_input(void);
48static void align(void);
49static void make_tree(void);
50static void get_tree(void);
51static Boolean user_mat(char *);
52static void clustal_out(FILE *);
53static void nbrf_out(FILE *);
54static void gcg_out(FILE *);
55static void phylip_out(FILE *);
56static Boolean open_alignment_output(char *);
57static void create_alignment_output(void);
58static void reset(void);
59static void profile_input(int);                 /* DES read a profile */
60static void format_options_menu(void);          /* format of alignment output */
61static void profile_align(void);                /* Align 2 alignments */
62
63/*
64*        Global variables
65*/
66
67char *lin1, *lin2;
68
69extern char *amino_acid_order;
70extern char *nucleic_acid_order;
71extern void *ckalloc(size_t);
72extern Boolean          percent,is_weight,dnaflag;
73extern int              wind_gap,ktup,window,signif;
74extern int               dna_wind_gap, dna_ktup, dna_window, dna_signif;
75extern int              prot_wind_gap,prot_ktup,prot_window,prot_signif;
76extern unsigned int boot_ran_seed;
77extern int      gap_open,      gap_extend;
78extern int  dna_gap_open,  dna_gap_extend;
79extern int prot_gap_open, prot_gap_extend;
80extern int boot_ntrials;                /* number of bootstrap trials */
81extern int              xover,big_pam;
82extern char *params;
83extern char             seqname[],treename[];
84extern char     *matptr,pam100mt[],pam250mt[],idmat[];
85extern int              nseqs,nblocks;
86extern int              weights[21][21];
87extern FILE *   tree;
88extern FILE *clustal_outfile, *gcg_outfile, *nbrf_outfile, *phylip_outfile;
89extern char **  names, **titles;
90extern int *seqlen_array;
91extern char **seq_array;
92extern char ntrials;            /* number of bootstrap trials (trees.c) */
93
94Boolean usemenu;
95char mtrxnam[FILENAMELEN+1];
96Boolean explicit_dnaflag;  /* Explicit setting of sequence type on comm.line*/
97Boolean output_clustal, output_nbrf, output_phylip, output_gcg;
98Boolean empty;
99Boolean profile1_empty, profile2_empty;   /* whether or not profiles   */
100int profile1_nseqs;            /* have been filled; the no. of seqs in prof 1*/
101Boolean tossgaps, kimura;
102int  matnum;
103static char clustal_outname[FILENAMELEN+1], gcg_outname[FILENAMELEN+1];
104static char  phylip_outname[FILENAMELEN+1],nbrf_outname[FILENAMELEN+1];
105
106static const char *pam_matrix_name[] = {
107                "PAM 100",
108                "PAM 250",
109                "Identity matrix",
110                "user defined"   };
111               
112char usermat[210];
113
114void init_amenu()
115{
116        matnum=2;
117        empty=TRUE;
118        explicit_dnaflag = FALSE;     
119    profile1_empty=TRUE;     /* DES */
120    profile2_empty=TRUE;     /* DES */
121        output_clustal = TRUE;
122        output_gcg     = FALSE;
123        output_phylip  = FALSE;
124        output_nbrf    = FALSE;
125
126    lin1 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );
127    lin2 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );
128}
129
130
131
132
133static Boolean check_param(const char *paramstr, const char *probe, char *arg, int *len)
134{
135        int i,j,k,paramlen;
136        char testarg[80];
137
138        paramlen = strlen(paramstr);
139        i = 1;
140        while(i < paramlen) {
141                for(j=i+1; j<=paramlen      &&
142                        paramstr[j] != '/'  &&
143                        paramstr[j] != EOS; j++) ;
144                for(k=i+1; k<j; k++)
145                        if(paramstr[k] == '=') {
146                                strncpy(testarg,&paramstr[i],k-i);
147                                if(!strncmp(testarg,probe,k-i)){
148                                        strncpy(arg,&paramstr[k+1],j-k-1);
149                                        arg[j-k-1] = EOS;
150                                        *len = strlen(arg);
151                                        return TRUE;
152                                }
153                        }
154                strncpy(testarg,&paramstr[i],j-i);
155                if(!strncmp(testarg,probe,j-i)) {
156                        *len = 0;
157                        arg[0] = EOS;
158                        return TRUE;
159                }
160                i = j + 1;
161        }
162        return FALSE;
163}
164
165
166
167 
168#if UNIX
169FILE *open_path(const char *fname)  /* to open in read-only file fname searching for
170                                 it through all path directories */
171{
172#define Mxdir 70
173        char dir[Mxdir+1], *path, *deb, *fin;
174        FILE *fich;
175        int lf, ltot;
176 
177        path=getenv("PATH");    /* get the list of path directories,
178                                        separated by :
179                                */
180        if (path == NULL ) return fopen(fname,"r");
181        lf=strlen(fname);
182        deb=path;
183        do
184                {
185                fin=strchr(deb,':');
186                if(fin!=NULL)
187                        { strncpy(dir,deb,fin-deb); ltot=fin-deb; }
188                else
189                        { strcpy(dir,deb); ltot=strlen(dir); }
190                /* now one directory is in string dir */
191                if( ltot + lf + 1 <= Mxdir)
192                        {
193                        dir[ltot]='/';
194                        strcpy(dir+ltot+1,fname); /* now dir is appended with fi
195   lename */
196                        if( (fich = fopen(dir,"r") ) != NULL) break;
197                        }
198                else fich = NULL;
199                deb=fin+1;
200                }
201        while (fin != NULL);
202        return fich;
203}
204#endif
205
206 
207
208static void get_help(int help_pointer)    /* Help procedure */
209{       
210        FILE *help_file;
211        int  i, number, nlines;
212        Boolean found_help;
213        char temp[MAXLINE+1];
214        const char *digits = "0123456789";
215        const char *help_marker    = ">>HELP<<";
216
217#if MSDOS
218        const char *help_file_name = "clustalv.hlp";
219#else
220        const char *help_file_name = "clustalv_help";
221#endif
222
223#if VMS
224        if((help_file=fopen(help_file_name,"r","rat=cr","rfm=var"))==NULL) {
225#else
226#if UNIX
227        if((help_file=open_path(help_file_name))==NULL) {
228#else
229        if((help_file=fopen(help_file_name,"r"))==NULL) {
230#endif
231#endif
232                error("Cannot open help file [%s]",help_file_name);
233                return;
234        }
235
236        nlines = 0;
237        number = -1;
238        found_help = FALSE;
239
240        while(TRUE) {
241                if(fgets(temp,MAXLINE+1,help_file) == NULL) {
242                        if(!found_help)
243                                error("No help found in help file");
244                        fclose(help_file);
245                        return;
246                }
247                if(strstr(temp,help_marker)) {
248                        for(i=0; i<MAXLINE; i++)
249                                if(strchr(digits, temp[i])) {
250                                        number = temp[i] - '0';
251                                        break;
252                                }
253                }
254                if(number == help_pointer) {
255                        found_help = TRUE;
256                        while(fgets(temp,MAXLINE+1,help_file)) {
257                                if(strstr(temp, help_marker)){
258                                        if(usemenu) {
259                                                fprintf(stdout,"\n");
260                                                getstr("Press [RETURN] to continue",lin2, MAXLINE+1);
261                                        }
262                                        fclose(help_file);
263                                        return;
264                                }
265                               fputs(temp,stdout);
266                               ++nlines;
267                               if(usemenu) {
268                                  if(nlines >= PAGE_LEN) {
269                                           fprintf(stdout,"\n");
270                                           getstr("Press [RETURN] to continue or  X  to stop",lin2, MAXLINE+1);
271                                           if(toupper(*lin2) == 'X') {
272                                                   fclose(help_file);
273                                                   return;
274                                           }
275                                           else
276                                                   nlines = 0;
277                                   }
278                               }
279                        }
280                        if(usemenu) {
281                                fprintf(stdout,"\n");
282                                getstr("Press [RETURN] to continue",lin2, MAXLINE+1);
283                        }
284                        fclose(help_file);
285                }
286        }
287
288}
289
290
291void parse_params()
292{
293        int i,len,temp;
294        char param_arg[80];
295        int param_arg_len;
296
297        Boolean do_align, do_tree, do_boot, do_profile, do_something;
298
299/* command line switches for PARAMETERS **************************/
300        static const char *fixedgapst    = "fixedgap";
301        static const char *floatgapst    = "floatgap";
302        static const char *kimurast      = "kimura";
303        static const char *ktuplest      = "ktuple";
304        static const char *matrixst      = "matrix";
305        static const char *outputst      = "output";
306        static const char *pairgapst     = "pairgap";
307        static const char *scorest       = "score";
308        static const char *seedst        = "seed";
309        static const char *topdiagsst    = "topdiags";
310        static const char *tossgapsst    = "tossgaps";
311        static const char *transitionsst = "transitions";
312        static const char *typest        = "type";
313        static const char *windowst      = "window";
314
315/* command line switches for DATA       **************************/
316        static const char *infilest   = "infile";
317        static const char *profile1st = "profile1";
318        static const char *profile2st = "profile2";
319
320/* command line switches for VERBS      **************************/
321        static const char *alignst     = "align";
322        static const char *bootstrapst = "bootstrap";
323        static const char *checkst     = "check"; /*  /check = /help */
324        static const char *helpst      = "help";
325        static const char *treest      = "tree";
326
327        fprintf(stdout,"\n\n\n");
328        fprintf(stdout," CLUSTAL V ... Multiple Sequence Alignments\n\n\n");
329
330        do_align = do_tree = do_boot = do_profile = do_something = FALSE;
331
332        *seqname=EOS;
333
334        len=strlen(params);
335        for(i=0;i<len;++i) params[i]=tolower(params[i]);
336
337        if(check_param(params, helpst,  param_arg, &param_arg_len) ||
338           check_param(params, checkst, param_arg, &param_arg_len) ) {
339                get_help(9);
340                exit(1);
341        }
342
343
344/*****************************************************************************/
345/*  Check to see if sequence type is explicitly stated..override ************/
346/* the automatic checking (DNA or Protein).   /type=d or /type=p *************/
347/*****************************************************************************/
348        if(check_param(params, typest, param_arg, &param_arg_len)) 
349                if(param_arg_len > 0) {
350                        if(param_arg[0] == 'p') {
351                                dnaflag = FALSE;
352                                explicit_dnaflag = TRUE;
353                                fprintf(stdout,
354                                "\nSequence type explicitly set to Protein\n");
355                        }
356                        else if(param_arg[0] == 'd') {
357                                fprintf(stdout,
358                                "\nSequence type explicitly set to DNA\n");
359                                dnaflag = TRUE;
360                                explicit_dnaflag = TRUE;
361                        }
362                        else
363                                fprintf(stdout,"\nUnknown sequence type %s\n",
364                                param_arg);
365                }
366
367
368/***************************************************************************
369*   check to see if 1st parameter does not start with '/' i.e. look for an *
370*   input file as first parameter.   The input file can also be specified  *
371*   by /infile=fname.                                                      *
372****************************************************************************/
373
374        for (i=0; params[i] != '/' && params[i] != EOS; i++) ;
375        if(i > 0) {
376                strncpy(seqname, &params[0], i);
377                seqname[i] = EOS;
378                nseqs = readseqs(1);
379                if(nseqs < 2) {
380                        fprintf(stderr,
381                        "\nNo. of seqs. read = %d. No alignment!\n",nseqs);
382                        exit(1);
383                }
384                for(i = 1; i<=nseqs; i++) 
385                        fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
386                        i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
387                empty = FALSE;
388                do_something = TRUE;
389        }
390
391/**************************************************/
392/*  Look for /infile=file.ext on the command line */
393/**************************************************/
394
395        if(check_param(params, infilest, param_arg, &param_arg_len)) {
396                if(param_arg_len == 0) {
397                        error("Bad sequence file name");
398                        exit(1);
399                }
400                strncpy(seqname, param_arg, param_arg_len);
401/*              seqname[param_arg_len-1] = EOS;  */
402                nseqs = readseqs(1);
403                if(nseqs < 2) {
404                        fprintf(stderr,
405                        "\nNo. of seqs. read = %d. No alignment!\n",nseqs);
406                        exit(1);
407                }
408                for(i = 1; i<=nseqs; i++) 
409                        fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
410                        i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
411                empty = FALSE;
412                do_something = TRUE;
413        }
414
415/*********************************************************/
416/* Look for /profile1=file.ext  AND  /profile2=file2.ext */
417/* You must give both file names OR neither.             */
418/*********************************************************/
419
420        if(check_param(params, profile1st, param_arg, &param_arg_len)) {
421                if(param_arg_len == 0) {
422                        error("Bad profile 1 file name");
423                        exit(1);
424                }
425                strncpy(seqname, param_arg, param_arg_len);
426/*              seqname[param_arg_len-1] = EOS; */
427                profile_input(1);
428                if(nseqs <= 0) 
429                        exit(1);
430        }
431
432        if(check_param(params, profile2st, param_arg, &param_arg_len)) {
433                if(param_arg_len == 0) {
434                        error("Bad profile 2 file name");
435                        exit(1);
436                }
437                if(profile1_empty) {
438                        error("Only 1 profile file (profile 2) specified.");
439                        exit(1);
440                }
441                strncpy(seqname, param_arg, param_arg_len);
442/*              seqname[param_arg_len-1] = EOS; */
443                profile_input(2);
444                if(nseqs > profile1_nseqs) 
445                        do_something = do_profile = TRUE;
446                else {
447                        error("No sequences read from profile 2");
448                        exit(1);
449                }
450        }
451
452/*************************************************************************/
453/* Look for /tree or /bootstrap or /align ********************************/
454/*************************************************************************/
455
456        if(check_param(params, treest, param_arg, &param_arg_len))
457                if(empty) {
458                        error("Cannot draw tree.  No input alignment file");
459                        exit(1);
460                }
461                else 
462                        do_tree = TRUE;
463
464        if(check_param(params, bootstrapst, param_arg, &param_arg_len))
465                if(empty) {
466                        error("Cannot bootstrap tree. No input alignment file");
467                        exit(1);
468                }
469                else {
470                        temp = 0;
471                        if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
472                        if(temp > 0)          boot_ntrials = temp;
473                        do_boot = TRUE;
474                }
475
476        if(check_param(params, alignst, param_arg, &param_arg_len))
477                if(empty) {
478                        error("Cannot align sequences.  No input file");
479                        exit(1);
480                }
481                else 
482                        do_align = TRUE;
483
484        if( (!do_tree) && (!do_boot) && (!empty) && (!do_profile) ) 
485                do_align = TRUE;
486
487        if(!do_something) {
488                error("No input file(s) specified");
489                exit(1);
490        }
491
492       
493        if(dnaflag) {
494                gap_open   = dna_gap_open;
495                gap_extend = dna_gap_extend;
496                ktup       = dna_ktup;
497                window     = dna_window;
498                signif     = dna_signif;
499                wind_gap   = dna_wind_gap;
500        }
501        else {
502                gap_open   = prot_gap_open;
503                gap_extend = prot_gap_extend;
504                ktup       = prot_ktup;
505                window     = prot_window;
506                signif     = prot_signif;
507                wind_gap   = prot_wind_gap;
508        }
509
510
511/****************************************************************************/
512/* look for parameters on command line  e.g. gap penalties, k-tuple etc.    */
513/****************************************************************************/
514
515/*** ? /kimura  */
516        if(check_param(params, kimurast, param_arg, &param_arg_len))
517                kimura = TRUE;
518
519
520/*** ? /tossgaps */
521        if(check_param(params, tossgapsst, param_arg, &param_arg_len))
522                tossgaps = TRUE;
523
524
525/*** ? /score=percent or /score=absolute */
526        if(check_param(params, scorest, param_arg, &param_arg_len))
527                if(param_arg_len > 0) {
528                        if(param_arg[0] == 'p')
529                                percent = TRUE;
530                        else if(param_arg[0] == 'a') 
531                                percent = FALSE;
532                        else
533                                fprintf(stdout,"\nUnknown SCORE type: %s\n",
534                                param_arg);
535                }
536
537
538/*** ? /transitions */
539        if(check_param(params, transitionsst, param_arg, &param_arg_len))
540                is_weight = FALSE;
541
542
543/*** ? /seed=n */
544        if(check_param(params, seedst, param_arg, &param_arg_len)) {
545                temp = 0;
546                if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
547                if(temp > 0) boot_ran_seed = temp;
548        fprintf(stdout,"\ntemp = %d; seed = %u;\n",temp,boot_ran_seed);
549        }
550
551
552/*** ? /output=PIR, GCG or PHYLIP  Only 1 can be switched on at a time */
553        if(check_param(params, outputst, param_arg, &param_arg_len))
554                if(param_arg_len > 0) {
555                        if(param_arg[0] == 'g')  {      /* GCG */
556                                output_gcg     = TRUE;
557                                output_clustal = FALSE;
558                        }
559                        else if(param_arg[0] == 'p') 
560                                if(param_arg[1] == 'i') {
561                                        output_nbrf    = TRUE;
562                                        output_clustal = FALSE;
563                                }
564                                else if(param_arg[1] == 'h') {
565                                        output_phylip  = TRUE;
566                                        output_clustal = FALSE;
567                                }
568                                else
569                                        fprintf(stdout,"\nUnknown OUTPUT type: %s\n",
570                                        param_arg);
571                        else
572                                fprintf(stdout,"\nUnknown OUTPUT type: %s\n",
573                                param_arg);
574                }
575
576
577/*** ? /ktuple=n */
578        if(check_param(params, ktuplest, param_arg, &param_arg_len)) {
579                temp = 0;
580                if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
581                if(temp > 0) {
582                        if(dnaflag) {
583                                if(temp <= 4) { 
584                                        ktup         = temp;
585                                        dna_ktup     = ktup;
586                                        wind_gap     = ktup + 4;
587                                        dna_wind_gap = wind_gap;
588                                }
589                        }
590                        else {
591                                if(temp <= 2) {
592                                        ktup          = temp;
593                                        prot_ktup     = ktup;
594                                        wind_gap      = ktup + 3;
595                                        prot_wind_gap = wind_gap;
596                                }
597                        }
598                }
599        }
600
601
602/*** ? /pairgap=n */
603        if(check_param(params, pairgapst, param_arg, &param_arg_len)) {
604                temp = 0;
605                if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
606                if(temp > 0)
607                        if(dnaflag) {
608                                if(temp > ktup) {
609                                        wind_gap     = temp;
610                                        dna_wind_gap = wind_gap;
611                                }
612                        }
613                        else {
614                                if(temp > ktup) {
615                                        wind_gap      = temp;
616                                        prot_wind_gap = wind_gap;
617                                }
618                        }
619        }
620       
621
622/*** ? /topdiags=n   */
623        if(check_param(params, topdiagsst, param_arg, &param_arg_len)) {
624                temp = 0;
625                if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
626                if(temp > 0)
627                        if(dnaflag) {
628                                if(temp > ktup) { 
629                                        signif       = temp;
630                                        dna_signif   = signif;
631                                }
632                        }
633                        else {
634                                if(temp > ktup) {
635                                        signif        = temp;
636                                        prot_signif   = signif;
637                                }
638                        }
639        }
640       
641
642/*** ? /window=n  */
643        if(check_param(params, windowst, param_arg, &param_arg_len)) {
644                temp = 0;
645                if(param_arg_len > 0) sscanf(param_arg,"%d",&temp);
646                if(temp > 0)
647                        if(dnaflag) {
648                                if(temp > ktup) { 
649                                        window       = temp;
650                                        dna_window   = window;
651                                }
652                        }
653                        else {
654                                if(temp > ktup) {
655                                        window        = temp;
656                                        prot_window   = window;
657                                }
658                        }
659        }
660
661
662/*** ? /matrix=pam100, or ID or file.ext (user's file)  */
663        if(check_param(params, matrixst, param_arg, &param_arg_len))
664                if(param_arg_len > 0) {
665                        if( !strcmp(param_arg, "pam100") )  {
666                                matptr = pam100mt;
667                                make_pamo(0);
668                                prot_gap_open   = 13;
669                                prot_gap_extend = 13;
670                                if(!dnaflag) gap_open   = prot_gap_open;
671                                if(!dnaflag) gap_extend = prot_gap_extend;
672                                matnum = 1;
673                        }
674                        else if( !strcmp(param_arg, "id") )  {
675                                matptr = idmat;
676                                make_pamo(0);
677                                matnum = 3;
678                        }
679                        else if(user_mat(param_arg))
680                                matnum = 4;
681                        else
682                                fprintf(stdout,"\nUnknown MATRIX type: %s\n",
683                                param_arg);
684                }
685
686
687/*** ? /fixedgap=n  */
688        if(check_param(params, fixedgapst, param_arg, &param_arg_len)) {
689                temp = 0;
690                if(param_arg_len > 0)
691                        sscanf(param_arg,"%d",&temp);
692                if(temp > 0)
693                        if(dnaflag) {
694                                        gap_open     = temp;
695                                        dna_gap_open = gap_open;
696                        }
697                        else {
698                                        gap_open      = temp;
699                                        prot_gap_open = gap_open;
700                        }
701        }
702
703
704/*** ? /floatgap=n   */
705        if(check_param(params, floatgapst, param_arg, &param_arg_len)) {
706                temp = 0;
707                if(param_arg_len > 0)
708                        sscanf(param_arg,"%d",&temp);
709                if(temp > 0)
710                        if(dnaflag) {
711                                        gap_extend      = temp;
712                                        dna_gap_extend  = gap_extend;
713                        }
714                        else {
715                                        gap_extend      = temp;
716                                        prot_gap_extend = gap_extend;
717                        }
718        }
719
720
721/****************************************************************************/
722/* Now do whatever has been requested ***************************************/
723/****************************************************************************/
724
725        if(do_profile)
726                profile_align();
727
728        if(do_align)
729                align();
730
731        if(do_tree)
732                phylogenetic_tree();
733
734        if(do_boot)
735                bootstrap_tree();
736
737        exit(1);
738
739/*******whew!***now*go*home****/
740}
741
742
743
744
745
746
747
748void main_menu()
749{
750        while(TRUE) {
751                fprintf(stdout,"\n\n\n");
752                fprintf(stdout," **************************************************************\n");
753                fprintf(stdout," ********* CLUSTAL V ... Multiple Sequence Alignments  ********\n");
754                fprintf(stdout," **************************************************************\n");
755                fprintf(stdout,"\n\n");
756               
757                fprintf(stdout,"     1. Sequence Input From Disc\n");
758                fprintf(stdout,"     2. Multiple Alignments\n");
759                fprintf(stdout,"     3. Profile Alignments\n");
760                fprintf(stdout,"     4. Phylogenetic trees\n");
761                fprintf(stdout,"\n");
762                fprintf(stdout,"     S. Execute a system command\n");
763                fprintf(stdout,"     H. HELP\n");
764                fprintf(stdout,"     X. EXIT (leave program)\n\n\n");
765               
766                getstr("Your choice",lin1, MAXLINE+1);
767
768                switch(toupper(*lin1)) {
769                        case '1': seq_input();
770                                break;
771                        case '2': multiple_align_menu();
772                                break;
773                        case '3': profile_align_menu();
774                                break;
775                        case '4': phylogenetic_tree_menu();
776                                break;
777                        case 'S': do_system();
778                                break;
779                        case '?':
780                        case 'H': get_help(1);
781                                break;
782                        case 'Q':
783                        case 'X': exit(0);
784                                break;
785                        default: fprintf(stderr,"\n\nUnrecognised Command\n\n");
786                                break;
787                }
788        }
789}
790
791
792
793
794
795
796
797
798
799static void multiple_align_menu()
800{
801    while(TRUE)
802    {
803        fprintf(stdout,"\n\n\n");
804        fprintf(stdout,"******Multiple*Alignment*Menu******\n");
805        fprintf(stdout,"\n\n");
806
807
808        fprintf(stdout,"    1.  Do complete multiple alignment now\n");
809        fprintf(stdout,"    2.  Produce dendrogram file only\n");
810        fprintf(stdout,"    3.  Use old dendrogram file\n");
811        fprintf(stdout,"    4.  Pairwise alignment parameters\n");
812        fprintf(stdout,"    5.  Multiple alignment parameters\n");
813        fprintf(stdout,"    6.  Output format options\n");
814        fprintf(stdout,"\n");
815        fprintf(stdout,"    S.  Execute a system command\n");
816        fprintf(stdout,"    H.  HELP\n");
817        fprintf(stdout,"    or press [RETURN] to go back to main menu\n\n\n");
818
819        getstr("Your choice",lin1, MAXLINE+1);
820        if(*lin1 == EOS) return;
821
822        switch(toupper(*lin1))
823        {
824        case '1': align();
825            break;
826        case '2': make_tree();
827            break;
828        case '3': get_tree();
829            break;
830        case '4': pair_menu();
831            break;
832        case '5': multi_menu();
833            break;
834        case '6': format_options_menu();
835            break;
836        case 'S': do_system();
837            break;
838        case '?':
839        case 'H': get_help(2);
840            break;
841        case 'Q':
842        case 'X': return;
843
844        default: fprintf(stderr,"\n\nUnrecognised Command\n\n");
845            break;
846        }
847    }
848}
849
850
851
852
853
854
855
856
857
858static void profile_align_menu()
859{
860    while(TRUE)
861    {
862        fprintf(stdout,"\n\n\n");
863        fprintf(stdout,"******Profile*Alignment*Menu******\n");
864        fprintf(stdout,"\n\n");
865
866        fprintf(stdout,"    1.  Input 1st. profile/sequence\n");
867        fprintf(stdout,"    2.  Input 2nd. profile/sequence\n");
868        fprintf(stdout,"    3.  Do alignment now\n");
869        fprintf(stdout,"    4.  Alignment parameters\n");
870        fprintf(stdout,"    5.  Output format options\n");
871        fprintf(stdout,"\n");
872        fprintf(stdout,"    S.  Execute a system command\n");
873        fprintf(stdout,"    H.  HELP\n");
874        fprintf(stdout,"    or press [RETURN] to go back to main menu\n\n\n");
875
876        getstr("Your choice",lin1, MAXLINE+1);
877        if(*lin1 == EOS) return;
878
879        switch(toupper(*lin1))
880        {
881        case '1': profile_input(1);      /* 1 => 1st profile */
882            break;
883        case '2': profile_input(2);      /* 2 => 2nd profile */
884            break;
885        case '3': profile_align();       /* align the 2 alignments now */
886            break;
887        case '4': multi_menu();
888            break;
889        case '5': format_options_menu();
890            break;
891        case 'S': do_system();
892            break;
893        case '?':
894        case 'H': get_help(6);
895            break;
896        case 'Q':
897        case 'X': return;
898
899        default: fprintf(stderr,"\n\nUnrecognised Command\n\n");
900            break;
901        }
902    }
903}
904
905
906
907
908
909
910
911
912
913
914
915static void phylogenetic_tree_menu()
916{
917    while(TRUE)
918    {
919        fprintf(stdout,"\n\n\n");
920        fprintf(stdout,"******Phylogenetic*tree*Menu******\n");
921        fprintf(stdout,"\n\n");
922
923        fprintf(stdout,"    1.  Input an alignment\n");
924        fprintf(stdout,"    2.  Exclude positions with gaps?        ");
925        if(tossgaps)
926                fprintf(stdout,"= ON\n");
927        else
928                fprintf(stdout,"= OFF\n");
929        fprintf(stdout,"    3.  Correct for multiple substitutions? ");
930        if(kimura)
931                fprintf(stdout,"= ON\n");
932        else
933                fprintf(stdout,"= OFF\n");
934        fprintf(stdout,"    4.  Draw tree now\n");
935        fprintf(stdout,"    5.  Bootstrap tree\n");
936        fprintf(stdout,"\n");
937        fprintf(stdout,"    S.  Execute a system command\n");
938        fprintf(stdout,"    H.  HELP\n");
939        fprintf(stdout,"    or press [RETURN] to go back to main menu\n\n\n");
940
941        getstr("Your choice",lin1, MAXLINE+1);
942        if(*lin1 == EOS) return;
943
944        switch(toupper(*lin1))
945        {
946                case '1': seq_input();
947                        break;
948                case '2': tossgaps ^= TRUE;
949                        break;
950                case '3': kimura ^= TRUE;;
951                        break;
952                case '4': phylogenetic_tree();
953                        break;
954                case '5': bootstrap_tree();
955                        break;
956                case 'S': do_system();
957                        break;
958                case '?':
959                case 'H': get_help(7);
960                        break;
961                case 'Q':
962                case 'X': return;
963
964                default: fprintf(stderr,"\n\nUnrecognised Command\n\n");
965                break;
966        }
967    }
968}
969
970
971
972
973
974
975static void format_options_menu()      /* format of alignment output */
976{       
977        char path[FILENAMELEN+1];
978
979        while(TRUE) {
980        fprintf(stdout,"\n\n\n");
981        fprintf(stdout," ********* Format of Alignment Output *********\n");
982        fprintf(stdout,"\n\n");
983        fprintf(stdout,"     1. Toggle CLUSTAL format output   =  %s\n",
984                                        (!output_clustal) ? "OFF" : "ON");
985        fprintf(stdout,"     2. Toggle NBRF/PIR format output  =  %s\n",
986                                        (!output_nbrf) ? "OFF" : "ON");
987        fprintf(stdout,"     3. Toggle GCG format output       =  %s\n",
988                                        (!output_gcg) ? "OFF" : "ON");
989        fprintf(stdout,"     4. Toggle PHYLIP format output    =  %s\n\n",
990                                        (!output_phylip) ? "OFF" : "ON");
991      if(empty)
992        fprintf(stdout,"\n");
993      else
994        fprintf(stdout,"     5. Create alignment output file(s) now?\n");
995        fprintf(stdout,"     H. HELP\n\n\n");   
996       
997        getstr("Enter number (or [RETURN] to exit)",lin2, MAXLINE+1);
998                if(*lin2 == EOS) return;
999               
1000                switch(toupper(*lin2)) {
1001                        case '1':
1002                                output_clustal ^= TRUE;
1003                                break;
1004                        case '2':
1005                                output_nbrf ^= TRUE;
1006                                break;
1007                        case '3':
1008                                output_gcg ^= TRUE;
1009                                break;
1010                        case '4':
1011                                output_phylip ^= TRUE;
1012                                break;
1013                        case '5':
1014                                if(empty) break;
1015                                get_path(seqname,path);
1016                                if(!open_alignment_output(path)) break;
1017                                create_alignment_output();
1018                                break;
1019                        case '?':
1020                        case 'H':
1021                                get_help(5);
1022                                break;
1023                        default:
1024                                fprintf(stderr,"\n\nUnrecognised Command\n\n");
1025                                break;
1026                }
1027        }
1028}
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041static void pair_menu()
1042{
1043        if(dnaflag) {
1044                ktup     = dna_ktup;
1045                window   = dna_window;
1046                signif   = dna_signif;
1047                wind_gap = dna_wind_gap;
1048        }
1049        else {
1050                ktup     = prot_ktup;
1051                window   = prot_window;
1052                signif   = prot_signif;
1053                wind_gap = prot_wind_gap;
1054        }
1055
1056        while(TRUE) {
1057                const char *lin3 = percent ? "Percentage" : "Absolute";
1058       
1059                fprintf(stdout,"\n\n\n");
1060                fprintf(stdout," ********* WILBUR/LIPMAN PAIRWISE ALIGNMENT PARAMETERS *********\n");
1061                fprintf(stdout,"\n\n");
1062
1063                fprintf(stdout,"     1. Toggle Scoring Method  :%s\n",lin3);
1064                fprintf(stdout,"     2. Gap Penalty            :%d\n",wind_gap);
1065                fprintf(stdout,"     3. K-tuple                :%d\n",ktup);
1066                fprintf(stdout,"     4. No. of top diagonals   :%d\n",signif);
1067                fprintf(stdout,"     5. Window size            :%d\n\n",window);
1068                fprintf(stdout,"     H. HELP\n\n\n");
1069               
1070                getstr("Enter number (or [RETURN] to exit)",lin2, MAXLINE+1);
1071                if( *lin2 == EOS) {
1072                        if(dnaflag) {
1073                                dna_ktup      = ktup;
1074                                dna_window    = window;
1075                                dna_signif    = signif;
1076                                dna_wind_gap  = wind_gap;
1077                        }
1078                        else {
1079                                prot_ktup     = ktup;
1080                                prot_window   = window;
1081                                prot_signif   = signif;
1082                                prot_wind_gap = wind_gap;
1083                        }
1084                        return;
1085                }
1086               
1087                switch(toupper(*lin2)) {
1088                        case '1':
1089                                percent ^= TRUE;
1090                                break;
1091                        case '2':
1092                                fprintf(stdout,"Gap Penalty Currently: %d\n",wind_gap);
1093                                wind_gap=getint("Enter number",1,500,wind_gap);
1094                                break;
1095                        case '3':
1096                                fprintf(stdout,"K-tuple Currently: %d\n",ktup);
1097                                ktup=getint("Enter number",1,4,ktup);
1098                                break;
1099                        case '4':
1100                                fprintf(stdout,"Top diagonals Currently: %d\n",signif);
1101                                signif=getint("Enter number",1,MAXLEN,signif);
1102                                break;
1103                        case '5':
1104                                fprintf(stdout,"Window size Currently: %d\n",window);
1105                                window=getint("Enter number",1,50,window);
1106                                break;
1107                        case '?':
1108                        case 'H':
1109                                get_help(3);
1110                                break;
1111                        default:
1112                                fprintf(stderr,"\n\nUnrecognised Command\n\n");
1113                                break;
1114                }
1115        }
1116}
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129static void multi_menu()
1130{
1131        if(dnaflag) {
1132                gap_open   = dna_gap_open;
1133                gap_extend = dna_gap_extend;
1134        }
1135        else {
1136                gap_open   = prot_gap_open;
1137                gap_extend = prot_gap_extend;
1138        }
1139
1140        while(TRUE) {
1141                const char *lin3 = is_weight ? "Weighted" :"Unweighted";
1142
1143                fprintf(stdout,"\n\n\n");
1144                fprintf(stdout," ********* MYERS/MILLER MULTIPLE ALIGNMENT PARAMETERS *********\n");
1145                fprintf(stdout,"\n\n");
1146               
1147                fprintf(stdout,"     1. Fixed Gap Penalty       :%d\n",gap_open);
1148                fprintf(stdout,"     2. Floating Gap Penalty    :%d\n",gap_extend);
1149                fprintf(stdout,"     3. Toggle Transitions (DNA):%s\n",lin3);
1150                fprintf(stdout,"     4. Protein weight matrix   :%s\n\n"
1151                                        ,pam_matrix_name[matnum-1]);
1152                fprintf(stdout,"     H. HELP\n\n\n");           
1153
1154                getstr("Enter number (or [RETURN] to exit)",lin2, MAXLINE+1);
1155
1156                if(*lin2 == EOS) {
1157                        if(dnaflag) {
1158                                dna_gap_open    = gap_open;
1159                                dna_gap_extend  = gap_extend;
1160                        }
1161                        else {
1162                                prot_gap_open   = gap_open;
1163                                prot_gap_extend = gap_extend;
1164                        }
1165                        return;
1166                }
1167               
1168                switch(toupper(*lin2)) {
1169                        case '1':
1170                        fprintf(stdout,"Fixed Gap Penalty Currently: %d\n",gap_open);
1171                                gap_open=getint("Enter number",1,100,gap_open);
1172                                break;
1173                        case '2':
1174                                fprintf(stdout,"Floating Gap Penalty Currently: %d\n",gap_extend);
1175                                gap_extend=getint("Enter number",1,100,gap_extend);
1176                                break;
1177                        case '3':
1178                                is_weight ^= TRUE;
1179                                break;
1180                        case '4':
1181                                read_matrix();
1182                                break;
1183                        case '?':
1184                        case 'H':
1185                                get_help(4);
1186                                break;
1187                        default:
1188                                fprintf(stderr,"\n\nUnrecognised Command\n\n");
1189                                break;
1190                }
1191        }
1192}
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204static void read_matrix()
1205{       static char userfile[FILENAMELEN+1];
1206       
1207        while(TRUE)
1208        {
1209                fprintf(stdout,"\n\n\n");
1210                fprintf(stdout," ********* PROTEIN WEIGHT MATRIX MENU *********\n");
1211                fprintf(stdout,"\n\n");
1212               
1213               
1214                fprintf(stdout,"     1. %s\n",pam_matrix_name[0]);
1215                fprintf(stdout,"     2. %s\n",pam_matrix_name[1]);
1216                fprintf(stdout,"     3. %s\n",pam_matrix_name[2]);
1217                fprintf(stdout,"     4. %s\n\n",pam_matrix_name[3]);
1218                fprintf(stdout,"     H. HELP\n\n");
1219                fprintf(stdout,
1220"     -- Current matrix is the %s ",pam_matrix_name[matnum-1]);
1221                if(matnum == 4) fprintf(stdout,"(file = %s)",userfile);
1222                fprintf(stdout,"--\n");
1223               
1224                       
1225                getstr("\n\nEnter number (or [RETURN] to exit)",lin2, MAXLINE+1);
1226                if(*lin2 == EOS) return;
1227
1228                switch(toupper(*lin2))  {
1229                        case '1':
1230                                matptr=pam100mt;
1231                                make_pamo(0);
1232                                prot_gap_open   = 13;
1233                                prot_gap_extend = 13;
1234                                matnum=1;
1235                                break;
1236                        case '2':
1237                                matptr=pam250mt;
1238                                make_pamo(0);
1239                                prot_gap_open   = 10;
1240                                prot_gap_extend = 10;
1241                                matnum=2;
1242                                break;
1243                        case '3':
1244                                matptr=idmat;
1245                                make_pamo(0);
1246                                prot_gap_open   = 10;
1247                                prot_gap_extend = 10;
1248                                matnum=3;
1249                                break;
1250                        case '4':
1251                                if(user_mat(userfile)) matnum=4;
1252                                break;
1253                        case '?':
1254                        case 'H':
1255                                get_help(8);
1256                                break;
1257                        default:
1258                                fprintf(stderr,"\n\nUnrecognised Command\n\n");
1259                                break;
1260                }
1261        }
1262}
1263
1264
1265static Boolean user_mat(char *str)
1266{
1267        int i,j,idx,val;
1268        FILE *infile;
1269       
1270        if(usemenu)
1271            getstr("Enter name of the matrix file",lin2, MAXLINE+1);
1272        else
1273                strcpy(lin2,str);
1274
1275        if(*lin2 == EOS) return FALSE;
1276       
1277        if((infile=fopen(lin2,"r"))==NULL) {
1278                error("Cannot find matrix file [%s]",lin2);
1279                return FALSE;
1280        }
1281
1282        strcpy(str,lin2);
1283        strcpy(mtrxnam,lin2);
1284       
1285        idx=0;
1286        for(i=0;i<20;++i) 
1287                for(j=0;j<=i;++j) {
1288                        if( fscanf(infile,"%d",&val) == EOF) {
1289                                error("Input matrix has too few values");
1290                                return FALSE;
1291                        }
1292                        usermat[idx++]=(char)val;
1293                }
1294
1295        fclose(infile);
1296        matptr=usermat;
1297        make_pamo(0);
1298        return TRUE;
1299}
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310static void seq_input()
1311{
1312        int i;
1313       
1314        if(usemenu) {
1315fprintf(stdout,"\n\nSequences should all be in 1 file.\n"); 
1316fprintf(stdout,
1317"\n3 formats accepted:  NBRF/PIR, EMBL/SwissProt, Pearson (Fasta).\n");
1318fprintf(stdout,
1319"\nGCG users should use TOPIR to convert their sequence files before use.\n\n\n");
1320        }
1321
1322       
1323       nseqs = readseqs(1);        /* DES   1 is the first seq to be read */
1324       if(nseqs < 0)               /* file could not be opened */
1325           { 
1326               nseqs = 0;
1327               empty = TRUE;
1328           }
1329       else if(nseqs == 0)         /* no sequences */
1330           {
1331               error("No sequences in file!  Bad format?");
1332               empty = TRUE;
1333           }
1334       else if(nseqs == 1)
1335           {
1336               error("Only one sequence in file!");
1337               empty = TRUE;
1338               nseqs = 0;
1339           }
1340       else 
1341           {
1342                fprintf(stdout,"\nSequences assumed to be %s \n\n",
1343                        dnaflag?"DNA":"PROTEIN");
1344                for(i=1; i<=nseqs; i++) {
1345                        fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
1346                        i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
1347                }       
1348                        if(dnaflag) {
1349                                gap_open   = dna_gap_open;
1350                                gap_extend = dna_gap_extend;
1351                        }
1352                        else {
1353                                gap_open   = prot_gap_open;
1354                                gap_extend = prot_gap_extend;
1355                        }
1356                        empty=FALSE;
1357           }
1358       
1359}
1360
1361
1362
1363
1364
1365
1366
1367static void profile_input(int profile_no)   /* DES  read a profile   */
1368{                                           /* profile_no is 1 or 2  */
1369        int local_nseqs, i;
1370       
1371        if(profile_no == 2 && profile1_empty) 
1372           {
1373             error("You must read in profile number 1 first");
1374             return;
1375           }
1376
1377
1378    if(profile_no == 1)     /* for the 1st profile */
1379      {
1380       local_nseqs = readseqs(1); /* (1) means 1st seq to be read = no. 1 */
1381       if(local_nseqs < 0)               /* file could not be opened */
1382               return;
1383       else if(local_nseqs == 0)         /* no sequences  */
1384           {
1385               error("No sequences in file!  Bad format?");
1386               return;
1387           }
1388       else 
1389           {                            /* success; found some seqs. */
1390                nseqs = profile1_nseqs = local_nseqs;
1391                fprintf(stdout,"\nNo. of seqs=%d\n",nseqs);
1392                profile1_empty=FALSE;
1393                profile2_empty=TRUE;
1394           }
1395      }
1396    else
1397      {                         /* first seq to be read = profile1_nseqs + 1 */
1398       local_nseqs = readseqs(profile1_nseqs+1); 
1399       if(local_nseqs < 0)               /* file could not be opened */
1400               profile2_empty = TRUE;
1401       else if(local_nseqs == 0)         /* no sequences */
1402           {
1403               error("No sequences in file!  Bad format?");
1404               profile2_empty = TRUE;
1405           }
1406       else 
1407           {
1408                fprintf(stdout,"\nNo. of seqs in profile=%d\n",local_nseqs);
1409                nseqs = profile1_nseqs + local_nseqs;
1410                fprintf(stdout,"\nTotal no. of seqs     =%d\n",nseqs);
1411                profile2_empty=FALSE;
1412                empty = FALSE;
1413           }
1414
1415      }
1416       
1417        fprintf(stdout,"\nSequences assumed to be %s \n\n",
1418                dnaflag?"DNA":"PROTEIN");
1419        for(i=profile2_empty?1:profile1_nseqs+1; i<=nseqs; i++) {
1420                fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
1421                   i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
1422        }       
1423        if(dnaflag) {
1424                gap_open   = dna_gap_open;
1425                gap_extend = dna_gap_extend;
1426        }
1427        else {
1428                gap_open   = prot_gap_open;
1429                gap_extend = prot_gap_extend;
1430        }
1431}
1432
1433
1434
1435
1436
1437
1438FILE *  open_output_file(const char *prompt, const char *path,
1439                         char *file_name, const char *file_extension)
1440 
1441{       static char temp[FILENAMELEN+1];
1442        static char local_prompt[MAXLINE];
1443        FILE * file_handle;
1444
1445        strcpy(file_name,path);
1446        strcat(file_name,file_extension);
1447        strcpy(local_prompt,prompt);
1448        strcat(local_prompt," [%s]: ");         
1449
1450        if(usemenu) {
1451                fprintf(stdout,local_prompt,file_name);
1452                gets_noOverflow(temp, FILENAMELEN+1);
1453                if(*temp != EOS) strcpy(file_name,temp);
1454        }
1455
1456#if VMS
1457        if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
1458#else
1459        if((file_handle=fopen(file_name,"w"))==NULL) {
1460#endif
1461                error("Cannot open output file [%s]",file_name);
1462                return NULL;
1463        }
1464        return file_handle;
1465}
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477static void align()
1478{ 
1479        char path[FILENAMELEN+1];
1480        int oldmax,oldneut;
1481       
1482        if(empty && usemenu) {
1483                error("No sequences in memory. Load sequences first.");
1484                return;
1485        }
1486
1487        get_path(seqname,path);
1488       
1489        if(!open_alignment_output(path)) return;
1490
1491        if((tree = open_output_file(
1492                        "\nEnter a name for the dendrogram file     ",
1493                        path,treename,"dnd"))==NULL) return;
1494       
1495        fclose(tree);
1496       
1497        if(dnaflag) {
1498                ktup     = dna_ktup;
1499                window   = dna_window;
1500                signif   = dna_signif;
1501                wind_gap = dna_wind_gap;
1502        }
1503        else {
1504                ktup     = prot_ktup;
1505                window   = prot_window;
1506                signif   = prot_signif;
1507                wind_gap = prot_wind_gap;
1508        }
1509
1510        reset();
1511       
1512        fprintf(stdout,"\nStart of Pairwise alignments\n");
1513        fprintf(stdout,"Aligning...\n");
1514        show_pair();
1515        upgma(nseqs);
1516       
1517        if((tree=fopen(treename,"r"))==NULL) {
1518                error("Cannot open file [%s]",treename);
1519                return;
1520        }
1521       
1522        oldneut=xover;
1523        oldmax=big_pam;
1524        myers(0);
1525        big_pam=oldmax;
1526        xover=oldneut;
1527        fclose(tree);
1528       
1529        fprintf(stdout,"\n\n\n");
1530        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1531       
1532        create_alignment_output();
1533}
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543static void make_tree()
1544{
1545        char path[FILENAMELEN+1],temp[FILENAMELEN+1];
1546       
1547        if(empty) {
1548                error("No sequences in memory. Load sequences first.");
1549                return;
1550        }
1551
1552        get_path(seqname,path);
1553
1554        strcpy(treename,path);
1555        strcat(treename,"dnd");
1556       
1557        fprintf(stdout,"\nEnter a name for the DENDROGRAM file [%s]: ",treename);
1558        gets_noOverflow(temp, FILENAMELEN+1);
1559        if(*temp != EOS)
1560                strcpy(treename,temp);
1561       
1562#if VMS
1563        if((tree=fopen(treename,"w","rat=cr","rfm=var"))==NULL) {
1564#else
1565        if((tree=fopen(treename,"w"))==NULL) {
1566#endif
1567                error("Cannot open file [%s]",treename);
1568                return;
1569        }
1570        fclose(tree);
1571       
1572        if(dnaflag) {
1573                ktup     = dna_ktup;
1574                window   = dna_window;
1575                signif   = dna_signif;
1576                wind_gap = dna_wind_gap;
1577        }
1578        else {
1579                ktup     = prot_ktup;
1580                window   = prot_window;
1581                signif   = prot_signif;
1582                wind_gap = prot_wind_gap;
1583        }
1584
1585        reset();
1586        fprintf(stdout,"\nStart of Pairwise alignments\n");
1587        fprintf(stdout,"Aligning...\n");
1588        show_pair();
1589        upgma(nseqs);
1590       
1591        fprintf(stdout,"\nDENDROGRAM file created [%s]\n",treename);
1592}
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602static void get_tree()
1603{
1604        char path[FILENAMELEN+1],temp[MAXLINE+1];
1605        int count,oldmax,oldneut;
1606       
1607        if(empty) {
1608                error("No sequences in memory. Load sequences first.");
1609                return;
1610        }
1611
1612        get_path(seqname,path);
1613       
1614        strcpy(treename,path);
1615        strcat(treename,"dnd");
1616       
1617        fprintf(stdout,"\nEnter a name for the DENDROGRAM file      [%s]:",treename);
1618        gets_noOverflow(temp, MAXLINE+1);
1619        if(*temp != EOS)
1620                strcpy(treename,temp);
1621       
1622        if((tree=fopen(treename,"r"))==NULL) {
1623                error("Cannot open file [%s]",treename);
1624                return;
1625        }
1626       
1627        count=0;
1628       
1629        while(fgets(temp,MAXLINE+1,tree)!=NULL) ++count;
1630       
1631        fclose(tree);
1632       
1633        if(++count != nseqs) {
1634                error("Dendrogram file is not consistent with the sequence data");
1635                return;
1636        }
1637       
1638        if(!open_alignment_output(path)) return;
1639
1640        reset();
1641       
1642        if((tree=fopen(treename,"r"))==NULL) {
1643                error("Cannot open file [%s]",treename);
1644                return;
1645        }
1646       
1647        oldmax=big_pam;
1648        oldneut=xover;
1649        myers(0);
1650        xover=oldneut;
1651        big_pam=oldmax;
1652        fclose(tree);
1653
1654        fprintf(stdout,"\n\n\n");
1655        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1656
1657        create_alignment_output();
1658}
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669static void profile_align()
1670{
1671        char path[FILENAMELEN+1];
1672        int oldmax,oldneut;
1673       
1674        if(profile2_empty) {
1675                error("No sequences in memory. Load sequences first.");
1676                return;
1677        }
1678
1679        get_path(seqname,path);
1680       
1681        if(!open_alignment_output(path)) return;
1682
1683        reset();
1684               
1685        oldmax=big_pam;
1686        oldneut=xover;
1687        myers(1);
1688        xover=oldneut;
1689        big_pam=oldmax;
1690
1691        fprintf(stdout,"\n\n\n");
1692        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1693
1694        create_alignment_output();
1695}
1696
1697
1698
1699
1700
1701
1702
1703
1704static void clustal_out(FILE *clusout)
1705{
1706        static char seq1[MAXLEN+1];
1707        char        temp[MAXLINE];
1708        int val,i,j,k,len;
1709        int chunks,ident,lv1,pos,ptr,copt,flag;
1710
1711
1712        fprintf(clusout,"CLUSTAL V multiple sequence alignment\n\n\n");
1713
1714        len=seqlen_array[1];
1715       
1716        chunks = len/LINELENGTH;
1717        if(len % LINELENGTH != 0)
1718                ++chunks;
1719               
1720        for(lv1=1;lv1<=chunks;++lv1) {
1721                pos = ((lv1-1)*LINELENGTH)+1;
1722                ptr = (len<pos+LINELENGTH-1) ? len : pos+LINELENGTH-1;
1723                for(i=1;i<=nseqs;++i) {
1724                        for(j=pos;j<=ptr;++j) {
1725                                val=seq_array[i][j];
1726                                if(val==0 || val>20)
1727                                        seq1[j]='X';
1728                                else
1729                                        if(val<0) seq1[j]='-';
1730                                else {
1731                                        if(dnaflag)
1732                                                seq1[j]=nucleic_acid_order[val];
1733                                        else
1734                                                seq1[j]=amino_acid_order[val];
1735                                }
1736                        }
1737                        strncpy(temp,&seq1[pos],ptr-pos+1);
1738                        temp[ptr-pos+1]=EOS;
1739                        fprintf(clusout,"%-15s %s\n",names[i],temp);
1740                }
1741       
1742                copt = ((nseqs*nseqs) - nseqs) / 2;
1743                for(i=pos;i<=ptr;++i) {
1744                        seq1[i]=' ';
1745                        ident=0;
1746                        for(j=1;j<=nseqs;++j)
1747                                if(seq_array[1][i] == seq_array[j][i])
1748                                        ++ident;
1749                        if(ident==nseqs)
1750                                seq1[i]='*';
1751                        else {
1752                                if(!dnaflag) {
1753                                        ident=flag=0;
1754                                        for(j=1;j<=nseqs;++j) {
1755                                                for(k=j+1;k<=nseqs;++k)
1756                                                        if(seq_array[j][i]>0 && seq_array[k][i]>0) {
1757                                                                if(weights[seq_array[j][i]][seq_array[k][i]]<xover)
1758                                                                        ++ident;
1759                                                                else {
1760                                                                        flag=TRUE;
1761                                                                        break;
1762                                                                }
1763                                                        }
1764                                                        else {
1765                                                                flag=TRUE;
1766                                                                break;
1767                                                        }
1768                                                        if(flag)
1769                                                                break;
1770                                        }
1771                                        if(flag)
1772                                                continue;
1773                                        if(ident==copt)
1774                                                seq1[i]='.';
1775                                }
1776                        }
1777                }
1778                strncpy(temp,&seq1[pos],ptr-pos+1);
1779                temp[ptr-pos+1]=EOS;
1780                fprintf(clusout,"                %s\n\n",temp);
1781                ++nblocks;
1782        }
1783       
1784/*
1785        fprintf(stdout,"\nEnd of Multiple Alignment\n\n");
1786        fprintf(clusout,"\nEnd of Multiple Alignment\n\n");
1787        fprintf(clusout,"\nKey:\n     Identity:     #\n     Conservative: ^\n\n");
1788
1789        fprintf(clusout,"\nParameters:\n\n");
1790        fprintf(clusout,"Ktup:         %d\nWilbur Gap:   %d\nCutoff:       %d\nDiagonal:     %d\n",
1791        ktup,wind_gap,signif,window);
1792        fprintf(clusout,"Fixed Gap:    %d\nFloating Gap: %d\n\n",gap_open,gap_extend);
1793
1794        fprintf(clusout,"Sequences were: ");
1795        fprintf(clusout,dnaflag ? "Nucleic Acid\n" : "Proteins\n");
1796        fprintf(clusout,percent ? "Percentage" : "Absolute");
1797        fprintf(clusout," Identities Were Scored\n");
1798        if(dnaflag) {
1799                fprintf(clusout,"Nucleotide Transitions: ");
1800                fprintf(clusout,is_weight ? "WEIGHTED\n" : "UNWEIGHTED\n");
1801        }
1802        fprintf(clusout,"Sequence Input file: %s\n",seqname);
1803        fprintf(clusout,"Matrix used: ");
1804        if(matptr==idmat)
1805                fprintf(clusout,"No penalty\n");
1806        else
1807                if(matptr==pam250mt)
1808                        fprintf(clusout,"PAM 250\n");
1809        else
1810                if(matptr==pam100mt)
1811                        fprintf(clusout,"PAM 100\n");
1812        else
1813                fprintf(clusout,"%s\n",mtrxnam);
1814*/
1815}
1816
1817
1818
1819
1820
1821
1822static void gcg_out(FILE *gcgout)
1823{
1824/*        static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1825/*      static char *nbases = "XACGT";  */
1826        char seq[MAXLEN+1], residue;
1827        int all_checks[MAXN+1];
1828        int i,j,k,len,val,chunks,block,pos1,pos2;       
1829        long grand_checksum;
1830       
1831        len = seqlen_array[1];
1832
1833        for(i=1; i<=nseqs; i++) {
1834                for(j=1; j<=len; j++) {
1835                        val = seq_array[i][j];
1836                        if(val == 0 || val > 20) 
1837                                residue = 'X';
1838                        else if(val < 0) 
1839                                residue = '.';
1840                        else {
1841                                if(dnaflag)
1842                                        residue = nucleic_acid_order[val];
1843                                else
1844                                        residue = amino_acid_order[val];
1845                        }
1846                        seq[j] = residue;
1847                }
1848                all_checks[i] = SeqGCGCheckSum(seq+1, len);
1849        }       
1850
1851        grand_checksum = 0;
1852        for(i=1; i<=nseqs; i++) grand_checksum += all_checks[i];
1853        grand_checksum = grand_checksum % 10000;
1854        fprintf(gcgout,"\n\n   MSF:%5d  Type: ",len);
1855        if(dnaflag)
1856                fprintf(gcgout,"N");
1857        else
1858                fprintf(gcgout,"P");
1859        fprintf(gcgout,"    Check:%6ld   .. \n\n", grand_checksum);
1860        for(i=1; i<=nseqs; i++) 
1861/*              for(j=0; j<MAXNAMES; j++)
1862                        if(names[i][j] == ' ')  names[i][j] = '_';     */
1863                fprintf(gcgout,
1864                        " Name: %-15s oo  Len:%5d  Check:%6d  Weight:  1.00\n",
1865                        names[i],len,all_checks[i]);
1866        fprintf(gcgout,"\n//\n"); 
1867
1868        chunks = len/GCG_LINELENGTH;
1869        if(len % GCG_LINELENGTH != 0) ++chunks;
1870
1871        for(block=1; block<=chunks; block++) {
1872                fprintf(gcgout,"\n\n");
1873                pos1 = ((block-1) * GCG_LINELENGTH) + 1;
1874                pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
1875                for(i=1; i<=nseqs; i++) {
1876                        fprintf(gcgout,"\n%-15s ",names[i]);
1877                        for(j=pos1, k=1; j<=pos2; j++, k++) {
1878                                val = seq_array[i][j];
1879                                if(val == 0 || val > 20) 
1880                                        residue = 'X';
1881                                else if(val < 0) 
1882                                        residue = '.';
1883                                else {
1884                                        if(dnaflag)
1885                                                residue = nucleic_acid_order[val];
1886                                        else
1887                                                residue = amino_acid_order[val];
1888                                }
1889                                fprintf(gcgout,"%c",residue);
1890                                if(j % 10 == 0) fprintf(gcgout," ");
1891                        }
1892                }
1893        }
1894        fprintf(gcgout,"\n\n");
1895}
1896
1897
1898
1899static void phylip_out(FILE *phyout)
1900{
1901/*      static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1902/*              static char *nbases = "XACGT";  */
1903        char residue;
1904        int i,j,k,len,val,chunks,block,pos1,pos2;       
1905       
1906        len = seqlen_array[1];
1907
1908
1909        chunks = len/GCG_LINELENGTH;
1910        if(len % GCG_LINELENGTH != 0) ++chunks;
1911
1912        fprintf(phyout,"%6d %6d",nseqs,len);
1913
1914        for(block=1; block<=chunks; block++) {
1915                pos1 = ((block-1) * GCG_LINELENGTH) + 1;
1916                pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
1917                for(i=1; i<=nseqs; i++) {
1918                        if(block == 1) 
1919                                fprintf(phyout,"\n%-10s ",names[i]);
1920                        else
1921                                fprintf(phyout,"\n           ");
1922                        for(j=pos1, k=1; j<=pos2; j++, k++) {
1923                                val = seq_array[i][j];
1924                                if(val == 0 || val > 20) 
1925                                        residue = 'X';
1926                                else if(val < 0) 
1927                                        residue = '-';
1928                                else {
1929                                        if(dnaflag)
1930                                                residue = nucleic_acid_order[val];
1931                                        else
1932                                                residue = amino_acid_order[val];
1933                                }
1934                                fprintf(phyout,"%c",residue);
1935                                if(j % 10 == 0) fprintf(phyout," ");
1936                        }
1937                }
1938                fprintf(phyout,"\n");
1939        }
1940}
1941
1942
1943
1944
1945
1946static void nbrf_out(FILE *nbout)
1947{
1948/*      static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1949/*              static char *nbases = "XACGT";  */
1950        char seq[MAXLEN+1], residue;
1951        int i,j,len,val;       
1952       
1953        len = seqlen_array[1];
1954
1955        for(i=1; i<=nseqs; i++) {
1956                fprintf(nbout, dnaflag ? ">DL;" : ">P1;");
1957                fprintf(nbout, "%s\n%s\n", names[i], titles[i]);
1958                for(j=1; j<=len; j++) {
1959                        val = seq_array[i][j];
1960                        if(val == 0 || val > 20) 
1961                                residue = 'X';
1962                        else if(val < 0) 
1963                                residue = '-';
1964                        else {
1965                                if(dnaflag)
1966                                        residue = nucleic_acid_order[val];
1967                                else
1968                                        residue = amino_acid_order[val];
1969                        }
1970                        seq[j] = residue;
1971                }
1972                for(j=1; j<=len; j++) {
1973                        fprintf(nbout,"%c",seq[j]);
1974                        if((j % LINELENGTH == 0) || (j == len)) 
1975                                fprintf(nbout,"\n");
1976                }
1977                fprintf(nbout,"*\n");
1978        }       
1979}
1980
1981
1982static Boolean open_alignment_output(char *path)
1983{
1984        if(output_clustal) 
1985                if((clustal_outfile = open_output_file(
1986                        "\nEnter a name for the CLUSTAL output file ",path,
1987                        clustal_outname,"aln"))==NULL) return FALSE;
1988        if(output_nbrf) 
1989                if((nbrf_outfile = open_output_file(
1990                        "\nEnter a name for the NBRF/PIR output file",path,
1991                        nbrf_outname,"pir"))==NULL) return FALSE;
1992        if(output_gcg) 
1993                if((gcg_outfile = open_output_file(
1994                        "\nEnter a name for the GCG output file     ",path,
1995                        gcg_outname,"msf"))==NULL) return FALSE;
1996        if(output_phylip) 
1997                if((phylip_outfile = open_output_file(
1998                        "\nEnter a name for the PHYLIP output file  ",path,
1999                        phylip_outname,"phy"))==NULL) return FALSE;
2000        return TRUE;
2001}
2002
2003
2004
2005
2006
2007static void create_alignment_output()
2008{
2009        if(output_clustal) {
2010                clustal_out(clustal_outfile);
2011                fclose(clustal_outfile);
2012                fprintf(stdout,"\nCLUSTAL-Alignment file created  [%s]\n",clustal_outname);
2013        }
2014        if(output_nbrf)  {
2015                nbrf_out(nbrf_outfile);
2016                fclose(nbrf_outfile);
2017                fprintf(stdout,"\nNBRF/PIR-Alignment file created [%s]\n",nbrf_outname);
2018        }
2019        if(output_gcg)  {
2020                gcg_out(gcg_outfile);
2021                fclose(gcg_outfile);
2022                fprintf(stdout,"\nGCG-Alignment file created      [%s]\n",gcg_outname);
2023        }
2024        if(output_phylip)  {
2025                phylip_out(phylip_outfile);
2026                fclose(phylip_outfile);
2027                fprintf(stdout,"\nPHYLIP-Alignment file created   [%s]\n",phylip_outname);
2028        }
2029}
2030
2031
2032
2033
2034
2035
2036static void reset()   /* remove gaps from older alignments (code = -1) */
2037{                     /* EXCEPT for gaps that were INPUT with the seqs.*/
2038        register int i,j,sl;                 /* which have  code = -2  */
2039       
2040        for(i=1;i<=nseqs;++i) {
2041                sl=0;
2042                for(j=1;j<=seqlen_array[i];++j) {
2043                        if(seq_array[i][j] == -1) continue;
2044                        ++sl;
2045                        seq_array[i][sl]=seq_array[i][j];
2046                }
2047                seqlen_array[i]=sl;
2048        }
2049}
2050
2051
Note: See TracBrowser for help on using the repository browser.