p2_complib2.c

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

	    tm_escore = t_rescore;
	    t_rbest = -BIGNUM;
	    t_rescore = FLT_MAX;
	  }
	  stats[nstats].escore = tm_escore;
	  stats[nstats++].score = tm_best;
	  tm_escore = FLT_MAX;
	  t_best = -BIGNUM;
	}
      }
      else if (pst.zsflag >= 0) {	/* nstats >= MAX_STATS, zsflag >=0 */
	if (!stats_done ) {
	  pst.n0 = qm_msp0->n0;
	  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;
	  kstats = nstats;
	  for (i=0; i<nbest; i++) {
	    bptr[i]->zscore = (*find_zp)(bptr[i]->rst.score[pst.score_ix],
					 bptr[i]->rst.escore,bptr[i]->seq->n1,
					 bptr[i]->rst.comp, m_msp0->pstat_void);
	  }
	}
#ifdef SAMP_STATS
	if (!m_msp0->escore_flg) {
	  jstats = nrand(kstats++);
	  if (jstats < MAX_STATS) {
	    stats[jstats].n1 = t_n1;	/* save the best score */
	    stats[jstats].comp =  k_comp;
	    stats[jstats].H = k_H;
	    if (pst.zsflag >=10) t_best = t_rbest;
	    stats[jstats].score = t_best;
	  }
	}
#endif
      }

      if (stats_done) {
	zscore=(*find_zp)(i_score,e_score,seq_p->n1,k_comp,
			  m_msp0->pstat_void);
	if (bestr[ires].frame == m_msg0.nitt1) {
	  addhistz((*find_zp)(tm_best,tm_escore,t_n1,k_comp,
			      m_msp0->pstat_void),
		   &(m_msp0->hist));
	  t_best = t_rbest = -BIGNUM;
	}

      }
      else zscore = (double) i_score;

      if (zscore > zbestcut) {
	if (nbest>=MAX_BEST) {
	  selectbestz(bptr, nbest-MAX_BEST/4-1, nbest); 
	  nbest -= MAX_BEST/4;
	  zbestcut = bptr[nbest-1]->zscore;
	  best_flag = 0;
	}
	/* if zbestcut == -BIGNUM, bptr[] has not been reinitialized */
	else if (best_flag) bptr[nbest]=&best[nbest];

	bptr[nbest]->seq = seq_p;
	memcpy(&(bptr[nbest]->rst), &(bestr[ires].rst), sizeof(struct rstruct));
	bptr[nbest]->zscore = zscore;
	bptr[nbest]->frame = bestr[ires].frame;
	bptr[nbest]->seq->seqnm = bestr[ires].seqnm;

	nbest++;
      }
    }	/* for loop */
    if (naa0 < nnodes-FIRSTNODE) continue;

#ifdef DEBUG
    fprintf(stderr," All results returned\n");
#endif

    info_gstring2p[0][0]= info_gstring2p[1][0]='\0';

    /* get info_gstring2,3 - algorithm/parameter description */
#ifdef PVM_SRC
    bufid = pvm_recv(pinums[FIRSTNODE],PARAMTYPE);
    pvm_upkbyte(info_gstring2p[0],MAX_STR,1);
    pvm_upkbyte(info_gstring2p[1],MAX_STR,1);
    pvm_upkbyte(info_gstring3,sizeof(info_gstring3),1);
    pvm_upkbyte(lib_range,MAX_SSTR,1);
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv(info_gstring2p[0],MAX_STR,MPI_BYTE,FIRSTNODE,PARAMTYPE,
	     MPI_COMM_WORLD,&mpi_status);
    MPI_Recv(info_gstring2p[1],MAX_STR,MPI_BYTE,FIRSTNODE,PARAMTYPE,
	     MPI_COMM_WORLD,&mpi_status);
    MPI_Recv(info_gstring3,sizeof(info_gstring3),MPI_BYTE,FIRSTNODE,PARAMTYPE,
	     MPI_COMM_WORLD,&mpi_status);
    MPI_Recv(lib_range,MAX_SSTR,MPI_BYTE,FIRSTNODE,PARAMTYPE,
	     MPI_COMM_WORLD,&mpi_status);
#endif

/* ********************** */
/* analyze the results    */
/* ********************** */
    
    if (!stats_done) {
      if (nbest < 20 || pst.zsflag < 0) {
	pst.zsflag_f = -1;
      }
      else {
	pst.n0 = qm_msp0->n0;
	pst.zsflag_f = process_hist(stats,nstats,*m_msp0,&pst,
				    &m_msp0->hist, &m_msp0->pstat_void,stats_done);

	if (pst.zdb_size <= 1) pst.zdb_size = ntt.entries;
	for (i=0; i<nbest; i++)
	  bptr[i]->zscore = (*find_zp)(bptr[i]->rst.score[pst.score_ix],
				       bptr[i]->rst.escore, bptr[i]->seq->n1,
				       bptr[i]->rst.comp, m_msp0->pstat_void);
      }
    }

    m_msp0->db.entries = ntt.entries;
    m_msp0->db.length = ntt.length;
    m_msp0->db.carry = ntt.carry;

    if (pst.zdb_size < 1) pst.zdb_size = ntt.entries;

    if (!qm_msp0->qshuffle) {
      last_stats(aa0p0, m_msp0->n0,
		 stats,nstats, bptr,nbest, *m_msp0, pst, 
		 &m_msp0->hist, &m_msp0->pstat_void);
    }
    else {
      last_stats(aa0p0, m_msp0->n0,
		 qstats,nqstats, bptr,nbest, *m_msp0, pst, 
		 &m_msp0->hist, &m_msp0->pstat_void);
    }

    if (m_msp0->last_calc_flg) {
      nbest = last_calc(bptr, nbest, *m_msp0, &pst, qm_msp0,
			m_msp0->pstat_void);
    }

    sortbeste(bptr,nbest);
    scale_scores(bptr,nbest,m_msp0->db,&pst,m_msp0->pstat_void);

    if (pst.zsflag >= 0 && bptr[0]->rst.escore >= m_msg0.e_cut) goto no_results;

    /*      else sortorder(bptr,nbest,wlsn,nnodes); */

