interface.c

经典生物信息学多序列比对工具clustalw

}
					}
				info("No. of seqs in profile=%d",(pint)local_nseqs);
                nseqs = profile1_nseqs + local_nseqs;
                info("Total no. of seqs     =%d",(pint)nseqs);
				profile2_empty=FALSE;
				empty = FALSE;
	   }

      }
	if (sec_struct_mask != NULL) sec_struct_mask=ckfree(sec_struct_mask);
	if (gap_penalty_mask != NULL) gap_penalty_mask=ckfree(gap_penalty_mask);
	if (ss_name != NULL) ss_name=ckfree(ss_name);

	if(local_nseqs<=0) return local_nseqs;
	
	info("Sequences assumed to be %s",
		dnaflag?"DNA":"PROTEIN");
	if (usemenu) fprintf(stdout,"\n\n");
        for(i=profile2_empty?1:profile1_nseqs+1; i<=nseqs; i++) {
                info("Sequence %d: %-*s   %6.d %s",
                   (pint)i,max_names,names[i],(pint)seqlen_array[i],dnaflag?"bp":"aa");
        }	
	if(dnaflag) {
		gap_open   = dna_gap_open;
		gap_extend = dna_gap_extend;
	}
	else {
		gap_open   = prot_gap_open;
		gap_extend = prot_gap_extend;
	}

	return nseqs;
}



static void calc_gap_penalty_mask(int prf_length, char *mask, char *gap_mask)
{
	int i,j;
	char *struct_mask;

	struct_mask = (char *)ckalloc((prf_length+1) * sizeof(char));
/*
    calculate the gap penalty mask from the secondary structures
*/
	i=0;
	while (i<prf_length) {
		if (tolower(mask[i]) == 'a' || mask[i] == '$') {
			for (j = -helix_end_plus; j<0; j++) {
				if ((i+j>=0) && (tolower(struct_mask[i+j]) != 'a')
				             && (tolower(struct_mask[i+j]) != 'b'))
					struct_mask[i+j] = 'a';
			}
			for (j = 0; j<helix_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'a'
				                    && mask[i+j] != '$')) break;
				struct_mask[i+j] = 'a';
			}
			i += j;
			while (tolower(mask[i]) == 'a'
				                    || mask[i] == '$') {
				if (i>=prf_length) break;
				if (mask[i] == '$') {
					struct_mask[i] = 'A';
					i++;
					break;
				}
				else struct_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<helix_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'a'
				                    || mask[i-j-1] == '$'))
					struct_mask[i-j-1] = 'a';
			}
			for (j = 0; j<helix_end_plus; j++) {
				if (i+j>=prf_length) break;
				struct_mask[i+j] = 'a';
			}
		}
	 	else if (tolower(mask[i]) == 'b' || mask[i] == '%') {
			for (j = -strand_end_plus; j<0; j++) {
				if ((i+j>=0) && (tolower(struct_mask[i+j]) != 'a')
				             && (tolower(struct_mask[i+j]) != 'b'))
					struct_mask[i+j] = 'b';
			}
			for (j = 0; j<strand_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'b'
				                    && mask[i+j] != '%')) break;
				struct_mask[i+j] = 'b';
			}
			i += j;
			while (tolower(mask[i]) == 'b'
				                    || mask[i] == '%') {
				if (i>=prf_length) break;
				if (mask[i] == '%') {
					struct_mask[i] = 'B';
					i++;
					break;
				}
				else struct_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<strand_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'b'
				                    || mask[i-j-1] == '%'))
				struct_mask[i-j-1] = 'b';
			}
			for (j = 0; j<strand_end_plus; j++) {
				if (i+j>=prf_length) break;
 				struct_mask[i+j] = 'b';
			}
		}
	else i++;
	}

	for(i=0;i<prf_length;i++) {
		switch (struct_mask[i]) {
			case 'A':
				gap_mask[i] = helix_penalty+'0';
				break;
			case 'a':
				gap_mask[i] = helix_end_penalty+'0';
				break;
			case 'B':
				gap_mask[i] = strand_penalty+'0';
				break;
			case 'b':
				gap_mask[i] = strand_end_penalty+'0';
				break;
			default:
				gap_mask[i] = loop_penalty+'0';
				break;
		}
	}

	struct_mask=ckfree(struct_mask);
	
}

