p2_complib2.c

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

	    (seq_p->sfnum[1]=sfnum[1])>0 &&
	    (seq_p->sfnum[2]=sfnum[2])>0 &&
	    (seq_p->sfnum[3]=sfnum[3])>0 &&
	    (seq_p->sfnum[4]=sfnum[4])>0 &&
	    (seq_p->sfnum[5]=sfnum[5])>0 &&
	    (seq_p->sfnum[6]=sfnum[6])>0 &&
	    (seq_p->sfnum[7]=sfnum[7])>0 &&
	    (seq_p->sfnum[8]=sfnum[8])>0 &&
	    (seq_p->sfnum[9]=sfnum[9])>0) ;
#endif
	m_seqnm++;
	nseq++;

	if (m_seqnm >= max_sql) {
	  max_sql += MAXSQL;
	  if ((ldes=(struct seq_record *)realloc(ldes,max_sql*sizeof(struct seq_record)))
	      ==NULL) {
	    fprintf(stderr," failure to realloc ldes(%d) %ld\n",
		    max_sql,max_sql*sizeof(struct seq_record));
	    s_abort("cannot allocate ldes","");
	    exit(1);
	  }
	}

	/* increment ptrs */
	aa1prev = aa1;

	aa1 += n1+1+SEQ_PAD;
	ntbuff += n1+1+SEQ_PAD;

	/* if the buffer is filled */
	if (nseq >= max_buf_cnt || ntbuff >= m_msg0.pbuf_siz - m_msg0.maxn) {
	  /* provide filled buffer to workers */
#ifdef PVM_SRC
	  pvm_initsend(PvmDataRaw);
	  pvm_pkint(&nseq,1,1);
	  pvm_pkint(&ntbuff,1,1);
	  pvm_pkint(n1_arr,nseq,1);
	  pvm_pkint(aa1i_arr,nseq,1);
	  pvm_pkint(m_seqnm_arr,nseq,1);
	  pvm_send(pinums[node],STARTTYPE4);

	  pvm_initsend(PvmDataRaw);
	  pvm_pkbyte((char *)aa1_buf,ntbuff,1);
	  pvm_send(pinums[node],STARTTYPE5);
#endif
#ifdef MPI_SRC
	  MPI_Send(&nseq,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
	  MPI_Send(&ntbuff,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
	  MPI_Send(n1_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
	  MPI_Send(aa1i_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
	  MPI_Send(m_seqnm_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);

	  MPI_Send(aa1_buf,ntbuff,MPI_BYTE,node,STARTTYPE5,MPI_COMM_WORLD);
#endif
	  nseq =  0;

	  aa1 = aa1_buf;
	  ntbuff = 0;
	  if (++node >= nnodes) node = FIRSTNODE;
	}

      loop1:
	if (lcont) {
	  memcpy(aa1,&aa1prev[n1-m_msg0.l_overlap],m_msg0.l_overlap);
	  aa1ptr = &aa1[m_msg0.l_overlap];
	  ocont = lcont;
	  maxt = m_msg0.maxt3;
	  loffset += n1 - m_msg0.l_overlap;
	}
	else {
	  if (ocont) *n1tot_cur = n1tot_v;
	  n1tot_v = 0;
	  n1tot_cur = NULL;

	  ocont = 0;
	  aa1ptr = aa1;
	  maxt = m_msg0.maxn;
	  loffset = 0l;
	}
      }
    }
  }	/* for (iln < nln) */

  if (nseq > 0) {
#ifdef PVM_SRC
    pvm_initsend(PvmDataRaw);
    pvm_pkint(&nseq,1,1);
    pvm_pkint(&ntbuff,1,1);
    pvm_pkint(n1_arr,nseq,1);
    pvm_pkint(aa1i_arr,nseq,1);
    pvm_pkint(m_seqnm_arr,nseq,1);
    pvm_send(pinums[node],STARTTYPE4);

    pvm_initsend(PvmDataRaw);
    pvm_pkbyte((char *)aa1_buf,ntbuff,1);
    pvm_send(pinums[node],STARTTYPE5);
#endif
#ifdef MPI_SRC
    MPI_Send(&nseq,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(&ntbuff,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(n1_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(aa1i_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(m_seqnm_arr,nseq,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);

    MPI_Send(aa1_buf,ntbuff,MPI_BYTE,node,STARTTYPE5,MPI_COMM_WORLD);
#endif
  }

  /*   fprintf(stderr," all sequences sent\n"); */

  if (ntt.entries <= 0) {
    s_abort("no reference library sequences found\n","");
  }

  zero = 0;
  for (node=FIRSTNODE; node < nnodes; node++) {
#ifdef PVM_SRC
    pvm_initsend(PvmDataRaw);
    pvm_pkint(&zero,1,1);
    pvm_pkint(&zero,1,1);
    pvm_pkint(n1_arr,1,1);
    pvm_pkint(aa1i_arr,1,1);
    pvm_pkint(m_seqnm_arr,1,1);
    pvm_send(pinums[node],STARTTYPE4);
#endif
#ifdef MPI_SRC
    MPI_Send(&zero,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(&zero,1,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(n1_arr,0,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(aa1i_arr,0,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
    MPI_Send(m_seqnm_arr,0,MPI_INT,node,STARTTYPE4, MPI_COMM_WORLD);
#endif
  }

  for (node = FIRSTNODE; node < nnodes; node++) {
#ifdef PVM_SRC
    bufid = pvm_recv(-1,STARTTYPE0);
    pvm_bufinfo(bufid,NULL,NULL,&tid);
    snode = tidtonode(tid);
    pvm_upkint(&lcnt,1,1);
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv(&lcnt,1,MPI_INT,MPI_ANY_SOURCE,STARTTYPE0,
	     MPI_COMM_WORLD,&mpi_status);
    snode= mpi_status.MPI_SOURCE;
#endif
    wlsn [snode-FIRSTNODE] = lcnt;
    fprintf(stderr," %d sequences at %d\n",lcnt,snode);
  }

  /* print out all descriptions */
  /*
  for (node = FIRSTNODE; node < nnodes; node++)
    for (lcnt = 0; lcnt < wlsn[node-FIRSTNODE]; lcnt ++)
      printf("%2d:%3d\t%s\n",node,lcnt,ldes[lcnt].bline);
  */

  /* Calculate cumulative totals and send to workers for a self search */

  clsn [0] = nclib= 0;
  for (node = FIRSTNODE; node < nnodes-1; node++) {
    /* clsn[] is for the next node */
    clsn[node-FIRSTNODE+1] = nclib += wlsn[node-FIRSTNODE];
  }

  if (m_msg0.self)
    for (node = FIRSTNODE; node < nnodes; node++) {
#ifdef PVM_SRC
      pvm_initsend(PvmDataRaw);
      pvm_pkint(&clsn[node-FIRSTNODE],1,1);
      pvm_send(pinums[node],STARTTYPE1);
#endif
#ifdef MPI_SRC
    MPI_Send(&clsn[node-FIRSTNODE],1,MPI_INT,node,STARTTYPE1,MPI_COMM_WORLD);
#endif
      fprintf(stderr,"sending lend: %d to worker %d\n",clsn[node-FIRSTNODE],node);
    }

  m_msg0.db.entries = m_msg1.db.entries = ntt.entries;
  m_msg0.db.length = m_msg1.db.length = ntt.length;
  m_msg0.db.carry = m_msg1.db.carry = ntt.carry;

  m_msg0.ldb.length = m_msg1.ldb.length = 0l;
  m_msg0.ldb.carry = m_msg1.ldb.carry = 0;
  m_msg0.ldb.entries = m_msg1.db.entries = 0;

  last_msg_b[0] = m_msg0.nbr_seq = m_msg1.nbr_seq = ntt.entries;

  qres_bufsize = BFR;
  /* if BFR is too big for this library, reduce it */
  while ( ntt.entries*(m_msg0.nitt1+1)/(2*nnodes) < qres_bufsize) {
    qres_bufsize /= 2;
    if ((qres_bufsize%(m_msg0.nitt1+1))!= 0) {
      qres_bufsize *= (m_msg0.nitt1+1);
      break;
    }
    if (qres_bufsize < 50) break;
  }
  last_msg_b[1] = qres_bufsize;

  fprintf(stderr," using BFR=%d/%d\n",qres_bufsize,BFR);

#ifdef PVM_SRC
  pvm_initsend(PvmDataRaw);
  pvm_pkint(last_msg_b,2,1);
  for (node=FIRSTNODE; node < nnodes; node++)
    pvm_send(pinums[node],STARTTYPE0);
#endif
#ifdef MPI_SRC
  for (node=FIRSTNODE; node < nnodes; node++)
    MPI_Send(last_msg_b,2,MPI_INT,node,STARTTYPE0,MPI_COMM_WORLD);
#endif  

  tscan = tprev = s_time();

/**************************************
  The logic of this section has been simplified to allow multistage
  comparison functions to be used and alignments to be generated.

  	send 1st query to workers
	get next query sequence from host (m_msp1)
    L1:	get results from next-1 search (m_msp0)
	sort the results of the next-1 search
	(possibly) do additional stages of search
	(possibly produce alignments for search
	send next query to workers (m_msp1)
	display result of next-1 search (m_msp0)
	get next query sequence from host (m_msp1)
	goto L1;

As a result of the interleaving, there must be two qm_msg structures,
one for the next-1 sequence (which is required for labeling the
output), and one for the next sequence (which is sent to the workers
while the results are being displayed.  qm_msp0 and qm_msp1 alternate
between these two structures.
***************************************/

/*
  qm_msp0 points to the older qm_msg 
  qm_msp1 points to the newer qm_msg
  the assignment below goes with curtype==ONETYPE
*/
  m_msp0 = &m_msg0;
  m_msp1 = &m_msg1;

  qm_msp0 = &qm_msg0;
  qm_msp1 = &qm_msg1;
  
  aa0p0 = aa00;		/* aa0p0 is the "old" sequence */
  aa0p1 = aa01;		/* aa0p1 is the "new" sequence */

  last_params(aa00,m_msp0->n0,m_msp0,&pst,qm_msp0);

  /* process_hist() is called here to get find_zp(), and some other
     structures initialized that would otherwise not be initialized
     because z-scores are not being calculated */

  if (m_msp0->escore_flg) {
     pst.zsflag_f = process_hist(stats,nstats,*m_msp0,&pst,
				 &m_msp0->hist,&m_msp0->pstat_void,0);
     if (pst.zdb_size <= 1) pst.zdb_size = ntt.entries;
     stats_done=1;
  }

  if (m_msp0->qshuffle && qstats==NULL) {
    if ((qstats =
	 (struct stat_str *)calloc(m_msg0.shuff_max+1,sizeof(struct stat_str)))==NULL)
      s_abort ("Cannot allocate qstats struct","");
  }
  nqstats = 0;

/* Send first query sequence to each worker */

  if (m_msg0.dfile[0] && (fdata=fopen(m_msg0.dfile,"w"))!=NULL)
    fprintf(fdata,"%3d>%-50s\n",qlib,qm_msp0->libstr);

#ifdef PVM_SRC
  pvm_initsend(PvmDataRaw);
  pvm_pkbyte((char *)qm_msp0,sizeof(qm_msg0),1);
  if (qm_msp0->n0 > 0) {
    pvm_pkbyte((char *)aa0p0,qm_msp0->n0+1+SEQ_PAD,1);
    if (m_msg0.q_ann_flg) pvm_pkbyte((char *)m_msp0->aa0a,qm_msp0->n0+1,1);
  }
  for (node = FIRSTNODE; node < nnodes; node++)
    pvm_send(pinums[node],MSEQTYPE);
#endif
#ifdef MPI_SRC
  for (node = FIRSTNODE; node < nnodes; node++) {
    MPI_Send(qm_msp0,sizeof(qm_msg0),MPI_BYTE,node,MSEQTYPE,MPI_COMM_WORLD);
    if (qm_msp0->n0 > 0) {
      MPI_Send(aa0p0,qm_msp0->n0+1+SEQ_PAD,MPI_BYTE,node,
	       MSEQTYPE1,MPI_COMM_WORLD);
      if (m_msg0.q_ann_flg) {
	if (m_msp0->aa0a == NULL) {
	  fprintf(stderr," m_msp0: %o/%oaa0a is null\n",m_msp0,m_msp0->aa0a);
	}
	MPI_Send(m_msp0->aa0a,qm_msp0->n0+1,MPI_BYTE,node, MSEQTYPE2,MPI_COMM_WORLD);
      }
    }
  }
#endif

  /* Get second query sequence (additional query sequences are read in
     the main loop */

  m_msp1->n0 = qm_msp1->n0 = 
    QGETLIB(aa0p1,MAXTST,q_bline, sizeof(q_bline),&qseek, &qlcont,q_file_p,&m_msp1->sq0off);
  strncpy(qm_msp1->libstr,q_bline,sizeof(qm_msg0.libstr)-20);
  qm_msp1->libstr[sizeof(qm_msg0.libstr)-21]='\0';
  if ((bp=strchr(qm_msp1->libstr,' '))!=NULL) *bp='\0';

  /* if annotations are included in sequence, remove them */
  if (m_msg0.q_ann_flg) {
    m_msp1->n0 = qm_msp1->n0 = 
      ann_scan(aa0p1,qm_msp1->n0,m_msp1,m_msp1->qdnaseq);
#ifdef DEBUG
    fprintf(stderr,"m_msp1->/aa0a is: %o/%o\n",m_msp1,m_msp1->aa0a);
#endif
  }

  if (qm_msp1->n0 > 0 && m_msg0.term_code && !qlcont && 
      m_msg0.qdnaseq == SEQT_PROT &&
      aa0p1[m_msp1->n0-1]!=m_msg0.term_code) {
    aa0p1[m_msp1->n0++]=m_msg0.term_code;
    aa0p1[m_msp1->n0]=0;
    qm_msp1->n0 = m_msp1->n0;
  }

  /* for ALTIVEC, must pad with 15 NULL's */
  if (m_msp1->n0 > 0) {
    for (i=0; i<SEQ_PAD+1; i++) {aa0p1[m_msp1->n0+i]=0;}
  }

  qm_msp1->slist = 0;
  qm_msp1->seqnm = qlib;

  last_params(aa0p1,m_msp1->n0,m_msp1,&pst,qm_msp1);

  sprintf(tmp_str," - %d %s", qm_msp1->n0, q_sqnam);
  if (strlen(qm_msp1->libstr) + strlen(tmp_str) >= sizeof(qm_msg0.libstr))
    qm_msp1->libstr[sizeof(qm_msg0.libstr)-strlen(tmp_str)-2] = '\0';
  strncat(qm_msp1->libstr,tmp_str,
	  sizeof(qm_msg0.libstr)-strlen(qm_msp1->libstr)-1);
  qm_msp1->libstr[sizeof(qm_msg0.libstr)-1]='\0';

  naa0 = 0;  /* reset node counter */

  /* sit in loop and collect results */
  nbest = nopt = 0;
  zbestcut = -BIGNUM;

  while (1) {

#ifdef PVM_SRC
    bufid = pvm_recv(-1,curtype);
    pvm_bufinfo(bufid,NULL,NULL,&tid);
    pvm_upkbyte((char *)&bestr[0],sizeof(struct comstr)*(qres_bufsize+1),1);
    snode = tidtonode(tid);
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv(bestr,sizeof(struct comstr)*(qres_bufsize+1),
	     MPI_BYTE,MPI_ANY_SOURCE,curtype,MPI_COMM_WORLD,&mpi_status);
    snode = mpi_status.MPI_SOURCE;
#endif

    nres = bestr[qres_bufsize].seqnm & ~FINISHED;

#ifdef DEBUG
    fprintf(stderr,"%d results from %d\n",nres,snode);
#endif

    if (bestr[qres_bufsize].seqnm&FINISHED) {	/* a worker is finished */

      /* get the ldb struct information */
#ifdef PVM_SRC
      bufid = pvm_recv(-1,LDBTYPE);
      pvm_bufinfo(bufid,NULL,NULL,&tid);
      pvm_upkbyte((char *)&lldb,sizeof(struct db_str),1);
      pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
      MPI_Recv(&lldb,sizeof(struct db_str),
	     MPI_BYTE,MPI_ANY_SOURCE,LDBTYPE,MPI_COMM_WORLD,&mpi_status);
      snode = mpi_status.MPI_SOURCE;
#endif
 /*       fprintf(stderr," node %d ldb: %d %ld [%d]\n",snode, lldb.entries, lldb.length, lldb.carry); */

      /* add it to the current ldb */
      m_msp0->ldb.entries += lldb.entries;
      m_msp0->ldb.carry += lldb.carry;
      m_msp0->ldb.length += lldb.length;
      if (lldb.length > LONG_MAX) {m_msp0->ldb.length -= LONG_MAX; m_msp0->ldb.carry++;}

      naa0++;

      /* fast_flag == 1 => send new sequences immediately */
      fast_flag = ((m_msp0->stages==1) && !(m_msp0->markx & MX_M9SUMM) &&
		   (m_msp0->ashow == 0) && (m_msp0->last_calc_flg==0));
      /* send a new query sequence if no more processing required */
      if (fast_flag) {
#ifdef PVM_SRC
	pvm_initsend(PvmDataRaw);
	pvm_pkbyte((char *)qm_msp1,sizeof(qm_msg1),1);
	if (qm_msp1->n0 != -1) {
	  pvm_pkbyte((char *)aa0p1,qm_msp1->n0+1+SEQ_PAD,1);
	  if (m_msp1->q_ann_flg) pvm_pkbyte((char *)m_msp1->aa0a,qm_msp1->n0+1,1);
	}
	pvm_send(tid,MSEQTYPE);
#endif
#ifdef MPI_SRC
	MPI_Send(qm_msp1,sizeof(qm_msg1),MPI_BYTE,snode,MSEQTYPE,MPI_COMM_WORLD);
	if (qm_msp1->n0 != -1) {
	  MPI_Send(aa0p1,qm_msp1->n0+1+SEQ_PAD,MPI_BYTE,snode,MSEQTYPE1,MPI_COMM_WORLD);
	  if (m_msp1->q_ann_flg)
	    MPI_Send(m_msp1->aa0a,qm_msp1->n0+1,MPI_BYTE,snode,MSEQTYPE2,MPI_COMM_WORLD);
	}
#endif
      }
    }

#ifdef DEBUG
    if (pst.debug_lib)
      fprintf(stderr," unpacking %d from %d; nbest %d\n",nres,snode,nbest);
#endif

    /* this section is now more complex because can get groups of
       sequence results; e.g. forward and reverse frame */

    t_best = t_rbest = t_qrbest = -BIGNUM;
    tm_escore = t_rescore = t_qrescore = FLT_MAX;
    for (ires = 0; ires < nres; ires++) {
      seq_p = &ldes[bestr[ires].m_seqnm];

      /* save raw results */
      if (fdata) {
	strncpy(tlibstr,seq_p->bline,10);
	if ((bp=strchr(tlibstr,' '))!=NULL) *bp='\0';
	fprintf(fdata,"%-10s\t%4d\t%d\t%4d\t%4d\t%4d\t%8ld\n",
		tlibstr,seq_p->n1,bestr[ires].frame,
		bestr[ires].rst.score[0],bestr[ires].rst.score[1],bestr[ires].rst.score[2],
		seq_p->lseek);
      }

      i_score = bestr[ires].rst.score[pst.score_ix];
      e_score = bestr[ires].rst.escore;
      k_comp = bestr[ires].rst.comp;
      k_H = bestr[ires].rst.H;

      t_n1 = seq_p->n1;
      if (i_score > t_best) {tm_best = t_best = i_score;}
      if (e_score < tm_escore) tm_escore = e_score;

      if (m_msp0->qshuffle) {
	if (bestr[ires].qr_score > t_qrbest) 
	  t_qrbest = bestr[ires].qr_score;
	if (bestr[ires].qr_escore < t_qrescore) 
	  t_qrescore = bestr[ires].qr_escore;

	if (bestr[ires].frame==m_msp0->nitt1 && 
	    nqstats < m_msp0->shuff_max ) {
	  qstats[nqstats].n1 = t_n1;	/* save the best score */
	  qstats[nqstats].comp =  k_comp;
	  qstats[nqstats].H = k_H;
	  qstats[nqstats].escore = t_qrescore;
	  qstats[nqstats++].score = t_qrbest;
	  t_qrbest = -BIGNUM;	/* reset t_qrbest, t_qrescore */
	  t_qrescore = FLT_MAX;
	}
      }

      if (pst.zsflag >= 10 && bestr[ires].r_score > t_rbest) {
	t_rbest = bestr[ires].r_score;
	t_rescore = bestr[ires].r_escore;
      }

      if (nstats < MAX_STATS) {
	if (bestr[ires].frame == m_msg0.nitt1) {
	  stats[nstats].n1 = t_n1;
	  stats[nstats].comp = k_comp;
	  stats[nstats].H = k_H;

	  if (pst.zsflag > 10) {
	    tm_best = t_rbest;