/* if more than one stage or markx==9, calculate opt scores or do alignment */
/* send results to workers as available */

    if (m_msg0.stages > 1 || m_msg0.markx & MX_M9SUMM) {

      /* to determine how many sequences to re-align (either for
	 do_opt() or calc_id() we need to modify m_msg.mshow to get
	 the correct number of alignments */

      if (m_msg0.mshow_flg != 1 && pst.zsflag >= 0) {
	for (i=0; i<nbest && bptr[i]->rst.escore< m_msg0.e_cut; i++) {}
	m_msg0.mshow = i;
      }

      /* allocate space for a_struct info */
      if (m_msg0.markx & MX_M9SUMM && m_msg0.mshow > 0) {
	if ((aln_d_base=(struct a_struct *)
	     calloc((size_t)m_msg0.mshow,sizeof(struct a_struct)))==NULL) {
	  fprintf(stderr," cannot allocate a_struct %d\n", m_msg0.mshow);
	  exit(1);
	}

	for (is = 0; is < m_msg0.mshow; is++ ) {
	  memcpy(&(bptr[is]->aln_d),&aln_d_base[is],sizeof(struct a_struct));
	}
      }

      do_stage2(bptr,m_msg0.mshow, *m_msp0, DO_OPT_FLG, qm_msp0);
    }

  no_results:
    tdone = s_time();
    tddone = time(NULL);

    /* changed from >> to >>> because qm_msp0->libstr is missing '>' */
    fprintf (outfd, "%3d>>>%s\n", qlib,qm_msp0->libstr);

    /* make certain that m_msp0->n0, libstr are current */
    m_msp0->n0 = qm_msp0->n0;
    /*    strncpy(m_msp0->libstr,qm_msp0->libstr,sizeof(m_msg0.libstr)); */

    prhist (outfd,*m_msp0,&pst,m_msp0->hist,nstats,m_msp0->ldb,
	    lib_range,info_gstring2p,info_hstring_p);

    if (bptr[0]->rst.escore < m_msg0.e_cut) {

      showbest (outfd, bptr, nbest, qlib, m_msp0,&pst,ntt,info_gstring2p);

      if (m_msg0.markx & (MX_AMAP+ MX_HTML + MX_M9SUMM) && !(m_msg0.markx & MX_M10FORM)) {
	fprintf(outfd,"\n>>>%s#%d %s%s, %d %s vs %s library\n",
		m_msg0.tname,qlib,qm_msp0->libstr,
		(m_msg0.revcomp ? "-":"\0"), qm_msp0->n0, m_msg0.sqnam,
		m_msg0.lname);
      }

      if (m_msg0.markx & MX_M10FORM) {
	if ((bp=strchr(qm_msp0->libstr,' '))!=NULL) *bp = '\0';
	fprintf(outfd,"\n>>>%s#%d %s%s, %d %s vs %s library\n",
		m_msg0.tname,qlib,qm_msp0->libstr,
		(m_msg0.revcomp ? "-":"\0"), qm_msp0->n0, m_msg0.sqnam,
		m_msg0.lname);
	if (bp!=NULL) *bp=' ';
	fprintf(outfd,"; pg_name: %s\n",argv[0]);
	fprintf(outfd,"; pg_ver: %s\n",mp_verstr);
	fprintf(outfd,"; pg_argv:");
	for (i=0; i<argc; i++)
	  fprintf(outfd," %s",argv[i]);
	fputc('\n',outfd);
	fputs(info_gstring3,outfd);
	fputs(info_hstring_p[0],outfd);
	fputs(info_hstring_p[1],outfd);
      }

      /* ashow is -1 if not set, -d 0 indicates no alignments, > 0 if set */
      /* if ashow is -1, m_msg.nshow (set by e_cut above) sets limit
	 in showalign */
      
      if (m_msp0->ashow != 0) {
	/* showalign needs m_msp->qtitle, so fill it in */
	strncpy(m_msp0->qtitle,qm_msp0->libstr,MAX_FN-1);
	m_msp0->qtitle[MAX_FN-1]='\0';
	showalign (outfd, bptr, nbest, qlib, *m_msp0, &pst, info_gstring2p);
      }
    }
    else {
      if (m_msg0.markx & (MX_M9SUMM + MX_M10FORM)) {
	fprintf(outfd,"\n>>>%s#%d %s%s, %d %s vs %s library\n",
		m_msg0.tname,qlib,qm_msp0->libstr,(m_msg0.revcomp ? "-":"\0"), qm_msg0.n0, m_msg0.sqnam,
		m_msg0.lname);
	fprintf(outfd,">>>!!! No sequences with E() < %f\n",m_msg0.e_cut);
      }
      else fprintf(outfd,"!! No sequences with E() < %f\n",m_msg0.e_cut);
    }

    if (! (m_msg0.markx & (MX_M9SUMM + MX_M10FORM))) {
      fprintf(outfd,"/** search time: ");
      ptime(outfd,tdone-tprev);
      fprintf(outfd," **/\n");
      tprev = tdone;
    }
    else if (m_msg0.markx & MX_M9SUMM) {
      if (aln_d_base != NULL) {
	free((void *)aln_d_base);
	aln_d_base = NULL;
      }
      fprintf(outfd,">>>***\n");
      fprintf(outfd,"/** %s **/\n",info_gstring2p[0]);
      fprintf(outfd,"/** %s **/\n",info_gstring2p[1]);
      fprintf(outfd,"/** %s **/\n",m_msp0->hist.stat_info);
      fprintf(outfd,">>><<<\n");
    }
    else if (m_msg0.markx & MX_M10FORM) {
      fprintf(outfd,">>><<<\n");
    }
    fflush(outfd);
    
/* *********************** */
/* end of analysis/display */
/* *********************** */


/* *********************** */
/* start the next search   */                                           
/* *********************** */

    if (fdata) {		/* label the results file */
      fprintf(fdata,"/** %s **/\n",info_gstring2p[0]);
      fprintf(fdata,"/** %s **/\n",info_gstring2p[1]);
      fprintf(fdata,"%3d>%-50s\n",qlib-1,qm_msp1->libstr);
      fflush(fdata);
    }
    
    if (m_msp1->escore_flg) {	/* re-initialize some stats stuff before search */
      pst.zsflag_f = process_hist(stats,nstats,*m_msp1,&pst,
				  &m_msp1->hist,&m_msp1->pstat_void,0);
      stats_done=1;
      if (pst.zdb_size <= 1) pst.zdb_size = ntt.entries;
    }
    else stats_done = 0;

    /* set up qstats if necessary - different queries have different qshuffle */
    if (m_msp1->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 = nstats = 0;

    /* send new qm_msp, sequence */
    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);
	}	  
      }
      for (node = FIRSTNODE; node < nnodes; node++)
	pvm_send(pinums[node],MSEQTYPE);
