Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: tags/initial/GDE/CLUSTALW/interface.c

Visit:

Last change on this file was 2, checked in by oldcode, 24 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 97.0 KB

Line
1	/* command line interface for Clustal W */
2	/* DES was here MARCH. 1994 */
3	/* DES was here SEPT. 1994 */
4	#include <stdio.h>
5	#include <string.h>
6	#include <ctype.h>
7	#include <stdlib.h>
8	#include <signal.h>
9	#include <setjmp.h>
10	#include "clustalw.h"
11	#include "param.h"
12
13	/*
14	* Prototypes
15	*/
16
17	#ifdef UNIX
18	FILE open_path(char );
19	#endif
20
21
22	static sint check_param(char inline1,char params[], char *param_arg[]);
23	static void set_optional_param(void);
24	static sint find_match(char probe, char list[], sint n);
25	static void show_aln(void);
26	static void create_parameter_output(void);
27	static void reset_align(void);
28	static void reset_prf1(void);
29	static void reset_prf2(void);
30	static void calc_gap_penalty_mask(int prf_length,char struct_mask,char gap_mask);
31	void print_sec_struct_mask(int prf_length,char mask,char struct_mask);
32	/*
33	* Global variables
34	*/
35
36	extern sint max_names;
37
38	extern Boolean interactive;
39
40	extern double **tmat;
41	extern float gap_open, gap_extend;
42	extern float dna_gap_open, dna_gap_extend;
43	extern float prot_gap_open, prot_gap_extend;
44	extern float pw_go_penalty, pw_ge_penalty;
45	extern float dna_pw_go_penalty, dna_pw_ge_penalty;
46	extern float prot_pw_go_penalty, prot_pw_ge_penalty;
47	extern char revision_level[];
48	extern sint wind_gap,ktup,window,signif;
49	extern sint dna_wind_gap, dna_ktup, dna_window, dna_signif;
50	extern sint prot_wind_gap,prot_ktup,prot_window,prot_signif;
51	extern sint boot_ntrials; /* number of bootstrap trials */
52	extern sint nseqs;
53	extern sint new_seq;
54	extern sint *seqlen_array;
55	extern sint divergence_cutoff;
56	extern sint debug;
57	extern Boolean neg_matrix;
58	extern Boolean quick_pairalign;
59	extern Boolean reset_alignments; /* DES */
60	extern sint gap_dist;
61	extern Boolean no_hyd_penalties, no_pref_penalties;
62	extern sint max_aa;
63	extern sint gap_pos1, gap_pos2;
64	extern sint max_aln_length;
65	extern sint *output_index, output_order;
66	extern sint profile_no;
67	extern short usermat[], pw_usermat[];
68	extern short aa_xref[], pw_aa_xref[];
69	extern short userdnamat[], pw_userdnamat[];
70	extern short dna_xref[], pw_dna_xref[];
71
72	extern Boolean lowercase; /* Flag for GDE output - set on comm. line*/
73	extern Boolean cl_seq_numbers;
74	extern Boolean output_clustal, output_nbrf, output_phylip, output_gcg, output_gde;
75	extern Boolean output_tree_clustal, output_tree_phylip, output_tree_distances;
76	extern sint bootstrap_format;
77	extern Boolean tossgaps, kimura;
78	extern Boolean percent;
79	extern Boolean explicit_dnaflag; /* Explicit setting of sequence type on comm.line*/
80	extern Boolean usemenu;
81	extern Boolean showaln, save_parameters;
82	extern Boolean dnaflag;
83	extern float transition_weight;
84	extern unsigned sint boot_ran_seed;
85
86
87	extern FILE *tree;
88	extern FILE clustal_outfile, gcg_outfile, nbrf_outfile, phylip_outfile;
89	extern FILE *gde_outfile;
90
91	extern char hyd_residues[];
92	extern char *amino_acid_codes;
93	extern char *paramstr;
94	extern char seqname[];
95
96	extern char **seq_array;
97	extern char names, titles;
98
99	extern char gap_penalty_mask1,gap_penalty_mask2;
100	extern char sec_struct_mask1,sec_struct_mask2;
101	extern sint struct_penalties,struct_penalties1,struct_penalties2;
102	extern sint output_struct_penalties;
103	extern Boolean use_ss1, use_ss2;
104	extern char ss_name1,ss_name2;
105
106	char *ss_name = NULL;
107	char *sec_struct_mask = NULL;
108	char *gap_penalty_mask = NULL;
109
110	char profile1_name[FILENAMELEN+1];
111	char profile2_name[FILENAMELEN+1];
112
113	Boolean empty;
114	Boolean profile1_empty, profile2_empty; /* whether or not profiles */
115
116	char outfile_name[FILENAMELEN+1]="";
117
118	static char clustal_outname[FILENAMELEN+1], gcg_outname[FILENAMELEN+1];
119	static char phylip_outname[FILENAMELEN+1],nbrf_outname[FILENAMELEN+1];
120	static char gde_outname[FILENAMELEN+1];
121
122	char clustal_tree_name[FILENAMELEN+1]="";
123	char dist_tree_name[FILENAMELEN+1]="";
124	char phylip_tree_name[FILENAMELEN+1]="";
125	char p1_tree_name[FILENAMELEN+1]="";
126	char p2_tree_name[FILENAMELEN+1]="";
127
128	static char *params[MAXARGS];
129	static char *param_arg[MAXARGS];
130
131	static char *cmd_line_type[] = {
132	" ",
133	"=n ",
134	"=f ",
135	"=string ",
136	"=filename ",
137	""};
138
139	static sint numparams;
140	static Boolean check_tree = TRUE;
141
142	sint profile1_nseqs; /* have been filled; the no. of seqs in prof 1*/
143	Boolean use_tree_file = FALSE,new_tree_file = FALSE;
144	Boolean use_tree1_file = FALSE, use_tree2_file = FALSE;
145	Boolean new_tree1_file = FALSE, new_tree2_file = FALSE;
146
147	static char *lin2;
148
149	char *dnamatrix_txt[] = {
150	"IUB",
151	"CLUSTALW(1.6)",
152	"user defined" };
153
154	char *matrix_txt[] = {
155	"BLOSUM series",
156	"PAM series",
157	"Gonnet series",
158	"Identity matrix",
159	"user defined" };
160
161	char *pw_matrix_txt[] = {
162	"BLOSUM30",
163	"PAM350",
164	"Gonnet250",
165	"Identity matrix",
166	"user defined" };
167
168	void init_interface(void)
169	{
170	empty=TRUE;
171
172	profile1_empty = TRUE; /* */
173	profile2_empty = TRUE; /* */
174
175	lin2 = (char )ckalloc( (MAXLINE+1) sizeof (char) );
176
177	}
178
179
180
181
182	static sint check_param(char inline1,char params[], char *param_arg[])
183	{
184
185	char *inptr;
186	char command_sep[5];
187	/*
188	#ifndef MAC
189	char strtok(char s1, const char *s2);
190	#endif
191	*/
192	sint len,i,j,k,s,n,match[MAXARGS];
193	Boolean name1 = FALSE;
194
195	command_sep[0] = COMMANDSEP;
196	command_sep[1] = '\0';
197
198	if (inline1[0] != COMMANDSEP) name1 = TRUE;
199
200	inptr=inline1;
201	for (i=0;i<MAXARGS;i++) {
202	if ((params[i]=strtok(inptr,command_sep))==NULL)
203	break;
204	inptr=NULL;
205	}
206	if (i==MAXARGS) {
207	fprintf(stdout,"Error: too many command line arguments\n");
208	return(-1);
209	}
210	/*
211	special case - first parameter is input filename
212	*/
213	s = 0;
214	if(paramstr[0] != COMMANDSEP) {
215	strcpy(seqname, params[0]);
216	/* JULIE
217	convert to lower case now
218	*/
219	#ifndef UNIX
220	for(k=0;k<(sint)strlen(params[0]);++k) seqname[k]=tolower(params[0][k]);
221	#else
222	for(k=0;k<(sint)strlen(params[0]);++k) seqname[k]=params[0][k];
223	#endif
224	s++;
225	}
226
227	n = i;
228	for (i=s;i<n;i++) {
229	param_arg[i] = NULL;
230	len = (sint)strlen(params[i]);
231	for(j=0; j<len; j++)
232	if(params[i][j] == '=') {
233	param_arg[i] = (char )ckalloc((len-j) sizeof(char));
234	strncpy(param_arg[i],&params[i][j+1],len-j-1);
235	params[i][j] = EOS;
236	/* JULIE
237	convert keywords to lower case now
238	*/
239	for(k=0;k<j;++k) params[i][k]=tolower(params[i][k]);
240	param_arg[i][len-j-1] = EOS;
241	break;
242	}
243	}
244
245	/*
246	for each parameter given on the command line, first search the list of recognised optional
247	parameters....
248	*/
249
250	for (i=0;i<n;i++) {
251	if ((i==0) && (name1 == TRUE)) continue;
252	j = 0;
253	match[i] = -1;
254	for(;;) {
255	if (cmd_line_para[j].str[0] == '\0') break;
256	if (!strcmp(params[i],cmd_line_para[j].str)) {
257	match[i] = j;
258	*cmd_line_para[match[i]].flag = i;
259	if ((cmd_line_para[match[i]].type != NOARG) &&
260	(param_arg[i] == NULL)) {
261	fprintf(stdout,
262	"Error: parameter required for /%s\n",params[i]);
263	exit(1);
264	}
265	/* JULIE
266	convert parameters to lower case now, unless the parameter is a filename
267	*/
268	#ifdef UNIX
269	else if (cmd_line_para[match[i]].type != FILARG
270	&& param_arg[i] != NULL)
271	#endif
272	if (param_arg[i]!=0)
273	{
274	for(k=0;k<strlen(param_arg[i]);++k)
275	param_arg[i][k]=tolower(param_arg[i][k]);
276	}
277	break;
278	}
279	j++;
280	}
281	}
282	/*
283	....then the list of recognised input files,....
284	*/
285	for (i=0;i<n;i++) {
286	if ((i==0) && (name1 == TRUE)) continue;
287	if (match[i] != -1) continue;
288	j = 0;
289	for(;;) {
290	if (cmd_line_file[j].str[0] == '\0') break;
291	if (!strcmp(params[i],cmd_line_file[j].str)) {
292	match[i] = j;
293	*cmd_line_file[match[i]].flag = i;
294	if ((cmd_line_file[match[i]].type != NOARG) &&
295	(param_arg[i] == NULL)) {
296	fprintf(stdout,
297	"Error: parameter required for /%s\n",params[i]);
298	exit(1);
299	}
300	break;
301	}
302	j++;
303	}
304	}
305	/*
306	....and finally the recognised verbs.
307	*/
308	for (i=0;i<n;i++) {
309	if ((i==0) && (name1 == TRUE)) continue;
310	if (match[i] != -1) continue;
311	j = 0;
312	for(;;) {
313	if (cmd_line_verb[j].str[0] == '\0') break;
314	if (!strcmp(params[i],cmd_line_verb[j].str)) {
315	match[i] = j;
316	*cmd_line_verb[match[i]].flag = i;
317	if ((cmd_line_verb[match[i]].type != NOARG) &&
318	(param_arg[i] == NULL)) {
319	fprintf(stdout,
320	"Error: parameter required for /%s\n",params[i]);
321	exit(1);
322	}
323	break;
324	}
325	j++;
326	}
327	}
328
329	/*
330	check for any unrecognised parameters.
331	*/
332	for (i=0;i<n;i++) {
333	if (match[i] == -1) {
334	fprintf(stdout,
335	"Error: unknown option %c%s\n",COMMANDSEP,params[i]);
336	exit(1);
337	}
338	}
339	return(n);
340	}
341
342	static void set_optional_param(void)
343	{
344	int i,temp;
345	int c;
346	float ftemp;
347
348	/****************************************************************************/
349	/* look for parameters on command line e.g. gap penalties, k-tuple etc. */
350	/****************************************************************************/
351
352	/*** ? /score=percent or /score=absolute */
353	if(setscore != -1)
354	if(strlen(param_arg[setscore]) > 0) {
355	temp = find_match(param_arg[setscore],score_arg,2);
356	if(temp == 0)
357	percent = TRUE;
358	else if(temp == 1)
359	percent = FALSE;
360	else
361	fprintf(stdout,"\nUnknown SCORE type: %s\n",
362	param_arg[setscore]);
363	}
364
365	/*** ? /seed=n */
366	if(setseed != -1) {
367	temp = 0;
368	if(strlen(param_arg[setseed]) > 0)
369	if (sscanf(param_arg[setseed],"%d",&temp)!=1) {
370	fprintf(stdout,"Bad option for /seed (must be integer)\n");
371	temp = 0;
372	}
373	if(temp > 0) boot_ran_seed = temp;
374	fprintf(stdout,"\ntemp = %d; seed = %u;\n",(pint)temp,boot_ran_seed);
375	}
376
377
378
379
380	/*** ? /output=PIR, GCG, GDE or PHYLIP */
381	if(setoutput != -1)
382	if(strlen(param_arg[setoutput]) > 0) {
383	temp = find_match(param_arg[setoutput],output_arg,4);
384	if (temp >= 0 && temp <= 3) {
385	output_clustal = FALSE;
386	output_gcg = FALSE;
387	output_phylip = FALSE;
388	output_nbrf = FALSE;
389	output_gde = FALSE;
390	}
391	switch (temp) {
392	case 0: /* GCG */
393	output_gcg = TRUE;
394	break;
395	case 1: /* GDE */
396	output_gde = TRUE;
397	break;
398	case 2: /* PIR */
399	output_nbrf = TRUE;
400	break;
401	case 3: /* PHYLIP */
402	output_phylip = TRUE;
403	break;
404	default:
405	fprintf(stdout,"\nUnknown OUTPUT type: %s\n",
406	param_arg[setoutput]);
407	}
408	}
409
410	/*** ? /outputtree=NJ or PHYLIP or DIST */
411	if(setoutputtree != -1)
412	if(strlen(param_arg[setoutputtree]) > 0) {
413	temp = find_match(param_arg[setoutputtree],outputtree_arg,3);
414	switch (temp) {
415	case 0: /* NJ */
416	output_tree_clustal = TRUE;
417	break;
418	case 1: /* PHYLIP */
419	output_tree_phylip = TRUE;
420	break;
421	case 2: /* DIST */
422	output_tree_distances = TRUE;
423	break;
424	default:
425	fprintf(stdout,"\nUnknown OUTPUT TREE type: %s\n",
426	param_arg[setoutputtree]);
427	}
428	}
429
430	/*** ? /profile (sets type of second input file to profile) */
431	if(setprofile != -1)
432	profile_type = PROFILE;
433
434	/*** ? /sequences (sets type of second input file to list of sequences) */
435	if(setsequences != -1)
436	profile_type = SEQUENCE;
437
438
439
440	/*** ? /ktuple=n */
441	if(setktuple != -1) {
442	temp = 0;
443	if(strlen(param_arg[setktuple]) > 0)
444	if (sscanf(param_arg[setktuple],"%d",&temp)!=1) {
445	fprintf(stdout,"Bad option for /ktuple (must be integer)\n");
446	temp = 0;
447	}
448	if(temp > 0) {
449	if(dnaflag) {
450	if(temp <= 4) {
451	ktup = temp;
452	dna_ktup = ktup;
453	wind_gap = ktup + 4;
454	dna_wind_gap = wind_gap;
455	}
456	}
457	else {
458	if(temp <= 2) {
459	ktup = temp;
460	prot_ktup = ktup;
461	wind_gap = ktup + 3;
462	prot_wind_gap = wind_gap;
463	}
464	}
465	}
466	}
467
468	/*** ? /pairgap=n */
469	if(setpairgap != -1) {
470	temp = 0;
471	if(strlen(param_arg[setpairgap]) > 0)
472	if (sscanf(param_arg[setpairgap],"%d",&temp)!=1) {
473	fprintf(stdout,"Bad option for /pairgap (must be integer)\n");
474	temp = 0;
475	}
476	if(temp > 0)
477	if(dnaflag) {
478	if(temp > ktup) {
479	wind_gap = temp;
480	dna_wind_gap = wind_gap;
481	}
482	}
483	else {
484	if(temp > ktup) {
485	wind_gap = temp;
486	prot_wind_gap = wind_gap;
487	}
488	}
489	}
490
491
492	/*** ? /topdiags=n */
493	if(settopdiags != -1) {
494	temp = 0;
495	if(strlen(param_arg[settopdiags]) > 0)
496	if (sscanf(param_arg[settopdiags],"%d",&temp)!=1) {
497	fprintf(stdout,"Bad option for /topdiags (must be integer)\n");
498	temp = 0;
499	}
500	if(temp > 0)
501	if(dnaflag) {
502	if(temp > ktup) {
503	signif = temp;
504	dna_signif = signif;
505	}
506	}
507	else {
508	if(temp > ktup) {
509	signif = temp;
510	prot_signif = signif;
511	}
512	}
513	}
514
515
516	/*** ? /window=n */
517	if(setwindow != -1) {
518	temp = 0;
519	if(strlen(param_arg[setwindow]) > 0)
520	if (sscanf(param_arg[setwindow],"%d",&temp)!=1) {
521	fprintf(stdout,"Bad option for /window (must be integer)\n");
522	temp = 0;
523	}
524	if(temp > 0)
525	if(dnaflag) {
526	if(temp > ktup) {
527	window = temp;
528	dna_window = window;
529	}
530	}
531	else {
532	if(temp > ktup) {
533	window = temp;
534	prot_window = window;
535	}
536	}
537	}
538
539	/*** ? /kimura */
540	if(setkimura != -1)
541	kimura = TRUE;
542
543	/*** ? /tossgaps */
544	if(settossgaps != -1)
545	tossgaps = TRUE;
546
547
548	/*** ? /negative */
549	if(setnegative != -1)
550	neg_matrix = TRUE;
551
552
553	/*** ? /pwmatrix=ID (user's file) */
554	if(setpwmatrix != -1)
555	if(strlen(param_arg[setpwmatrix]) > 0) {
556	if (strcmp(param_arg[setpwmatrix],"blosum")==0) {
557	strcpy(pw_mtrxname, param_arg[setpwmatrix]);
558	pw_matnum = 1;
559	}
560	else if (strcmp(param_arg[setpwmatrix],"pam")==0) {
561	strcpy(pw_mtrxname, param_arg[setpwmatrix]);
562	pw_matnum = 2;
563	}
564	else if (strcmp(param_arg[setpwmatrix],"gonnet")==0) {
565	strcpy(pw_mtrxname, param_arg[setpwmatrix]);
566	pw_matnum = 3;
567	}
568	else if (strcmp(param_arg[setpwmatrix],"id")==0) {
569	strcpy(pw_mtrxname, param_arg[setpwmatrix]);
570	pw_matnum = 4;
571	}
572	else {
573	if(user_mat(param_arg[setpwmatrix], pw_usermat, pw_aa_xref))
574	{
575	strcpy(pw_mtrxname,param_arg[setpwmatrix]);
576	pw_matnum=5;
577	}
578	else exit(1);
579	}
580
581	}
582
583
584	/*** ? /matrix=ID (user's file) */
585	if(setmatrix != -1)
586	if(strlen(param_arg[setmatrix]) > 0) {
587	if (strcmp(param_arg[setmatrix],"blosum")==0) {
588	strcpy(mtrxname, param_arg[setmatrix]);
589	matnum = 1;
590	}
591	else if (strcmp(param_arg[setmatrix],"pam")==0) {
592	strcpy(mtrxname, param_arg[setmatrix]);
593	matnum = 2;
594	}
595	else if (strcmp(param_arg[setmatrix],"gonnet")==0) {
596	strcpy(mtrxname, param_arg[setmatrix]);
597	matnum = 3;
598	}
599	else if (strcmp(param_arg[setmatrix],"id")==0) {
600	strcpy(mtrxname, param_arg[setmatrix]);
601	matnum = 4;
602	}
603	else {
604	if(user_mat(param_arg[setmatrix], usermat, aa_xref))
605	{
606	strcpy(mtrxname,param_arg[setmatrix]);
607	matnum=5;
608	}
609	else exit(1);
610	}
611
612	}
613
614
615	/*** ? /pwdnamatrix=ID (user's file) */
616	if(setpwdnamatrix != -1)
617	if(strlen(param_arg[setpwdnamatrix]) > 0) {
618	if (strcmp(param_arg[setpwdnamatrix],"iub")==0) {
619	strcpy(pw_dnamtrxname, param_arg[setpwdnamatrix]);
620	pw_dnamatnum = 1;
621	}
622	else if (strcmp(param_arg[setpwdnamatrix],"clustalw")==0) {
623	strcpy(pw_dnamtrxname, param_arg[setpwdnamatrix]);
624	pw_dnamatnum = 2;
625	}
626	else {
627	if(user_mat(param_arg[setpwdnamatrix], pw_userdnamat, pw_dna_xref))
628	{
629	strcpy(pw_dnamtrxname,param_arg[setpwdnamatrix]);
630	pw_dnamatnum=3;
631	}
632	else exit(1);
633	}
634
635	}
636
637
638	/*** ? /matrix=ID (user's file) */
639	if(setdnamatrix != -1)
640	if(strlen(param_arg[setdnamatrix]) > 0) {
641	if (strcmp(param_arg[setdnamatrix],"iub")==0) {
642	strcpy(dnamtrxname, param_arg[setdnamatrix]);
643	dnamatnum = 1;
644	}
645	else if (strcmp(param_arg[setdnamatrix],"clustalw")==0) {
646	strcpy(dnamtrxname, param_arg[setdnamatrix]);
647	dnamatnum = 2;
648	}
649	else {
650	if(user_mat(param_arg[setdnamatrix], userdnamat, dna_xref))
651	{
652	strcpy(dnamtrxname,param_arg[setdnamatrix]);
653	dnamatnum=3;
654	}
655	else exit(1);
656	}
657
658	}
659	/*** ? /maxdiv= n */
660	if(setmaxdiv != -1) {
661	temp = 0;
662	if(strlen(param_arg[setmaxdiv]) > 0)
663	if (sscanf(param_arg[setmaxdiv],"%d",&temp)!=1) {
664	fprintf(stdout,"Bad option for /maxdiv (must be integer)\n");
665	temp = 0;
666	}
667	if (temp >= 0)
668	divergence_cutoff = temp;
669	}
670
671	/*** ? /gapdist= n */
672	if(setgapdist != -1) {
673	temp = 0;
674	if(strlen(param_arg[setgapdist]) > 0)
675	if (sscanf(param_arg[setgapdist],"%d",&temp)!=1) {
676	fprintf(stdout,"Bad option for /gapdist (must be integer)\n");
677	temp = 0;
678	}
679	if (temp >= 0)
680	gap_dist = temp;
681	}
682
683	/*** ? /debug= n */
684	if(setdebug != -1) {
685	temp = 0;
686	if(strlen(param_arg[setdebug]) > 0)
687	if (sscanf(param_arg[setdebug],"%d",&temp)!=1) {
688	fprintf(stdout,"Bad option for /debug (must be integer)\n");
689	temp = 0;
690	}
691	if (temp >= 0)
692	debug = temp;
693	}
694
695	/*** ? /outfile= (user's file) */
696	if(setoutfile != -1)
697	if(strlen(param_arg[setoutfile]) > 0) {
698	strcpy(outfile_name, param_arg[setoutfile]);
699	}
700
701	/*** ? /case= lower/upper */
702	if(setcase != -1)
703	if(strlen(param_arg[setcase]) > 0) {
704	temp = find_match(param_arg[setcase],case_arg,2);
705	if(temp == 0) {
706	lowercase = TRUE;
707	}
708	else if(temp == 1) {
709	lowercase = FALSE;
710	}
711	else
712	fprintf(stdout,"\nUnknown case %s\n",
713	param_arg[setcase]);
714	}
715
716	/*** ? /seqnos=off/on */
717	if(setseqno != -1)
718	if(strlen(param_arg[setseqno]) > 0) {
719	temp = find_match(param_arg[setseqno],seqno_arg,2);
720	if(temp == 0) {
721	cl_seq_numbers = FALSE;
722	}
723	else if(temp == 1) {
724	cl_seq_numbers = TRUE;
725	}
726	else
727	fprintf(stdout,"\nUnknown SEQNO option %s\n",
728	param_arg[setseqno]);
729	}
730
731
732	/*** ? /gapopen=n */
733	if(setgapopen != -1) {
734	ftemp = 0.0;
735	if(strlen(param_arg[setgapopen]) > 0)
736	if (sscanf(param_arg[setgapopen],"%f",&ftemp)!=1) {
737	fprintf(stdout,"Bad option for /gapopen (must be real number)\n");
738	ftemp = 0.0;
739	}
740	if(ftemp >= 0.0)
741	if(dnaflag) {
742	gap_open = ftemp;
743	dna_gap_open = gap_open;
744	}
745	else {
746	gap_open = ftemp;
747	prot_gap_open = gap_open;
748	}
749	}
750
751
752	/*** ? /gapext=n */
753	if(setgapext != -1) {
754	ftemp = 0.0;
755	if(strlen(param_arg[setgapext]) > 0)
756	if (sscanf(param_arg[setgapext],"%f",&ftemp)!=1) {
757	fprintf(stdout,"Bad option for /gapext (must be real number)\n");
758	ftemp = 0.0;
759	}
760	if(ftemp >= 0)
761	if(dnaflag) {
762	gap_extend = ftemp;
763	dna_gap_extend = gap_extend;
764	}
765	else {
766	gap_extend = ftemp;
767	prot_gap_extend = gap_extend;
768	}
769	}
770
771	/*** ? /transweight=n*/
772	if(settransweight != -1) {
773	ftemp = 0.0;
774	if(strlen(param_arg[settransweight]) > 0)
775	if (sscanf(param_arg[settransweight],"%f",&ftemp)!=1) {
776	fprintf(stdout,"Bad option for /transweight (must be real number)\n");
777	ftemp = 0.0;
778	}
779	transition_weight=temp;
780	}
781
782	/*** ? /pwgapopen=n */
783	if(setpwgapopen != -1) {
784	ftemp = 0.0;
785	if(strlen(param_arg[setpwgapopen]) > 0)
786	if (sscanf(param_arg[setpwgapopen],"%f",&ftemp)!=1) {
787	fprintf(stdout,"Bad option for /pwgapopen (must be real number)\n");
788	ftemp = 0.0;
789	}
790	if(ftemp >= 0.0)
791	if(dnaflag) {
792	pw_go_penalty = ftemp;
793	dna_pw_go_penalty = pw_go_penalty;
794	}
795	else {
796	pw_go_penalty = ftemp;
797	prot_pw_go_penalty = pw_go_penalty;
798	}
799	}
800
801
802	/*** ? /gapext=n */
803	if(setpwgapext != -1) {
804	ftemp = 0.0;
805	if(strlen(param_arg[setpwgapext]) > 0)
806	if (sscanf(param_arg[setpwgapext],"%f",&ftemp)!=1) {
807	fprintf(stdout,"Bad option for /pwgapext (must be real number)\n");
808	ftemp = 0.0;
809	}
810	if(ftemp >= 0)
811	if(dnaflag) {
812	pw_ge_penalty = ftemp;
813	dna_pw_ge_penalty = pw_ge_penalty;
814	}
815	else {
816	pw_ge_penalty = ftemp;
817	prot_pw_ge_penalty = pw_ge_penalty;
818	}
819	}
820
821
822
823	/*** ? /outorder=n */
824	if(setoutorder != -1) {
825	if(strlen(param_arg[setoutorder]) > 0)
826	temp = find_match(param_arg[setoutorder],outorder_arg,2);
827	if(temp == 0) {
828	output_order = INPUT;
829	}
830	else if(temp == 1) {
831	output_order = ALIGNED;
832	}
833	else
834	fprintf(stdout,"\nUnknown OUTPUT ORDER type %s\n",
835	param_arg[setoutorder]);
836	}
837
838	/*** ? /bootlabels=n */
839	if(setbootlabels != -1) {
840	if(strlen(param_arg[setbootlabels]) > 0)
841	temp = find_match(param_arg[setbootlabels],bootlabels_arg,2);
842	if(temp == 0) {
843	bootstrap_format = BS_NODE_LABELS;
844	}
845	else if(temp == 1) {
846	bootstrap_format = BS_BRANCH_LABELS;
847	}
848	else
849	fprintf(stdout,"\nUnknown bootlabels type %s\n",
850	param_arg[setoutorder]);
851	}
852
853	/*** ? /endgaps */
854	if(setuseendgaps != -1)
855	use_endgaps = FALSE;
856
857	/*** ? /nopgap */
858	if(setnopgap != -1)
859	no_pref_penalties = TRUE;
860
861	/*** ? /nohgap */
862	if(setnohgap != -1)
863	no_hyd_penalties = TRUE;
864
865	/*** ? /hgapresidues="string" */
866	if(sethgapres != -1)
867	if(strlen(param_arg[sethgapres]) > 0) {
868	for (i=0;i<strlen(hyd_residues) && i<26;i++) {
869	c = param_arg[sethgapres][i];
870	if (isalpha(c))
871	hyd_residues[i] = (char)toupper(c);
872	else
873	break;
874	}
875	}
876
877
878	/*** ? /nosecstr1 */
879	if(setsecstr1 != -1)
880	use_ss1 = FALSE;
881
882	/*** ? /nosecstr2 */
883	if(setsecstr2 != -1)
884	use_ss2 = FALSE;
885
886	/*** ? /secstroutput */
887	if(setsecstroutput != -1)
888	if(strlen(param_arg[setsecstroutput]) > 0) {
889	temp = find_match(param_arg[setsecstroutput],outputsecstr_arg,4);
890	if(temp >= 0 && temp <= 3)
891	output_struct_penalties = temp;
892	else
893	fprintf(stdout,"\nUnknown case %s\n",
894	param_arg[setsecstroutput]);
895	}
896
897
898	/*** ? /helixgap= n */
899	if(sethelixgap != -1) {
900	temp = 0;
901	if(strlen(param_arg[sethelixgap]) > 0)
902	if (sscanf(param_arg[sethelixgap],"%d",&temp)!=1) {
903	fprintf(stdout,"Bad option for /helixgap (must be integer)\n");
904	temp = 0;
905	}
906	if (temp >= 1 && temp <= 9)
907	helix_penalty = temp;
908	}
909
910	/*** ? /strandgap= n */
911	if(setstrandgap != -1) {
912	temp = 0;
913	if(strlen(param_arg[setstrandgap]) > 0)
914	if (sscanf(param_arg[setstrandgap],"%d",&temp)!=1) {
915	fprintf(stdout,"Bad option for /strandgap (must be integer)\n");
916	temp = 0;
917	}
918	if (temp >= 1 && temp <= 9)
919	strand_penalty = temp;
920	}
921
922	/*** ? /loopgap= n */
923	if(setloopgap != -1) {
924	temp = 0;
925	if(strlen(param_arg[setloopgap]) > 0)
926	if (sscanf(param_arg[setloopgap],"%d",&temp)!=1) {
927	fprintf(stdout,"Bad option for /loopgap (must be integer)\n");
928	temp = 0;
929	}
930	if (temp >= 1 && temp <= 9)
931	loop_penalty = temp;
932	}
933
934	/*** ? /terminalgap= n */
935	if(setterminalgap != -1) {
936	temp = 0;
937	if(strlen(param_arg[setterminalgap]) > 0)
938	if (sscanf(param_arg[setterminalgap],"%d",&temp)!=1) {
939	fprintf(stdout,"Bad option for /terminalgap (must be integer)\n");
940	temp = 0;
941	}
942	if (temp >= 1 && temp <= 9) {
943	helix_end_penalty = temp;
944	strand_end_penalty = temp;
945	}
946	}
947
948	/*** ? /helixendin= n */
949	if(sethelixendin != -1) {
950	temp = 0;
951	if(strlen(param_arg[sethelixendin]) > 0)
952	if (sscanf(param_arg[sethelixendin],"%d",&temp)!=1) {
953	fprintf(stdout,"Bad option for /helixendin (must be integer)\n");
954	temp = 0;
955	}
956	if (temp >= 0 && temp <= 3)
957	helix_end_minus = temp;
958	}
959
960	/*** ? /helixendout= n */
961	if(sethelixendout != -1) {
962	temp = 0;
963	if(strlen(param_arg[sethelixendout]) > 0)
964	if (sscanf(param_arg[sethelixendout],"%d",&temp)!=1) {
965	fprintf(stdout,"Bad option for /helixendout (must be integer)\n");
966	temp = 0;
967	}
968	if (temp >= 0 && temp <= 3)
969	helix_end_plus = temp;
970	}
971
972	/*** ? /strandendin= n */
973	if(setstrandendin != -1) {
974	temp = 0;
975	if(strlen(param_arg[setstrandendin]) > 0)
976	if (sscanf(param_arg[setstrandendin],"%d",&temp)!=1) {
977	fprintf(stdout,"Bad option for /strandendin (must be integer)\n");
978	temp = 0;
979	}
980	if (temp >= 0 && temp <= 3)
981	strand_end_minus = temp;
982	}
983
984	/*** ? /strandendout= n */
985	if(setstrandendout != -1) {
986	temp = 0;
987	if(strlen(param_arg[setstrandendout]) > 0)
988	if (sscanf(param_arg[setstrandendout],"%d",&temp)!=1) {
989	fprintf(stdout,"Bad option for /strandendout (must be integer)\n");
990	temp = 0;
991	}
992	if (temp >= 0 && temp <= 3)
993	strand_end_plus = temp;
994	}
995
996	}
997
998	#ifdef UNIX
999	FILE open_path(char fname) /* to open in read-only file fname searching for
1000	it through all path directories */
1001	{
1002	#define Mxdir 70
1003	char dir[Mxdir+1], path, deb, *fin;
1004	FILE *fich;
1005	sint lf, ltot;
1006
1007	path=getenv("PATH"); /* get the list of path directories,
1008	separated by :
1009	*/
1010	if (path == NULL ) return fopen(fname,"r");
1011	lf=(sint)strlen(fname);
1012	deb=path;
1013	do
1014	{
1015	fin=strchr(deb,':');
1016	if(fin!=NULL)
1017	{ strncpy(dir,deb,fin-deb); ltot=fin-deb; }
1018	else
1019	{ strcpy(dir,deb); ltot=(sint)strlen(dir); }
1020	/* now one directory is in string dir */
1021	if( ltot + lf + 1 <= Mxdir)
1022	{
1023	dir[ltot]='/';
1024	strcpy(dir+ltot+1,fname); /* now dir is appended with fi
1025	lename */
1026	if( (fich = fopen(dir,"r") ) != NULL) break;
1027	}
1028	else fich = NULL;
1029	deb=fin+1;
1030	}
1031	while (fin != NULL);
1032	return fich;
1033	}
1034	#endif
1035
1036
1037	void get_help(char help_pointer) /* Help procedure */
1038	{
1039	FILE *help_file;
1040	sint i, number, nlines;
1041	Boolean found_help;
1042	char temp[MAXLINE+1];
1043	char token = '\0';
1044	char *digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1045	char *help_marker = ">>HELP";
1046
1047	extern char *help_file_name;
1048
1049	#ifdef VMS
1050	if((help_file=fopen(help_file_name,"r","rat=cr","rfm=var"))==NULL) {
1051	error("Cannot open help file [%s]",help_file_name);
1052	return;
1053	}
1054	#else
1055
1056	#ifdef UNIX
1057	if((help_file=open_path(help_file_name))==NULL) {
1058	if((help_file=fopen(help_file_name,"r"))==NULL) {
1059	error("Cannot open help file [%s]",help_file_name);
1060	return;
1061	}
1062	}
1063
1064	#else
1065	if((help_file=fopen(help_file_name,"r"))==NULL) {
1066	error("Cannot open help file [%s]",help_file_name);
1067	return;
1068	}
1069	#endif
1070
1071	#endif
1072	/* error("Cannot open help file [%s]",help_file_name);
1073	return;
1074	}
1075	*/
1076	nlines = 0;
1077	number = -1;
1078	found_help = FALSE;
1079
1080	while(TRUE) {
1081	if(fgets(temp,MAXLINE+1,help_file) == NULL) {
1082	if(!found_help)
1083	error("No help found in help file");
1084	fclose(help_file);
1085	return;
1086	}
1087	if(strstr(temp,help_marker)) {
1088	token = ' ';
1089	for(i=strlen(help_marker); i<8; i++)
1090	if(strchr(digits, temp[i])) {
1091	token = temp[i];
1092	break;
1093	}
1094	}
1095	if(token == help_pointer) {
1096	found_help = TRUE;
1097	while(fgets(temp,MAXLINE+1,help_file)) {
1098	if(strstr(temp, help_marker)){
1099	if(usemenu) {
1100	fprintf(stdout,"\n");
1101	getstr("Press [RETURN] to continue",lin2);
1102	}
1103	fclose(help_file);
1104	return;
1105	}
1106	fputs(temp,stdout);
1107	++nlines;
1108	if(usemenu) {
1109	if(nlines >= PAGE_LEN) {
1110	fprintf(stdout,"\n");
1111	getstr("Press [RETURN] to continue or X to stop",lin2);
1112	if(toupper(*lin2) == 'X') {
1113	fclose(help_file);
1114	return;
1115	}
1116	else
1117	nlines = 0;
1118	}
1119	}
1120	}
1121	if(usemenu) {
1122	fprintf(stdout,"\n");
1123	getstr("Press [RETURN] to continue",lin2);
1124	}
1125	fclose(help_file);
1126	}
1127	}
1128	}
1129
1130	static void show_aln(void) /* Alignment screen display procedure */
1131	{
1132	FILE *file;
1133	sint nlines;
1134	char temp[MAXLINE+1];
1135	char file_name[FILENAMELEN+1];
1136
1137	if(output_clustal) strcpy(file_name,clustal_outname);
1138	else if(output_nbrf) strcpy(file_name,nbrf_outname);
1139	else if(output_gcg) strcpy(file_name,gcg_outname);
1140	else if(output_phylip) strcpy(file_name,phylip_outname);
1141	else if(output_gde) strcpy(file_name,gde_outname);
1142
1143	#ifdef VMS
1144	if((file=fopen(file_name,"r","rat=cr","rfm=var"))==NULL) {
1145	#else
1146	if((file=fopen(file_name,"r"))==NULL) {
1147	#endif
1148	error("Cannot open file [%s]",file_name);
1149	return;
1150	}
1151
1152	fprintf(stdout,"\n\n");
1153	nlines = 0;
1154
1155	while(fgets(temp,MAXLINE+1,file)) {
1156	fputs(temp,stdout);
1157	++nlines;
1158	if(nlines >= PAGE_LEN) {
1159	fprintf(stdout,"\n");
1160	getstr("Press [RETURN] to continue or X to stop",lin2);
1161	if(toupper(*lin2) == 'X') {
1162	fclose(file);
1163	return;
1164	}
1165	else
1166	nlines = 0;
1167	}
1168	}
1169	fclose(file);
1170	fprintf(stdout,"\n");
1171	getstr("Press [RETURN] to continue",lin2);
1172	}
1173
1174
1175	void parse_params(Boolean xmenus)
1176	{
1177	sint i,j,len,temp;
1178	char path[FILENAMELEN];
1179
1180
1181	Boolean do_align, do_convert, do_align_only, do_tree_only, do_tree, do_boot, do_profile, do_something;
1182
1183	if (!xmenus)
1184	{
1185	fprintf(stdout,"\n\n\n");
1186	fprintf(stdout," CLUSTAL %s Multiple Sequence Alignments\n\n\n",revision_level);
1187	}
1188
1189	do_align = do_convert = do_align_only = do_tree_only = do_tree = do_boot = do_profile = do_something = FALSE;
1190
1191	*seqname=EOS;
1192
1193	len=(sint)strlen(paramstr);
1194	/* JULIE
1195	Stop converting command line to lower case - unix, mac, pc are case sensitive
1196	for(i=0;i<len;++i) paramstr[i]=tolower(paramstr[i]);
1197	*/
1198
1199	numparams = check_param(paramstr, params, param_arg);
1200	if (numparams <0) exit(1);
1201
1202	if(sethelp != -1) {
1203	get_help('9');
1204	exit(1);
1205	}
1206
1207	if(setoptions != -1) {
1208	fprintf(stdout,"clustalw option list:-\n");
1209	for (i=0;cmd_line_verb[i].str[0] != '\0';i++) {
1210	fprintf(stdout,"\t\t%c%s%s",COMMANDSEP,cmd_line_verb[i].str,cmd_line_type[cmd_line_verb[i].type]);
1211	if (cmd_line_verb[i].type == OPTARG) {
1212	if (cmd_line_verb[i].arg[0][0] != '\0')
1213	fprintf(stdout,"=%s",cmd_line_verb[i].arg[0]);
1214	for (j=1;cmd_line_verb[i].arg[j][0] != '\0';j++)
1215	fprintf(stdout," OR %s",cmd_line_verb[i].arg[j]);
1216	}
1217	fprintf(stdout,"\n");
1218	}
1219	for (i=0;cmd_line_file[i].str[0] != '\0';i++) {
1220	fprintf(stdout,"\t\t%c%s%s",COMMANDSEP,cmd_line_file[i].str,cmd_line_type[cmd_line_file[i].type]);
1221	if (cmd_line_file[i].type == OPTARG) {
1222	if (cmd_line_file[i].arg[0][0] != '\0')
1223	fprintf(stdout,"=%s",cmd_line_file[i].arg[0]);
1224	for (j=1;cmd_line_file[i].arg[j][0] != '\0';j++)
1225	fprintf(stdout," OR %s",cmd_line_file[i].arg[j]);
1226	}
1227	fprintf(stdout,"\n");
1228	}
1229	for (i=0;cmd_line_para[i].str[0] != '\0';i++) {
1230	fprintf(stdout,"\t\t%c%s%s",COMMANDSEP,cmd_line_para[i].str,cmd_line_type[cmd_line_para[i].type]);
1231	if (cmd_line_para[i].type == OPTARG) {
1232	if (cmd_line_para[i].arg[0][0] != '\0')
1233	fprintf(stdout,"=%s",cmd_line_para[i].arg[0]);
1234	for (j=1;cmd_line_para[i].arg[j][0] != '\0';j++)
1235	fprintf(stdout," OR %s",cmd_line_para[i].arg[j]);
1236	}
1237	fprintf(stdout,"\n");
1238	}
1239	exit(1);
1240	}
1241
1242
1243	/*****************************************************************************/
1244	/* Check to see if sequence type is explicitely stated..override ************/
1245	/* the automatic checking (DNA or Protein). /type=d or /type=p *************/
1246	/*****************************************************************************/
1247	if(settype != -1)
1248	if(strlen(param_arg[settype])>0) {
1249	temp = find_match(param_arg[settype],type_arg,2);
1250	if(temp == 0) {
1251	dnaflag = FALSE;
1252	explicit_dnaflag = TRUE;
1253	info("Sequence type explicitly set to Protein");
1254	}
1255	else if(temp == 1) {
1256	info("Sequence type explicitly set to DNA");
1257	dnaflag = TRUE;
1258	explicit_dnaflag = TRUE;
1259	}
1260	else
1261	fprintf(stdout,"\nUnknown sequence type %s\n",
1262	param_arg[settype]);
1263	}
1264
1265
1266	/***************************************************************************
1267	* check to see if 1st parameter does not start with '/' i.e. look for an *
1268	* input file as first parameter. The input file can also be specified *
1269	* by /infile=fname. *
1270	****************************************************************************/
1271	/* JULIE - moved to check_param()
1272	if(paramstr[0] != '/') {
1273	strcpy(seqname, params[0]);
1274	}
1275	*/
1276
1277	/**************************************************/
1278	/* Look for /infile=file.ext on the command line */
1279	/**************************************************/
1280
1281	if(setinfile != -1) {
1282	if(strlen(param_arg[setinfile]) <= 0) {
1283	error("Bad sequence file name");
1284	exit(1);
1285	}
1286	strcpy(seqname, param_arg[setinfile]);
1287	}
1288
1289	if(*seqname != EOS) {
1290	profile_no = 0;
1291	nseqs = readseqs((sint)1);
1292	if(nseqs < 2) {
1293	fprintf(stdout,
1294	"\nNo. of seqs. read = %d. No alignment!\n",(pint)nseqs);
1295	exit(1);
1296	}
1297	for(i = 1; i<=nseqs; i++)
1298	info("Sequence %d: %-*s %6.d %s",
1299	(pint)i,max_names,names[i],(pint)seqlen_array[i],dnaflag?"bp":"aa");
1300	empty = FALSE;
1301	do_something = TRUE;
1302	}
1303
1304	set_optional_param();
1305
1306	/*********************************************************/
1307	/* Look for /profile1=file.ext AND /profile2=file2.ext */
1308	/* You must give both file names OR neither. */
1309	/*********************************************************/
1310
1311	if(setprofile1 != -1) {
1312	if(strlen(param_arg[setprofile1]) <= 0) {
1313	error("Bad profile 1 file name");
1314	exit(1);
1315	}
1316	strcpy(seqname, param_arg[setprofile1]);
1317	profile_no = 1;
1318	profile_input();
1319	if(nseqs <= 0)
1320	exit(1);
1321	strcpy(profile1_name,seqname);
1322	}
1323
1324	if(setprofile2 != -1) {
1325	if(strlen(param_arg[setprofile2]) <= 0) {
1326	error("Bad profile 2 file name");
1327	exit(1);
1328	}
1329	if(profile1_empty) {
1330	error("Only 1 profile file (profile 2) specified.");
1331	exit(1);
1332	}
1333	strcpy(seqname, param_arg[setprofile2]);
1334	profile_no = 2;
1335	profile_input();
1336	if(nseqs > profile1_nseqs)
1337	do_something = do_profile = TRUE;
1338	else {
1339	error("No sequences read from profile 2");
1340	exit(1);
1341	}
1342	strcpy(profile2_name,seqname);
1343	}
1344
1345	/*************************************************************************/
1346	/* Look for /tree or /bootstrap or /align or /usetree ******************/
1347	/*************************************************************************/
1348
1349	if (setbatch != -1)
1350	interactive=FALSE;
1351
1352	if (setinteractive != -1)
1353	interactive=TRUE;
1354
1355	if (interactive) {
1356	settree = -1;
1357	setbootstrap = -1;
1358	setalign = -1;
1359	setusetree = -1;
1360	setusetree1 = -1;
1361	setusetree2 = -1;
1362	setnewtree = -1;
1363	setconvert = -1;
1364	}
1365
1366	if(settree != -1 )
1367	if(empty) {
1368	error("Cannot draw tree. No input alignment file");
1369	exit(1);
1370	}
1371	else
1372	do_tree = TRUE;
1373
1374	if(setbootstrap != -1)
1375	if(empty) {
1376	error("Cannot bootstrap tree. No input alignment file");
1377	exit(1);
1378	}
1379	else {
1380	temp = 0;
1381	if(param_arg[setbootstrap] != NULL)
1382	if (sscanf(param_arg[setbootstrap],"%d",&temp)!=1) {
1383	fprintf(stdout,"Bad option for /bootstrap (must be integer)\n");
1384	temp = 0;
1385	};
1386	if(temp > 0) boot_ntrials = temp;
1387	do_boot = TRUE;
1388	}
1389
1390	if(setalign != -1)
1391	if(empty) {
1392	error("Cannot align sequences. No input file");
1393	exit(1);
1394	}
1395	else
1396	do_align = TRUE;
1397
1398	if(setconvert != -1)
1399	if(empty) {
1400	error("Cannot convert sequences. No input file");
1401	exit(1);
1402	}
1403	else
1404	do_convert = TRUE;
1405
1406	if(setusetree != -1)
1407	if(empty) {
1408	error("Cannot align sequences. No input file");
1409	exit(1);
1410	}
1411	else {
1412	if(strlen(param_arg[setusetree]) == 0) {
1413	error("Cannot align sequences. No tree file specified");
1414	exit(1);
1415	}
1416	else {
1417	strcpy(phylip_tree_name, param_arg[setusetree]);
1418	}
1419	use_tree_file = TRUE;
1420	do_align_only = TRUE;
1421	}
1422
1423	if(setnewtree != -1)
1424	if(empty) {
1425	error("Cannot align sequences. No input file");
1426	exit(1);
1427	}
1428	else {
1429	if(strlen(param_arg[setnewtree]) == 0) {
1430	error("Cannot align sequences. No tree file specified");
1431	exit(1);
1432	}
1433	else {
1434	strcpy(phylip_tree_name, param_arg[setnewtree]);
1435	}
1436	new_tree_file = TRUE;
1437	do_tree_only = TRUE;
1438	}
1439
1440	if(setusetree1 != -1)
1441	if(profile1_empty) {
1442	error("Cannot align profiles. No input file");
1443	exit(1);
1444	}
1445	else if(profile_type == SEQUENCE) {
1446	error("Invalid option /usetree1.");
1447	exit(1);
1448	}
1449	else {
1450	if(strlen(param_arg[setusetree1]) == 0) {
1451	error("Cannot align profiles. No tree file specified");
1452	exit(1);
1453	}
1454	else {
1455	strcpy(p1_tree_name, param_arg[setusetree1]);
1456	}
1457	use_tree1_file = TRUE;
1458	do_align_only = TRUE;
1459	}
1460
1461	if(setnewtree1 != -1)
1462	if(profile1_empty) {
1463	error("Cannot align profiles. No input file");
1464	exit(1);
1465	}
1466	else if(profile_type == SEQUENCE) {
1467	error("Invalid option /newtree1.");
1468	exit(1);
1469	}
1470	else {
1471	if(strlen(param_arg[setnewtree1]) == 0) {
1472	error("Cannot align profiles. No tree file specified");
1473	exit(1);
1474	}
1475	else {
1476	strcpy(p1_tree_name, param_arg[setnewtree1]);
1477	}
1478	new_tree1_file = TRUE;
1479	}
1480
1481	if(setusetree2 != -1)
1482	if(profile2_empty) {
1483	error("Cannot align profiles. No input file");
1484	exit(1);
1485	}
1486	else if(profile_type == SEQUENCE) {
1487	error("Invalid option /usetree2.");
1488	exit(1);
1489	}
1490	else {
1491	if(strlen(param_arg[setusetree2]) == 0) {
1492	error("Cannot align profiles. No tree file specified");
1493	exit(1);
1494	}
1495	else {
1496	strcpy(p2_tree_name, param_arg[setusetree2]);
1497	}
1498	use_tree2_file = TRUE;
1499	do_align_only = TRUE;
1500	}
1501
1502	if(setnewtree2 != -1)
1503	if(profile2_empty) {
1504	error("Cannot align profiles. No input file");
1505	exit(1);
1506	}
1507	else if(profile_type == SEQUENCE) {
1508	error("Invalid option /newtree2.");
1509	exit(1);
1510	}
1511	else {
1512	if(strlen(param_arg[setnewtree2]) == 0) {
1513	error("Cannot align profiles. No tree file specified");
1514	exit(1);
1515	}
1516	else {
1517	strcpy(p2_tree_name, param_arg[setnewtree2]);
1518	}
1519	new_tree2_file = TRUE;
1520	}
1521
1522
1523	if( (!do_tree) && (!do_boot) && (!empty) && (!do_profile) && (!do_align_only) && (!do_tree_only) && (!do_convert))
1524	do_align = TRUE;
1525
1526	/*** ? /quicktree */
1527	if(setquicktree != -1)
1528	quick_pairalign = TRUE;
1529
1530	if(dnaflag) {
1531	gap_open = dna_gap_open;
1532	gap_extend = dna_gap_extend;
1533	pw_go_penalty = dna_pw_go_penalty;
1534	pw_ge_penalty = dna_pw_ge_penalty;
1535	ktup = dna_ktup;
1536	window = dna_window;
1537	signif = dna_signif;
1538	wind_gap = dna_wind_gap;
1539
1540	}
1541	else {
1542	gap_open = prot_gap_open;
1543	gap_extend = prot_gap_extend;
1544	pw_go_penalty = prot_pw_go_penalty;
1545	pw_ge_penalty = prot_pw_ge_penalty;
1546	ktup = prot_ktup;
1547	window = prot_window;
1548	signif = prot_signif;
1549	wind_gap = prot_wind_gap;
1550	}
1551
1552	if(interactive) {
1553	if (!xmenus) usemenu = TRUE;
1554	return;
1555	}
1556
1557
1558	if(!do_something) {
1559	error("No input file(s) specified");
1560	exit(1);
1561	}
1562
1563
1564
1565
1566	/****************************************************************************/
1567	/* Now do whatever has been requested ***************************************/
1568	/****************************************************************************/
1569
1570
1571	if(do_profile) {
1572	if (profile_type == PROFILE) profile_align(p1_tree_name,p2_tree_name);
1573	else new_sequence_align(phylip_tree_name);
1574	}
1575
1576	else if(do_align)
1577	align(phylip_tree_name);
1578
1579	else if(do_convert) {
1580	get_path(seqname,path);
1581	if(!open_alignment_output(path)) exit(1);
1582	create_alignment_output(1,nseqs);
1583	}
1584
1585	else if (do_align_only)
1586	get_tree(phylip_tree_name);
1587
1588	else if(do_tree_only)
1589	make_tree(phylip_tree_name);
1590
1591	else if(do_tree)
1592	phylogenetic_tree(phylip_tree_name,clustal_tree_name,dist_tree_name);
1593
1594	else if(do_boot)
1595	bootstrap_tree(phylip_tree_name,clustal_tree_name);
1596
1597	fprintf(stdout,"\n");
1598	exit(1);
1599
1600	/*****whew!nowgohome***/
1601	}
1602
1603
1604	Boolean user_mat(char str, short mat, short *xref)
1605	{
1606	sint maxres;
1607
1608	FILE *infile;
1609
1610	if(usemenu)
1611	getstr("Enter name of the matrix file",lin2);
1612	else
1613	strcpy(lin2,str);
1614
1615	if(*lin2 == EOS) return FALSE;
1616
1617	if((infile=fopen(lin2,"r"))==NULL) {
1618	error("Cannot find matrix file [%s]",lin2);
1619	return FALSE;
1620	}
1621
1622	strcpy(str, lin2);
1623
1624	maxres = read_user_matrix(str, mat, xref);
1625	if (maxres <= 0) return FALSE;
1626
1627	return TRUE;
1628	}
1629
1630
1631
1632
1633
1634
1635	sint seq_input(Boolean append)
1636	{
1637	sint i;
1638	sint local_nseqs;
1639
1640	if(usemenu) {
1641	fprintf(stdout,"\n\nSequences should all be in 1 file.\n");
1642	fprintf(stdout,"\n7 formats accepted: \n");
1643	fprintf(stdout,
1644	"NBRF/PIR, EMBL/SwissProt, Pearson (Fasta), GDE, Clustal, GCG/MSF, RSF.\n\n\n");
1645	/*fprintf(stdout,
1646	"\nGCG users should use TOPIR to convert their sequence files before use.\n\n\n");*/
1647	}
1648
1649	if (append)
1650	local_nseqs = readseqs(nseqs+(sint)1);
1651	else
1652	local_nseqs = readseqs((sint)1); /* 1 is the first seq to be read */
1653	if(local_nseqs < 0) /* file could not be opened */
1654	{
1655	return local_nseqs;
1656	}
1657	else if(local_nseqs == 0) /* no sequences */
1658	{
1659	error("No sequences in file! Bad format?");
1660	return local_nseqs;
1661	}
1662	else
1663	{
1664	struct_penalties1 = struct_penalties2 = NONE;
1665	if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
1666	if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
1667	if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
1668	if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
1669	if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
1670	if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
1671
1672	if(append) nseqs+=local_nseqs;
1673	else nseqs=local_nseqs;
1674	info("Sequences assumed to be %s",
1675	dnaflag?"DNA":"PROTEIN");
1676	if (usemenu) {
1677	fprintf(stdout,"\n\n");
1678	for(i=1; i<=nseqs; i++) {
1679	/* DES fprintf(stdout,"%s: = ",names[i]); */
1680	info("Sequence %d: %-*s %6.d %s",
1681	(pint)i,max_names,names[i],(pint)seqlen_array[i],dnaflag?"bp":"aa");
1682	}
1683	}
1684	if(dnaflag) {
1685	gap_open = dna_gap_open;
1686	gap_extend = dna_gap_extend;
1687	}
1688	else {
1689	gap_open = prot_gap_open;
1690	gap_extend = prot_gap_extend;
1691	}
1692	empty=FALSE;
1693	}
1694	return local_nseqs;
1695	}
1696
1697
1698
1699
1700
1701
1702
1703	sint profile_input(void) /* read a profile */
1704	{ /* profile_no is 1 or 2 */
1705	sint local_nseqs, i;
1706
1707	if(profile_no == 2 && profile1_empty)
1708	{
1709	error("You must read in profile number 1 first");
1710	return 0;
1711	}
1712
1713	if(profile_no == 1) /* for the 1st profile */
1714	{
1715	local_nseqs = readseqs((sint)1); /* (1) means 1st seq to be read = no. 1 */
1716	if(local_nseqs == 0) /* no sequences */
1717	{
1718	error("No sequences in file! Bad format?");
1719	}
1720	else if (local_nseqs > 0)
1721	{ /* success; found some seqs. */
1722	struct_penalties1 = NONE;
1723	if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
1724	if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
1725	if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
1726	if (struct_penalties != NONE) /* feature table / mask in alignment */
1727	{
1728	struct_penalties1 = struct_penalties;
1729	if (struct_penalties == SECST) {
1730	sec_struct_mask1 = (char )ckalloc((max_aln_length) sizeof (char));
1731	for (i=0;i<max_aln_length;i++)
1732	sec_struct_mask1[i] = sec_struct_mask[i];
1733	}
1734	gap_penalty_mask1 = (char )ckalloc((max_aln_length) sizeof (char));
1735	for (i=0;i<max_aln_length;i++)
1736	gap_penalty_mask1[i] = gap_penalty_mask[i];
1737	ss_name1 = (char )ckalloc( (MAXNAMES+1) sizeof (char));
1738
1739	strcpy(ss_name1,ss_name);
1740	if (debug>0) {
1741	for (i=0;i<seqlen_array[1];i++)
1742	fprintf(stdout,"%c",gap_penalty_mask1[i]);
1743	fprintf(stdout,"\n");
1744	}
1745	}
1746	nseqs = profile1_nseqs = local_nseqs;
1747	info("No. of seqs=%d",(pint)nseqs);
1748	profile1_empty=FALSE;
1749	profile2_empty=TRUE;
1750	}
1751	}
1752	else
1753	{ /* first seq to be read = profile1_nseqs + 1 */
1754	local_nseqs = readseqs(profile1_nseqs+(sint)1);
1755	if(local_nseqs == 0) /* no sequences */
1756	{
1757	error("No sequences in file! Bad format?");
1758	}
1759	else if(local_nseqs > 0)
1760	{
1761	struct_penalties2 = NONE;
1762	if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
1763	if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
1764	if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
1765	if (struct_penalties != NONE) /* feature table / mask in alignment */
1766	{
1767	struct_penalties2 = struct_penalties;
1768	if (struct_penalties == SECST) {
1769	sec_struct_mask2 = (char )ckalloc((max_aln_length) sizeof (char));
1770	for (i=0;i<max_aln_length;i++)
1771	sec_struct_mask2[i] = sec_struct_mask[i];
1772	}
1773	gap_penalty_mask2 = (char )ckalloc((max_aln_length) sizeof (char));
1774	for (i=0;i<max_aln_length;i++)
1775	gap_penalty_mask2[i] = gap_penalty_mask[i];
1776	ss_name2 = (char )ckalloc( (MAXNAMES+1) sizeof (char));
1777	strcpy(ss_name2,ss_name);
1778	if (debug>0) {
1779	for (i=0;i<seqlen_array[profile1_nseqs+1];i++)
1780	fprintf(stdout,"%c",gap_penalty_mask2[i]);
1781	fprintf(stdout,"\n");
1782	}
1783	}
1784	info("No. of seqs in profile=%d",(pint)local_nseqs);
1785	nseqs = profile1_nseqs + local_nseqs;
1786	info("Total no. of seqs =%d",(pint)nseqs);
1787	profile2_empty=FALSE;
1788	empty = FALSE;
1789	}
1790
1791	}
1792	if (sec_struct_mask != NULL) sec_struct_mask=ckfree(sec_struct_mask);
1793	if (gap_penalty_mask != NULL) gap_penalty_mask=ckfree(gap_penalty_mask);
1794	if (ss_name != NULL) ss_name=ckfree(ss_name);
1795
1796	if(local_nseqs<=0) return local_nseqs;
1797
1798	info("Sequences assumed to be %s",
1799	dnaflag?"DNA":"PROTEIN");
1800	if (usemenu) fprintf(stdout,"\n\n");
1801	for(i=profile2_empty?1:profile1_nseqs+1; i<=nseqs; i++) {
1802	info("Sequence %d: %-*s %6.d %s",
1803	(pint)i,max_names,names[i],(pint)seqlen_array[i],dnaflag?"bp":"aa");
1804	}
1805	if(dnaflag) {
1806	gap_open = dna_gap_open;
1807	gap_extend = dna_gap_extend;
1808	}
1809	else {
1810	gap_open = prot_gap_open;
1811	gap_extend = prot_gap_extend;
1812	}
1813
1814	return nseqs;
1815	}
1816
1817
1818
1819	static void calc_gap_penalty_mask(int prf_length, char mask, char gap_mask)
1820	{
1821	int i,j;
1822	char *struct_mask;
1823
1824	struct_mask = (char )ckalloc((prf_length+1) sizeof(char));
1825	/*
1826	calculate the gap penalty mask from the secondary structures
1827	*/
1828	i=0;
1829	while (i<prf_length) {
1830	if (tolower(mask[i]) == 'a' \|\| mask[i] == '$') {
1831	for (j = -helix_end_plus; j<0; j++) {
1832	if ((i+j>=0) && (tolower(struct_mask[i+j]) != 'a')
1833	&& (tolower(struct_mask[i+j]) != 'b'))
1834	struct_mask[i+j] = 'a';
1835	}
1836	for (j = 0; j<helix_end_minus; j++) {
1837	if (i+j>=prf_length \|\| (tolower(mask[i+j]) != 'a'
1838	&& mask[i+j] != '$')) break;
1839	struct_mask[i+j] = 'a';
1840	}
1841	i += j;
1842	while (tolower(mask[i]) == 'a'
1843	\|\| mask[i] == '$') {
1844	if (i>=prf_length) break;
1845	if (mask[i] == '$') {
1846	struct_mask[i] = 'A';
1847	i++;
1848	break;
1849	}
1850	else struct_mask[i] = mask[i];
1851	i++;
1852	}
1853	for (j = 0; j<helix_end_minus; j++) {
1854	if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'a'
1855	\|\| mask[i-j-1] == '$'))
1856	struct_mask[i-j-1] = 'a';
1857	}
1858	for (j = 0; j<helix_end_plus; j++) {
1859	if (i+j>=prf_length) break;
1860	struct_mask[i+j] = 'a';
1861	}
1862	}
1863	else if (tolower(mask[i]) == 'b' \|\| mask[i] == '%') {
1864	for (j = -strand_end_plus; j<0; j++) {
1865	if ((i+j>=0) && (tolower(struct_mask[i+j]) != 'a')
1866	&& (tolower(struct_mask[i+j]) != 'b'))
1867	struct_mask[i+j] = 'b';
1868	}
1869	for (j = 0; j<strand_end_minus; j++) {
1870	if (i+j>=prf_length \|\| (tolower(mask[i+j]) != 'b'
1871	&& mask[i+j] != '%')) break;
1872	struct_mask[i+j] = 'b';
1873	}
1874	i += j;
1875	while (tolower(mask[i]) == 'b'
1876	\|\| mask[i] == '%') {
1877	if (i>=prf_length) break;
1878	if (mask[i] == '%') {
1879	struct_mask[i] = 'B';
1880	i++;
1881	break;
1882	}
1883	else struct_mask[i] = mask[i];
1884	i++;
1885	}
1886	for (j = 0; j<strand_end_minus; j++) {
1887	if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'b'
1888	\|\| mask[i-j-1] == '%'))
1889	struct_mask[i-j-1] = 'b';
1890	}
1891	for (j = 0; j<strand_end_plus; j++) {
1892	if (i+j>=prf_length) break;
1893	struct_mask[i+j] = 'b';
1894	}
1895	}
1896	else i++;
1897	}
1898
1899	for(i=0;i<prf_length;i++) {
1900	switch (struct_mask[i]) {
1901	case 'A':
1902	gap_mask[i] = helix_penalty+'0';
1903	break;
1904	case 'a':
1905	gap_mask[i] = helix_end_penalty+'0';
1906	break;
1907	case 'B':
1908	gap_mask[i] = strand_penalty+'0';
1909	break;
1910	case 'b':
1911	gap_mask[i] = strand_end_penalty+'0';
1912	break;
1913	default:
1914	gap_mask[i] = loop_penalty+'0';
1915	break;
1916	}
1917	}
1918
1919	struct_mask=ckfree(struct_mask);
1920
1921	}
1922
1923	void print_sec_struct_mask(int prf_length, char mask, char struct_mask)
1924	{
1925	int i,j;
1926
1927	/*
1928	calculate the gap penalty mask from the secondary structures
1929	*/
1930	i=0;
1931	while (i<prf_length) {
1932	if (tolower(mask[i]) == 'a' \|\| mask[i] == '$') {
1933	for (j = 0; j<helix_end_minus; j++) {
1934	if (i+j>=prf_length \|\| (tolower(mask[i+j]) != 'a'
1935	&& mask[i+j] != '$')) break;
1936	struct_mask[i+j] = 'a';
1937	}
1938	i += j;
1939	while (tolower(mask[i]) == 'a'
1940	\|\| mask[i] == '$') {
1941	if (i>=prf_length) break;
1942	if (mask[i] == '$') {
1943	struct_mask[i] = 'A';
1944	i++;
1945	break;
1946	}
1947	else struct_mask[i] = mask[i];
1948	i++;
1949	}
1950	for (j = 0; j<helix_end_minus; j++) {
1951	if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'a'
1952	\|\| mask[i-j-1] == '$'))
1953	struct_mask[i-j-1] = 'a';
1954	}
1955	}
1956	else if (tolower(mask[i]) == 'b' \|\| mask[i] == '%') {
1957	for (j = 0; j<strand_end_minus; j++) {
1958	if (i+j>=prf_length \|\| (tolower(mask[i+j]) != 'b'
1959	&& mask[i+j] != '%')) break;
1960	struct_mask[i+j] = 'b';
1961	}
1962	i += j;
1963	while (tolower(mask[i]) == 'b'
1964	\|\| mask[i] == '%') {
1965	if (i>=prf_length) break;
1966	if (mask[i] == '%') {
1967	struct_mask[i] = 'B';
1968	i++;
1969	break;
1970	}
1971	else struct_mask[i] = mask[i];
1972	i++;
1973	}
1974	for (j = 0; j<strand_end_minus; j++) {
1975	if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'b'
1976	\|\| mask[i-j-1] == '%'))
1977	struct_mask[i-j-1] = 'b';
1978	}
1979	}
1980	else i++;
1981	}
1982	}
1983
1984
1985
1986	FILE * open_output_file(char prompt, char path,
1987	char file_name, char file_extension)
1988
1989	{ static char temp[FILENAMELEN+1];
1990	static char local_prompt[MAXLINE];
1991	FILE * file_handle;
1992
1993	/* if (*file_name == EOS) {
1994	*/ strcpy(file_name,path);
1995	strcat(file_name,file_extension);
1996	/* }
1997	*/
1998	if(strcmp(file_name,seqname)==0) {
1999	warning("Output file name is the same as input file.");
2000	if (usemenu) {
2001	strcpy(local_prompt,"\n\nEnter new name to avoid overwriting ");
2002	strcat(local_prompt," [%s]: ");
2003	fprintf(stdout,local_prompt,file_name);
2004	gets(temp);
2005	if(*temp != EOS) strcpy(file_name,temp);
2006	}
2007	}
2008	else if (usemenu) {
2009	strcpy(local_prompt,prompt);
2010	strcat(local_prompt," [%s]: ");
2011	fprintf(stdout,local_prompt,file_name);
2012	gets(temp);
2013	if(*temp != EOS) strcpy(file_name,temp);
2014	}
2015
2016	#ifdef VMS
2017	if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
2018	#else
2019	if((file_handle=fopen(file_name,"w"))==NULL) {
2020	#endif
2021	error("Cannot open output file [%s]",file_name);
2022	return NULL;
2023	}
2024	return file_handle;
2025	}
2026
2027
2028
2029	FILE * open_explicit_file(char *file_name)
2030	{
2031	FILE * file_handle;
2032
2033	if (*file_name == EOS) {
2034	error("Bad output file [%s]",file_name);
2035	return NULL;
2036	}
2037	#ifdef VMS
2038	if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
2039	#else
2040	if((file_handle=fopen(file_name,"w"))==NULL) {
2041	#endif
2042	error("Cannot open output file [%s]",file_name);
2043	return NULL;
2044	}
2045	return file_handle;
2046	}
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056	void align(char *phylip_name)
2057	{
2058	char path[FILENAMELEN+1];
2059	FILE *tree;
2060	sint count;
2061
2062	if(empty && usemenu) {
2063	error("No sequences in memory. Load sequences first.");
2064	return;
2065	}
2066
2067	struct_penalties1 = struct_penalties2 = NONE;
2068	if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
2069	if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
2070	if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
2071	if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
2072	if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
2073	if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
2074
2075
2076	get_path(seqname,path);
2077	/* DES DEBUG
2078	fprintf(stdout,"\n\n Seqname = %s \n Path = %s \n\n",seqname,path);
2079	*/
2080	if(usemenu \|\| !interactive) {
2081	if(!open_alignment_output(path)) return;
2082	}
2083
2084	if (nseqs > 3) {
2085
2086	get_path(seqname,path);
2087	if (phylip_name[0]!=EOS) {
2088	if((tree = open_explicit_file(
2089	phylip_name))==NULL) return;
2090	}
2091	else {
2092	if((tree = open_output_file(
2093	"\nEnter name for new GUIDE TREE file ",path,
2094	phylip_name,"dnd")) == NULL) return;
2095	}
2096	}
2097
2098	if (save_parameters) create_parameter_output();
2099
2100	if(reset_alignments) reset_align();
2101
2102	info("Start of Pairwise alignments");
2103	info("Aligning...");
2104	if(dnaflag) {
2105	gap_open = dna_gap_open;
2106	gap_extend = dna_gap_extend;
2107	pw_go_penalty = dna_pw_go_penalty;
2108	pw_ge_penalty = dna_pw_ge_penalty;
2109	ktup = dna_ktup;
2110	window = dna_window;
2111	signif = dna_signif;
2112	wind_gap = dna_wind_gap;
2113
2114	}
2115	else {
2116	gap_open = prot_gap_open;
2117	gap_extend = prot_gap_extend;
2118	pw_go_penalty = prot_pw_go_penalty;
2119	pw_ge_penalty = prot_pw_ge_penalty;
2120	ktup = prot_ktup;
2121	window = prot_window;
2122	signif = prot_signif;
2123	wind_gap = prot_wind_gap;
2124
2125	}
2126
2127	if (quick_pairalign)
2128	show_pair();
2129	else
2130	pairalign((sint)0,nseqs,(sint)0,nseqs);
2131
2132	if (nseqs > 3) {
2133
2134	guide_tree(tree,1,nseqs);
2135	info("Guide tree file created: [%s]",
2136	phylip_name);
2137	}
2138
2139
2140	count = malign((sint)0,phylip_name);
2141
2142	if (count <= 0) return;
2143
2144	if (usemenu) fprintf(stdout,"\n\n\n");
2145
2146	create_alignment_output(1,nseqs);
2147	if (showaln && usemenu) show_aln();
2148	phylip_name[0]=EOS;
2149	}
2150
2151
2152
2153
2154
2155	void new_sequence_align(char *phylip_name)
2156	{
2157	char path[FILENAMELEN+1];
2158	char tree_name[FILENAMELEN+1],temp[MAXLINE+1];
2159	Boolean use_tree;
2160	FILE *tree;
2161	sint i,j,count;
2162	float dscore;
2163	Boolean save_ss2;
2164
2165	if(profile1_empty && usemenu) {
2166	error("No profile in memory. Input 1st profile first.");
2167	return;
2168	}
2169
2170	if(profile2_empty && usemenu) {
2171	error("No sequences in memory. Input sequences first.");
2172	return;
2173	}
2174
2175	get_path(profile2_name,path);
2176
2177	if(usemenu \|\| !interactive) {
2178	if(!open_alignment_output(path)) return;
2179	}
2180
2181	new_seq = profile1_nseqs+1;
2182
2183	/* check for secondary structure information for list of sequences */
2184
2185	save_ss2 = use_ss2;
2186	if (struct_penalties2 != NONE && use_ss2 == TRUE && (nseqs - profile1_nseqs >
2187	1)) {
2188	if (struct_penalties2 == SECST)
2189	warning("Warning: ignoring secondary structure for a list of sequences");
2190	else if (struct_penalties2 == GMASK)
2191	warning("Warning: ignoring gap penalty mask for a list of sequences");
2192	use_ss2 = FALSE;
2193	}
2194
2195	for (i=1;i<=new_seq;i++) {
2196	for (j=i+1;j<=new_seq;j++) {
2197	dscore = countid(i,j);
2198	tmat[i][j] = ((double)100.0 - (double)dscore)/(double)100.0;
2199	tmat[j][i] = tmat[i][j];
2200	}
2201	}
2202
2203	tree_name[0] = EOS;
2204	use_tree = FALSE;
2205	if (nseqs > 3) {
2206	if (check_tree && usemenu) {
2207	strcpy(tree_name,path);
2208	strcat(tree_name,"dnd");
2209	#ifdef VMS
2210	if((tree=fopen(tree_name,"r","rat=cr","rfm=var"))!=NULL) {
2211	#else
2212	if((tree=fopen(tree_name,"r"))!=NULL) {
2213	#endif
2214	if (usemenu)
2215	fprintf(stdout,"\nUse the existing GUIDE TREE file, %s (y/n) ? [y]: ",
2216	tree_name);
2217	gets(temp);
2218	if(temp != 'n' && temp != 'N') {
2219	strcpy(phylip_name,tree_name);
2220	use_tree = TRUE;
2221	}
2222	fclose(tree);
2223	}
2224	}
2225	else if (!usemenu && use_tree_file) {
2226	use_tree = TRUE;
2227	}
2228	}
2229
2230	if (save_parameters) create_parameter_output();
2231
2232	if(reset_alignments) {
2233	reset_prf1();
2234	reset_prf2();
2235	}
2236	else fix_gaps();
2237
2238	if (struct_penalties1 == SECST)
2239
2240	calc_gap_penalty_mask(seqlen_array[1],sec_struct_mask1,gap_penalty_mask1);
2241
2242	if (struct_penalties2 == SECST)
2243
2244	calc_gap_penalty_mask(seqlen_array[profile1_nseqs+1],sec_struct_mask2,gap_penalty_mask2);
2245
2246
2247	/* create the new tree file, if necessary */
2248
2249	if (use_tree == FALSE) {
2250
2251	if (nseqs > 3) {
2252	get_path(profile2_name,path);
2253	if (phylip_name[0]!=EOS) {
2254	if((tree = open_explicit_file(
2255	phylip_name))==NULL) return;
2256	}
2257	else {
2258	if((tree = open_output_file(
2259	"\nEnter name for new GUIDE TREE file ",path,
2260	phylip_name,"dnd")) == NULL) return;
2261	}
2262	}
2263	info("Start of Pairwise alignments");
2264	info("Aligning...");
2265	if(dnaflag) {
2266	gap_open = dna_gap_open;
2267	gap_extend = dna_gap_extend;
2268	pw_go_penalty = dna_pw_go_penalty;
2269	pw_ge_penalty = dna_pw_ge_penalty;
2270	ktup = dna_ktup;
2271	window = dna_window;
2272	signif = dna_signif;
2273	wind_gap = dna_wind_gap;
2274
2275	}
2276	else {
2277	gap_open = prot_gap_open;
2278	gap_extend = prot_gap_extend;
2279	pw_go_penalty = prot_pw_go_penalty;
2280	pw_ge_penalty = prot_pw_ge_penalty;
2281	ktup = prot_ktup;
2282	window = prot_window;
2283	signif = prot_signif;
2284	wind_gap = prot_wind_gap;
2285
2286	}
2287
2288	if (quick_pairalign)
2289	show_pair();
2290	else
2291	pairalign((sint)0,nseqs,new_seq-2,nseqs);
2292
2293	if (nseqs > 3) {
2294	guide_tree(tree,1,nseqs);
2295	info("Guide tree file created: [%s]",
2296	phylip_name);
2297	}
2298	}
2299
2300	if (new_tree_file) return;
2301
2302	count = seqalign(new_seq-2,phylip_name);
2303
2304	use_ss2 = save_ss2;
2305
2306	if (count <= 0) return;
2307
2308	if (usemenu) fprintf(stdout,"\n\n\n");
2309
2310	create_alignment_output(1,nseqs);
2311	if (showaln && usemenu) show_aln();
2312	phylip_name[0]=EOS;
2313
2314	}
2315
2316
2317
2318
2319
2320	void make_tree(char *phylip_name)
2321	{
2322	char path[FILENAMELEN+1];
2323	FILE *tree;
2324
2325	if(empty) {
2326	error("No sequences in memory. Load sequences first.");
2327	return;
2328	}
2329
2330	struct_penalties1 = struct_penalties2 = NONE;
2331	if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
2332	if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
2333	if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
2334	if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
2335	if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
2336	if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
2337
2338	if(reset_alignments) reset_align();
2339
2340	get_path(seqname,path);
2341
2342	if (nseqs <= 3) {
2343	error("Less than 3 sequences in memory. Phylogenetic tree cannot be built.");
2344	return;
2345	}
2346
2347	if (save_parameters) create_parameter_output();
2348
2349	info("Start of Pairwise alignments");
2350	info("Aligning...");
2351	if(dnaflag) {
2352	gap_open = dna_gap_open;
2353	gap_extend = dna_gap_extend;
2354	pw_go_penalty = dna_pw_go_penalty;
2355	pw_ge_penalty = dna_pw_ge_penalty;
2356	ktup = dna_ktup;
2357	window = dna_window;
2358	signif = dna_signif;
2359	wind_gap = dna_wind_gap;
2360
2361	}
2362	else {
2363	gap_open = prot_gap_open;
2364	gap_extend = prot_gap_extend;
2365	pw_go_penalty = prot_pw_go_penalty;
2366	pw_ge_penalty = prot_pw_ge_penalty;
2367	ktup = prot_ktup;
2368	window = prot_window;
2369	signif = prot_signif;
2370	wind_gap = prot_wind_gap;
2371
2372
2373	}
2374
2375	if (quick_pairalign)
2376	show_pair();
2377	else
2378	pairalign((sint)0,nseqs,(sint)0,nseqs);
2379
2380	if (nseqs > 3) {
2381	get_path(seqname,path);
2382	if (phylip_name[0]!=EOS) {
2383	if((tree = open_explicit_file(
2384	phylip_name))==NULL) return;
2385	}
2386	else {
2387	if((tree = open_output_file(
2388	"\nEnter name for new GUIDE TREE file ",path,
2389	phylip_name,"dnd")) == NULL) return;
2390	}
2391
2392	guide_tree(tree,1,nseqs);
2393	info("Guide tree file created: [%s]",
2394	phylip_name);
2395	}
2396
2397	if(reset_alignments) reset_align();
2398
2399	phylip_name[0]=EOS;
2400	}
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410	void get_tree(char *phylip_name)
2411	{
2412	char path[FILENAMELEN+1],temp[MAXLINE+1];
2413	sint count;
2414
2415	if(empty) {
2416	error("No sequences in memory. Load sequences first.");
2417	return;
2418	}
2419	struct_penalties1 = struct_penalties2 = NONE;
2420	if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
2421	if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
2422	if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
2423	if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
2424	if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
2425	if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
2426
2427
2428	get_path(seqname,path);
2429
2430	if(usemenu \|\| !interactive) {
2431	if(!open_alignment_output(path)) return;
2432	}
2433
2434	if(reset_alignments) reset_align();
2435
2436	get_path(seqname,path);
2437
2438	if (nseqs > 3) {
2439
2440	if(usemenu) {
2441	strcpy(phylip_name,path);
2442	strcat(phylip_name,"dnd");
2443
2444	fprintf(stdout,"\nEnter a name for the guide tree file [%s]: ",
2445	phylip_name);
2446	gets(temp);
2447	if(*temp != EOS)
2448	strcpy(phylip_name,temp);
2449	}
2450
2451	if(usemenu \|\| !interactive) {
2452	#ifdef VMS
2453	if((tree=fopen(phylip_name,"r","rat=cr","rfm=var"))==NULL) {
2454	#else
2455	if((tree=fopen(phylip_name,"r"))==NULL) {
2456	#endif
2457	error("Cannot open tree file [%s]",phylip_name);
2458	return;
2459	}
2460	}
2461	}
2462	else {
2463	info("Start of Pairwise alignments");
2464	info("Aligning...");
2465	if(dnaflag) {
2466	gap_open = dna_gap_open;
2467	gap_extend = dna_gap_extend;
2468	pw_go_penalty = dna_pw_go_penalty;
2469	pw_ge_penalty = dna_pw_ge_penalty;
2470	ktup = dna_ktup;
2471	window = dna_window;
2472	signif = dna_signif;
2473	wind_gap = dna_wind_gap;
2474
2475	}
2476	else {
2477	gap_open = prot_gap_open;
2478	gap_extend = prot_gap_extend;
2479	pw_go_penalty = prot_pw_go_penalty;
2480	pw_ge_penalty = prot_pw_ge_penalty;
2481	ktup = prot_ktup;
2482	window = prot_window;
2483	signif = prot_signif;
2484	wind_gap = prot_wind_gap;
2485
2486	}
2487
2488	if (quick_pairalign)
2489	show_pair();
2490	else
2491	pairalign((sint)0,nseqs,(sint)0,nseqs);
2492	}
2493
2494	if (save_parameters) create_parameter_output();
2495
2496	count = malign(0,phylip_name);
2497	if (count <= 0) return;
2498
2499	if (usemenu) fprintf(stdout,"\n\n\n");
2500
2501	create_alignment_output(1,nseqs);
2502	if (showaln && usemenu) show_aln();
2503
2504	phylip_name[0]=EOS;
2505	}
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
2516	void profile_align(char p1_tree_name,char p2_tree_name)
2517	{
2518	char path[FILENAMELEN+1];
2519	char tree_name[FILENAMELEN+1];
2520	char temp[MAXLINE+1];
2521	Boolean use_tree1,use_tree2;
2522	FILE *tree;
2523	sint count,i,j,dscore;
2524
2525	if(profile1_empty \|\| profile2_empty) {
2526	error("No sequences in memory. Load sequences first.");
2527	return;
2528	}
2529
2530	get_path(profile1_name,path);
2531
2532	if(usemenu \|\| !interactive) {
2533	if(!open_alignment_output(path)) return;
2534	}
2535
2536	if(reset_alignments) {
2537	reset_prf1();
2538	reset_prf2();
2539	}
2540	else fix_gaps();
2541
2542	tree_name[0] = EOS;
2543	use_tree1 = FALSE;
2544	if (profile1_nseqs > 3) {
2545	if (check_tree && usemenu) {
2546	strcpy(tree_name,path);
2547	strcat(tree_name,"dnd");
2548	#ifdef VMS
2549	if((tree=fopen(tree_name,"r","rat=cr","rfm=var"))!=NULL) {
2550	#else
2551	if((tree=fopen(tree_name,"r"))!=NULL) {
2552	#endif
2553	fprintf(stdout,"\nUse the existing GUIDE TREE file for Profile 1, %s (y/n) ? [y]: ",
2554	tree_name);
2555	gets(temp);
2556	if(temp != 'n' && temp != 'N') {
2557	strcpy(p1_tree_name,tree_name);
2558	use_tree1 = TRUE;
2559	}
2560	fclose(tree);
2561	}
2562	}
2563	else if (!usemenu && use_tree1_file) {
2564	use_tree1 = TRUE;
2565	}
2566	}
2567	tree_name[0] = EOS;
2568	use_tree2 = FALSE;
2569	get_path(profile2_name,path);
2570	if (nseqs-profile1_nseqs > 3) {
2571	if (check_tree && usemenu) {
2572	strcpy(tree_name,path);
2573	strcat(tree_name,"dnd");
2574	#ifdef VMS
2575	if((tree=fopen(tree_name,"r","rat=cr","rfm=var"))!=NULL) {
2576	#else
2577	if((tree=fopen(tree_name,"r"))!=NULL) {
2578	#endif
2579	fprintf(stdout,"\nUse the existing GUIDE TREE file for Profile 2, %s (y/n) ? [y]: ",
2580	tree_name);
2581	gets(temp);
2582	if(temp != 'n' && temp != 'N') {
2583	strcpy(p2_tree_name,tree_name);
2584	use_tree2 = TRUE;
2585	}
2586	fclose(tree);
2587	}
2588	}
2589	else if (!usemenu && use_tree2_file) {
2590	use_tree2 = TRUE;
2591	}
2592	}
2593
2594	if (save_parameters) create_parameter_output();
2595
2596	if (struct_penalties1 == SECST)
2597
2598	calc_gap_penalty_mask(seqlen_array[1],sec_struct_mask1,gap_penalty_mask1);
2599
2600	if (struct_penalties2 == SECST)
2601
2602	calc_gap_penalty_mask(seqlen_array[profile1_nseqs+1],sec_struct_mask2,gap_penalty_mask2);
2603
2604	if (use_tree1 == FALSE)
2605	if (profile1_nseqs > 3) {
2606	for (i=1;i<=profile1_nseqs;i++) {
2607	for (j=i+1;j<=profile1_nseqs;j++) {
2608	dscore = countid(i,j);
2609	tmat[i][j] = (100.0 - dscore)/100.0;
2610	tmat[j][i] = tmat[i][j];
2611	}
2612	}
2613	get_path(profile1_name,path);
2614	if (p1_tree_name[0]!=EOS) {
2615	if((tree = open_explicit_file(p1_tree_name))==NULL) return;
2616	}
2617	else {
2618	if((tree = open_output_file(
2619	"\nEnter name for new GUIDE TREE file for profile 1 ",path,
2620	p1_tree_name,"dnd")) == NULL) return;
2621	}
2622
2623	guide_tree(tree,1,profile1_nseqs);
2624	info("Guide tree file created: [%s]",
2625	p1_tree_name);
2626	}
2627	if (use_tree2 == FALSE)
2628	if(nseqs-profile1_nseqs > 3) {
2629	for (i=1;i<=nseqs-profile1_nseqs;i++) {
2630	for (j=i+1;j<=nseqs-profile1_nseqs;j++) {
2631	dscore = countid(i+profile1_nseqs,j+profile1_nseqs);
2632	tmat[i][j] = (100.0 - dscore)/100.0;
2633	tmat[j][i] = tmat[i][j];
2634	}
2635	}
2636	if (p2_tree_name[0]!=EOS) {
2637	if((tree = open_explicit_file(p2_tree_name))==NULL) return;
2638	}
2639	else {
2640	get_path(profile2_name,path);
2641	if((tree = open_output_file(
2642	"\nEnter name for new GUIDE TREE file for profile 2 ",path,
2643	p2_tree_name,"dnd")) == NULL) return;
2644	}
2645	guide_tree(tree,profile1_nseqs+1,nseqs-profile1_nseqs);
2646	info("Guide tree file created: [%s]",
2647	p2_tree_name);
2648	}
2649
2650	if (new_tree1_file \|\| new_tree2_file) return;
2651
2652	/* do an initial alignment to get the pairwise identities between the two
2653	profiles - used to set parameters for the final alignment */
2654	count = palign1();
2655	if (count == 0) return;
2656
2657	reset_prf1();
2658	reset_prf2();
2659
2660	count = palign2(p1_tree_name,p2_tree_name);
2661
2662	if (count == 0) return;
2663
2664	if(usemenu) fprintf(stdout,"\n\n\n");
2665
2666	create_alignment_output(1,nseqs);
2667	if (showaln && usemenu) show_aln();
2668
2669	p1_tree_name[0]=EOS;
2670	p2_tree_name[0]=EOS;
2671	}
2672
2673
2674
2675
2676
2677	void clustal_out(FILE *clusout, sint fres, sint len, sint fseq, sint lseq)
2678	{
2679	static char *seq1;
2680	static sint *seq_no;
2681	static sint *print_seq_no;
2682	char ss_mask1, ss_mask2;
2683	char temp[MAXLINE];
2684	char c;
2685	sint val;
2686	sint ii,lv1,catident1[NUMRES],catident2[NUMRES],ident,chunks;
2687	sint i,j,k,l;
2688	sint pos,ptr;
2689	sint line_length;
2690
2691	/*
2692	stop doing this ...... opens duplicate files in VMS DES
2693	fclose(clusout);
2694	if ((clusout=fopen(clustal_outname,"w")) == NULL)
2695	{
2696	fprintf(stdout,"Error opening %s\n",clustal_outfile);
2697	return;
2698	}
2699	*/
2700	seq_no = (sint )ckalloc((nseqs+1) sizeof(sint));
2701	print_seq_no = (sint )ckalloc((nseqs+1) sizeof(sint));
2702	for (i=fseq;i<=lseq;i++)
2703	print_seq_no[i] = seq_no[i] = 0;
2704
2705	seq1 = (char )ckalloc((max_aln_length+1) sizeof(char));
2706
2707	if (struct_penalties1 == SECST && use_ss1 == TRUE) {
2708	ss_mask1 = (char )ckalloc((seqlen_array[1]+10) sizeof(char));
2709	for (i=0;i<seqlen_array[1];i++)
2710	ss_mask1[i] = sec_struct_mask1[i];
2711	print_sec_struct_mask(seqlen_array[1],sec_struct_mask1,ss_mask1);
2712	}
2713	if (struct_penalties2 == SECST && use_ss2 == TRUE) {
2714	ss_mask2 = (char )ckalloc((seqlen_array[profile1_nseqs+1]+10) sizeof(char));
2715	for (i=0;i<seqlen_array[profile1_nseqs+1];i++)
2716	ss_mask2[i] = sec_struct_mask2[i];
2717	print_sec_struct_mask(seqlen_array[profile1_nseqs+1],sec_struct_mask2,ss_mask2);
2718	}
2719
2720	fprintf(clusout,"CLUSTAL %s multiple sequence alignment\n\n",
2721	revision_level);
2722
2723	/* decide the line length for this alignment - maximum is LINELENGTH */
2724	line_length=PAGEWIDTH-max_names;
2725	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
2726	if (line_length > LINELENGTH) line_length=LINELENGTH;
2727
2728	chunks = len/line_length;
2729	if(len % line_length != 0)
2730	++chunks;
2731
2732	for(lv1=1;lv1<=chunks;++lv1) {
2733	pos = ((lv1-1)*line_length)+1;
2734	ptr = (len<pos+line_length-1) ? len : pos+line_length-1;
2735
2736	fprintf(clusout,"\n");
2737
2738	if (output_struct_penalties == 0 \|\| output_struct_penalties == 2) {
2739	if (struct_penalties1 == SECST && use_ss1 == TRUE) {
2740	for(i=pos;i<=ptr;++i) {
2741	val=ss_mask1[i+fres-2];
2742	if (val == gap_pos1)
2743	temp[i-pos]='-';
2744	else
2745	temp[i-pos]=val;
2746	}
2747	temp[ptr-pos+1]=EOS;
2748	fprintf(clusout,"!SS_%-*s %s\n",max_names,ss_name1,temp);
2749	}
2750	}
2751	if (output_struct_penalties == 1 \|\| output_struct_penalties == 2) {
2752	if (struct_penalties1 != NONE && use_ss1 == TRUE) {
2753	for(i=pos;i<=ptr;++i) {
2754	val=gap_penalty_mask1[i+fres-2];
2755	if (val == gap_pos1)
2756	temp[i-pos]='-';
2757	else
2758	temp[i-pos]=val;
2759	}
2760	temp[ptr-pos+1]=EOS;
2761	fprintf(clusout,"!GM_%-*s %s\n",max_names,ss_name1,temp);
2762	}
2763	}
2764	if (output_struct_penalties == 0 \|\| output_struct_penalties == 2) {
2765	if (struct_penalties2 == SECST && use_ss2 == TRUE) {
2766	for(i=pos;i<=ptr;++i) {
2767	val=ss_mask2[i+fres-2];
2768	if (val == gap_pos1)
2769	temp[i-pos]='-';
2770	else
2771	temp[i-pos]=val;
2772	}
2773	temp[ptr-pos+1]=EOS;
2774	fprintf(clusout,"!SS_%-*s %s\n",max_names,ss_name2,temp);
2775	}
2776	}
2777	if (output_struct_penalties == 1 \|\| output_struct_penalties == 2) {
2778	if (struct_penalties2 != NONE && use_ss2 == TRUE) {
2779	for(i=pos;i<=ptr;++i) {
2780	val=gap_penalty_mask2[i+fres-2];
2781	if (val == gap_pos1)
2782	temp[i-pos]='-';
2783	else
2784	temp[i-pos]=val;
2785	}
2786	temp[ptr-pos+1]=EOS;
2787	fprintf(clusout,"!GM_%-*s %s\n",max_names,ss_name2,temp);
2788	}
2789	}
2790
2791	for(ii=fseq;ii<=lseq;++ii) {
2792	i=output_index[ii];
2793	print_seq_no[i] = 0;
2794	for(j=pos;j<=ptr;++j) {
2795	if (j+fres-1<=seqlen_array[i])
2796	val = seq_array[i][j+fres-1];
2797	else val = -3;
2798	if((val == -3) \|\| (val == 253)) break;
2799	else if((val < 0) \|\| (val > max_aa))
2800	seq1[j]='-';
2801	else {
2802	seq1[j]=amino_acid_codes[val];
2803	seq_no[i]++;
2804	print_seq_no[i]=1;
2805	}
2806	}
2807	for(;j<=ptr;++j) seq1[j]='-';
2808	strncpy(temp,&seq1[pos],ptr-pos+1);
2809	temp[ptr-pos+1]=EOS;
2810	fprintf(clusout,"%-*s %s",max_names+5,names[i],temp);
2811	if (cl_seq_numbers && print_seq_no[i])
2812	fprintf(clusout," %d",seq_no[i]);
2813	fprintf(clusout,"\n");
2814	}
2815
2816	for(i=pos;i<=ptr;++i) {
2817	seq1[i]=' ';
2818	ident=0;
2819	for(j=1;res_cat1[j-1]!=NULL;j++) catident1[j-1] = 0;
2820	for(j=1;res_cat2[j-1]!=NULL;j++) catident2[j-1] = 0;
2821	for(j=fseq;j<=lseq;++j) {
2822	if((seq_array[fseq][i] >=0) &&
2823	(seq_array[fseq][i] <= max_aa)) {
2824	if(seq_array[fseq][i] == seq_array[j][i])
2825	++ident;
2826	for(k=1;res_cat1[k-1]!=NULL;k++) {
2827	for(l=0;(c=res_cat1[k-1][l]);l++) {
2828	if (amino_acid_codes[seq_array[j][i]]==c)
2829	{
2830	catident1[k-1]++;
2831	break;
2832	}
2833	}
2834	}
2835	for(k=1;res_cat2[k-1]!=NULL;k++) {
2836	for(l=0;(c=res_cat2[k-1][l]);l++) {
2837	if (amino_acid_codes[seq_array[j][i]]==c)
2838	{
2839	catident2[k-1]++;
2840	break;
2841	}
2842	}
2843	}
2844	}
2845	}
2846	if(ident==lseq-fseq+1)
2847	seq1[i]='*';
2848	else if (!dnaflag) {
2849	for(k=1;res_cat1[k-1]!=NULL;k++) {
2850	if (catident1[k-1]==lseq-fseq+1) {
2851	seq1[i]=':';
2852	break;
2853	}
2854	}
2855	if(seq1[i]==' ')
2856	for(k=1;res_cat2[k-1]!=NULL;k++) {
2857	if (catident2[k-1]==lseq-fseq+1) {
2858	seq1[i]='.';
2859	break;
2860	}
2861	}
2862	}
2863	}
2864	strncpy(temp,&seq1[pos],ptr-pos+1);
2865	temp[ptr-pos+1]=EOS;
2866	for(k=0;k<max_names+6;k++) fprintf(clusout," ");
2867	fprintf(clusout,"%s\n",temp);
2868	}
2869
2870
2871	seq1=ckfree((void *)seq1);
2872	if (struct_penalties1 == SECST && use_ss1 == TRUE) ckfree(ss_mask1);
2873	if (struct_penalties2 == SECST && use_ss2 == TRUE) ckfree(ss_mask2);
2874	/* DES ckfree(output_index); */
2875
2876	}
2877
2878
2879
2880
2881
2882
2883	void gcg_out(FILE *gcgout, sint fres, sint len, sint fseq, sint lseq)
2884	{
2885	/* static char aacids = "XCSTPAGNDEQHRKMILVFYW";/
2886	/* static char nbases = "XACGT"; /
2887	char *seq, residue;
2888	sint val;
2889	sint *all_checks;
2890	sint i,ii,chunks,block;
2891	sint j,k,pos1,pos2;
2892	long grand_checksum;
2893
2894	seq = (char )ckalloc((max_aln_length+1) sizeof(char));
2895	all_checks = (sint )ckalloc((lseq+1) sizeof(sint));
2896
2897	for(i=fseq; i<=lseq; i++) {
2898	for(j=fres; j<=fres+len-1; j++) {
2899	val = seq_array[i][j];
2900	if((val == -3) \|\| (val == 253)) break;
2901	else if((val < 0) \|\| (val > max_aa))
2902	residue = '.';
2903	else {
2904	residue = amino_acid_codes[val];
2905	}
2906	seq[j-fres+1] = residue;
2907	}
2908	all_checks[i] = SeqGCGCheckSum(seq+1, (int)len);
2909	}
2910
2911	grand_checksum = 0;
2912	for(i=1; i<=nseqs; i++) grand_checksum += all_checks[output_index[i]];
2913	grand_checksum = grand_checksum % 10000;
2914	fprintf(gcgout,"PileUp\n\n");
2915	fprintf(gcgout,"\n\n MSF:%5d Type: ",(pint)len);
2916	if(dnaflag)
2917	fprintf(gcgout,"N");
2918	else
2919	fprintf(gcgout,"P");
2920	fprintf(gcgout," Check:%6ld .. \n\n", (long)grand_checksum);
2921	for(ii=fseq; ii<=lseq; ii++) {
2922	i = output_index[ii];
2923	/* for(j=0; j<max_names; j++)
2924	if(names[i][j] == ' ') names[i][j] = '_'; */
2925	fprintf(gcgout,
2926	" Name: %s oo Len:%5d Check:%6ld Weight: 1.00\n",
2927	names[i],(pint)len,(long)all_checks[i]);
2928	}
2929	fprintf(gcgout,"\n//\n");
2930
2931	chunks = len/GCG_LINELENGTH;
2932	if(len % GCG_LINELENGTH != 0) ++chunks;
2933
2934	for(block=1; block<=chunks; block++) {
2935	fprintf(gcgout,"\n\n");
2936	pos1 = ((block-1) * GCG_LINELENGTH) + 1;
2937	pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
2938	for(ii=fseq; ii<=lseq; ii++) {
2939	i = output_index[ii];
2940	fprintf(gcgout,"\n%-*s ",max_names+5,names[i]);
2941	for(j=pos1, k=1; j<=pos2; j++, k++) {
2942	/*
2943	JULIE -
2944	check for sint sequences - pad out with '.' characters to end of alignment
2945	*/
2946	if (j+fres-1<=seqlen_array[i])
2947	val = seq_array[i][j+fres-1];
2948	else val = -3;
2949	if((val == -3) \|\| (val == 253))
2950	residue = '.';
2951	else if((val < 0) \|\| (val > max_aa))
2952	residue = '.';
2953	else {
2954	residue = amino_acid_codes[val];
2955	}
2956	fprintf(gcgout,"%c",residue);
2957	if(j % 10 == 0) fprintf(gcgout," ");
2958	}
2959	}
2960	}
2961	/* DES ckfree(output_index); */
2962
2963	seq=ckfree((void *)seq);
2964	all_checks=ckfree((void *)all_checks);
2965
2966	fprintf(gcgout,"\n\n");
2967	}
2968
2969
2970
2971	void phylip_out(FILE *phyout, sint fres, sint len, sint fseq, sint lseq)
2972	{
2973	/* static char aacids = "XCSTPAGNDEQHRKMILVFYW";/
2974	/* static char nbases = "XACGT"; /
2975	char residue;
2976	sint val;
2977	sint i,ii,chunks,block;
2978	sint j,k,pos1,pos2;
2979	sint name_len;
2980	Boolean warn;
2981	char **snames;
2982
2983	snames=(char *)ckalloc((lseq-fseq+2)sizeof(char *));
2984	name_len=0;
2985	for(i=fseq; i<=lseq; i++) {
2986	snames[i]=(char )ckalloc((11)sizeof(char));
2987	ii=strlen(names[i]);
2988	strncpy(snames[i],names[i],10);
2989	if(name_len<ii) name_len=ii;
2990	}
2991	if(name_len>10) {
2992	warn=FALSE;
2993	for(i=fseq; i<=lseq; i++) {
2994	for(j=i+1;j<=lseq;j++) {
2995	if (strcmp(snames[i],snames[j]) == 0)
2996	warn=TRUE;
2997	}
2998	}
2999	if(warn)
3000	warning("Truncating sequence names to 10 characters for PHYLIP output.\n"
3001	"Names in the PHYLIP format file are NOT unambiguous.");
3002	else
3003	warning("Truncating sequence names to 10 characters for PHYLIP output.");
3004	}
3005
3006
3007	chunks = len/GCG_LINELENGTH;
3008	if(len % GCG_LINELENGTH != 0) ++chunks;
3009
3010	fprintf(phyout,"%6d %6d",(pint)nseqs,(pint)len);
3011
3012	for(block=1; block<=chunks; block++) {
3013	pos1 = ((block-1) * GCG_LINELENGTH)+1;
3014	pos2 = (len<pos1+GCG_LINELENGTH-1)? len : pos1+GCG_LINELENGTH-1;
3015	for(ii=fseq; ii<=lseq; ii++) {
3016	i = output_index[ii];
3017	if(block == 1) {
3018	fprintf(phyout,"\n%-10s ",snames[i]);
3019	}
3020	else
3021	fprintf(phyout,"\n ");
3022	for(j=pos1, k=1; j<=pos2; j++, k++) {
3023	if (j+fres-1<=seqlen_array[i])
3024	val = seq_array[i][j+fres-1];
3025	else val = -3;
3026	if((val == -3) \|\| (val == 253))
3027	break;
3028	else if((val < 0) \|\| (val > max_aa))
3029	residue = '-';
3030	else {
3031	residue = amino_acid_codes[val];
3032	}
3033	fprintf(phyout,"%c",residue);
3034	if(j % 10 == 0) fprintf(phyout," ");
3035	}
3036	}
3037	fprintf(phyout,"\n");
3038	}
3039	/* DES ckfree(output_index); */
3040
3041	for(i=fseq;i<=lseq;i++)
3042	ckfree(snames[i]);
3043	ckfree(snames);
3044	}
3045
3046
3047
3048
3049
3050	void nbrf_out(FILE *nbout, sint fres, sint len, sint fseq, sint lseq)
3051	{
3052	/* static char aacids = "XCSTPAGNDEQHRKMILVFYW";/
3053	/* static char nbases = "XACGT"; /
3054	char *seq, residue;
3055	sint val;
3056	sint i,ii;
3057	sint j,slen;
3058	sint line_length;
3059
3060	seq = (char )ckalloc((max_aln_length+1) sizeof(char));
3061
3062	/* decide the line length for this alignment - maximum is LINELENGTH */
3063	line_length=PAGEWIDTH-max_names;
3064	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
3065	if (line_length > LINELENGTH) line_length=LINELENGTH;
3066
3067	for(ii=fseq; ii<=lseq; ii++) {
3068	i = output_index[ii];
3069	fprintf(nbout, dnaflag ? ">DL;" : ">P1;");
3070	fprintf(nbout, "%s\n%s\n", names[i], titles[i]);
3071	slen = 0;
3072	for(j=fres; j<fres+len; j++) {
3073	val = seq_array[i][j];
3074	if((val == -3) \|\| (val == 253))
3075	break;
3076	else if((val < 0) \|\| (val > max_aa))
3077	residue = '-';
3078	else {
3079	residue = amino_acid_codes[val];
3080	}
3081	seq[j-fres] = residue;
3082	slen++;
3083	}
3084	for(j=1; j<=slen; j++) {
3085	fprintf(nbout,"%c",seq[j-1]);
3086	if((j % line_length == 0) \|\| (j == slen))
3087	fprintf(nbout,"\n");
3088	}
3089	fprintf(nbout,"*\n");
3090	}
3091	/* DES ckfree(output_index); */
3092
3093	seq=ckfree((void *)seq);
3094	}
3095
3096
3097	void gde_out(FILE *gdeout, sint fres, sint len, sint fseq, sint lseq)
3098	{
3099	/* static char aacids = "XCSTPAGNDEQHRKMILVFYW";/
3100	/* static char nbases = "XACGT"; /
3101	char *seq, residue;
3102	sint val;
3103	char ss_mask1, ss_mask2;
3104	sint i,ii;
3105	sint j,slen;
3106	sint line_length;
3107
3108	seq = (char )ckalloc((max_aln_length+1) sizeof(char));
3109
3110	/* decide the line length for this alignment - maximum is LINELENGTH */
3111	line_length=PAGEWIDTH-max_names;
3112	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
3113	if (line_length > LINELENGTH) line_length=LINELENGTH;
3114
3115	if (struct_penalties1 == SECST && use_ss1 == TRUE) {
3116	ss_mask1 = (char )ckalloc((seqlen_array[1]+10) sizeof(char));
3117	for (i=0;i<seqlen_array[1];i++)
3118	ss_mask1[i] = sec_struct_mask1[i];
3119	print_sec_struct_mask(seqlen_array[1],sec_struct_mask1,ss_mask1);
3120	}
3121	if (struct_penalties2 == SECST && use_ss2 == TRUE) {
3122	ss_mask2 = (char )ckalloc((seqlen_array[profile1_nseqs+1]+10)
3123	sizeof(char));
3124	for (i=0;i<seqlen_array[profile1_nseqs+1];i++)
3125	ss_mask2[i] = sec_struct_mask2[i];
3126	print_sec_struct_mask(seqlen_array[profile1_nseqs+1],sec_struct_mask2,ss_mask2);
3127
3128
3129	}
3130
3131
3132	for(ii=fseq; ii<=lseq; ii++) {
3133	i = output_index[ii];
3134	fprintf(gdeout, dnaflag ? "#" : "%%");
3135	fprintf(gdeout, "%s\n", names[i]);
3136	slen = 0;
3137	for(j=fres; j<fres+len; j++) {
3138	val = seq_array[i][j];
3139	if((val == -3) \|\| (val == 253))
3140	break;
3141	else if((val < 0) \|\| (val > max_aa))
3142	residue = '-';
3143	else {
3144	residue = amino_acid_codes[val];
3145	}
3146	if (lowercase)
3147	seq[j-fres] = (char)tolower((int)residue);
3148	else
3149	seq[j-fres] = residue;
3150	slen++;
3151	}
3152	for(j=1; j<=slen; j++) {
3153	fprintf(gdeout,"%c",seq[j-1]);
3154	if((j % line_length == 0) \|\| (j == slen))
3155	fprintf(gdeout,"\n");
3156	}
3157	}
3158	/* DES ckfree(output_index); */
3159
3160	if (output_struct_penalties == 0 \|\| output_struct_penalties == 2) {
3161	if (struct_penalties1 == SECST && use_ss1 == TRUE) {
3162	fprintf(gdeout,"\"SS_%-*s\n",max_names,ss_name1);
3163	for(i=fres; i<fres+len; i++) {
3164	val=ss_mask1[i-1];
3165	if (val == gap_pos1)
3166	seq[i-fres]='-';
3167	else
3168	seq[i-fres]=val;
3169	}
3170	seq[i-fres]=EOS;
3171	for(i=1; i<=len; i++) {
3172	fprintf(gdeout,"%c",seq[i-1]);
3173	if((i % line_length == 0) \|\| (i == len))
3174	fprintf(gdeout,"\n");
3175	}
3176	}
3177
3178	if (struct_penalties2 == SECST && use_ss2 == TRUE) {
3179	fprintf(gdeout,"\"SS_%-*s\n",max_names,ss_name2);
3180	for(i=fres; i<fres+len; i++) {
3181	val=ss_mask2[i-1];
3182	if (val == gap_pos1)
3183	seq[i-fres]='-';
3184	else
3185	seq[i-fres]=val;
3186	}
3187	seq[i]=EOS;
3188	for(i=1; i<=len; i++) {
3189	fprintf(gdeout,"%c",seq[i-1]);
3190	if((i % line_length == 0) \|\| (i == len))
3191	fprintf(gdeout,"\n");
3192	}
3193	}
3194	}
3195	if (output_struct_penalties == 1 \|\| output_struct_penalties == 2) {
3196	if (struct_penalties1 != NONE && use_ss1 == TRUE) {
3197	fprintf(gdeout,"\"GM_%-*s\n",max_names,ss_name1);
3198	for(i=fres; i<fres+len; i++) {
3199	val=gap_penalty_mask1[i-1];
3200	if (val == gap_pos1)
3201	seq[i-fres]='-';
3202	else
3203	seq[i-fres]=val;
3204	}
3205	seq[i]=EOS;
3206	for(i=1; i<=len; i++) {
3207	fprintf(gdeout,"%c",seq[i-1]);
3208	if((i % line_length == 0) \|\| (i == len))
3209	fprintf(gdeout,"\n");
3210	}
3211	}
3212
3213	if (struct_penalties2 != NONE && use_ss2 == TRUE) {
3214	fprintf(gdeout,"\"GM_%-*s\n",max_names,ss_name2);
3215	for(i=fres; i<fres+len; i++) {
3216	val=gap_penalty_mask2[i-1];
3217	if (val == gap_pos1)
3218	seq[i-fres]='-';
3219	else
3220	seq[i-fres]=val;
3221	}
3222	seq[i]=EOS;
3223	for(i=1; i<=len; i++) {
3224	fprintf(gdeout,"%c",seq[i-1]);
3225	if((i % line_length == 0) \|\| (i == len))
3226	fprintf(gdeout,"\n");
3227	}
3228	}
3229	}
3230
3231	if (struct_penalties1 == SECST && use_ss1 == TRUE) ckfree(ss_mask1);
3232	if (struct_penalties2 == SECST && use_ss2 == TRUE) ckfree(ss_mask2);
3233
3234	seq=ckfree((void *)seq);
3235	}
3236
3237
3238	Boolean open_alignment_output(char *path)
3239	{
3240
3241	if(!output_clustal && !output_nbrf && !output_gcg &&
3242	!output_phylip && !output_gde) {
3243	error("You must select an alignment output format");
3244	return FALSE;
3245	}
3246
3247	if(output_clustal)
3248	if (outfile_name[0]!=EOS) {
3249	strcpy(clustal_outname,outfile_name);
3250	if((clustal_outfile = open_explicit_file(
3251	clustal_outname))==NULL) return FALSE;
3252	}
3253	else {
3254	/* DES DEBUG
3255	fprintf(stdout,"\n\n path = %s\n clustal_outname = %s\n\n",
3256	path,clustal_outname);
3257	*/
3258	if((clustal_outfile = open_output_file(
3259	"\nEnter a name for the CLUSTAL output file ",path,
3260	clustal_outname,"aln"))==NULL) return FALSE;
3261	/* DES DEBUG
3262	fprintf(stdout,"\n\n path = %s\n clustal_outname = %s\n\n",
3263	path,clustal_outname);
3264	*/
3265	}
3266	if(output_nbrf)
3267	if (outfile_name[0]!=EOS) {
3268	strcpy(nbrf_outname,outfile_name);
3269	if((nbrf_outfile = open_explicit_file(
3270	nbrf_outname))==NULL) return FALSE;
3271	}
3272	else
3273	if((nbrf_outfile = open_output_file(
3274	"\nEnter a name for the NBRF/PIR output file",path,
3275	nbrf_outname,"pir"))==NULL) return FALSE;
3276	if(output_gcg)
3277	if (outfile_name[0]!=EOS) {
3278	strcpy(gcg_outname,outfile_name);
3279	if((gcg_outfile = open_explicit_file(
3280	gcg_outname))==NULL) return FALSE;
3281	}
3282	else
3283	if((gcg_outfile = open_output_file(
3284	"\nEnter a name for the GCG output file ",path,
3285	gcg_outname,"msf"))==NULL) return FALSE;
3286	if(output_phylip)
3287	if (outfile_name[0]!=EOS) {
3288	strcpy(phylip_outname,outfile_name);
3289	if((phylip_outfile = open_explicit_file(
3290	phylip_outname))==NULL) return FALSE;
3291	}
3292	else
3293	if((phylip_outfile = open_output_file(
3294	"\nEnter a name for the PHYLIP output file ",path,
3295	phylip_outname,"phy"))==NULL) return FALSE;
3296	if(output_gde)
3297	if (outfile_name[0]!=EOS) {
3298	strcpy(gde_outname,outfile_name);
3299	if((gde_outfile = open_explicit_file(
3300	gde_outname))==NULL) return FALSE;
3301	}
3302	else
3303	if((gde_outfile = open_output_file(
3304	"\nEnter a name for the GDE output file ",path,
3305	gde_outname,"gde"))==NULL) return FALSE;
3306	return TRUE;
3307	}
3308
3309
3310
3311
3312
3313	void create_alignment_output(sint fseq, sint lseq)
3314	{
3315	sint i,length;
3316
3317	length=0;
3318	for (i=fseq;i<=lseq;i++)
3319	if (length < seqlen_array[i])
3320	length = seqlen_array[i];
3321	if (usemenu) info("Consensus length = %d",(pint)length);
3322
3323	if(output_clustal) {
3324	clustal_out(clustal_outfile, 1, length, fseq, lseq);
3325	fclose(clustal_outfile);
3326	info("CLUSTAL-Alignment file created [%s]",clustal_outname);
3327	}
3328	if(output_nbrf) {
3329	nbrf_out(nbrf_outfile, 1, length, fseq, lseq);
3330	fclose(nbrf_outfile);
3331	info("NBRF/PIR-Alignment file created [%s]",nbrf_outname);
3332	}
3333	if(output_gcg) {
3334	gcg_out(gcg_outfile, 1, length, fseq, lseq);
3335	fclose(gcg_outfile);
3336	info("GCG-Alignment file created [%s]",gcg_outname);
3337	}
3338	if(output_phylip) {
3339	phylip_out(phylip_outfile, 1, length, fseq, lseq);
3340	fclose(phylip_outfile);
3341	info("PHYLIP-Alignment file created [%s]",phylip_outname);
3342	}
3343	if(output_gde) {
3344	gde_out(gde_outfile, 1, length, fseq, lseq);
3345	fclose(gde_outfile);
3346	info("GDE-Alignment file created [%s]",gde_outname);
3347	}
3348	}
3349
3350
3351
3352
3353
3354
3355	static void reset_align(void) /* remove gaps from older alignments (code = gap_pos1) */
3356	{ /* EXCEPT for gaps that were INPUT with the seqs.*/
3357	register sint sl; /* which have code = gap_pos2 */
3358	sint i,j;
3359
3360	for(i=1;i<=nseqs;++i) {
3361	sl=0;
3362	for(j=1;j<=seqlen_array[i];++j) {
3363	if(seq_array[i][j] == gap_pos1) continue;
3364	++sl;
3365	seq_array[i][sl]=seq_array[i][j];
3366	}
3367	seqlen_array[i]=sl;
3368	}
3369
3370
3371	}
3372
3373
3374
3375	static void reset_prf1(void) /* remove gaps from older alignments (code = gap_pos1) */
3376	{ /* EXCEPT for gaps that were INPUT with the seqs.*/
3377	register sint sl; /* which have code = gap_pos2 */
3378	sint i,j;
3379
3380	if (struct_penalties1 != NONE) {
3381	sl=0;
3382	for (j=0;j<seqlen_array[1];++j) {
3383	if (gap_penalty_mask1[j] == gap_pos1) continue;
3384	gap_penalty_mask1[sl]=gap_penalty_mask1[j];
3385	++sl;
3386	}
3387	}
3388
3389	if (struct_penalties1 == SECST) {
3390	sl=0;
3391	for (j=0;j<seqlen_array[1];++j) {
3392	if (sec_struct_mask1[j] == gap_pos1) continue;
3393	sec_struct_mask1[sl]=sec_struct_mask1[j];
3394	++sl;
3395	}
3396	}
3397
3398	for(i=1;i<=profile1_nseqs;++i) {
3399	sl=0;
3400	for(j=1;j<=seqlen_array[i];++j) {
3401	if(seq_array[i][j] == gap_pos1) continue;
3402	++sl;
3403	seq_array[i][sl]=seq_array[i][j];
3404	}
3405	seqlen_array[i]=sl;
3406	}
3407
3408
3409	}
3410
3411
3412
3413	static void reset_prf2(void) /* remove gaps from older alignments (code = gap_pos1) */
3414	{ /* EXCEPT for gaps that were INPUT with the seqs.*/
3415	register sint sl; /* which have code = gap_pos2 */
3416	sint i,j;
3417
3418	if (struct_penalties2 != NONE) {
3419	sl=0;
3420	for (j=0;j<seqlen_array[profile1_nseqs+1];++j) {
3421	if (gap_penalty_mask2[j] == gap_pos1) continue;
3422	gap_penalty_mask2[sl]=gap_penalty_mask2[j];
3423	++sl;
3424	}
3425	}
3426
3427	if (struct_penalties2 == SECST) {
3428	sl=0;
3429	for (j=0;j<seqlen_array[profile1_nseqs+1];++j) {
3430	if (sec_struct_mask2[j] == gap_pos1) continue;
3431	sec_struct_mask2[sl]=sec_struct_mask2[j];
3432	++sl;
3433	}
3434	}
3435
3436	for(i=profile1_nseqs+1;i<=nseqs;++i) {
3437	sl=0;
3438	for(j=1;j<=seqlen_array[i];++j) {
3439	if(seq_array[i][j] == gap_pos1) continue;
3440	++sl;
3441	seq_array[i][sl]=seq_array[i][j];
3442	}
3443	seqlen_array[i]=sl;
3444	}
3445
3446
3447	}
3448
3449
3450
3451	void fix_gaps(void) /* fix gaps introduced in older alignments (code = gap_pos1) */
3452	{
3453	sint i,j;
3454
3455	if (struct_penalties1 != NONE) {
3456	for (j=0;j<seqlen_array[1];++j) {
3457	if (gap_penalty_mask1[j] == gap_pos1)
3458	gap_penalty_mask1[j]=gap_pos2;
3459	}
3460	}
3461
3462	if (struct_penalties1 == SECST) {
3463	for (j=0;j<seqlen_array[1];++j) {
3464	if (sec_struct_mask1[j] == gap_pos1)
3465	sec_struct_mask1[j]=gap_pos2;
3466	}
3467	}
3468
3469	for(i=1;i<=nseqs;++i) {
3470	for(j=1;j<=seqlen_array[i];++j) {
3471	if(seq_array[i][j] == gap_pos1)
3472	seq_array[i][j]=gap_pos2;
3473	}
3474	}
3475	}
3476
3477	static sint find_match(char probe, char list[], sint n)
3478	{
3479	sint i,j,len;
3480	sint count,match=0;
3481
3482	len = (sint)strlen(probe);
3483	for (i=0;i<len;i++) {
3484	count = 0;
3485	for (j=0;j<n;j++) {
3486	if (probe[i] == list[j][i]) {
3487	match = j;
3488	count++;
3489	}
3490	}
3491	if (count == 0) return((sint)-1);
3492	if (count == 1) return(match);
3493	}
3494	return((sint)-1);
3495	}
3496
3497	static void create_parameter_output(void)
3498	{
3499	char parname[FILENAMELEN+1], temp[FILENAMELEN+1];
3500	char path[FILENAMELEN+1];
3501	FILE *parout;
3502
3503	get_path(seqname,path);
3504	strcpy(parname,path);
3505	strcat(parname,"par");
3506
3507	if(usemenu) {
3508	fprintf(stdout,"\nEnter a name for the parameter output file [%s]: ",
3509	parname);
3510	gets(temp);
3511	if(*temp != EOS)
3512	strcpy(parname,temp);
3513	}
3514
3515	/* create a file with execute permissions first */
3516	remove(parname);
3517	/*
3518	fd = creat(parname, 0777);
3519	close(fd);
3520	*/
3521
3522	if((parout = open_explicit_file(parname))==NULL) return;
3523
3524	fprintf(parout,"clustalw \\\n");
3525	if (!empty && profile1_empty) fprintf(parout,"/infile=%s \\\n",seqname);
3526	if (!profile1_empty) fprintf(parout,"/profile1=%s\\\n",profile1_name);
3527	if (!profile2_empty) fprintf(parout,"/profile2=%s \\\n",profile2_name);
3528	if (dnaflag == TRUE) fprintf(parout,"/type=dna \\\n");
3529	else fprintf(parout,"/type=protein \\\n");
3530
3531	if (quick_pairalign) {
3532	fprintf(parout,"/quicktree \\\n");
3533	fprintf(parout,"/ktuple=%d \\\n",(pint)ktup);
3534	fprintf(parout,"/window=%d \\\n",(pint)window);
3535	fprintf(parout,"/pairgap=%d \\\n",(pint)wind_gap);
3536	fprintf(parout,"/topdiags=%d \\\n",(pint)signif);
3537	if (percent) fprintf(parout,"/score=percent \\\n");
3538	else fprintf(parout,"/score=absolute \\\n");
3539	}
3540	else {
3541	if (!dnaflag) {
3542	fprintf(parout,"/pwmatrix=%s \\\n",pw_mtrxname);
3543	fprintf(parout,"/pwgapopen=%.2f \\\n",prot_pw_go_penalty);
3544	fprintf(parout,"/pwgapext=%.2f \\\n",prot_pw_ge_penalty);
3545	}
3546	else {
3547	fprintf(parout,"/pwgapopen=%.2f \\\n",pw_go_penalty);
3548	fprintf(parout,"/pwgapext=%.2f \\\n",pw_ge_penalty);
3549	}
3550	}
3551
3552	if (!dnaflag) {
3553	fprintf(parout,"/matrix=%s \\\n",mtrxname);
3554	fprintf(parout,"/gapopen=%.2f \\\n",prot_gap_open);
3555	fprintf(parout,"/gapext=%.2f \\\n",prot_gap_extend);
3556	}
3557	else {
3558	fprintf(parout,"/gapopen=%.2f \\\n",dna_gap_open);
3559	fprintf(parout,"/gapext=%.2f \\\n",dna_gap_extend);
3560	}
3561
3562	fprintf(parout,"/maxdiv=%d \\\n",(pint)divergence_cutoff);
3563	if (!use_endgaps) fprintf(parout,"/endgaps \\\n");
3564
3565	if (!dnaflag) {
3566	if (neg_matrix) fprintf(parout,"/negative \\\n");
3567	if (no_pref_penalties) fprintf(parout,"/nopgap \\\n");
3568	if (no_hyd_penalties) fprintf(parout,"/nohgap \\\n");
3569	fprintf(parout,"/hgapresidues=%s \\\n",hyd_residues);
3570	fprintf(parout,"/gapdist=%d \\\n",(pint)gap_dist);
3571	}
3572	else {
3573	fprintf(parout,"/transweight=%.2f \\\n",transition_weight);
3574	}
3575
3576	if (output_gcg) fprintf(parout,"/output=gcg \\\n");
3577	else if (output_gde) fprintf(parout,"/output=gde \\\n");
3578	else if (output_nbrf) fprintf(parout,"/output=pir \\\n");
3579	else if (output_phylip) fprintf(parout,"/output=phylip \\\n");
3580	if (outfile_name[0]!=EOS) fprintf(parout,"/outfile=%s \\\n",outfile_name);
3581	if (output_order==ALIGNED) fprintf(parout,"/outorder=aligned \\\n");
3582	else fprintf(parout,"/outorder=input \\\n");
3583	if (output_gde)
3584	if (lowercase) fprintf(parout,"/case=lower \\\n");
3585	else fprintf(parout,"/case=upper \\\n");
3586
3587
3588	fprintf(parout,"/interactive\n");
3589
3590	fclose(parout);
3591	/*
3592	if (kimura) fprintf(parout,"/kimura \\\n");
3593	if (tossgaps) fprintf(parout,"/tossgaps \\\n");
3594	fprintf(parout,"/seed=%d \\\n",(pint)boot_ran_seed);
3595	fprintf(parout,"/bootstrap=%d \\\n",(pint)boot_ntrials);
3596	*/
3597	}
3598

Note: See TracBrowser for help on using the repository browser.

Download in other formats: