comp_lib2.c

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

    if (m_msg.qframe == 2) free(aa0[1]-1);

    if (have_f_str) {
      if (f_str[1]!=f_str[0]) {
	close_work (aa0[1], m_msg.n0, &pst, &f_str[1]);
      }
      close_work (aa0[0], m_msg.n0, &pst, &f_str[0]);
      have_f_str = 0;
#ifndef COMP_THR
      if (m_msg.qshuffle) close_work (aa0s, m_msg.n0, &pst, &qf_str);
#endif
      if (pst.pam_pssm) {
	free_pam2p(pst.pam2p[0]);
	free_pam2p(pst.pam2p[1]);
      }
    }

    if (aa1shuff_b != NULL) {
      free(aa1shuff_b);
      aa1shuff_b = NULL;
    }

    for (iln=0; iln < m_msg.nln; iln++) {
      if (m_msg.lb_mfd[iln]!=NULL) {
	closelib(m_msg.lb_mfd[iln]);
      }
    }

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

    ttdisp += tdone-tscan;

    /* reset pst parameters to original */
    pst.zsflag = m_msg.zsflag;
    pst.n1_low = m_msg.n1_low;
    pst.n1_high = m_msg.n1_high;

    /*     maxn = m_msg.max_tot; */

    m_msg.n0 = 
      QGETLIB (aa0[0], MAXTST, m_msg.qtitle, sizeof(m_msg.qtitle),
	       &qseek, &qlcont,q_file_p,&m_msg.q_off);
    if (m_msg.n0 <= 0) break;
    if ((bp=strchr(m_msg.qtitle,' '))!=NULL) *bp='\0';
    strncpy(info_qlabel, m_msg.qtitle,sizeof(info_qlabel));
    if (bp != NULL) *bp=' ';
    info_qlabel[sizeof(info_qlabel)-1]='\0';

#ifdef SUPERFAMNUM
    m_msg.nqsfnum = nsfnum;
    for (i=0; i <= nsfnum & i<10; i++) m_msg.qsfnum[i] = sfnum[i];
    m_msg.nqsfnum_n = nsfnum_n;
    for (i=0; i <= nsfnum_n & i<10; i++) m_msg.qsfnum_n[i] = sfnum_n[i];
#endif

    if (m_msg.q_ann_flg) {
      m_msg.n0 = ann_scan(aa0[0],m_msg.n0,&m_msg,m_msg.qdnaseq);
    }

    if (m_msg.term_code && m_msg.qdnaseq==SEQT_PROT &&
	aa0[0][m_msg.n0-1]!=m_msg.term_code) {
      aa0[0][m_msg.n0++]=m_msg.term_code;
      aa0[0][m_msg.n0]=0;
    }

#if defined(SW_ALTIVEC) || defined(SW_SSE2)
    /* for ALTIVEC, must pad with 16 NULL's */
    for (id=0; id<SEQ_PAD; id++) {aa0[0][m_msg.n0+id]=0;}
#endif

  }	/* end of while(1) for multiple queries */

  if (m_msg.markx & MX_M10FORM)
    fprintf(outfd,">>><<<\n");

  tdone = s_time();
  if ( m_msg.markx & MX_HTML) fputs("<p><pre>\n",outfd); 
  if (m_msg.std_output) {
    printsum(outfd, m_msg.db);
  }
  if ( m_msg.markx & MX_HTML) fputs("</pre>\n",outfd);
#ifdef HTML_HEAD
  if (m_msg.markx & MX_HTML) fprintf(outfd,"</body>\n</html>\n");
#endif
  if (m_msg.std_output && outfd!=stdout) printsum(stdout,m_msg.db);

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

void
printsum(FILE *fd, struct db_str ntt)
{
  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 */
  fprintf(fd, "\n\n%ld residues in %ld query   sequences\n", qtt.length, 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, "%.0f residues in %ld library sequences\n", db_tt, ntt.entries);
  }

#ifndef COMP_THR
  fprintf(fd," Scomplib [%s]\n start: %s done: %s\n",mp_verstr,tstr1,tstr2);
#else
  fprintf(fd," Tcomplib [%s] (%d proc)\n start: %s done: %s\n", mp_verstr,
    max_workers,tstr1,tstr2);
#endif
#ifndef COMP_MLIB
  fprintf(fd," Scan time: ");
  ptime(fd, tscan - tprev);
  fprintf (fd," Display time: ");
  ptime (fd, tdone - tscan);
#else
  fprintf(fd," Total Scan time: ");
  ptime(fd, ttscan);
  fprintf (fd," Total Display time: ");
  ptime (fd, ttdisp);
#endif
  fprintf (fd,"\n");
  fprintf (fd, "\nFunction used was %s [%s]\n", prog_func,verstr);
}

void fsigint()
{
  struct db_str db;

  db.entries = db.length = db.carry = 0;
  tdone = s_time();
  tddone = time(NULL);

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

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

  exit(1);
}

/* save the sequence meta-data for a sequence if we need to re-use the
   rbuf/wbuf buffer */

void preserve_seq(struct buf2_str *lib_buf2_p,
		  struct seq_record *best_seqs,
		  struct beststr *best) {
  struct seq_record *dest_seq_p, *saved_seq_p;
  struct beststr *next_bbp;

  saved_seq_p = lib_buf2_p->best_save->seq;
  dest_seq_p = &best_seqs[lib_buf2_p->best_save - best];
  lib_buf2_p->best_save->seq = dest_seq_p;

  for (next_bbp = lib_buf2_p->best_save->bbp_link;
       (next_bbp != NULL) && (next_bbp->seq == saved_seq_p)
	 && (next_bbp->n1 == saved_seq_p->n1);
       next_bbp = next_bbp->bbp_link) {
    next_bbp->seq = dest_seq_p;
  }

  memcpy(dest_seq_p,lib_buf2_p->seq,sizeof(struct seq_record));
  dest_seq_p->aa1b = NULL;
}

/* in the current version (fasta_35_01) save_best is used by both
   threaded and unthreaded versions */

#define COPY_RST_P(d,s) 		\
{ d->rst.score[0] = s->rst.score[0];	\
  d->rst.score[1] = s->rst.score[1];	\
  d->rst.score[2] = s->rst.score[2];	\
  d->rst.comp = s->rst.comp;		\
  d->rst.H = s->rst.H;			\
  d->rst.escore = s->rst.escore;	\
  d->rst.segnum = s->rst.segnum;	\
  d->rst.seglen = s->rst.seglen;	\
}

void
save_best(struct buf_head *lib_bhead_p, struct mngmsg m_msg, struct pstruct *ppst, 
	  struct db_str *ldb, FILE *fdata, int *qsfnum, struct hist_str *histp,
	  void **pstat_voidp
	  )
{
  double zscore;
  int i_score;
  struct beststr *bbp;
  struct buf2_str *rbuf_p, *lib_buf2_p;
  int i, t_best, t_rbest, t_qrbest, tm_best, t_n1, sc_ix;
  double e_score, tm_escore, t_rescore, t_qrescore;
  int buf2_cnt;
  int jstats;

  sc_ix = ppst->score_ix;

  lib_buf2_p = lib_bhead_p->buf2;
  buf2_cnt = lib_bhead_p->buf2_cnt;
  
  t_best = t_rbest = t_qrbest = -BIGNUM;
  tm_escore = t_rescore = t_qrescore = FLT_MAX;

  while (buf2_cnt--) { /* count down the number of results */
    rbuf_p = lib_buf2_p++;	/* step through the results buffer */

    if (rbuf_p->rst.score[0] == -BIGNUM) continue;

    i_score = rbuf_p->rst.score[sc_ix];
    e_score = rbuf_p->rst.escore;

    /* need to look for frame 0 if TFASTA, then save stats at frame 6 */
    if (fdata) {
      fprintf(fdata,
	      "%-12s %5d %6d %d %.5f %.5f %4d %4d %4d %g %d %d %8lld\n",
	      rbuf_p->seq->libstr,
	      0,
	      rbuf_p->seq->n1,rbuf_p->frame,rbuf_p->rst.comp,rbuf_p->rst.H,
	      rbuf_p->rst.score[0],rbuf_p->rst.score[1],rbuf_p->rst.score[2],
	      rbuf_p->rst.escore, rbuf_p->rst.segnum, rbuf_p->rst.seglen,
	      rbuf_p->seq->lseek);
    }

    t_n1 = rbuf_p->seq->n1;
    if (i_score > t_best) tm_best = t_best = i_score;
    if (e_score < tm_escore) tm_escore = e_score;

    if (m_msg.qshuffle) {
      if (rbuf_p->qr_score > t_qrbest)
	t_qrbest = rbuf_p->qr_score;
      if (rbuf_p->qr_escore < t_qrescore)
	t_qrescore = rbuf_p->qr_escore;
      
      if (rbuf_p->frame == m_msg.nitt1 && nqstats < m_msg.shuff_max) {
	qstats[nqstats].n1 = rbuf_p->seq->n1;	/* save the best score */
	qstats[nqstats].comp =  rbuf_p->rst.comp;
	qstats[nqstats].H = rbuf_p->rst.H;
	qstats[nqstats].escore = t_qrescore;
	qstats[nqstats++].score = t_qrbest;
	t_qrbest = -BIGNUM;	/* reset t_qrbest, t_qrescore */
	t_qrescore = FLT_MAX;
      }
    }

    if (ppst->zsflag >= 10 && rbuf_p->r_score > t_rbest) {
      t_rbest = rbuf_p->r_score;
      t_rescore = rbuf_p->r_escore;
    }

    /* statistics done for best score of set */

    if (rbuf_p->frame == m_msg.nitt1) {
      ldb->entries++;
      ldb->length += t_n1;
      if (ldb->length > LONG_MAX) {
	ldb->length -= LONG_MAX; ldb->carry++;
      }

      if (nstats < MAX_STATS ) {
	stats[nstats].n1 = t_n1;
	stats[nstats].comp = rbuf_p->rst.comp;
	stats[nstats].H = rbuf_p->rst.H;
	if (ppst->zsflag >= 10) {
	  tm_best = t_rbest;
	  tm_escore = t_rescore;
	  t_rbest = -BIGNUM;
	  t_rescore = FLT_MAX;
	}
	stats[nstats].escore  = tm_escore;
	stats[nstats++].score = tm_best;
	t_best = -BIGNUM;
	tm_escore = FLT_MAX;
      }
      else if (ppst->zsflag > 0) {
	if (!stats_done) {
	  ppst->zsflag_f = process_hist(stats,nstats,m_msg, ppst,
				      histp, pstat_voidp,0);
	  kstats = nstats;
	  stats_done = 1;
	  for (i=0; i<MAX_BEST; i++) {
	    bestp_arr[i]->zscore = 
	      (*find_zp)(bestp_arr[i]->rst.score[ppst->score_ix],
			 bestp_arr[i]->rst.escore, bestp_arr[i]->seq->n1,
			 bestp_arr[i]->rst.comp, *pstat_voidp);
	  }
	}
#ifdef SAMP_STATS
	else {
	  if (!m_msg.escore_flg) {
	    jstats = nrand(++kstats);
	    if (jstats < MAX_STATS) {
	      stats[jstats].n1 = t_n1;
	      stats[jstats].comp = rbuf_p->rst.comp;
	      stats[jstats].H = rbuf_p->rst.H;
	      if (ppst->zsflag >= 10) {
		tm_best = t_rbest;
	      }
	      stats[jstats].score = tm_best;
	    }
	  }
	}
#endif
      }
    }

    /* calculate a z-score if possible, for histogram, and save threshold */
    if (stats_done) {

      zscore=(*find_zp)(i_score, e_score, rbuf_p->seq->n1,(double)rbuf_p->rst.comp,
			*pstat_voidp);

      if (rbuf_p->frame == m_msg.nitt1) {
	addhistz((*find_zp)(t_best, tm_escore, rbuf_p->seq->n1, (double) rbuf_p->rst.comp,
			    *pstat_voidp), histp);
	t_best = t_rbest = -BIGNUM;
	tm_escore = t_rescore = FLT_MAX;
      }
    }
    else zscore = (double) i_score;

    if (zscore > zbestcut) {
      if (nbest >= MAX_BEST) {
	bestfull = nbest-MAX_BEST/4;
	selectbestz(bestp_arr,bestfull-1,nbest);
	zbestcut = bestp_arr[bestfull-1]->zscore;
	nbest = bestfull;
      }
      bbp = bestp_arr[nbest++];

      COPY_RST_P(bbp, rbuf_p);

      bbp->seq = rbuf_p->seq;
      bbp->n1 = rbuf_p->seq->n1;
      if (rbuf_p->best_save) {	/* have we saved this seq before ? */
	/* yes - link back to the previous save, if it points to the
	   same seq_p */
	if (rbuf_p->best_save->seq == rbuf_p->seq) {
	  bbp->bbp_link = rbuf_p->best_save;
	}
	else {	/* break the link */
	  bbp->bbp_link = NULL;
	}
      }
      rbuf_p->best_save = bbp;	/* note where this seq is saved */	

      lib_bhead_p->have_best_save = 1;

      bbp->zscore = zscore;
      bbp->frame = rbuf_p->frame;
    }
  }
}

