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

Last change on this file was 13845, checked in by westram, 10 years ago
  • remove executable flag
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.6 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
10/*
11*       Prototypes
12*/
13
14extern void     getstr(char *,char *);
15extern double   getreal(char *,double,double);
16extern int              getint(char *,int,int,int);
17extern void             do_system(void);
18extern void             make_pamo(int);
19extern int          readseqs(int);
20extern void             get_path(char *,char *);
21extern void             show_pair(void);
22extern void             upgma(int);
23extern void             myers(int);
24extern void     phylogenetic_tree(void);
25extern void         bootstrap_tree(void);
26extern void             error(char *,...);
27extern int              SeqGCGCheckSum(char *, int);
28
29
30void init_amenu(void);
31void parse_params(void);
32void main_menu(void);
33FILE * open_output_file(char *, char *, char *, char *);
34#if UNIX
35FILE * open_path(char *);
36#endif
37
38static Boolean check_param(char *, char *, char *, int *);
39static void get_help(int);                      /* Help procedure */
40static void pair_menu(void);
41static void multi_menu(void);
42static void multiple_align_menu(void);          /* multiple alignments menu */
43static void profile_align_menu(void);           /* profile       "      "   */
44static void phylogenetic_tree_menu(void);       /* NJ trees/distances menu  */
45static void read_matrix(void);
46static void seq_input(void);
47static void align(void);
48static void make_tree(void);
49static void get_tree(void);
50static Boolean user_mat(char *);
51static void clustal_out(FILE *);
52static void nbrf_out(FILE *);
53static void gcg_out(FILE *);
54static void phylip_out(FILE *);
55static Boolean open_alignment_output(char *);
56static void create_alignment_output(void);
57static void reset(void);
58static void profile_input(int);                 /* DES read a profile */
59static void format_options_menu(void);          /* format of alignment output */
60static void profile_align(void);                /* Align 2 alignments */
61
62/*
63*        Global variables
64*/
65
66char *lin1, *lin2, *lin3;
67
68extern char *amino_acid_order;
69extern char *nucleic_acid_order;
70extern void *ckalloc(size_t);
71extern Boolean          percent,is_weight,dnaflag;
72extern int              wind_gap,ktup,window,signif;
73extern int               dna_wind_gap, dna_ktup, dna_window, dna_signif;
74extern int              prot_wind_gap,prot_ktup,prot_window,prot_signif;
75extern unsigned int boot_ran_seed;
76extern int      gap_open,      gap_extend;
77extern int  dna_gap_open,  dna_gap_extend;
78extern int prot_gap_open, prot_gap_extend;
79extern int boot_ntrials;                /* number of bootstrap trials */
80extern int              xover,big_pam;
81extern char *params;
82extern char             seqname[],treename[];
83extern char     *matptr,pam100mt[],pam250mt[],idmat[];
84extern int              nseqs,nblocks;
85extern int              weights[21][21];
86extern FILE *   tree;
87extern FILE *clustal_outfile, *gcg_outfile, *nbrf_outfile, *phylip_outfile;
88extern char **  names, **titles;
89extern int *seqlen_array;
90extern char **seq_array;
91extern char ntrials;            /* number of bootstrap trials (trees.c) */
92
93Boolean usemenu;
94char mtrxnam[FILENAMELEN+1];
95Boolean explicit_dnaflag;  /* Explicit setting of sequence type on comm.line*/
96Boolean output_clustal, output_nbrf, output_phylip, output_gcg;
97Boolean empty;
98Boolean profile1_empty, profile2_empty;   /* whether or not profiles   */
99int profile1_nseqs;            /* have been filled; the no. of seqs in prof 1*/
100Boolean tossgaps, kimura;
101int  matnum;
102static char clustal_outname[FILENAMELEN+1], gcg_outname[FILENAMELEN+1];
103static char  phylip_outname[FILENAMELEN+1],nbrf_outname[FILENAMELEN+1];
104
105static char *pam_matrix_name[] = {
106                "PAM 100",
107                "PAM 250",
108                "Identity matrix",
109                "user defined"   };
110               
111char usermat[210];
112
113void init_amenu()
114{
115        matnum=2;
116        empty=TRUE;
117        explicit_dnaflag = FALSE;     
118    profile1_empty=TRUE;     /* DES */
119    profile2_empty=TRUE;     /* DES */
120        output_clustal = TRUE;
121        output_gcg     = FALSE;
122        output_phylip  = FALSE;
123        output_nbrf    = FALSE;
124
125    lin1 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );
126    lin2 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );
127    lin3 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );
128}
129
130
131
132
133static Boolean check_param(char *paramstr, 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(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        char *digits = "0123456789";
215        char *help_marker    = ">>HELP<<";
216
217#if MSDOS
218        char *help_file_name = "clustalv.hlp";
219#else
220        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);
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);
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);
283                        }
284                        fclose(help_file);
285                }
286        }
287
288}
289
290
291void parse_params()
292{
293        int i,j,k,len,lenp,temp;
294        char probe[80], 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 char *fixedgapst         = "fixedgap";
301        static char *floatgapst         = "floatgap";
302        static char *kimurast           = "kimura";
303        static char *ktuplest           = "ktuple";
304        static char *matrixst           = "matrix";
305        static char *outputst           = "output";
306        static char *pairgapst          = "pairgap";
307        static char *scorest            = "score";
308        static char *seedst             = "seed";
309        static char *topdiagsst         = "topdiags";
310        static char *tossgapsst         = "tossgaps";
311        static char *transitionsst      = "transitions";
312        static char *typest             = "type";
313        static char *windowst           = "window";
314
315/* command line switches for DATA       **************************/
316        static char *infilest           = "infile";
317        static char *profile1st         = "profile1";
318        static char *profile2st         = "profile2";
319
320/* command line switches for VERBS      **************************/
321        static char *alignst            = "align";
322        static char *bootstrapst        = "bootstrap";
323        static char *checkst            = "check";        /*  /check = /help */
324        static char *helpst             = "help";
325        static 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);
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);
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);
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);
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);
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                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);
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                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);
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);
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,nv,pos,idx,val;
1268        FILE *infile;
1269       
1270        if(usemenu)
1271                getstr("Enter name of the matrix file",lin2);
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        char c;
1313        int i;
1314       
1315        if(usemenu) {
1316fprintf(stdout,"\n\nSequences should all be in 1 file.\n"); 
1317fprintf(stdout,
1318"\n3 formats accepted:  NBRF/PIR, EMBL/SwissProt, Pearson (Fasta).\n");
1319fprintf(stdout,
1320"\nGCG users should use TOPIR to convert their sequence files before use.\n\n\n");
1321        }
1322
1323       
1324       nseqs = readseqs(1);        /* DES   1 is the first seq to be read */
1325       if(nseqs < 0)               /* file could not be opened */
1326           { 
1327               nseqs = 0;
1328               empty = TRUE;
1329           }
1330       else if(nseqs == 0)         /* no sequences */
1331           {
1332               error("No sequences in file!  Bad format?");
1333               empty = TRUE;
1334           }
1335       else if(nseqs == 1)
1336           {
1337               error("Only one sequence in file!");
1338               empty = TRUE;
1339               nseqs = 0;
1340           }
1341       else 
1342           {
1343                fprintf(stdout,"\nSequences assumed to be %s \n\n",
1344                        dnaflag?"DNA":"PROTEIN");
1345                for(i=1; i<=nseqs; i++) {
1346                        fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
1347                        i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
1348                }       
1349                        if(dnaflag) {
1350                                gap_open   = dna_gap_open;
1351                                gap_extend = dna_gap_extend;
1352                        }
1353                        else {
1354                                gap_open   = prot_gap_open;
1355                                gap_extend = prot_gap_extend;
1356                        }
1357                        empty=FALSE;
1358           }
1359       
1360}
1361
1362
1363
1364
1365
1366
1367
1368static void profile_input(int profile_no)   /* DES  read a profile   */
1369{                                           /* profile_no is 1 or 2  */
1370        char c;
1371        int local_nseqs, i;
1372       
1373        if(profile_no == 2 && profile1_empty) 
1374           {
1375             error("You must read in profile number 1 first");
1376             return;
1377           }
1378
1379
1380    if(profile_no == 1)     /* for the 1st profile */
1381      {
1382       local_nseqs = readseqs(1); /* (1) means 1st seq to be read = no. 1 */
1383       if(local_nseqs < 0)               /* file could not be opened */
1384               return;
1385       else if(local_nseqs == 0)         /* no sequences  */
1386           {
1387               error("No sequences in file!  Bad format?");
1388               return;
1389           }
1390       else 
1391           {                            /* success; found some seqs. */
1392                nseqs = profile1_nseqs = local_nseqs;
1393                fprintf(stdout,"\nNo. of seqs=%d\n",nseqs);
1394                profile1_empty=FALSE;
1395                profile2_empty=TRUE;
1396           }
1397      }
1398    else
1399      {                         /* first seq to be read = profile1_nseqs + 1 */
1400       local_nseqs = readseqs(profile1_nseqs+1); 
1401       if(local_nseqs < 0)               /* file could not be opened */
1402               profile2_empty = TRUE;
1403       else if(local_nseqs == 0)         /* no sequences */
1404           {
1405               error("No sequences in file!  Bad format?");
1406               profile2_empty = TRUE;
1407           }
1408       else 
1409           {
1410                fprintf(stdout,"\nNo. of seqs in profile=%d\n",local_nseqs);
1411                nseqs = profile1_nseqs + local_nseqs;
1412                fprintf(stdout,"\nTotal no. of seqs     =%d\n",nseqs);
1413                profile2_empty=FALSE;
1414                empty = FALSE;
1415           }
1416
1417      }
1418       
1419        fprintf(stdout,"\nSequences assumed to be %s \n\n",
1420                dnaflag?"DNA":"PROTEIN");
1421        for(i=profile2_empty?1:profile1_nseqs+1; i<=nseqs; i++) {
1422                fprintf(stdout,"Sequence %d: %-*.s   %6.d %s\n",
1423                   i,MAXNAMES,names[i],seqlen_array[i],dnaflag?"bp":"aa");
1424        }       
1425        if(dnaflag) {
1426                gap_open   = dna_gap_open;
1427                gap_extend = dna_gap_extend;
1428        }
1429        else {
1430                gap_open   = prot_gap_open;
1431                gap_extend = prot_gap_extend;
1432        }
1433}
1434
1435
1436
1437
1438
1439
1440FILE *  open_output_file(char *prompt,      char *path, 
1441                                char *file_name,   char *file_extension)
1442 
1443{       static char temp[FILENAMELEN+1];
1444        static char local_prompt[MAXLINE];
1445        FILE * file_handle;
1446
1447        strcpy(file_name,path);
1448        strcat(file_name,file_extension);
1449        strcpy(local_prompt,prompt);
1450        strcat(local_prompt," [%s]: ");         
1451
1452        if(usemenu) {
1453                fprintf(stdout,local_prompt,file_name);
1454                gets(temp);
1455                if(*temp != EOS) strcpy(file_name,temp);
1456        }
1457
1458#if VMS
1459        if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
1460#else
1461        if((file_handle=fopen(file_name,"w"))==NULL) {
1462#endif
1463                error("Cannot open output file [%s]",file_name);
1464                return NULL;
1465        }
1466        return file_handle;
1467}
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479static void align()
1480{ 
1481        char path[FILENAMELEN+1],temp[FILENAMELEN+1];
1482        int oldmax,oldneut;
1483       
1484        if(empty && usemenu) {
1485                error("No sequences in memory. Load sequences first.");
1486                return;
1487        }
1488
1489        get_path(seqname,path);
1490       
1491        if(!open_alignment_output(path)) return;
1492
1493        if((tree = open_output_file(
1494                        "\nEnter a name for the dendrogram file     ",
1495                        path,treename,"dnd"))==NULL) return;
1496       
1497        fclose(tree);
1498       
1499        if(dnaflag) {
1500                ktup     = dna_ktup;
1501                window   = dna_window;
1502                signif   = dna_signif;
1503                wind_gap = dna_wind_gap;
1504        }
1505        else {
1506                ktup     = prot_ktup;
1507                window   = prot_window;
1508                signif   = prot_signif;
1509                wind_gap = prot_wind_gap;
1510        }
1511
1512        reset();
1513       
1514        fprintf(stdout,"\nStart of Pairwise alignments\n");
1515        fprintf(stdout,"Aligning...\n");
1516        show_pair();
1517        upgma(nseqs);
1518       
1519        if((tree=fopen(treename,"r"))==NULL) {
1520                error("Cannot open file [%s]",treename);
1521                return;
1522        }
1523       
1524        oldneut=xover;
1525        oldmax=big_pam;
1526        myers(0);
1527        big_pam=oldmax;
1528        xover=oldneut;
1529        fclose(tree);
1530       
1531        fprintf(stdout,"\n\n\n");
1532        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1533       
1534        create_alignment_output();
1535}
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545static void make_tree()
1546{
1547        char path[FILENAMELEN+1],temp[FILENAMELEN+1];
1548       
1549        if(empty) {
1550                error("No sequences in memory. Load sequences first.");
1551                return;
1552        }
1553
1554        get_path(seqname,path);
1555
1556        strcpy(treename,path);
1557        strcat(treename,"dnd");
1558       
1559        fprintf(stdout,"\nEnter a name for the DENDROGRAM file [%s]: ",treename);
1560        gets(temp);
1561        if(*temp != EOS)
1562                strcpy(treename,temp);
1563       
1564#if VMS
1565        if((tree=fopen(treename,"w","rat=cr","rfm=var"))==NULL) {
1566#else
1567        if((tree=fopen(treename,"w"))==NULL) {
1568#endif
1569                error("Cannot open file [%s]",treename);
1570                return;
1571        }
1572        fclose(tree);
1573       
1574        if(dnaflag) {
1575                ktup     = dna_ktup;
1576                window   = dna_window;
1577                signif   = dna_signif;
1578                wind_gap = dna_wind_gap;
1579        }
1580        else {
1581                ktup     = prot_ktup;
1582                window   = prot_window;
1583                signif   = prot_signif;
1584                wind_gap = prot_wind_gap;
1585        }
1586
1587        reset();
1588        fprintf(stdout,"\nStart of Pairwise alignments\n");
1589        fprintf(stdout,"Aligning...\n");
1590        show_pair();
1591        upgma(nseqs);
1592       
1593        fprintf(stdout,"\nDENDROGRAM file created [%s]\n",treename);
1594}
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604static void get_tree()
1605{
1606        char path[FILENAMELEN+1],temp[MAXLINE+1];
1607        int count,oldmax,oldneut;
1608       
1609        if(empty) {
1610                error("No sequences in memory. Load sequences first.");
1611                return;
1612        }
1613
1614        get_path(seqname,path);
1615       
1616        strcpy(treename,path);
1617        strcat(treename,"dnd");
1618       
1619        fprintf(stdout,"\nEnter a name for the DENDROGRAM file      [%s]:",treename);
1620        gets(temp);
1621        if(*temp != EOS)
1622                strcpy(treename,temp);
1623       
1624        if((tree=fopen(treename,"r"))==NULL) {
1625                error("Cannot open file [%s]",treename);
1626                return;
1627        }
1628       
1629        count=0;
1630       
1631        while(fgets(temp,MAXLINE+1,tree)!=NULL) ++count;
1632       
1633        fclose(tree);
1634       
1635        if(++count != nseqs) {
1636                error("Dendrogram file is not consistent with the sequence data");
1637                return;
1638        }
1639       
1640        if(!open_alignment_output(path)) return;
1641
1642        reset();
1643       
1644        if((tree=fopen(treename,"r"))==NULL) {
1645                error("Cannot open file [%s]",treename);
1646                return;
1647        }
1648       
1649        oldmax=big_pam;
1650        oldneut=xover;
1651        myers(0);
1652        xover=oldneut;
1653        big_pam=oldmax;
1654        fclose(tree);
1655
1656        fprintf(stdout,"\n\n\n");
1657        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1658
1659        create_alignment_output();
1660}
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671static void profile_align()
1672{
1673        char path[FILENAMELEN+1],temp[MAXLINE+1];
1674        int count,oldmax,oldneut;
1675       
1676        if(profile2_empty) {
1677                error("No sequences in memory. Load sequences first.");
1678                return;
1679        }
1680
1681        get_path(seqname,path);
1682       
1683        if(!open_alignment_output(path)) return;
1684
1685        reset();
1686               
1687        oldmax=big_pam;
1688        oldneut=xover;
1689        myers(1);
1690        xover=oldneut;
1691        big_pam=oldmax;
1692
1693        fprintf(stdout,"\n\n\n");
1694        fprintf(stdout,"Consensus length = %d\n",seqlen_array[1]);
1695
1696        create_alignment_output();
1697}
1698
1699
1700
1701
1702
1703
1704
1705
1706static void clustal_out(FILE *clusout)
1707{
1708        static char seq1[MAXLEN+1];
1709        char        temp[MAXLINE];
1710        int val,i,j,k,a,b,len;
1711        int chunks,ident,lv1,pos,ptr,copt,flag;
1712
1713
1714        fprintf(clusout,"CLUSTAL V multiple sequence alignment\n\n\n");
1715
1716        len=seqlen_array[1];
1717       
1718        chunks = len/LINELENGTH;
1719        if(len % LINELENGTH != 0)
1720                ++chunks;
1721               
1722        for(lv1=1;lv1<=chunks;++lv1) {
1723                pos = ((lv1-1)*LINELENGTH)+1;
1724                ptr = (len<pos+LINELENGTH-1) ? len : pos+LINELENGTH-1;
1725                for(i=1;i<=nseqs;++i) {
1726                        for(j=pos;j<=ptr;++j) {
1727                                val=seq_array[i][j];
1728                                if(val==0 || val>20)
1729                                        seq1[j]='X';
1730                                else
1731                                        if(val<0) seq1[j]='-';
1732                                else {
1733                                        if(dnaflag)
1734                                                seq1[j]=nucleic_acid_order[val];
1735                                        else
1736                                                seq1[j]=amino_acid_order[val];
1737                                }
1738                        }
1739                        strncpy(temp,&seq1[pos],ptr-pos+1);
1740                        temp[ptr-pos+1]=EOS;
1741                        fprintf(clusout,"%-15s %s\n",names[i],temp);
1742                }
1743       
1744                copt = ((nseqs*nseqs) - nseqs) / 2;
1745                for(i=pos;i<=ptr;++i) {
1746                        seq1[i]=' ';
1747                        ident=0;
1748                        for(j=1;j<=nseqs;++j)
1749                                if(seq_array[1][i] == seq_array[j][i])
1750                                        ++ident;
1751                        if(ident==nseqs)
1752                                seq1[i]='*';
1753                        else {
1754                                if(!dnaflag) {
1755                                        ident=flag=0;
1756                                        for(j=1;j<=nseqs;++j) {
1757                                                for(k=j+1;k<=nseqs;++k)
1758                                                        if(seq_array[j][i]>0 && seq_array[k][i]>0) {
1759                                                                if(weights[seq_array[j][i]][seq_array[k][i]]<xover)
1760                                                                        ++ident;
1761                                                                else {
1762                                                                        flag=TRUE;
1763                                                                        break;
1764                                                                }
1765                                                        }
1766                                                        else {
1767                                                                flag=TRUE;
1768                                                                break;
1769                                                        }
1770                                                        if(flag)
1771                                                                break;
1772                                        }
1773                                        if(flag)
1774                                                continue;
1775                                        if(ident==copt)
1776                                                seq1[i]='.';
1777                                }
1778                        }
1779                }
1780                strncpy(temp,&seq1[pos],ptr-pos+1);
1781                temp[ptr-pos+1]=EOS;
1782                fprintf(clusout,"                %s\n\n",temp);
1783                ++nblocks;
1784        }
1785       
1786/*
1787        fprintf(stdout,"\nEnd of Multiple Alignment\n\n");
1788        fprintf(clusout,"\nEnd of Multiple Alignment\n\n");
1789        fprintf(clusout,"\nKey:\n     Identity:     #\n     Conservative: ^\n\n");
1790
1791        fprintf(clusout,"\nParameters:\n\n");
1792        fprintf(clusout,"Ktup:         %d\nWilbur Gap:   %d\nCutoff:       %d\nDiagonal:     %d\n",
1793        ktup,wind_gap,signif,window);
1794        fprintf(clusout,"Fixed Gap:    %d\nFloating Gap: %d\n\n",gap_open,gap_extend);
1795
1796        fprintf(clusout,"Sequences were: ");
1797        fprintf(clusout,dnaflag ? "Nucleic Acid\n" : "Proteins\n");
1798        fprintf(clusout,percent ? "Percentage" : "Absolute");
1799        fprintf(clusout," Identities Were Scored\n");
1800        if(dnaflag) {
1801                fprintf(clusout,"Nucleotide Transitions: ");
1802                fprintf(clusout,is_weight ? "WEIGHTED\n" : "UNWEIGHTED\n");
1803        }
1804        fprintf(clusout,"Sequence Input file: %s\n",seqname);
1805        fprintf(clusout,"Matrix used: ");
1806        if(matptr==idmat)
1807                fprintf(clusout,"No penalty\n");
1808        else
1809                if(matptr==pam250mt)
1810                        fprintf(clusout,"PAM 250\n");
1811        else
1812                if(matptr==pam100mt)
1813                        fprintf(clusout,"PAM 100\n");
1814        else
1815                fprintf(clusout,"%s\n",mtrxnam);
1816*/
1817}
1818
1819
1820
1821
1822
1823
1824static void gcg_out(FILE *gcgout)
1825{
1826/*        static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1827/*      static char *nbases = "XACGT";  */
1828        char seq[MAXLEN+1], residue;
1829        int all_checks[MAXN+1];
1830        int i,j,k,len,val,check,chunks,block,pos1,pos2; 
1831        long grand_checksum;
1832       
1833        len = seqlen_array[1];
1834
1835        for(i=1; i<=nseqs; i++) {
1836                for(j=1; j<=len; j++) {
1837                        val = seq_array[i][j];
1838                        if(val == 0 || val > 20) 
1839                                residue = 'X';
1840                        else if(val < 0) 
1841                                residue = '.';
1842                        else {
1843                                if(dnaflag)
1844                                        residue = nucleic_acid_order[val];
1845                                else
1846                                        residue = amino_acid_order[val];
1847                        }
1848                        seq[j] = residue;
1849                }
1850                all_checks[i] = SeqGCGCheckSum(seq+1, len);
1851        }       
1852
1853        grand_checksum = 0;
1854        for(i=1; i<=nseqs; i++) grand_checksum += all_checks[i];
1855        grand_checksum = grand_checksum % 10000;
1856        fprintf(gcgout,"\n\n   MSF:%5d  Type: ",len);
1857        if(dnaflag)
1858                fprintf(gcgout,"N");
1859        else
1860                fprintf(gcgout,"P");
1861        fprintf(gcgout,"    Check:%6ld   .. \n\n", grand_checksum);
1862        for(i=1; i<=nseqs; i++) 
1863/*              for(j=0; j<MAXNAMES; j++)
1864                        if(names[i][j] == ' ')  names[i][j] = '_';     */
1865                fprintf(gcgout,
1866                        " Name: %-15s oo  Len:%5d  Check:%6d  Weight:  1.00\n",
1867                        names[i],len,all_checks[i]);
1868        fprintf(gcgout,"\n//\n"); 
1869
1870        chunks = len/GCG_LINELENGTH;
1871        if(len % GCG_LINELENGTH != 0) ++chunks;
1872
1873        for(block=1; block<=chunks; block++) {
1874                fprintf(gcgout,"\n\n");
1875                pos1 = ((block-1) * GCG_LINELENGTH) + 1;
1876                pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
1877                for(i=1; i<=nseqs; i++) {
1878                        fprintf(gcgout,"\n%-15s ",names[i]);
1879                        for(j=pos1, k=1; j<=pos2; j++, k++) {
1880                                val = seq_array[i][j];
1881                                if(val == 0 || val > 20) 
1882                                        residue = 'X';
1883                                else if(val < 0) 
1884                                        residue = '.';
1885                                else {
1886                                        if(dnaflag)
1887                                                residue = nucleic_acid_order[val];
1888                                        else
1889                                                residue = amino_acid_order[val];
1890                                }
1891                                fprintf(gcgout,"%c",residue);
1892                                if(j % 10 == 0) fprintf(gcgout," ");
1893                        }
1894                }
1895        }
1896        fprintf(gcgout,"\n\n");
1897}
1898
1899
1900
1901static void phylip_out(FILE *phyout)
1902{
1903/*      static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1904/*              static char *nbases = "XACGT";  */
1905        char residue;
1906        int i,j,k,len,val,chunks,block,pos1,pos2;       
1907       
1908        len = seqlen_array[1];
1909
1910
1911        chunks = len/GCG_LINELENGTH;
1912        if(len % GCG_LINELENGTH != 0) ++chunks;
1913
1914        fprintf(phyout,"%6d %6d",nseqs,len);
1915
1916        for(block=1; block<=chunks; block++) {
1917                pos1 = ((block-1) * GCG_LINELENGTH) + 1;
1918                pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
1919                for(i=1; i<=nseqs; i++) {
1920                        if(block == 1) 
1921                                fprintf(phyout,"\n%-10s ",names[i]);
1922                        else
1923                                fprintf(phyout,"\n           ");
1924                        for(j=pos1, k=1; j<=pos2; j++, k++) {
1925                                val = seq_array[i][j];
1926                                if(val == 0 || val > 20) 
1927                                        residue = 'X';
1928                                else if(val < 0) 
1929                                        residue = '-';
1930                                else {
1931                                        if(dnaflag)
1932                                                residue = nucleic_acid_order[val];
1933                                        else
1934                                                residue = amino_acid_order[val];
1935                                }
1936                                fprintf(phyout,"%c",residue);
1937                                if(j % 10 == 0) fprintf(phyout," ");
1938                        }
1939                }
1940                fprintf(phyout,"\n");
1941        }
1942}
1943
1944
1945
1946
1947
1948static void nbrf_out(FILE *nbout)
1949{
1950/*      static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
1951/*              static char *nbases = "XACGT";  */
1952        char seq[MAXLEN+1], residue;
1953        int i,j,len,val;       
1954       
1955        len = seqlen_array[1];
1956
1957        for(i=1; i<=nseqs; i++) {
1958                fprintf(nbout, dnaflag ? ">DL;" : ">P1;");
1959                fprintf(nbout, "%s\n%s\n", names[i], titles[i]);
1960                for(j=1; j<=len; j++) {
1961                        val = seq_array[i][j];
1962                        if(val == 0 || val > 20) 
1963                                residue = 'X';
1964                        else if(val < 0) 
1965                                residue = '-';
1966                        else {
1967                                if(dnaflag)
1968                                        residue = nucleic_acid_order[val];
1969                                else
1970                                        residue = amino_acid_order[val];
1971                        }
1972                        seq[j] = residue;
1973                }
1974                for(j=1; j<=len; j++) {
1975                        fprintf(nbout,"%c",seq[j]);
1976                        if((j % LINELENGTH == 0) || (j == len)) 
1977                                fprintf(nbout,"\n");
1978                }
1979                fprintf(nbout,"*\n");
1980        }       
1981}
1982
1983
1984static Boolean open_alignment_output(char *path)
1985{
1986        if(output_clustal) 
1987                if((clustal_outfile = open_output_file(
1988                        "\nEnter a name for the CLUSTAL output file ",path,
1989                        clustal_outname,"aln"))==NULL) return FALSE;
1990        if(output_nbrf) 
1991                if((nbrf_outfile = open_output_file(
1992                        "\nEnter a name for the NBRF/PIR output file",path,
1993                        nbrf_outname,"pir"))==NULL) return FALSE;
1994        if(output_gcg) 
1995                if((gcg_outfile = open_output_file(
1996                        "\nEnter a name for the GCG output file     ",path,
1997                        gcg_outname,"msf"))==NULL) return FALSE;
1998        if(output_phylip) 
1999                if((phylip_outfile = open_output_file(
2000                        "\nEnter a name for the PHYLIP output file  ",path,
2001                        phylip_outname,"phy"))==NULL) return FALSE;
2002        return TRUE;
2003}
2004
2005
2006
2007
2008
2009static void create_alignment_output()
2010{
2011        if(output_clustal) {
2012                clustal_out(clustal_outfile);
2013                fclose(clustal_outfile);
2014                fprintf(stdout,"\nCLUSTAL-Alignment file created  [%s]\n",clustal_outname);
2015        }
2016        if(output_nbrf)  {
2017                nbrf_out(nbrf_outfile);
2018                fclose(nbrf_outfile);
2019                fprintf(stdout,"\nNBRF/PIR-Alignment file created [%s]\n",nbrf_outname);
2020        }
2021        if(output_gcg)  {
2022                gcg_out(gcg_outfile);
2023                fclose(gcg_outfile);
2024                fprintf(stdout,"\nGCG-Alignment file created      [%s]\n",gcg_outname);
2025        }
2026        if(output_phylip)  {
2027                phylip_out(phylip_outfile);
2028                fclose(phylip_outfile);
2029                fprintf(stdout,"\nPHYLIP-Alignment file created   [%s]\n",phylip_outname);
2030        }
2031}
2032
2033
2034
2035
2036
2037
2038static void reset()   /* remove gaps from older alignments (code = -1) */
2039{                     /* EXCEPT for gaps that were INPUT with the seqs.*/
2040        register int i,j,sl;                 /* which have  code = -2  */
2041       
2042        for(i=1;i<=nseqs;++i) {
2043                sl=0;
2044                for(j=1;j<=seqlen_array[i];++j) {
2045                        if(seq_array[i][j] == -1) continue;
2046                        ++sl;
2047                        seq_array[i][sl]=seq_array[i][j];
2048                }
2049                seqlen_array[i]=sl;
2050        }
2051}
2052
2053
Note: See TracBrowser for help on using the repository browser.