Context Navigation

source: tags/arb-6.0/GDE/CLUSTALW/interface.c

Visit:

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

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

Download in other formats: