p2_workcomp2.c

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

#endif
    }
  }
  /* all done - lcnt has the total number of library sequences */

#ifdef DEBUG
  if (lcnt > 0)
    for (i=0; i<10; i++) {
      for (j=0; j<10; j++) libstr[j]=pst.sq[seqpt[i].aa1[j]];
      libstr[10]='\0';
      fprintf(stderr,"[%d] n1: %d seqnm: %d aa1: %s\n",
	      worker,seqpt[i].n1,seqpt[i].m_seqnm,libstr);
    }
#endif

  /* send back the number of descriptions received */

#ifdef PVM_SRC
  pvm_initsend(PvmDataRaw);
  pvm_pkint(&lcnt,1,1);
  pvm_send(hosttid,STARTTYPE0);
#endif
#ifdef MPI_SRC
/*  p4_dprintf(" have %d descriptions to send\n",lcnt); */
  MPI_Send(&lcnt,1,MPI_INT,hosttid,STARTTYPE0,MPI_COMM_WORLD);
#endif  

/*****************************************************************/
/* Library reads are finished, get ready to do searches          */
/*****************************************************************/

  /* get last set of numbers */
#ifdef PVM_SRC
  bufid = pvm_recv(hosttid,STARTTYPE0);
  pvm_upkint(last_msg_b,2,1);
  pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
  MPI_Recv(last_msg_b, 2, MPI_INT, hosttid, STARTTYPE0, MPI_COMM_WORLD,
	     &mpi_status);
#endif
  m_msg.nbr_seq = last_msg_b[0];
  qres_bufsize = last_msg_b[1];

#ifdef DEBUG
#ifdef PVM_SRC
  fprintf(stderr,"[%d] have nbr_seq %d qres_bufsize %d\n",worker,
	     m_msg.nbr_seq, qres_bufsize);
#endif
#ifdef MPI_SRC
  /*  p4_dprintf("[%d] have nbr_seq %d qres_bufsize %d\n",worker,
	     m_msg.nbr_seq, qres_bufsize);
  */;
#endif
#endif  
  /* If self search, receive sequence numbering data */
  if (m_msg.self) {
#ifdef PVM_SRC
    bufid = pvm_recv(hosttid,STARTTYPE1);
    pvm_upkint(&lend,1,1);
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv(&lend,1,MPI_INT,hosttid,STARTTYPE1,MPI_COMM_WORLD,&mpi_status);
#endif
  }
  
  /* allocate space for a_res flag array */

  if ((walign_done[0] = (int *)calloc(lcnt,sizeof(int)))==NULL) {
    w_abort("cannot allocate walign_done","");
  }
  walign_cnt[0]=0;

  if ((walign_done[1] = (int *)calloc(lcnt,sizeof(int)))==NULL) {
    w_abort("cannot allocate walign_done","");
  }
  walign_cnt[1]=0;

  /* was commented in for only FASTX/TFASTX, but do it always to
     simplify */
  aainit(pst.tr_type, pst.debug_lib);
  pst.maxlen = m_msg.maxn;

