comp_lib2.c

序列对齐 Compare a protein sequence to a protein sequence database or a DNA sequence to a DNA sequenc

/* copyright (c) 1996, 1997, 1998, 1999, 2002  William R. Pearson and the
   U. of Virginia */

/* $Name: fa35_03_06 $ - $Id: comp_lib2.c,v 1.25 2008/02/19 08:50:13 wrp Exp $ */

/*
 * Jan 17, 2007 - remove #ifdef PRSS - begin better statistics in place
 * for small libraries, related libraries
 *
 * Concurrent read version
 *
 *	Feb 20, 1998 modifications for prss3
 *
 *	December, 1998 - DNA searches are now down with forward and reverse
 *			 strands
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>

#include <limits.h>
#include <float.h>
#include <math.h>

#ifdef UNIX
#include <unistd.h>
#include <sys/types.h>
#include <signal.h>
#endif

#include "defs.h"
#include "mm_file.h"

#include "mw.h"			/* defines beststr */
#include "structs.h"		/* mngmsg, libstruct */
#include "param.h"		/* pstruct, thr_str, buf_head, rstruct */
#include "thr_bufs.h"
#include "drop_func.h"

#define XTERNAL
#include "uascii.h"

char *mp_verstr="35.03";

/********************************/
/* global variable declarations */
/********************************/

extern int max_workers;

#ifdef SUPERFAMNUM
int nsfnum;
int sfnum[10];
extern int sfn_cmp(int *q, int *s);
int nsfnum_n;
int sfnum_n[10];
#endif

/********************************/
/* extern variable declarations */
/********************************/
extern char *prog_func;		/* function label */
extern char *verstr, *iprompt0, *iprompt1, *iprompt2, *refstr;

/********************************/
/*extern function declarations  */
/********************************/
struct lmf_str *openlib(char *, int, int *, int, struct lmf_str *);

int re_getlib(unsigned char *, int, int, int, int, int, long *, long *, 
	      struct lmf_str *m_fptr);

void set_shuffle(struct mngmsg m_msg);
void closelib(struct lmf_str *m_fptr);

void irand(int);
int nrand(int);

extern int ann_scan(unsigned char *, int, struct mngmsg *, int );
extern int scanseq(unsigned char *seq, int n, char *str);
extern void re_ascii(int *qascii, int *sascii);
extern int recode(unsigned char *seq, int n, int *qascii, int nsq);
extern void revcomp(unsigned char *seq, int n, int *c_nt);

extern void init_ascii(int is_ext, int *sascii, int is_dna);
extern void qshuffle(unsigned char *aa0, int n0, int nm0);
extern void free_pam2p(int **);

/* initialize environment (doinit.c) */
extern void initenv (int argc, char **argv, struct mngmsg *m_msg,
		     struct pstruct *ppst, unsigned char **aa0);

/* print timing information */
extern void ptime (FILE *, time_t);

#ifdef COMP_MLIB 
#define QGETLIB (q_file_p->getlib)
#endif

#define GETLIB (m_file_p->getlib)

void
save_best(struct buf_head *lib_buf, struct mngmsg, struct pstruct *ppst,
	  struct db_str *, FILE *fdata, int *, struct hist_str *, void **);

void
save_shuf(struct buf_head *lib_buf, int nitt, int shuff_max);

void preserve_seq(struct buf2_str *, struct seq_record *, struct beststr *);
void
buf_do_work(unsigned char **aa0, int n0, struct buf_head *lib_bhead_p,
	    struct pstruct *ppst, void **f_str);
void
buf_qshuf_work(unsigned char *aa0s, int n0, struct buf_head *lib_bhead_p,
	       struct pstruct *ppst, void *qf_str, int score_ix);
void
buf_shuf_work(unsigned char **aa0, int n0, unsigned char *aa1s, int zs_win,
	      struct buf_head *lib_bhead_p, struct pstruct *ppst,
	      void **f_str, int score_ix);

/* statistics functions */
extern int
process_hist(struct stat_str *sptr, int nstat, 
	     struct mngmsg m_msg,
	     struct pstruct *ppst,
	     struct hist_str *hist, void **, int);
extern void addhistz(double, struct hist_str *); /* scaleswn.c */
void selectbestz(struct beststr **, int, int );
extern double (*find_zp)(int score, double escore, int length, double comp,void *);

void last_stats(const unsigned char *, int, 
		struct stat_str *sptr, int nstats,
		struct beststr **bestp_arr, int nbest,
		struct mngmsg m_msg, struct pstruct *ppst, 
		struct hist_str *histp, void *);

int last_calc( unsigned char **aa0, unsigned char *aa1, int maxn,
	       struct beststr **bestp_arr, int nbest,
	       struct mngmsg m_msg, struct pstruct *ppst, 
	       void **f_str, void *rs_str);

void scale_scores(struct beststr **bestp_arr, int nbest,
		  struct db_str,struct pstruct *ppst, void *);

extern int shuffle(unsigned char *, unsigned char *, int);
extern int wshuffle(unsigned char *, unsigned char *, int, int);

extern void set_db_size(int, struct db_str *, struct hist_str *);

/* display functions */
extern void
showbest (FILE *fp, unsigned char **aa0, unsigned char *aa1, int maxn,
	  struct beststr **bestp_arr, int nbest,
	  int qlib, struct mngmsg *m_msg,struct pstruct *ppst,
	  struct db_str db, char **gstring2p, void **f_str);

extern void
showalign (FILE *fp, unsigned char **aa0, unsigned char *aa1, int maxn,
	   struct beststr **bestp_arr, int nbest,
	   int qlib, struct mngmsg m_msg,struct pstruct *ppst,
	   char **gstring2p, void **f_str);

/* misc functions */
void h_init(struct pstruct *, struct mngmsg *, char *);		/* doinit.c */
void last_init(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
void last_params(unsigned char *, int, struct mngmsg *, struct pstruct *);

void s_abort(char *, char *);		/* compacc.c */

/* initfa/sw.c */
void resetp(struct mngmsg *, struct pstruct *); 

void gettitle(char *, char *, int);	/* nxgetaa.c */
void libchoice(char *lname, int, struct mngmsg *); /* lib_sel.c */
void libselect(char *lname, struct mngmsg *);	/* lib_sel.c */
void query_parm(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
void selectbestz(struct beststr **, int, int);

/* compacc.c */
void prhist(FILE *, struct mngmsg, struct pstruct *, 
	    struct hist_str hist, int nstats, struct db_str, char *, char **, char **);
void printsum(FILE *, struct db_str db);
int reset_maxn(struct mngmsg *, int);	/* set m_msg.maxt, maxn from maxl */

FILE *outfd;			/* Output file */

/* this information is global for fsigint() */
extern time_t s_time();			/* fetches time */
time_t tstart, tscan, tprev, tdone;	/* Timing */
#ifdef COMP_MLIB
time_t ttscan, ttdisp;
#endif
time_t tdstart, tddone;

static struct db_str qtt = {0l, 0l, 0};

#ifdef COMP_THR
/***************************************/
/* thread global variable declarations */
/***************************************/

/* functions for getting/sending buffers to threads (thr_sub.c) */
extern void init_thr(int , struct thr_str *, struct mngmsg, struct pstruct *,
		     unsigned char *, int);
extern void start_thr(void);
extern void get_rbuf(struct buf_head **lib_buf, int max_wor_buf);
extern void put_rbuf(struct buf_head *lib_buf, int max_work_buf);
extern void wait_rbuf(int max_work_buf);
extern void rbuf_done(int nthreads);
extern void put_rbuf_done(int nthreads, struct buf_head *lib_buf, 
			  int max_work_buf);
#undef XTERNAL
#include "thr.h"
#endif

struct buf_head *lib_buf2_list;

/* these variables must be global for comp_thr.c so that savebest()
   can use them */
static struct beststr **bestp_arr;	/* array of pointers */
static int nbest;	/* number of best scores */

static struct stat_str *stats; /* array of scores for statistics from real
			     (or shuffled) sequences*/
static struct stat_str *qstats;	/* array of scores for shuffled query stats */
  struct stat_str *rstats;	/* array of scores from shuffled library */

  /* these variables are global so they can be set both by the main()
     program and savebest() in threaded mode.
  */
static int nstats, nqstats, nrstats, kstats;
static double zbestcut;		/* cut off for best z-score */
static int bestfull;		/* index for selectbest() */
static int stats_done=0;	/* flag for z-value processing */

void fsigint();

int
main (int argc, char *argv[]) 
{
  unsigned char *aa0[6], *aa0s, *aa1, *aa1ptr, *aa1a = NULL;
  unsigned char *aa1save;	/* aa1shuff and aa1save must be distinct */
  unsigned char *aa1shuff, *aa1shuff_b=NULL;	/* for new unthreaded version */
  int n1;
  int j;
  int shuff_mult, n1lib_req;
  struct beststr *bbp;

  int *n1tot_ptr=NULL, *n1tot_cur;
  int n1tot_cnt=0;
  int n1tot_v, aa1_loff;

  long loffset, l_off;	/* loffset is the coordinate of first residue
			   when lcont > 0; l_off is not used in the
			   main loop, only in showbest and showalign */
  struct a_res_str *next_ares, *old_ares; /* used to free-up old a_res */

  /* status/parameter information */
  char info_lib_range[MAX_FN];
  char *info_lib_range_p;
  char info_pgm_abbr[MAX_SSTR];
  char info_qlabel[MAX_FN];
  char *info_gstring2p[2];
  char info_gstring3[MAX_STR];
  char *info_hstring_p[2];
  char l_bline[MAX_SSTR];
#ifdef COMP_MLIB
  char q_bline[MAX_STR];
  fseek_t qseek;
  int qlib;
  struct lmf_str *q_file_p;
  int sstart, sstop, is;
#endif
  int lstart, lstop;
  int id;
  struct lmf_str *m_file_p;

  int t_best, t_rbest, t_qrbest;	/* best score of two/six frames */
  double t_escore, t_rescore, t_qrescore; /* best evalues of two/six frames */
  double db_tt;
  int i_score;

  struct pstruct pst;
  void *f_str[6], *qf_str;	/* different f_str[]'s for forward,reverse */
  int have_f_str=0;

  /* move buf_head, buf2_str outside COMP_THR to be used consistently for everything */
  struct buf_head *lib_bhead_p, *t_lib_bhead_p;

  struct buf2_str *lib_buf2_p;
  struct seq_record *current_seq_p;
  int buf2_cnt;

  long ntbuff;			/* current length (in residues/bytes of buffer */
  int max_buf2_cnt;		/* number of entries per buffer for threads */
  int max_do_cnt;		/* number of entries to use for shuffles/opts-aligns */
  int ave_seq_len, buf_siz;
  int max_work_buf;
  int filled_worker_bufs;
#ifdef COMP_THR
  int empty_reader_bufs;
  int t_reader_buf_readp;
  struct thr_str *work_info;
#endif

  struct mngmsg m_msg;		/* Message from host to manager */
  int iln, itt;			/* index into library names */
  char rline[MAX_FN];
  char argv_line[MAX_STR];
  int t_quiet;

  struct rstruct  rst;		/* results structure */
  struct rstruct  rrst;		/* results structure for shuffle*/
  int i;

  FILE *fdata=NULL;		/* file for full results */
  fseek_t lseek;		/* seek into library of current sequence */
  char libstr[MAX_UID];		/* required because getlib() does not return
				   lcont > 0 */
  struct beststr *best;		/* array of best scores */

  /* save sequence meta info for sequences that are not currently available */
  struct seq_record *best_seqs;

  int n_libstr;			/* length of libstr */
  int jstats;
  int leng;			/* leng is length of the descriptive line */
  int maxn;			/* size of the library sequence examined */
  int maxl;			/* size of library buffer */
  int qlcont;			/* continued query sequence */
  int lcont, ocont, maxt;	/* continued sequence */
  int igncnt=0;			/* count for ignoring sequences warning */
  double zscore;		/* tmp value */
  char *bp;			/* general purpose string ptr */
  
  /* this is necessary because of an SGI Irix 64 issue */
  info_gstring2p[0] = calloc(MAX_STR,sizeof(char));
  info_gstring2p[1] = calloc(MAX_STR,sizeof(char));
  info_hstring_p[0] = calloc(MAX_STR,sizeof(char));
  info_hstring_p[1] = calloc(MAX_STR,sizeof(char));

  /* Initialization */

#if defined(UNIX)
  m_msg.quiet= !isatty(1);
#else
  m_msg.quiet = 0;
#endif

#ifdef PGM_DOC
  argv_line[0]='\0';
  for (i=0; i<argc; i++) {
    strncat(argv_line," ",sizeof(argv_line)-strlen(argv_line)-1);
    if (strchr(argv[i],' ')) {
      strncat(argv_line,"\"",sizeof(argv_line)-strlen(argv_line)-1);
      strncat(argv_line,argv[i],sizeof(argv_line)-strlen(argv_line)-1);
      strncat(argv_line,"\"",sizeof(argv_line)-strlen(argv_line)-1);
    }
    else {
      strncat(argv_line,argv[i],sizeof(argv_line)-strlen(argv_line)-1);
    }
  }
  argv_line[sizeof(argv_line)-1]='\0';

  fprintf(stdout, "#%s\n",argv_line);
#endif

  /* first initialization routine - nothing is known */
  h_init(&pst, &m_msg, info_pgm_abbr);
  
  m_msg.db.length = m_msg.ldb.length = qtt.length = 0l;
  m_msg.db.entries = m_msg.db.carry = 
    m_msg.ldb.entries = m_msg.ldb.carry = qtt.entries = qtt.carry = 0;
  m_msg.pstat_void = NULL;
  m_msg.hist.entries = 0;

  for (iln=0; iln<MAX_LF; iln++) m_msg.lb_mfd[iln]=NULL;

  f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] = f_str[0] = NULL;

  aa0[0] = NULL;
  /* second initialiation - get commmand line arguments */
  initenv (argc, argv, &m_msg, &pst,&aa0[0]);

  if (m_msg.markx & MX_M11OUT) {
    fprintf(stdout, "#:lav\n\nd {\n   \"%s\"\n}",argv_line+1);
  }

#ifdef COMP_THR
  /* now have max_workers - allocate work_info[] */
  if (max_workers >= MAX_WORKERS) max_workers = MAX_WORKERS;
  if ((work_info=
       (struct thr_str *)calloc(max_workers,sizeof(struct thr_str)))==NULL) {
    fprintf(stderr, " cannot allocate work_info[%d]\n",max_workers);
    exit(1);
  }
#else
  max_workers = 1;
#endif

  tstart = tscan = s_time();
  tdstart = time(NULL);

  /* Allocate space for the query and library sequences */
  /* pad aa0[] with an extra SEQ_PAD chars for ALTIVEC padding */
  if (aa0[0]==NULL) {
    if ((aa0[0] = (unsigned char *)malloc((m_msg.max_tot+1+SEQ_PAD)*sizeof(unsigned char)))
	== NULL)
      s_abort ("Unable to allocate query sequence", "");
    *aa0[0]=0;
    aa0[0]++;
  }
  aa0[5]=aa0[4]=aa0[3]=aa0[2]=aa0[1]=aa0[0];

  if ((aa1save = (unsigned char *)malloc((m_msg.max_tot+1)*sizeof (char))) == NULL) {
    s_abort ("Unable to allocate library overlap", "");
  }
  *aa1save=0;
  aa1save++;

  if (m_msg.markx & MX_HTML) {
#ifdef HTML_HEAD    
    fprintf(stdout,"<html>\n<head>\n<title>%s Results</title>\n</head>\n<body>\n",prog_func);
#endif
    fprintf(stdout,"<pre>\n");
  }

  if (m_msg.std_output) {
    fprintf(stdout,"%s\n",iprompt0);
    fprintf(stdout," %s%s\n",verstr,refstr);
  }
  if (m_msg.markx & MX_HTML) fputs("</pre>\n",stdout);

  /* Query library */
  if (m_msg.tname[0] == '\0') {