void print_sec_struct_mask(int prf_length, char *mask, char *struct_mask)
{
	int i,j;

/*
    calculate the gap penalty mask from the secondary structures
*/
	i=0;
	while (i<prf_length) {
		if (tolower(mask[i]) == 'a' || mask[i] == '$') {
			for (j = 0; j<helix_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'a'
				                    && mask[i+j] != '$')) break;
				struct_mask[i+j] = 'a';
			}
			i += j;
			while (tolower(mask[i]) == 'a'
				                    || mask[i] == '$') {
				if (i>=prf_length) break;
				if (mask[i] == '$') {
					struct_mask[i] = 'A';
					i++;
					break;
				}
				else struct_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<helix_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'a'
				                    || mask[i-j-1] == '$'))
					struct_mask[i-j-1] = 'a';
			}
		}
	 	else if (tolower(mask[i]) == 'b' || mask[i] == '%') {
			for (j = 0; j<strand_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'b'
				                    && mask[i+j] != '%')) break;
				struct_mask[i+j] = 'b';
			}
			i += j;
			while (tolower(mask[i]) == 'b'
				                    || mask[i] == '%') {
				if (i>=prf_length) break;
				if (mask[i] == '%') {
					struct_mask[i] = 'B';
					i++;
					break;
				}
				else struct_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<strand_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'b'
				                    || mask[i-j-1] == '%'))
				struct_mask[i-j-1] = 'b';
			}
		}
	else i++;
	}
}



FILE *  open_output_file(char *prompt,      char *path, 
				char *file_name,   char *file_extension)
 
{	static char temp[FILENAMELEN+1];
	static char local_prompt[MAXLINE];
	FILE * file_handle;

/*	if (*file_name == EOS) {
*/		strcpy(file_name,path);
		strcat(file_name,file_extension);
/*	}
*/
	if(strcmp(file_name,seqname)==0) {
		warning("Output file name is the same as input file.");
		if (usemenu) {
			strcpy(local_prompt,"\n\nEnter new name to avoid overwriting ");
			strcat(local_prompt," [%s]: ");          
			fprintf(stdout,local_prompt,file_name);
			gets(temp);
			if(*temp != EOS) strcpy(file_name,temp);
		}
	}
	else if (usemenu) {
		strcpy(local_prompt,prompt);
		strcat(local_prompt," [%s]: ");          
		fprintf(stdout,local_prompt,file_name);
		gets(temp);
		if(*temp != EOS) strcpy(file_name,temp);
	}

#ifdef VMS
	if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
#else
	if((file_handle=fopen(file_name,"w"))==NULL) {
#endif
		error("Cannot open output file [%s]",file_name);
		return NULL;
	}
	return file_handle;
}



FILE *  open_explicit_file(char *file_name)
{ 
	FILE * file_handle;

	if (*file_name == EOS) {
		error("Bad output file [%s]",file_name);
		return NULL;
	}
#ifdef VMS
	if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
#else
	if((file_handle=fopen(file_name,"w"))==NULL) {
#endif
		error("Cannot open output file [%s]",file_name);
		return NULL;
	}
	return file_handle;
}


/* Ramu void */

void align(char *phylip_name)
{ 
	char path[FILENAMELEN+1];
	FILE *tree;
	sint count;
	
	if(empty && usemenu) {
		error("No sequences in memory. Load sequences first.");
		return;
	}

	   struct_penalties1 = struct_penalties2 = NONE;
	   if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
	   if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
	   if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
	   if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
	   if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
	   if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);


        get_path(seqname,path);
/* DES DEBUG 
	fprintf(stdout,"\n\n Seqname = %s  \n Path = %s \n\n",seqname,path);
*/
	if(usemenu || !interactive) {
        	if(!open_alignment_output(path)) return;
	}

	if (nseqs >= 2) {

        	get_path(seqname,path);
        	if (phylip_name[0]!=EOS) {
                	if((tree = open_explicit_file(
                	phylip_name))==NULL) return;
        	}
        	else {
                 	if((tree = open_output_file(
                	"\nEnter name for new GUIDE TREE           file  ",path,
                	phylip_name,"dnd")) == NULL) return;
        	}
	}

	if (save_parameters) create_parameter_output();

	if(reset_alignments_new || reset_alignments_all) reset_align();

        info("Start of Pairwise alignments");
        info("Aligning...");
        if(dnaflag) {
                gap_open   = dna_gap_open;
                gap_extend = dna_gap_extend;
                pw_go_penalty  = dna_pw_go_penalty;
                pw_ge_penalty  = dna_pw_ge_penalty;
                ktup       = dna_ktup;
                window     = dna_window;
                signif     = dna_signif;
                wind_gap   = dna_wind_gap;

        }
        else {
                gap_open   = prot_gap_open;
                gap_extend = prot_gap_extend;
                pw_go_penalty  = prot_pw_go_penalty;
                pw_ge_penalty  = prot_pw_ge_penalty;
                ktup       = prot_ktup;
                window     = prot_window;
                signif     = prot_signif;
                wind_gap   = prot_wind_gap;

        }

        if (quick_pairalign)
           show_pair((sint)0,nseqs,(sint)0,nseqs);
        else
           pairalign((sint)0,nseqs,(sint)0,nseqs);

	if (nseqs >= 2) {

		guide_tree(tree,1,nseqs);
		info("Guide tree        file created:   [%s]",
                phylip_name);
	}

	
	count = malign((sint)0,phylip_name);
	
	if (count <= 0) return;

	if (usemenu) fprintf(stdout,"\n\n\n");
	
	create_alignment_output(1,nseqs);
        if (showaln && usemenu) show_aln();
	phylip_name[0]=EOS;
	return ;
}