void
save_shuf(struct buf_head *lib_bhead_p, int nitt1, int shuff_max)
{
  struct buf2_str *rbuf_p, *lib_buf2_p;
  int t_rbest;
  double t_rescore;
  int buf2_cnt, jstats;
  static int kstats=0;

  lib_buf2_p = lib_bhead_p->buf2;
  buf2_cnt = lib_bhead_p->buf2_cnt;
  
  t_rbest = -BIGNUM;

  while (buf2_cnt--) { /* count down the number of results */
    rbuf_p = lib_buf2_p++;	/* step through the results buffer */

    if (rbuf_p->r_score > t_rbest) {
      t_rbest = rbuf_p->r_score;
      t_rescore = rbuf_p->r_escore;
    }

    /* statistics done for best score of set */

    if (rbuf_p->frame == nitt1) {
      if (nrstats < shuff_max ) { kstats = jstats = nrstats++; }
      else {	/* randomly replace */
	jstats = nrand(++kstats);
	if (jstats >= shuff_max) goto done;
      }

      rstats[jstats].n1 = rbuf_p->seq->n1;
      rstats[jstats].comp = rbuf_p->rst.comp;
      rstats[jstats].H = rbuf_p->rst.H;
      rstats[jstats].escore  = t_rescore;
      rstats[jstats].score = t_rbest;
done:
      t_rbest = -BIGNUM;
    }
  }
}


void
buf_do_work(unsigned char **aa0,  int n0,
	    struct buf_head *lib_bhead_p, 
	    struct pstruct *ppst, void **f_str) {
  
  int buf2_cnt;
  struct buf2_str *lib_buf2_p;

  buf2_cnt = lib_bhead_p->buf2_cnt;
  lib_buf2_p = lib_bhead_p->buf2;

  while (buf2_cnt--) {

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

    if (lib_buf2_p->seq->n1 < ppst->n1_low ||
	lib_buf2_p->seq->n1 > ppst->n1_high ) {
      lib_buf2_p++;
      continue;
    }

    /*
    if (lib_buf2_p->seq->aa1b[0] == '\0') {
      fprintf(stderr,"invalid aa1 NULL at: %s frame: %d\n",
	      lib_buf2_p->seq->libstr, lib_buf2_p->frame);
    }
    */