interface.c

是有关基因比对的经典算法的实现。这对于初学计算生物学的人是非常重要的算法。

/* command line interface for Clustal W  */
/* DES was here MARCH. 1994 */
/* DES was here SEPT.  1994 */
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <signal.h>
#include <setjmp.h>
#include "clustalw.h"
#include "param.h"

/*
*	Prototypes
*/

#ifdef UNIX
FILE    *open_path(char *);
#endif


static sint check_param(char **args,char *params[], char *param_arg[]);
static void set_optional_param(void);
static sint find_match(char *probe, char *list[], sint n);
static void show_aln(void);
static void create_parameter_output(void);
static void reset_align(void);
static void reset_prf1(void);
static void reset_prf2(void);
static void calc_gap_penalty_mask(int prf_length,char *struct_mask,char *gap_mask);
void print_sec_struct_mask(int prf_length,char *mask,char *struct_mask);
/*
*	 Global variables
*/

extern sint max_names;

extern Boolean interactive;

extern double  **tmat;
extern float    gap_open,      gap_extend;
extern float  	dna_gap_open,  dna_gap_extend;
extern float 	prot_gap_open, prot_gap_extend;
extern float    pw_go_penalty,      pw_ge_penalty;
extern float  	dna_pw_go_penalty,  dna_pw_ge_penalty;
extern float 	prot_pw_go_penalty, prot_pw_ge_penalty;
extern char 	revision_level[];
extern sint    wind_gap,ktup,window,signif;
extern sint    dna_wind_gap, dna_ktup, dna_window, dna_signif;
extern sint    prot_wind_gap,prot_ktup,prot_window,prot_signif;
extern sint	boot_ntrials;		/* number of bootstrap trials */
extern sint	nseqs;
extern sint	new_seq;
extern sint 	*seqlen_array;
extern sint 	divergence_cutoff;
extern sint 	debug;
extern Boolean 	no_weights;
extern Boolean 	neg_matrix;
extern Boolean  quick_pairalign;
extern Boolean	reset_alignments_new;		/* DES */
extern Boolean	reset_alignments_all;		/* DES */
extern sint 	gap_dist;
extern Boolean 	no_hyd_penalties, no_pref_penalties;
extern sint 	max_aa;
extern sint 	gap_pos1, gap_pos2;
extern sint  	max_aln_length;
extern sint 	*output_index, output_order;
extern sint profile_no;
extern short 	usermat[], pw_usermat[];
extern short 	aa_xref[], pw_aa_xref[];
extern short 	userdnamat[], pw_userdnamat[];
extern short 	dna_xref[], pw_dna_xref[];
extern sint	*seq_weight;

extern Boolean 	lowercase; /* Flag for GDE output - set on comm. line*/
extern Boolean 	cl_seq_numbers;
extern Boolean 	output_clustal, output_nbrf, output_phylip, output_gcg, output_gde, output_nexus;
extern Boolean 	output_tree_clustal, output_tree_phylip, output_tree_distances, output_tree_nexus;
extern sint     bootstrap_format;
extern Boolean 	tossgaps, kimura;
extern Boolean  percent;
extern Boolean 	explicit_dnaflag;  /* Explicit setting of sequence type on comm.line*/
extern Boolean 	usemenu;
extern Boolean 	showaln, save_parameters;
extern Boolean	dnaflag;
extern float	transition_weight;
extern unsigned sint boot_ran_seed;


extern FILE 	*tree;
extern FILE 	*clustal_outfile, *gcg_outfile, *nbrf_outfile, *phylip_outfile, *nexus_outfile;
extern FILE 	*gde_outfile;

extern char 	hyd_residues[];
extern char 	*amino_acid_codes;
extern char 	**args;
extern char	seqname[];

extern char 	**seq_array;
extern char 	**names, **titles;

extern char *gap_penalty_mask1,*gap_penalty_mask2;
extern char *sec_struct_mask1,*sec_struct_mask2;
extern sint struct_penalties,struct_penalties1,struct_penalties2;
extern sint output_struct_penalties;
extern Boolean use_ss1, use_ss2;
extern char *ss_name1,*ss_name2;


char *ss_name = NULL;
char *sec_struct_mask = NULL;
char *gap_penalty_mask = NULL;

char  	profile1_name[FILENAMELEN+1];
char  	profile2_name[FILENAMELEN+1];

Boolean empty;
Boolean profile1_empty, profile2_empty;   /* whether or not profiles   */

char 	outfile_name[FILENAMELEN+1]="";

static char 	clustal_outname[FILENAMELEN+1], gcg_outname[FILENAMELEN+1];
static char  	phylip_outname[FILENAMELEN+1],nbrf_outname[FILENAMELEN+1];
static char  	gde_outname[FILENAMELEN+1],nexus_outname[FILENAMELEN+1];

