interface.c

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

			exit(1);
		}
		strcpy(seqname, param_arg[setprofile2]);
		profile_no = 2;
		profile_input();
		if(nseqs > profile1_nseqs) 
			do_something = do_profile = TRUE;
		else {
			if(nseqs<0) cl_error_code=2;
			else if(nseqs==0) cl_error_code=3;
			error("No sequences read from profile 2");
			exit(cl_error_code);
		}
		strcpy(profile2_name,seqname);
	}

/*************************************************************************/
/* Look for /tree or /bootstrap or /align or /usetree ******************/
/*************************************************************************/

	if (setbatch != -1)
		interactive=FALSE;

	if (setinteractive != -1)
		interactive=TRUE;

	if (interactive) {
		settree = -1;
		setbootstrap = -1;
		setalign = -1;
		setusetree = -1;
		setusetree1 = -1;
		setusetree2 = -1;
		setnewtree = -1;
		setconvert = -1;
	}

	if(settree != -1 )
		if(empty) {
			error("Cannot draw tree.  No input alignment file");
			exit(1);
		}
		else 
			do_tree = TRUE;

	if(setbootstrap != -1)
		if(empty) {
			error("Cannot bootstrap tree. No input alignment file");
			exit(1);
		}
		else {
			temp = 0;
			if(param_arg[setbootstrap] != NULL)
				 if (sscanf(param_arg[setbootstrap],"%d",&temp)!=1) {
                         fprintf(stdout,"Bad option for /bootstrap (must be integer)\n");
                         temp = 0;
                    };
			if(temp > 0)          boot_ntrials = temp;
			do_boot = TRUE;
		}

	if(setalign != -1)
		if(empty) {
			error("Cannot align sequences.  No input file");
			exit(1);
		}
		else 
			do_align = TRUE;

	if(setconvert != -1)
		if(empty) {
			error("Cannot convert sequences.  No input file");
			exit(1);
		}
		else 
			do_convert = TRUE;
 
	if(setusetree != -1)
		if(empty) {
			error("Cannot align sequences.  No input file");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setusetree]) == 0) {
				error("Cannot align sequences.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(phylip_tree_name, param_arg[setusetree]);
		        }
		        use_tree_file = TRUE;
		        do_align_only = TRUE;
		}

	if(setnewtree != -1)
		if(empty) {
			error("Cannot align sequences.  No input file");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setnewtree]) == 0) {
				error("Cannot align sequences.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(phylip_tree_name, param_arg[setnewtree]);
		        }
		    new_tree_file = TRUE;
			do_tree_only = TRUE;
		}
 
	if(setusetree1 != -1)
		if(profile1_empty) {
			error("Cannot align profiles.  No input file");
			exit(1);
		}
		else if(profile_type == SEQUENCE) {
			error("Invalid option /usetree1.");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setusetree1]) == 0) {
				error("Cannot align profiles.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(p1_tree_name, param_arg[setusetree1]);
		        }
		        use_tree1_file = TRUE;
		        do_align_only = TRUE;
		}

	if(setnewtree1 != -1)
		if(profile1_empty) {
			error("Cannot align profiles.  No input file");
			exit(1);
		}
		else if(profile_type == SEQUENCE) {
			error("Invalid option /newtree1.");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setnewtree1]) == 0) {
				error("Cannot align profiles.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(p1_tree_name, param_arg[setnewtree1]);
		        }
		    new_tree1_file = TRUE;
		}
 
	if(setusetree2 != -1)
		if(profile2_empty) {
			error("Cannot align profiles.  No input file");
			exit(1);
		}
		else if(profile_type == SEQUENCE) {
			error("Invalid option /usetree2.");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setusetree2]) == 0) {
				error("Cannot align profiles.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(p2_tree_name, param_arg[setusetree2]);
		        }
		        use_tree2_file = TRUE;
		        do_align_only = TRUE;
		}

	if(setnewtree2 != -1)
		if(profile2_empty) {
			error("Cannot align profiles.  No input file");
			exit(1);
		}
		else if(profile_type == SEQUENCE) {
			error("Invalid option /newtree2.");
			exit(1);
		}
		else  {
		        if(strlen(param_arg[setnewtree2]) == 0) {
				error("Cannot align profiles.  No tree file specified");
				exit(1);
		        }
                        else {
			        strcpy(p2_tree_name, param_arg[setnewtree2]);
		        }
		    new_tree2_file = TRUE;
		}
 

	if( (!do_tree) && (!do_boot) && (!empty) && (!do_profile) && (!do_align_only) && (!do_tree_only) && (!do_convert)) 
		do_align = TRUE;

