p2_workcomp2.c

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

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

/* $Name: fa35_03_06 $ - $Id: p2_workcomp2.c,v 1.5 2007/12/13 14:38:20 wrp Exp $ */

/* This version is modifed to read all files, query and database,
   through the manager process. Workers will now receive their
   database from the manager, rather than reading it themselves.  This
   cuts down considerably on NFS traffic, simplifies searches of
   multiple files, and allows use of clusters of slave nodes that do
   not have NFS access */

/* September, 1994 - this version has been modified to do two kinds of
   searches, a general library search, or list of library sequences search.
   The latter would be used to generate optimized scores for fasta and
   to produce alignments */

/* modified July, 2002, to provide query shuffle */

/* modified October, 2005, to support struct a_res_str a_res -
   coordinates of alignment in aa0[], aa1[].  Future modifications
   will cause do_walign to be run only once - subsequent calls for
   seqc[0,1] can be filled using a_res, by adding a_res to the
   struct sqs2 array.

   19-March-2006 - modifications to call do_walign() only once, and
   use the resulting a_res structure for subsequent calls to calc_id,
   calcons, calcons_a, have been implemented.  Also, the -V option is
   now valid with the parallel programs.

   31-May-2006 - some functions (e.g. dropfs and dropff do not store
   complete information in a_res - thus they cannot use this shortcut
   (yet).

*/



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

#include <limits.h>

#ifdef PVM_SRC
#include "pvm3.h"
#endif

#ifdef MPI_SRC
#include "mpi.h"
#endif

/*
#define PvmDataDefault 0
#define PvmTaskDefault 0
*/
#include "msg.h"
#include "defs.h"
#include "param.h"
#include "w_mw.h"
#include "structs.h"

#ifdef MPI_SRC
#define XTERNAL
#endif
#include "upam.h"
#include "uascii.h"

#ifdef PVM_SRC
int worker, mytid;
int nnodes, pinums[MAXNOD];
#endif

#include "drop_func.h"

extern void alloc_pam (int d1, int d2, struct pstruct *ppst); /* allocate ppst->pam12,pam12x */
extern int **alloc_pam2p (int len, int nsq); 
extern void w_init ();
extern void irand(int);
extern void revcomp(unsigned char *, int, int *);

extern void initseq(char **seqc0, char **seqc0a, char **seqc1, char **seqca, int seqsiz);
extern void freeseq(char **seqc0, char **seqc0a, char **seqc1, char **seqca);

void send_bestr(int, int, struct comstr *, int, int);
void send_bestr2(int, struct comstr2 *, int);
void send_struct(int hosttid, int msg_type, void *void_p, size_t slen);
void send_code(int, char *, int);

extern void update_param(struct qmng_str *qm_msg, struct mngmsg *m_msg,
			 struct pstruct *ppst);
extern int shuffle(unsigned char *, unsigned char *, int);
extern int wshuffle(unsigned char *, unsigned char *, int, int, int *);

extern char err_str[];

/* local function declarations */
void free_ares(struct sqs2 *, int itt, int *, int walign_cnt, int worker);



void w_abort (p, p1)
char *p, *p1;
{
    fprintf (stderr, " %s %s\n", p, p1);
#ifdef PVM_SRC
    pvm_exit();
    exit (1);
#endif
#ifdef MPI_SRC
  MPI_Abort(MPI_COMM_WORLD,1);
#endif
}

