Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

interface.c

Visit:

Last change on this file was 19480, checked in by westram, 2 years ago
reintegrates 'gde' into 'trunk' adds: log:branches/gde@19460:19479
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.6 KB

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

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

Context Navigation

source: trunk/GDE/CLUSTALW/interface.c

Download in other formats: