interface.c

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

/* command line interface for Clustal W  */
/* DES was here MARCH. 1994 */
/* DES was here SEPT.  1994 */
/* Fixed memory allocation bug in check_param() . Alan Bleasby Dec 2002 */

#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


char *nameonly(char *s) ;

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 seqRange; /*Ramu */

extern Boolean 	output_clustal, output_nbrf, output_phylip, output_gcg, output_gde, output_nexus, output_fasta;
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     *fasta_outfile; /* Ramu */
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];
static char     fasta_outname[FILENAMELEN+1];  /* Ramu */
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]="";

char pim_name[FILENAMELEN+1]=""; /* Ramu */

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;
		sint ajb;

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



	params[0]=(char *)ckalloc((strlen(args[0])+1)*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])+1)*sizeof(char));
		ajb=0;
		for(j=0;j<strlen(args[i])-1;j++)
			if(isprint(args[i][j+1])) params[i][ajb++]=args[i][j+1];
		params[i][ajb]='\0';
        }

        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,6);
			if (temp >= 0 && temp <= 5) {
				output_clustal = FALSE;
				output_gcg     = FALSE;
				output_phylip  = FALSE;
				output_nbrf    = FALSE;
				output_gde     = FALSE;
				output_nexus   = FALSE;
				output_fasta   = 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: /* NEXUS */
					output_nexus   = TRUE;
					break;
				case 5: /* NEXUS */
					output_fasta   = 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;
  
  /*** ? /sequences (sets type of second input file to list of sequences)  */
  if(setsequences != -1)
    profile_type = SEQUENCE;
  
  
  
  /*** ? /ktuple=n */
  if(setktuple != -1) {
    temp = 0;
    if(strlen(param_arg[setktuple]) > 0)
      if (sscanf(param_arg[setktuple],"%d",&temp)!=1) {
	fprintf(stdout,"Bad option for /ktuple (must be integer)\n");
	temp = 0;
      }
    if(temp > 0) {
      if(dnaflag) {
	if(temp <= 4) {