#ifdef PVM_SRC
main ()
#endif
#ifdef MPI_SRC
void
workcomp(int worker)
#endif
{
  unsigned char *aa0[6], *aa1s, *aa0s;	/* Query and library sequences */
  struct mngmsg m_msg;	/* start message from manager to worker 1 */
  struct qmng_str qm_msg; /* updated for each query */
  int last_msg_b[10];	/* last set of numbers */
  struct sqs2 *seqpt; 		/* sequence pointers for chunk */
  int seqbuf_n,seqbuf_s;	/* number of sequences, length of sequences */
  int max_sql;			/* maximum number of sequences/node */
  int *n1_arr;			/* array of sequence lengths in buffer */
  int *m_seqnm_arr;		/* array of sequence numbers in buffer */
  int *aa1i_arr;		/* array of offsets into the buffer */
  unsigned char *seq_buf;	/* space for sequence data */
  int ntx;
  int nsq;			/* effective alphabet size */
  long curtype = ONETYPE;	/* current send message type */
  int ieven=0;			/* flag for window shuffle */
  int cur_n0;
  int n1, n1over;		/* length of query, library sequences */
  struct comstr bestr[BFR+1];	/* best structures */
  struct comstr2 bestr2[BFR2+1];	/* best structures */
  struct a_struct aln, *aln_dp;
  struct db_str ldb;
  int qres_bufsize;		/* results buffer size */
  int bestcnt = 0;		/* how many best structures are full */
  char info_gstring2[2][MAX_STR];		/* parameter string for manager */
  char *info_gstring2p[2] = {info_gstring2[0], info_gstring2[1]};		/* parameter string for manager */
  char info_gstring3[MAX_STR];		/* parameter string for manager */
  struct pstruct pst;		/* parameter structure */
  struct rstruct rst, qrst, rrst;	/* results structure */
  void *f_str[6], *qf_str;
  int sw_score;
  int lcnt, count, seqnm;	/* counters */
  int *walign_done[2], walign_cnt[2];	/* index of current valid a_res in seqpt */
  int have_walign;
  int *tres;			/* allocated storage for seqpt[].a_res[].res */
  int lend;			/*global library sequence number information */
  int lsn;			/* library sequence number */
  struct stage2_str *liblist=NULL;	/* list of sequences to search */
  int i, j;		/* my turn to send sequence descriptions */
  char lib_range[MAX_SSTR];
  char libstr[21];
  char errstr[128];
  int itt=0;
  int bufid;
  char *seqc0, *seqc0a, *seqc1, *seqca;
  char *seqc, *seqc_buff;
  int seqc_buff_cnt, seqc_buff_len, seqc_flag;
  int maxc, lc, nc, nident, ngap, aln_code_n;
  float percent, gpercent;
  int old_shuffle=0;	/* did a qshuffle last time */
  int hosttid=0;
  char worker_str[5];

#ifdef MPI_SRC
  MPI_Status mpi_status;
#endif
  
#ifdef PVM_SRC
  mytid = pvm_mytid();
  hosttid = pvm_parent();
#endif

  w_init();	/* sets up default sascii, hsq, sq */

  /* Allocate space for the query sequence */
  if ((aa0[0] = (unsigned char *) malloc ((MAXTST+2+SEQ_PAD)*sizeof (char))) == NULL) {
    w_abort ("Unable to allocate sequence array[0] - exiting!","");
  }
  *aa0[0]='\0';
  aa0[0]++;

  /* initial messages set up various parameter structures:

     STARTTYPE0: &nnodes
     		 pinums
		 &m_msg

     STARTTYPE1  &pst

     STARTTYPE2	 pam12
     STARTTYPE3	 pam12x
  */

#ifdef PVM_SRC
#ifdef ROUTE_DIRECT
  pvm_setopt(PvmRoute,PvmRouteDirect);
#endif
  /* get number of nodes, pinums */
  bufid = pvm_recv(hosttid,STARTTYPE0);
  pvm_upkint(&nnodes,1,1);
  pvm_upkint(pinums,nnodes,1);
  pvm_upkbyte((char *)&m_msg,(int)sizeof(m_msg),1);
  worker = tidtonode(mytid);
  pvm_freebuf(bufid);
#endif

  sprintf(worker_str,"@%d",worker);

#ifdef MPI_SRC
  MPI_Recv(&m_msg,sizeof(m_msg),MPI_BYTE,hosttid,STARTTYPE0,MPI_COMM_WORLD,
	   &mpi_status);
#endif

  /* the aln structure needs some information from m_msg0.aln */
  memcpy(&aln,&m_msg.aln,sizeof(struct a_struct));

  /*
  fprintf(stderr,"d1: %d d2: %d\n",m_msg.pamd1,m_msg.pamd2);
  */

  /* get pst params */
#ifdef PVM_SRC
  bufid = pvm_recv(hosttid,STARTTYPE1);
  pvm_upkbyte((char *)&pst,(int)sizeof(pst),1);
  /* 31t nsq = pst.nsq; */
  pvm_upkbyte((char *)pascii,(int)sizeof(aascii),1);
  pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
  MPI_Recv(&pst,(int)sizeof(pst),MPI_BYTE,hosttid,STARTTYPE1,MPI_COMM_WORLD,
	   &mpi_status);

  MPI_Recv(pascii,(int)sizeof(aascii)/sizeof(int),MPI_INT,hosttid,STARTTYPE1,MPI_COMM_WORLD,
	   &mpi_status);
#endif

  if (pst.ext_sq_set) { nsq = pst.nsqx;}
  else { nsq = pst.nsq;}

  aa0[5] = aa0[4] = aa0[3] = aa0[2] = aa0[1] = aa0[0];
  if (m_msg.qframe == 2) {
    if ((aa0[1]=(unsigned char *)malloc((MAXTST+2)*sizeof (char)))==NULL)
      w_abort ("Unable to allocate sequence[1] array - exiting!","");
    *aa0[1]='\0';
    aa0[1]++;
  }

  if ((aa1s=(unsigned char *)malloc((m_msg.max_tot+1)*sizeof (char)))==NULL)
      w_abort ("Unable to allocate shuffled library sequence", "");
  *aa1s=0;
  aa1s++;

  irand(0);	/* necessary for shuffled sequences */

  /* this function allocates pam12, pam12x
     assigns pst.pam[0][0]=pam12, pst.pam[1][0] = pam12x
     and sets up the correct pst.pam[0][0][0] pointers */

  alloc_pam(m_msg.pamd1,m_msg.pamd2,&pst);

#ifdef PVM_SRC
  bufid = pvm_recv(hosttid,STARTTYPE2);
  pvm_upkint(pam12,m_msg.pamd1*m_msg.pamd2,1);
  pvm_freebuf(bufid);

  bufid = pvm_recv(hosttid,STARTTYPE3);
  pvm_upkint(pam12x,m_msg.pamd1*m_msg.pamd2,1);
  pvm_freebuf(bufid);
#endif  

#ifdef DEBUG
  if (worker==FIRSTNODE) {
    fprintf(stderr,"ext?: %d\tnsq: %d\tnsqx: %d\n",pst.ext_sq_set,pst.nsq, pst.nsqx);
    for (i=1; i<5; i++) {
      for (j=1; j <= i; j++) fprintf(stderr," %c,%c:%2d",pst.sq[i],pst.sq[j],pst.pam2[0][i][j]);
      fprintf(stderr,"\n");
    }
    for (i=pst.nsq+1; i<pst.nsq+5; i++) {
      for (j=pst.nsq+1; j <= i; j++) fprintf(stderr," %c,%c:%2d",pst.sqx[i],pst.sqx[j],pst.pam2[0][i][j]);
      fprintf(stderr,"\n");
    }

    for (i=1; i<5; i++) {
      for (j=1; j <= i; j++) fprintf(stderr," %c,%c:%2d",pst.sqx[i],pst.sqx[j],pst.pam2[1][i][j]);
      fprintf(stderr,"\n");
    }
    for (i=pst.nsq+1; i<pst.nsq+5; i++) {
      for (j=pst.nsq+1; j <= i; j++) fprintf(stderr," %c,%c:%2d",pst.sqx[i],pst.sqx[j],pst.pam2[1][i][j]);
      fprintf(stderr,"\n");
    }
  }
#endif

#ifdef MPI_SRC
  MPI_Recv(pam12,m_msg.pamd1*m_msg.pamd2,MPI_INT,hosttid,STARTTYPE2,
	   MPI_COMM_WORLD,&mpi_status);

  MPI_Recv(pam12x,m_msg.pamd1*m_msg.pamd2,MPI_INT,hosttid,STARTTYPE3,
	   MPI_COMM_WORLD,&mpi_status);
#endif

/*
  We have the PAM matrices - get the library sequences
*/

  /* Allocate space for the sequences */
  max_sql = MAXSQL/2;

  if ((seqpt=(struct sqs2 *)calloc(max_sql,sizeof(struct sqs2)))==NULL)
    w_abort("cannot allocate seqpt(sqs2)","");

  if ((n1_arr=(int *)calloc(m_msg.pbuf_siz+1,sizeof(int)))==NULL)
    w_abort("cannot allocate n1_arr","");

  if ((aa1i_arr=(int *)calloc(m_msg.pbuf_siz+1,sizeof(int)))==NULL)
    w_abort("cannot allocate n1_arr","");

  if ((m_seqnm_arr=(int *)calloc(m_msg.pbuf_siz+1,sizeof(int)))==NULL)
    w_abort("cannot allocate m_seqnm_arr","");

/*****************************************************************/
/* This section gets all the database sequences from the manager */
/*****************************************************************/

  lcnt = 0;
  while (1) {
#ifdef PVM_SRC
    /* get the number of sequences, sequence lengths */
    bufid = pvm_recv(hosttid,STARTTYPE4);
    pvm_upkint(&seqbuf_n,1,1);	/* number of sequences */
    pvm_upkint(&seqbuf_s,1,1);	/* size of sequence buffer */
    pvm_upkint(n1_arr,seqbuf_n,1);	/* length of each sequence in buffer */
    pvm_upkint(aa1i_arr,seqbuf_n,1);    /* indexes for each sequence */
    pvm_upkint(m_seqnm_arr,seqbuf_n,1);	/* number of each library sequence */
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv(&seqbuf_n,1,MPI_INT,hosttid,STARTTYPE4,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(&seqbuf_s,1,MPI_INT,hosttid,STARTTYPE4,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(n1_arr,seqbuf_n,MPI_INT,hosttid,STARTTYPE4,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(aa1i_arr,seqbuf_n,MPI_INT,hosttid,STARTTYPE4,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(m_seqnm_arr,seqbuf_n,MPI_INT,hosttid,STARTTYPE4,MPI_COMM_WORLD,
	     &mpi_status);
#endif

    if (seqbuf_n <= 0) break;
#ifdef DEBUG    
    /*
    fprintf(stderr,"[%d] seqbuf_n: %d seqbuf_s: %d\n",
	    worker,seqbuf_n,seqbuf_s);
    fprintf(stderr,"[%d] lcnt: %d n1: %d seqnm %d\n",
	    worker,0,n1_arr[0],m_seqnm_arr[0]);
    fprintf(stderr,"[%d] lcnt: %d n1: %d seqnm %d\n",
	    worker,1,n1_arr[1],m_seqnm_arr[1]);
    */
#endif

    /* allocate space for sequences */
    if ((seq_buf = (unsigned char *)calloc((size_t)seqbuf_s+1,sizeof(char)))
	==NULL) {
      w_abort("cannot allocate tmp_seq","");
    }
    seq_buf++; /* leave a '\0' at the start */

    /* get the sequence buffer */
#ifdef PVM_SRC
    bufid = pvm_recv(hosttid,STARTTYPE5);
    pvm_upkbyte((char *)seq_buf,seqbuf_s,1);
    pvm_freebuf(bufid);
#endif  
#ifdef MPI_SRC
    MPI_Recv(seq_buf,seqbuf_s,MPI_BYTE,hosttid,STARTTYPE5,MPI_COMM_WORLD,
	     &mpi_status);
#endif

    /* now we have everything  - update the pointers */
    if (lcnt+seqbuf_n >= max_sql) {
      max_sql += max(MAXSQL/2,seqbuf_n);
      if ((seqpt=(struct sqs2 *)realloc(seqpt,max_sql*sizeof(struct sqs2)))
	  ==NULL)
	w_abort("cannot allocate seqpt(sqs2)","");
    }

    /* convert from offsets to pointers into buffer */
    /* ntx = 0; */
    for (i=0; i<seqbuf_n; i++,lcnt++) {
      seqpt[lcnt].n1 = n1_arr[i];
      seqpt[lcnt].m_seqnm = m_seqnm_arr[i];
      seqpt[lcnt].aa1 = &seq_buf[aa1i_arr[i]];
      /*      ntx += n1_arr[i]+1 + SEQ_PAD */

#ifdef DEBUG
      /* must have null's at both ends of sequence */
      if (seqpt[lcnt].aa1[-1]!= '\0') {
	fprintf(stderr,"Missing null at start: %d %d\n",
		lcnt,seqpt[lcnt].aa1[-1]);
	seqpt[lcnt].aa1[-1]='\0';
      }
      if (seqpt[lcnt].aa1[seqpt[lcnt].n1]!= '\0') {
	fprintf(stderr,"Missing null at end: %d %d\n",
		lcnt,seqpt[lcnt].aa1[seqpt[lcnt].n1]);
	seqpt[lcnt].aa1[seqpt[lcnt].n1]='\0';
      }