char     clustal_tree_name[FILENAMELEN+1]="";
char     dist_tree_name[FILENAMELEN+1]="";
char 	phylip_tree_name[FILENAMELEN+1]="";
char 	nexus_tree_name[FILENAMELEN+1]="";
char 	p1_tree_name[FILENAMELEN+1]="";
char 	p2_tree_name[FILENAMELEN+1]="";

static char *params[MAXARGS];
static char *param_arg[MAXARGS];

static char *cmd_line_type[] = {
                " ",
                "=n ",
                "=f ",
                "=string ",
                "=filename ",
                ""};

static sint numparams;
static Boolean check_tree = TRUE;

sint 	profile1_nseqs;	/* have been filled; the no. of seqs in prof 1*/
Boolean use_tree_file = FALSE,new_tree_file = FALSE;
Boolean use_tree1_file = FALSE, use_tree2_file = FALSE;
Boolean new_tree1_file = FALSE, new_tree2_file = FALSE;

static char *lin2;

MatMenu dnamatrix_menu = {3,
                "IUB","iub",
                "CLUSTALW(1.6)","clustalw",
                "User defined",""
		};

MatMenu matrix_menu = {5,
                "BLOSUM series","blosum",
                "PAM series","pam",
                "Gonnet series","gonnet",
                "Identity matrix","id",
                "User defined",""
		};
 
MatMenu pw_matrix_menu = {5,
                "BLOSUM 30","blosum",
                "PAM 350","pam",
                "Gonnet 250","gonnet",
                "Identity matrix","id",
                "User defined",""
		};

void init_interface(void)
{
	empty=TRUE;
     
	profile1_empty = TRUE;     /*  */
	profile2_empty = TRUE;     /*  */

        lin2 = (char *)ckalloc( (MAXLINE+1) * sizeof (char) );

}




static sint check_param(char **args,char *params[], char *param_arg[])
{

/*
#ifndef MAC
        char *strtok(char *s1, const char *s2);
#endif
*/
        sint     len,i,j,k,s,n,match[MAXARGS];
		Boolean 	name1 = FALSE;


	if(args[0]==NULL) return;

	params[0]=(char *)ckalloc(strlen(args[0])*sizeof(char));
	if (args[0][0]!=COMMANDSEP)
	{
		name1 = TRUE;
		strcpy(params[0],args[0]);
	}
	else
		strcpy(params[0],&args[0][1]);

        for (i=1;i<MAXARGS;i++) {
		if(args[i]==NULL) break;
		params[i]=(char *)ckalloc(strlen(args[i])*sizeof(char));
		for(j=0;j<strlen(args[i])-1;j++)
			if(isprint(args[i][j+1])) params[i][j]=args[i][j+1];
        }
        if (i==MAXARGS) {
		fprintf(stdout,"Error: too many command line arguments\n");
 		return(-1);
	}
/*
    special case - first parameter is input filename
*/
	s = 0;
	if(name1 == TRUE) {
		strcpy(seqname, params[0]);
/*  JULIE
    convert to lower case now
*/
#ifndef UNIX
		for(k=0;k<(sint)strlen(params[0]);++k) seqname[k]=tolower(params[0][k]);
#else
		for(k=0;k<(sint)strlen(params[0]);++k) seqname[k]=params[0][k];
#endif 
		s++;
	}
	
	n = i;
    for (i=s;i<n;i++) {
        param_arg[i] = NULL;
		len = (sint)strlen(params[i]);
		for(j=0; j<len; j++)
			if(params[i][j] == '=') {
				param_arg[i] = (char *)ckalloc((len-j) * sizeof(char));
				strncpy(param_arg[i],&params[i][j+1],len-j-1);
				params[i][j] = EOS;
/*  JULIE
    convert keywords to lower case now
*/
				for(k=0;k<j;++k) params[i][k]=tolower(params[i][k]);
				param_arg[i][len-j-1] = EOS;
				break;
			}
	}
	
/*
    for each parameter given on the command line, first search the list of recognised optional 
    parameters....
*/

    for (i=0;i<n;i++) {
		if ((i==0) && (name1 == TRUE)) continue;
		j = 0;
		match[i] = -1;
		for(;;) {
			if (cmd_line_para[j].str[0] == '\0') break;
			if (!strcmp(params[i],cmd_line_para[j].str)) {
				match[i] = j;
				*cmd_line_para[match[i]].flag = i;
				if ((cmd_line_para[match[i]].type != NOARG) &&
                                    (param_arg[i] == NULL)) {
					fprintf(stdout,
                       				 "Error: parameter required for /%s\n",params[i]);
					exit(1);
				}
/*  JULIE
    convert parameters to lower case now, unless the parameter is a filename
*/
#ifdef UNIX
				else if (cmd_line_para[match[i]].type != FILARG
				         && param_arg[i] != NULL)
#endif 
				if (param_arg[i]!=0)
					{
					for(k=0;k<strlen(param_arg[i]);++k)
param_arg[i][k]=tolower(param_arg[i][k]);
					}
				break;
			}
			j++;
		}
	}
/*
    ....then the list of recognised input files,.... 
*/
    for (i=0;i<n;i++) {
		if ((i==0) && (name1 == TRUE)) continue;
		if (match[i] != -1) continue;
		j = 0;
		for(;;) {
			if (cmd_line_file[j].str[0] == '\0') break;
			if (!strcmp(params[i],cmd_line_file[j].str)) {
				match[i] = j;
				*cmd_line_file[match[i]].flag = i;
				if ((cmd_line_file[match[i]].type != NOARG) &&
                                    (param_arg[i] == NULL)) {
					fprintf(stdout,
                       				 "Error: parameter required for /%s\n",params[i]);
					exit(1);
				}
				break;
			}
			j++;
		}
	}
/*
	....and finally the recognised verbs. 
*/
    for (i=0;i<n;i++) {
		if ((i==0) && (name1 == TRUE)) continue;
		if (match[i] != -1) continue;
		j = 0;
		for(;;) {
			if (cmd_line_verb[j].str[0] == '\0') break;
			if (!strcmp(params[i],cmd_line_verb[j].str)) {
				match[i] = j;
				*cmd_line_verb[match[i]].flag = i;
				if ((cmd_line_verb[match[i]].type != NOARG) &&
                                    (param_arg[i] == NULL)) {
					fprintf(stdout,
                       				 "Error: parameter required for /%s\n",params[i]);
					exit(1);
				}
				break;
			}
			j++;
		}
	}

/*
	check for any unrecognised parameters.
*/
    for (i=0;i<n;i++) {
		if (match[i] == -1) {
			fprintf(stdout,
                        "Error: unknown option %c%s\n",COMMANDSEP,params[i]);
			exit(1);
		}
	}
        return(n);
}