/*** ? /quicktree  */
        if(setquicktree != -1)
		quick_pairalign = TRUE;

	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(interactive) {
		if (!xmenus) usemenu = TRUE;
		return;
	}


	if(!do_something) {
		error("No input file(s) specified");
		exit(1);
	}




/****************************************************************************/
/* Now do whatever has been requested ***************************************/
/****************************************************************************/


	if(do_profile) {
		if (profile_type == PROFILE) profile_align(p1_tree_name,p2_tree_name);
		else new_sequence_align(phylip_tree_name);
	}

	else if(do_align)
		align(phylip_tree_name);

        else if(do_convert) {
                get_path(seqname,path);
                if(!open_alignment_output(path)) exit(1);
                create_alignment_output(1,nseqs);
        }

        else if (do_align_only)
                get_tree(phylip_tree_name);

	else if(do_tree_only)
		make_tree(phylip_tree_name);

	else if(do_tree)
		phylogenetic_tree(phylip_tree_name,clustal_tree_name,dist_tree_name,nexus_tree_name,pim_name);

	else if(do_boot)
		bootstrap_tree(phylip_tree_name,clustal_tree_name,nexus_tree_name);

	fprintf(stdout,"\n");
	exit(0);

/*******whew!***now*go*home****/
}


Boolean user_mat(char *str, short *mat, short *xref)
{
        sint maxres;

        FILE *infile;

        if(usemenu)
                getstr("Enter name of the matrix file",lin2);
        else
                strcpy(lin2,str);

        if(*lin2 == EOS) return FALSE;

        if((infile=fopen(lin2,"r"))==NULL) {
                error("Cannot find matrix file [%s]",lin2);
                return FALSE;
        }

	strcpy(str, lin2);

	maxres = read_user_matrix(str, mat, xref);
        if (maxres <= 0) return FALSE;

	return TRUE;
}

Boolean user_mat_series(char *str, short *mat, short *xref)
{
        sint maxres;

        FILE *infile;

        if(usemenu)
                getstr("Enter name of the matrix file",lin2);
        else
                strcpy(lin2,str);

        if(*lin2 == EOS) return FALSE;

        if((infile=fopen(lin2,"r"))==NULL) {
                error("Cannot find matrix file [%s]",lin2);
                return FALSE;
        }

	strcpy(str, lin2);

	maxres = read_matrix_series(str, mat, xref);
        if (maxres <= 0) return FALSE;

	return TRUE;
}