void new_sequence_align(char *phylip_name)
{ 
	char path[FILENAMELEN+1];
	char tree_name[FILENAMELEN+1],temp[MAXLINE+1];
	Boolean use_tree;
	FILE *tree;
	sint i,j,count;
	float dscore;
	Boolean save_ss2;
	
	if(profile1_empty && usemenu) {
		error("No profile in memory. Input 1st profile first.");
		return;
	}

	if(profile2_empty && usemenu) {
		error("No sequences in memory. Input sequences first.");
		return;
	}

        get_path(profile2_name,path);

        if(usemenu || !interactive) {
        	if(!open_alignment_output(path)) return;
	}

	new_seq = profile1_nseqs+1;

/* check for secondary structure information for list of sequences */

	save_ss2 = use_ss2;
	if (struct_penalties2 != NONE && use_ss2 == TRUE && (nseqs - profile1_nseqs >
1)) {
		if (struct_penalties2 == SECST) 
			warning("Warning: ignoring secondary structure for a list of sequences");
		else if (struct_penalties2 == GMASK)
			warning("Warning: ignoring gap penalty mask for a list of sequences");
		use_ss2 = FALSE;
	}

	for (i=1;i<=new_seq;i++) {
     		for (j=i+1;j<=new_seq;j++) {
       			dscore = countid(i,j);
       			tmat[i][j] = ((double)100.0 - (double)dscore)/(double)100.0;
       			tmat[j][i] = tmat[i][j];
     		}
   	}

	tree_name[0] = EOS;
	use_tree = FALSE;
	if (nseqs >= 2) {
		if (check_tree && usemenu) {
			strcpy(tree_name,path);
			strcat(tree_name,"dnd");
#ifdef VMS
        	if((tree=fopen(tree_name,"r","rat=cr","rfm=var"))!=NULL) {
#else
        	if((tree=fopen(tree_name,"r"))!=NULL) {
#endif
		if (usemenu)
            	fprintf(stdout,"\nUse the existing GUIDE TREE file,  %s  (y/n) ? [y]: ",
                                           tree_name);
                gets(temp);
                if(*temp != 'n' && *temp != 'N') {
                    strcpy(phylip_name,tree_name);
                    use_tree = TRUE;
                }
                fclose(tree);
        	}
		}
		else if (!usemenu && use_tree_file) {
			use_tree = TRUE;
		}
	}
	
	if (save_parameters) create_parameter_output();

	if(reset_alignments_new || reset_alignments_all) {
/*
		reset_prf1();
*/
		reset_prf2();
	}
	else fix_gaps();

	if (struct_penalties1 == SECST)

		calc_gap_penalty_mask(seqlen_array[1],sec_struct_mask1,gap_penalty_mask1);

	if (struct_penalties2 == SECST)

calc_gap_penalty_mask(seqlen_array[profile1_nseqs+1],sec_struct_mask2,gap_penalty_mask2);


/* create the new tree file, if necessary */

	if (use_tree == FALSE) {

		if (nseqs >= 2) {
        		get_path(profile2_name,path);
        		if (phylip_name[0]!=EOS) {
                		if((tree = open_explicit_file(
                		phylip_name))==NULL) return;
        		}
        		else {
                 		if((tree = open_output_file(
                		"\nEnter name for new GUIDE TREE           file  ",path,
                		phylip_name,"dnd")) == NULL) return;
        		}
		}
        info("Start of Pairwise alignments");
        info("Aligning...");
        if(dnaflag) {
                gap_open   = dna_gap_open;
                gap_extend = dna_gap_extend;
                pw_go_