/*****************************************************************/
/* Main search loop, which calles do_work() repeatedly           */
/*****************************************************************/

  cur_n0 = 0;
  while (1)  {
/*
#ifdef DEBUG
#ifdef PVM_SRC
    fprintf(stderr," W: %d waiting MSEQTYPE\n",worker);
#endif
#ifdef MPI_SRC
    p4_dprintf(" W: %d waiting MSEQTYPE\n",worker);
#endif
#endif
*/

/*****************************************************************/
/* Wait for a query sequence from the manager                    */
/*****************************************************************/

#ifdef PVM_SRC
    bufid = pvm_recv(hosttid,MSEQTYPE);
    pvm_upkbyte((char *)&qm_msg,sizeof(qm_msg),1);
#endif
#ifdef MPI_SRC
    MPI_Recv(&qm_msg,sizeof(struct mngmsg),MPI_BYTE,hosttid,MSEQTYPE,
	     MPI_COMM_WORLD,&mpi_status);
#endif
#ifdef DEBUG
    fprintf(stderr,"[%d] have MSEQTYPE n0: %d s_func: %d slist: %d qf: %d\n",
	    worker,qm_msg.n0,qm_msg.s_func,qm_msg.slist,qm_msg.qshuffle);
#endif

/*****************************************************************/
/* New query sequence indicated by qm_msg.slist=0                */
/*****************************************************************/

    if (qm_msg.n0 > 0 && qm_msg.slist == 0) {

      if (cur_n0 > 0) {

/*****************************************************************/
/*    free everything associated with previous search            */
/*****************************************************************/

 	close_work (aa0[0], cur_n0, &pst, &f_str[0]);  
	free_ares(seqpt, 0, walign_done[0], walign_cnt[0], worker);
	walign_cnt[0] = 0;
	if (m_msg.q_ann_flg) free(m_msg.aa0a);


	if (m_msg.qframe == 2) {
	  close_work(aa0[1], cur_n0, &pst, &f_str[1]);
	  free_ares(seqpt, 1, walign_done[1], walign_cnt[1], worker);
	  walign_cnt[1] = 0;
	}
	if (old_shuffle) {
	  close_work(aa0s,cur_n0, &pst, &qf_str);
	  aa0s--;
	  free(aa0s);
	  old_shuffle = 0;
	}
	if (pst.pam_pssm) {
	  free_pam2p(pst.pam2p[0]);
	  free_pam2p(pst.pam2p[1]);
	}
      }

/*****************************************************************/
/*    Start allocating things for the next search                */
/*****************************************************************/

      pst.pam_pssm = qm_msg.pam_pssm;
      cur_n0 = qm_msg.n0;
      if (m_msg.q_ann_flg) {
	if ((m_msg.aa0a = calloc(qm_msg.n0+1,sizeof(char)))==NULL) {
	  w_abort(" cannot allocate aa0a");
	}
      }

/*****************************************************************/
/*    Get the next query sequence                                */
/*****************************************************************/

#ifdef PVM_SRC
      pvm_upkbyte((char *)aa0[0],qm_msg.n0+1+SEQ_PAD,1);
      if (m_msg.q_ann_flg) {
	pvm_upkbyte((char *)m_msg.aa0a,qm_msg.n0+1,1);
      }
#endif
#ifdef MPI_SRC
      MPI_Recv(aa0[0],qm_msg.n0+1+SEQ_PAD,MPI_BYTE,hosttid,
	       MSEQTYPE1,MPI_COMM_WORLD, &mpi_status);
      if (m_msg.q_ann_flg) {
	MPI_Recv(m_msg.aa0a,qm_msg.n0+1,MPI_BYTE,hosttid,
		 MSEQTYPE2,MPI_COMM_WORLD, &mpi_status);
      }
#endif

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

      /* This discovers most reasons for core dumps */
      if (pst.debug_lib)
	for (j=0; j<qm_msg.n0; j++)
	  if (aa0[0][j]>pst.nsq) {
	    fprintf(stderr,
		    "seq: %s residue[%d/%d] %d range (%d)\n",
		    qm_msg.libstr,j,qm_msg.n0,aa0[0][j],pst.nsq);
	    aa0[0][j]=0;
	    qm_msg.n0=j-1;
	    break;
	  }
#endif
      update_params(&qm_msg,&m_msg,&pst);
    }

/*****************************************************************/
/*    End of free()'s/ initialization for new sequence           */
/*****************************************************************/

#ifdef PVM_SRC
    pvm_freebuf(bufid);
#endif

    if (qm_msg.n0 == -1) {

/*****************************************************************/
/*   All done with searches                                      */
/*****************************************************************/
/*    printf(" %d: got n0 == -1\n",worker); */
      break;
    }

    /*    p4_dprintf(" W:%d n0:%d slist:%d s_func:%d (%d)\n",worker,qm_msg.n0,qm_msg.slist,qm_msg.s_func,qres_bufsize); */

/*****************************************************************/
/*   if qm_msg.slist > 0, search specific sequences, to be sent  */
/*****************************************************************/

    if (qm_msg.slist > 0) {  /* list search, not library search */
      if (liblist != NULL) free(liblist);

      /* get the list of sequences */
      if ((liblist=(struct stage2_str *)
	   calloc(qm_msg.slist,sizeof(struct stage2_str)))==NULL) {
	  sprintf(errstr,"sequence list %d",qm_msg.slist);
	  w_abort (errstr, "");
	}

#ifdef PVM_SRC
      bufid = pvm_recv(hosttid,LISTTYPE);
      pvm_upkbyte((char *)liblist,qm_msg.slist*sizeof(struct stage2_str),1);
      pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
      MPI_Recv(liblist,qm_msg.slist*sizeof(struct stage2_str),MPI_BYTE,
	       hosttid,LISTTYPE,MPI_COMM_WORLD, &mpi_status);
#endif
    }

/*****************************************************************/
/* have list of sequences to be compared/aligned                 */
/*****************************************************************/

    /* Initial stuff */
    if (qm_msg.slist == 0) {
/*****************************************************************/
/* New query - set up matrices and init_work()                   */
/*****************************************************************/
#ifdef DEBUG
/*
      fprintf(stderr,"n1: %d\t",qm_msg.n0);
      for (i=0; i<10; i++) fprintf(stderr,"%c",nt[aa0[0][i]]);
      fprintf(stderr,"\n");
*/
#endif
      if (pst.pam_pssm) {
	pst.pam2p[0] = alloc_pam2p(qm_msg.n0,nsq);
	pst.pam2p[1] = alloc_pam2p(qm_msg.n0,nsq);
      }

      init_work (aa0[0], qm_msg.n0, &pst, &f_str[0]);
      f_str[5]=f_str[4]=f_str[3]=f_str[2]=f_str[1]=f_str[0];

      if (pst.n1_low > 0 && pst.n1_high < BIGNUM) 
	sprintf(lib_range," (range: %d-%d)",pst.n1_low,pst.n1_high);
      else if (pst.n1_low > 0) 
	sprintf(lib_range," (range: >%d)",pst.n1_low);
      else if (pst.n1_high < BIGNUM)
	sprintf(lib_range," (range: <%d)",pst.n1_high);
      else
	lib_range[0]='\0';

      if (qm_msg.qshuffle) {
	if ((aa0s=(unsigned char *)malloc((qm_msg.n0+2)*sizeof (char)))==NULL)
	  w_abort ("Unable to allocate aa0s array - exiting!","");
	*aa0s='\0';
	aa0s++;

	memcpy(aa0s,aa0[0],qm_msg.n0+1);
	qshuffle(aa0s,qm_msg.n0,qm_msg.nm0);
#ifdef DEBUG
	fprintf(stderr,"[%d] shuffle: %d\n",worker,qm_msg.n0);
	fputs("   ",stderr);
	for (i=0; i<5; i++) {fprintf(stderr,"%c",pst.sq[aa0s[i]]);}
	fputc('\n',stderr);
#endif

	init_work (aa0s, qm_msg.n0, &pst, &qf_str);
	old_shuffle=1;
      }

      if (m_msg.qframe == 2) {
	memcpy(aa0[1],aa0[0],qm_msg.n0+1);
	revcomp(aa0[1],qm_msg.n0,&pst.c_nt[0]);
	init_work (aa0[1], qm_msg.n0, &pst, &f_str[1]);
      }
#ifdef DEBUG
/*
	fprintf(stderr,"[%d] init_work qf: %d nf: %d\n",worker,m_msg.qframe,m_msg.nframe);
*/
#endif
    }

/*****************************************************************/
/* Finished with initialization,                                 */
/* start doing comparisons or alignments                         */
/*****************************************************************/

    bestcnt = 0;

    ldb.length = 0l;
    ldb.carry = ldb.entries = 0;

    if (qm_msg.slist == 0) {	/* library search */

/*****************************************************************/
/* Start library search                                          */
/*****************************************************************/

      for (count=0; count < lcnt; count++) {

	if (seqpt[count].n1 < pst.n1_low ||
	    seqpt[count].n1 > pst.n1_high) continue;

	ldb.entries++;
	ldb.length += seqpt[count].n1;
	if (ldb.length > LONG_MAX) {
	  ldb.length -= LONG_MAX; ldb.carry++;
	}

	for (itt=m_msg.revcomp; itt<=m_msg.nitt1; itt++) {

	  rst.score[0] = rst.score[1] = rst.score[2] = -BIGNUM;

	  if (m_msg.self) {
	    lsn = lend + count;
	    if ((qm_msg.seqnm > lsn) && (((qm_msg.seqnm + lsn) % 2) != 0)) {
	      do_work (aa0[itt], qm_msg.n0,seqpt[count].aa1, seqpt[count].n1,
		       itt, &pst, f_str[itt], 0, &rst);
	    }
	    else if ((qm_msg.seqnm <= lsn) && (((qm_msg.seqnm+lsn)%2) == 0)) {
	      do_work (aa0[itt], qm_msg.n0, seqpt[count].aa1, seqpt[count].n1, 
		       itt, &pst, f_str[itt], 0, &rst);
	    }
	    else continue;
	  }
	  else {
	    do_work (aa0[itt], qm_msg.n0, seqpt[count].aa1, seqpt[count].n1,
		     itt, &pst, f_str[itt], 0, &rst);
	    if (qm_msg.qshuffle) {
	      do_work (aa0s, qm_msg.n0, seqpt[count].aa1, seqpt[count].n1,
		     itt, &pst, qf_str, 1, &qrst);
	    }
	  }
#ifdef DEBUG
/*
	  if (count < 10 || (count % 200 == 199)) {
	    fprintf(stderr,"[node %d] itt:%d/%d (%d) %3d %3d %3d - %d/%d\n",
		    worker,itt,m_msg.nitt1,count,
		    rst.score[0],rst.score[1],rst.score[2],
		    seqpt[count].m_seqnm,seqpt[count].n1);
	  }
*/
#endif
	  sw_score = -1;

	  bestr[bestcnt].seqnm  = count;
	  bestr[bestcnt].m_seqnm  = seqpt[count].m_seqnm;
	  memcpy(&(bestr[bestcnt].rst), &rst, sizeof(struct rstruct));
	  bestr[bestcnt].frame = itt;

	  bestr[bestcnt].qr_score = qrst.score[pst.score_ix];
	  bestr[bestcnt].qr_escore = qrst.escore;

	  if (pst.zsflag >= 10) {
	    if (pst.zs_win > 0) 
	      wshuffle(seqpt[count].aa1, aa1s,seqpt[count].n1,pst.zs_win,&ieven);
	    else 
	      shuffle(seqpt[count].aa1, aa1s,seqpt[count].n1);

	    do_work(aa0[itt],qm_msg.n0,aa1s,seqpt[count].n1,itt, &pst, 
		    f_str[itt], 0, &rst);
	    bestr[bestcnt].r_score = rst.score[pst.score_ix];
	  }

	  bestcnt++;
	  if (bestcnt >= qres_bufsize) {
#ifdef DEBUG
	    fprintf(stderr," worker: %d sending %d results\n",worker,qres_bufsize);
#endif
	    send_bestr(hosttid,curtype,bestr,qres_bufsize,bestcnt);
	    bestcnt = 0;
	  }
	}
      }	/* END - for count loop */
      send_bestr(hosttid, curtype, bestr,qres_bufsize, (bestcnt | FINISHED));
      send_struct(hosttid, LDBTYPE, &ldb, sizeof(ldb));
    }

/*****************************************************************/
/* End of library search section                                 */
/*****************************************************************/