static void set_optional_param(void)
{
	int i,temp;
	int c;
	float ftemp;
	char tstr[100];

/****************************************************************************/
/* look for parameters on command line  e.g. gap penalties, k-tuple etc.    */
/****************************************************************************/

/*** ? /score=percent or /score=absolute */
        if(setscore != -1)
                if(strlen(param_arg[setscore]) > 0) {
     			temp = find_match(param_arg[setscore],score_arg,2);
                        if(temp == 0)
                                percent = TRUE;
                        else if(temp == 1)
                                percent = FALSE;
                        else
                                fprintf(stdout,"\nUnknown SCORE type: %s\n",
                                param_arg[setscore]);
                }

/*** ? /seed=n */
        if(setseed != -1) {
                temp = 0;
                if(strlen(param_arg[setseed]) > 0)
                         if (sscanf(param_arg[setseed],"%d",&temp)!=1) {
                         	fprintf(stdout,"Bad option for /seed (must be integer)\n");
                         	temp = 0;
                         }
                if(temp > 0) boot_ran_seed = temp;
        fprintf(stdout,"\ntemp = %d; seed = %u;\n",(pint)temp,boot_ran_seed);
        }




/*** ? /output=PIR, GCG, GDE or PHYLIP */
		if(setoutput != -1)
		if(strlen(param_arg[setoutput]) > 0) {
			temp = find_match(param_arg[setoutput],output_arg,5);
			if (temp >= 0 && temp <= 3) {
				output_clustal = FALSE;
				output_gcg     = FALSE;
				output_phylip  = FALSE;
				output_nbrf    = FALSE;
				output_gde     = FALSE;
				output_nexus   = FALSE;
			}
			switch (temp) {
				case 0: /* GCG */
					output_gcg     = TRUE;
					break;
				case 1: /* GDE */
					output_gde     = TRUE;
					break;
				case 2: /* PIR */
					output_nbrf    = TRUE;
					break;
				case 3: /* PHYLIP */
					output_phylip  = TRUE;
					break;
				case 4: /* PHYLIP */
					output_nexus   = TRUE;
					break;
				default:
					fprintf(stdout,"\nUnknown OUTPUT type: %s\n",
					param_arg[setoutput]);
			}
		}

/*** ? /outputtree=NJ or PHYLIP or DIST or NEXUS */
	if(setoutputtree != -1)
		if(strlen(param_arg[setoutputtree]) > 0) {
			temp = find_match(param_arg[setoutputtree],outputtree_arg,4);
			switch (temp) {
				case 0: /* NJ */
					output_tree_clustal = TRUE;
					break;
				case 1: /* PHYLIP */
					output_tree_phylip  = TRUE;
					break;
				case 2: /* DIST */
					output_tree_distances = TRUE;
					break;
				case 3: /* NEXUS */
					output_tree_nexus = TRUE;
					break;
				default:
					fprintf(stdout,"\nUnknown OUTPUT TREE type: %s\n",
					param_arg[setoutputtree]);
			}
		}

/*** ? /profile (sets type of second input file to profile) */
	if(setprofile != -1)
		profile_type = PROFILE;