sint seq_input(Boolean append)
{
        sint i;
	sint local_nseqs;

	if(usemenu) {
fprintf(stdout,"\n\nSequences should all be in 1 file.\n"); 
fprintf(stdout,"\n7 formats accepted: \n");
fprintf(stdout,
"NBRF/PIR, EMBL/SwissProt, Pearson (Fasta), GDE, Clustal, GCG/MSF, RSF.\n\n\n");
/*fprintf(stdout,
"\nGCG users should use TOPIR to convert their sequence files before use.\n\n\n");*/
	}

       if (append)
          local_nseqs = readseqs(nseqs+(sint)1);
       else
          local_nseqs = readseqs((sint)1);  /*  1 is the first seq to be read */
       if(local_nseqs < 0)               /* file could not be opened */
           { 
		return local_nseqs;
           }
       else if(local_nseqs == 0)         /* no sequences */
           {
	       error("No sequences in file!  Bad format?");
               return local_nseqs;
           }
       else 
           {
	   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);
	   
		if(append) nseqs+=local_nseqs;
		else nseqs=local_nseqs;
		info("Sequences assumed to be %s",
			dnaflag?"DNA":"PROTEIN");
		if (usemenu) {
			fprintf(stdout,"\n\n");
                	for(i=1; i<=nseqs; i++) {
/* DES                         fprintf(stdout,"%s: = ",names[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;
			}
			empty=FALSE;
	   }
	return local_nseqs;	
}







sint profile_input(void)   /* read a profile   */
{                                           /* profile_no is 1 or 2  */
        sint local_nseqs, i;
	
        if(profile_no == 2 && profile1_empty) 
           {
             error("You must read in profile number 1 first");
             return 0;
           }

    if(profile_no == 1)     /* for the 1st profile */
      {
       local_nseqs = readseqs((sint)1); /* (1) means 1st seq to be read = no. 1 */
       if(local_nseqs < 0)               /* file could not be opened */
           { 
		return local_nseqs;
           }
       else if(local_nseqs == 0)         /* no sequences  */
           {
	       error("No sequences in file!  Bad format?");
		return local_nseqs;
           }
       else if (local_nseqs > 0)
           { 				/* success; found some seqs. */
		struct_penalties1 = NONE;
		if (sec_struct_mask1 != NULL) sec_struct_mask1=ckfree(sec_struct_mask1);
		if (gap_penalty_mask1 != NULL) gap_penalty_mask1=ckfree(gap_penalty_mask1);
		if (ss_name1 != NULL) ss_name1=ckfree(ss_name1);
                if (struct_penalties != NONE) /* feature table / mask in alignment */
                	{
					struct_penalties1 = struct_penalties;
					if (struct_penalties == SECST) {
						sec_struct_mask1 = (char *)ckalloc((max_aln_length) * sizeof (char));
						for (i=0;i<max_aln_length;i++)
							sec_struct_mask1[i] = sec_struct_mask[i];
					}
					gap_penalty_mask1 = (char *)ckalloc((max_aln_length) * sizeof (char));
					for (i=0;i<max_aln_length;i++)
						gap_penalty_mask1[i] = gap_penalty_mask[i];
        				ss_name1 = (char *)ckalloc( (MAXNAMES+1) * sizeof (char));

					strcpy(ss_name1,ss_name);
if (debug>0) {
for (i=0;i<seqlen_array[1];i++)
	fprintf(stdout,"%c",gap_penalty_mask1[i]);
fprintf(stdout,"\n");
}
					}
                nseqs = profile1_nseqs = local_nseqs;
				info("No. of seqs=%d",(pint)nseqs);
				profile1_empty=FALSE;
				profile2_empty=TRUE;
	   }
      }
    else
      {			        /* first seq to be read = profile1_nseqs + 1 */
       local_nseqs = readseqs(profile1_nseqs+(sint)1); 
       if(local_nseqs < 0)               /* file could not be opened */
           { 
		return local_nseqs;
           }
       else if(local_nseqs == 0)         /* no sequences */
           {
	       error("No sequences in file!  Bad format?");
		return local_nseqs;
           }
       else if(local_nseqs > 0)
           {
		struct_penalties2 = NONE;
		if (sec_struct_mask2 != NULL) sec_struct_mask2=ckfree(sec_struct_mask2);
		if (gap_penalty_mask2 != NULL) gap_penalty_mask2=ckfree(gap_penalty_mask2);
		if (ss_name2 != NULL) ss_name2=ckfree(ss_name2);
                if (struct_penalties != NONE) /* feature table / mask in alignment */
                	{
					struct_penalties2 = struct_penalties;
					if (struct_penalties == SECST) {
						sec_struct_mask2 = (char *)ckalloc((max_aln_length) * sizeof (char));
						for (i=0;i<max_aln_length;i++)
							sec_struct_mask2[i] = sec_struct_mask[i];
					}
					gap_penalty_mask2 = (char *)ckalloc((max_aln_length) * sizeof (char));
					for (i=0;i<max_aln_length;i++)
						gap_penalty_mask2[i] = gap_penalty_mask[i];
        				ss_name2 = (char *)ckalloc( (MAXNAMES+1) * sizeof (char));
					strcpy(ss_name2,ss_name);
if (debug>0) {
for (i=0;i<seqlen_array[profile1_nseqs+1];i++)
	fprintf(stdout,"%c",gap_penalty_mask2[i]);
fprintf(stdout,"\n");