#endif
#ifdef MPI_SRC
      for (node=FIRSTNODE; node < nnodes; node++) {
	MPI_Send(qm_msp1,sizeof(qm_msg1),MPI_BYTE,node,MSEQTYPE,
		 MPI_COMM_WORLD);
	if (qm_msp1->n0 != -1) {
	  MPI_Send(aa0p1,qm_msp1->n0+1+SEQ_PAD,MPI_BYTE,node,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
    }

    qlib++;
    if (qm_msp1->n0 != -1) {
      qtt.entries++;
      qtt.length += qm_msp1->n0;
    }
    else goto done;
    
/* ******************************** */
/* flip m_msg, qm_msg, aa0 pointers */
/* ******************************** */

    naa0 = 0;
    best_flag = 1;
    nbest = nopt = 0;
    zbestcut = -BIGNUM;
    m_msp0->ldb.entries = m_msp0->ldb.carry = 0;
    m_msp0->ldb.length = 0l;

    if (curtype == ONETYPE) {
      curtype = TWOTYPE;
      qm_msp0 = &qm_msg1;
      qm_msp1 = &qm_msg0;
      m_msp0 = &m_msg1;
      m_msp1 = &m_msg0;
      aa0p0 = aa01;
      aa0p1 = aa00;
    }
    else  {
      curtype = ONETYPE;
      qm_msp0 = &qm_msg0;
      qm_msp1 = &qm_msg1;
      m_msp0 = &m_msg0;
      m_msp1 = &m_msg1;
      aa0p0 = aa00;
      aa0p1 = aa01;
    }


/* **********************************************************/
/* all library sequences are done get next library sequence */
/* **********************************************************/

    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 ((qlib+1) >= m_msg0.ql_stop) { qm_msp1->n0 = m_msp1->n0 = -1;}

    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_msg0.n0 > 0) {
       for (i=0; i<SEQ_PAD+1; i++) {aa00[m_msg0.n0+i]=0;}
    }

    qm_msp1->slist = 0;
    /*
    leng = strlen (qm_msp1->libstr);
    sprintf (&(qm_msp1->libstr[leng]), " %d %s", qm_msp1->n0, q_sqnam);
    */
    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';

    qm_msp1->seqnm = qlib;

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

  }	    /* while loop */
  
  /* ************************** */
  /* end of query library while */
  /* ************************** */

 done:
  tdone = s_time();
  if (m_msg0.markx & (MX_M9SUMM + MX_M10FORM)) fputs(">>>///\n",outfd);
  printsum(outfd);
  if (outfd!=stdout) printsum(stdout);
  printsum(stderr);
#ifdef PVM_SRC
  pvm_exit();
#endif
#ifdef MPI_SRC
  MPI_Finalize();
#endif

  exit(0);
}   /* End of main program */

void
printsum(FILE *fd)
{
  double db_tt;
  char tstr1[26], tstr2[26];

  strncpy(tstr1,ctime(&tdstart),sizeof(tstr1));
  strncpy(tstr2,ctime(&tddone),sizeof(tstr1));
  tstr1[24]=tstr2[24]='\0';

  /* Print timing to output file as well */
  if (qtt.carry==0) {
    fprintf(fd, "\n%ld residues in %d query   sequences\n", qtt.length, qtt.entries);
  }
  else {
    db_tt = (double)qtt.carry*(double)LONG_MAX + (double)qtt.length;
    fprintf(fd, "\n%.0g residues in %d query   sequences\n", db_tt, qtt.entries);
  }

  if (ntt.carry==0) {
    fprintf(fd, "%ld residues in %ld library sequences\n", ntt.length, ntt.entries);
  }
  else {
    db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
    fprintf(fd, "%.6f residues in %ld library sequences\n", db_tt, ntt.entries);
  }

  fprintf(fd," %d processors (%d workers) were used\n",
	  nnodes+-FIRSTNODE+1,nnodes-FIRSTNODE);
  fprintf(fd," Pvcomplib [%s]\n start: %s done: %s\n",mp_verstr,tstr1,tstr2);
  fprintf(fd," Loading time: ");
  ptime(fd, tscan - tstart);
  fprintf (fd," Scan time: ");
  ptime (fd, tdone - tscan);
  fprintf (fd,"\n");
  fprintf (fd, "\nFunction used was %s [%s]\n", prog_func,verstr);
}

void fsigint()
{
  int i;

  tdone = s_time();
  tddone = time(NULL);

  if (outfd!=stdout) fprintf(outfd,"/*** interrupted ***/\n");
  fprintf(stderr,"/*** interrupted ***/\n");

  printsum(stdout);
  if (outfd!=stdout) printsum(outfd);

#ifdef PVM_SRC
  for (i=FIRSTNODE; i<nnodes; i++) pvm_kill(pinums[i]);
  pvm_exit();
#endif
#ifdef MPI_SRC
  MPI_Abort(MPI_COMM_WORLD,1);
  MPI_Finalize();
#endif
  exit(1);
}