comp_lib2.c

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

    if (nbest > 0 && nbest < m_msg.shuff_max) {

      /* (1) get the sequences into a buffer - the sequence
             information is currently in the bestp_arr - find out how
             many we have, and how many we will need - the number to shuffle */

      /* figure out how much space we need */
      n1lib_req = 0;
      for (i = 0; i < nbest; i++) {
	if (bestp_arr[i]->seq->aa1b == NULL) {
	  n1lib_req += bestp_arr[i]->seq->n1 + 2;
	}
      }

#ifndef COMP_THR
      if (n1lib_req >= maxn) { /* we need new space, aa1shuff is too small */
	if ((aa1shuff = aa1shuff_b = 
	     (unsigned char *)realloc(aa1shuff_b, n1lib_req*sizeof(char)))==NULL) {
	  fprintf(stderr," *** cannot realloc aa1shuff[%d]\n",n1lib_req);
	  exit(1);
	}
	*aa1shuff = '\0';
	aa1shuff++;
      }
#else
      if ((aa1shuff = aa1shuff_b = 
	   (unsigned char *)calloc(n1lib_req,sizeof(char)))==NULL) {
	fprintf(stderr," *** cannot realloc aa1shuff[%d]\n",n1lib_req);
	exit(1);
      }
      *aa1shuff = '\0';
      aa1shuff++;
#endif

      shuff_mult = (m_msg.shuff_max / nbest) + 1;
      
#ifdef COMP_THR
      max_do_cnt = min(max_buf2_cnt,m_msg.shuff_max / (2 * max_workers));
      /* we don't have a left over one, so we need one */
      if (empty_reader_bufs == 0) {
	get_rbuf(&lib_bhead_p,max_work_buf);
	empty_reader_bufs++;
      }
#else
      max_do_cnt = max_buf2_cnt;
      lib_bhead_p = lib_buf2_list;	/* equivalent to un-threaded get_rbuf() */
#endif
      lib_bhead_p->buf2_cnt = 0;
      lib_bhead_p->have_results = 0;
      lib_bhead_p->buf2_type=BUF2_DOSHUF;
      lib_buf2_p = lib_bhead_p->buf2;

      for (i = 0; i < nbest; i++) {
	bbp = bestp_arr[i];
	if (bbp->seq->aa1b == NULL) {
	  /* get the sequence */
	  (bbp->seq->m_file_p->ranlib)(l_bline, sizeof(l_bline),
				       bbp->seq->lseek,bbp->seq->libstr,bbp->seq->m_file_p);
	  n1 = re_getlib(aa1save,maxn,m_msg.maxt3,
			 m_msg.l_overlap,bbp->seq->cont,m_msg.term_code,
			 &loffset,&l_off,bbp->seq->m_file_p);
	  memcpy(aa1shuff, aa1save, n1+1);
	  bbp->seq->aa1b = aa1shuff;
	  aa1shuff += n1 + 1;
	}

	for (j = 0; j < shuff_mult; j++ ) {
	  for (itt = m_msg.revcomp; itt <= m_msg.nitt1; itt++) {
	    /* this invalidates lib_buf2_p->seq */
	    lib_buf2_p->seq = bbp->seq;
	    lib_buf2_p->frame = itt;
	    lib_buf2_p++;		/* point to next buf2 */
	    lib_bhead_p->buf2_cnt++;
	  } /* for (itt .. */

	  if (lib_bhead_p->buf2_cnt >= max_do_cnt) {
	    /* (2) send sequences for shuffling */
#ifdef COMP_THR		/* if COMP_THR - fill and empty buffers */
	    /* provide empty buffer to workers */
	    lib_bhead_p->have_data = 1;
	    put_rbuf(lib_bhead_p,max_work_buf);
	    empty_reader_bufs--;

	    /* get an empty buffer to fill */
	    get_rbuf(&lib_bhead_p,max_work_buf);
	    empty_reader_bufs++;
#else			/* just do the searches */
	    if (lib_bhead_p->buf2_type & BUF2_DOSHUF) {
	      buf_shuf_work(aa0,m_msg.n0, aa1save, pst.zs_win,
			    lib_bhead_p, &pst, f_str, pst.score_ix);
	    }
#endif
	    /* (3) save results in the rstats structure */
	    if (lib_bhead_p->buf2_cnt > 0 && lib_bhead_p->have_results) {
	      save_shuf(lib_bhead_p,m_msg.nitt1,m_msg.shuff_max);
	    }

	    lib_bhead_p->buf2_cnt = 0;
	    lib_bhead_p->have_results = 0;
	    lib_bhead_p->buf2_type=BUF2_DOSHUF;
	    lib_buf2_p = lib_bhead_p->buf2;
	  }
	}
      }			/* done with bestp_arr[] */

#ifdef COMP_THR		/* if COMP_THR - fill and empty buffers */
      /* check last buffers for any results */
      if (lib_bhead_p->buf2_cnt > 0) {
	put_rbuf(lib_bhead_p,max_work_buf);
	empty_reader_bufs--;
      }
    
      /* wait for the threads to finish */

      wait_rbuf(max_work_buf - empty_reader_bufs);
      /*
      fprintf(stderr, " num_reader[%d]-empty[%d]: %d\tnrstats: %d\n",
	      num_reader_bufs,empty_reader_bufs,
	      num_reader_bufs-empty_reader_bufs, nrstats);
      */

    for (i=0, t_reader_buf_readp = reader_buf_readp-1;
	 i < num_reader_bufs - empty_reader_bufs ; i++, t_reader_buf_readp--) {
      if (t_reader_buf_readp < 0) t_reader_buf_readp = max_work_buf - 1;
      t_lib_bhead_p = reader_buf[t_reader_buf_readp];

      /*
      fprintf(stderr," buf2_cnt[%d]: %d [data/res:%d/%d]\n",
	      t_reader_buf_readp,t_lib_bhead_p->buf2_cnt,
      	      t_lib_bhead_p->have_data, t_lib_bhead_p->have_results);
      */

      if (t_lib_bhead_p->buf2_cnt > 0 && t_lib_bhead_p->have_results) {
	  save_shuf(t_lib_bhead_p,m_msg.nitt1, m_msg.shuff_max);
	  t_lib_bhead_p->buf2_cnt = t_lib_bhead_p->have_results = 0;
	}

      /*
	fprintf(stderr," nrstats: %d\n",nrstats);
      */
      }
#else	/* just do the searches */
      /* aa1save is used for shuffles, not aa1shuf, because aa1shuf
	 has library sequences */
      buf_shuf_work(aa0,m_msg.n0, aa1save, pst.zs_win,
		    lib_bhead_p, &pst, f_str, pst.score_ix);

      save_shuf(lib_bhead_p,m_msg.nitt1,m_msg.shuff_max);
      lib_bhead_p->buf2_cnt = lib_bhead_p->have_results = 0;
#endif
      /* (4) analyze rstats */
      if (pst.zsflag < 10) pst.zsflag += 10;
      pst.zsflag_f = process_hist(rstats,nrstats,m_msg, &pst,&m_msg.hist,
				  &m_msg.pstat_void,0);
    }

    if (!pst.zdb_size_set) pst.zdb_size = m_msg.ldb.entries;

#ifdef COMP_THR
    /* before I call last_calc/showbest/showalign, I need init_work() to
       get an f_str. This duplicates some code above, which is used in
       the non-threaded version
    */

    if (!have_f_str) {
      init_work(aa0[0],m_msg.n0,&pst,&f_str[0]);
      have_f_str = 1;
      f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] =  f_str[0];

      if (m_msg.qframe == 2) {
	if ((aa0[1]=(unsigned char *)calloc((size_t)m_msg.n0+2,
					    sizeof(unsigned char)))==NULL) {
	  fprintf(stderr," cannot allocate aa0[1][%d] for alignments\n",
		  m_msg.n0+2);
	}
	*aa0[1]='\0';
	aa0[1]++;
	memcpy(aa0[1],aa0[0],m_msg.n0+1);
	revcomp(aa0[1],m_msg.n0,&pst.c_nt[0]);
	init_work(aa0[1],m_msg.n0,&pst,&f_str[1]);
      }
    }
#endif

    /* now we have one set of scaled scores for in bestp_arr  -
       for FASTS/F, we need to do some additional processing */

    if (!m_msg.qshuffle) {
      last_stats(aa0[0], m_msg.n0, stats,nstats, bestp_arr,nbest,
		 m_msg, &pst, &m_msg.hist, &m_msg.pstat_void);
    }
    else {
      last_stats(aa0[0], m_msg.n0,
		 qstats,nqstats, bestp_arr,nbest, m_msg, &pst, 
		 &m_msg.hist, &m_msg.pstat_void);
    }

    /* here is a contradiction: if pst.zsflag < 0, then m_msg.pstat_void
       should be NULL; if it is not, then process_hist() has been called */
    if (pst.zsflag < 0 && m_msg.pstat_void != NULL) pst.zsflag = 1;

    if (m_msg.last_calc_flg) {
      /* last_calc may need coefficients from last_stats() */
      nbest = last_calc(aa0, aa1save, maxn, bestp_arr, nbest, m_msg, &pst,
			f_str, m_msg.pstat_void);
    }

    /* in addition to scaling scores, this sorts bestp_arr[nbest] */
    scale_scores(bestp_arr,nbest,m_msg.db, &pst,m_msg.pstat_void);

    get_param(&pst, info_gstring2p,info_gstring3);

    /* **************************************************************** */
    /* label Library: output                                            */
    /* **************************************************************** */

    if (m_msg.std_output) {
      if (m_msg.db.carry==0) {
	fprintf(stdout, " %7ld residues in %5ld sequences\n", m_msg.db.length, m_msg.db.entries);
      }
      else {
	db_tt = (double)m_msg.db.carry*(double)LONG_MAX + (double)m_msg.db.length;
	fprintf(stdout, " %.0f residues in %5ld library sequences\n", db_tt, m_msg.db.entries);
      }

      prhist (stdout, m_msg, &pst, m_msg.hist, nstats, m_msg.ldb, info_lib_range_p,
	      info_gstring2p, info_hstring_p);

      tscan = s_time();
      printf (" Scan time: ");
      ptime(stdout,tscan-tprev);
      printf ("\n");
    }
#ifdef COMP_MLIB
    ttscan += tscan-tprev;
#endif

  l3:
    if (!m_msg.quiet) {
      printf("Enter filename for results [%s]: ", m_msg.outfile);
      fflush(stdout);
    }

    rline[0]='\0';
    if (!m_msg.quiet && fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
    if ((bp=strchr(rline,'\n'))!=NULL) *bp = '\0';
    if (rline[0]!='\0') strncpy(m_msg.outfile,rline,sizeof(m_msg.outfile));
    if (m_msg.outfile[0]!='\0') {
      if ((outfd=fopen(m_msg.outfile,"w"))==NULL) {
	fprintf(stderr," could not open %s\n",m_msg.outfile);
	if (!m_msg.quiet) goto l3;
	else goto l4;
      }

#ifdef PGM_DOC
      fputs(argv_line,outfd);
      fputc('\n',outfd);
#endif  
      fputs(iprompt0,outfd);
      fprintf(outfd," %s%s\n",verstr,refstr);

      fprintf(outfd,"Query: %s%s, %d %s\n",
	      info_qlabel, (m_msg.revcomp ? "-" : "\0"), m_msg.n0,
	      m_msg.sqnam);

#ifdef COMP_MLIB
      fprintf(outfd,"%3d>>>%s - %d %s%s\n", qlib,
	   m_msg.qtitle, m_msg.n0, m_msg.sqnam,
	   (m_msg.revcomp ? " (reverse complement)" : rline));
#else
      fprintf(outfd,"%.50s: %d %s%s\n %s\n",
	   info_qlabel, m_msg.n0, m_msg.sqnam,
	   (m_msg.revcomp ? " (reverse complement)" : rline));
#endif
      fprintf(outfd, "Library: %.60s%s", m_msg.ltitle,info_lib_range);
      if (m_msg.db.carry==0) {
	fprintf(outfd, " %7ld residues in %5ld sequences\n", m_msg.db.length, m_msg.db.entries);
      }
      else {
	db_tt = (double)m_msg.db.carry*(double)LONG_MAX + (double)m_msg.db.length;
	fprintf(outfd, " %.0f residues in %5ld library sequences\n", db_tt, m_msg.db.entries);
      }

      prhist(outfd,m_msg, &pst,m_msg.hist, nstats, m_msg.db, info_lib_range,
	     info_gstring2p, info_hstring_p);
    }

  l4:   
    if (m_msg.markx & MX_HTML) {
      fputs("</pre>\n<p>\n<hr>\n<p>\n",outfd);
    }

    /* code from p2_complib.c to pre-calculate -m 9 alignment info -
       requires -q with -m 9 */

    if (m_msg.quiet || m_msg.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_msg.mshow_flg != 1 && pst.zsflag >= 0) {
	for (i=0; i<nbest && bestp_arr[i]->rst.escore< m_msg.e_cut; i++) {}
	m_msg.mshow = i;
      }

      if (m_msg.mshow <= 0) { /* no results to display */
	fprintf(outfd,"!! No sequences with E() < %f\n",m_msg.e_cut);
	m_msg.nshow = 0;
	goto end_l;
      }
    }

    if (m_msg.do_showbest) {
      showbest (stdout, aa0, aa1save, maxn, bestp_arr, nbest, qtt.entries, &m_msg, &pst,
		m_msg.db, info_gstring2p, f_str);

      if (outfd != stdout) {
	t_quiet = m_msg.quiet;
	m_msg.quiet = -1;	/* should guarantee 1..nbest shown */
	showbest (outfd, aa0, aa1save, maxn, bestp_arr, nbest, qtt.entries, &m_msg, &pst,
		  m_msg.db, info_gstring2p, f_str);
	m_msg.quiet = t_quiet;
      }
    }

    if (m_msg.nshow > 0) {
      rline[0]='N';
      if (!m_msg.quiet){
	printf(" Display alignments also? (y/n) [n] "); fflush(stdout);
	if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
      }
      else rline[0]='Y';

      if (toupper((int)rline[0])=='Y') {
	if (!m_msg.quiet && m_msg.do_showbest) {
	  printf(" number of alignments [%d]? ",m_msg.nshow);
	  fflush(stdout);
	  if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
	  if (rline[0]!=0) sscanf(rline,"%d",&m_msg.nshow);
	  m_msg.ashow=m_msg.nshow;
	}

	if (m_msg.markx & (MX_AMAP+ MX_HTML + MX_M9SUMM) && !(m_msg.markx & MX_M10FORM)) {
	  fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
		  info_qlabel,(m_msg.revcomp ? "_rev":"\0"), m_msg.n0,
		  m_msg.sqnam,m_msg.lname);
	}

	if (m_msg.markx & MX_M10FORM) {
	  fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
		  info_qlabel,(m_msg.revcomp ? "-":"\0"), m_msg.n0, m_msg.sqnam,
		  m_msg.lname);
	  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);
	}

	showalign (outfd, aa0, aa1save, maxn, bestp_arr, nbest, qtt.entries,
		   m_msg, &pst, info_gstring2p, f_str);

	fflush(outfd);
      }
    }

  end_l:

    if (fdata) {
      fprintf(fdata,"/** Algorithm : %s  **/\n",info_gstring2p[0]);
      fprintf(fdata,"/** Parameters : %s  **/\n",info_gstring2p[1]);
      fprintf(fdata,"%3ld%-50s\n",qtt.entries-1,m_msg.qtitle);
      fflush(fdata);
    }
    
#if defined(COMP_THR)
    rbuf_done(max_workers);

    for (i=0; i<max_work_buf; i++) {
      lib_buf2_list[i].buf2_cnt=
	lib_buf2_list[i].have_results=
	lib_buf2_list[i].have_best_save = 0;
    }

    num_worker_bufs = 0;
    num_reader_bufs = max_work_buf;
    reader_done = 0;
    reader_wait = 1;
    worker_buf_workp = 0;
    worker_buf_readp = 0;
    reader_buf_workp = 0;
    reader_buf_readp = 0;

    start_thread = 1;	/* stop thread from starting again */
    
#endif
    /* clean up best_seqs */
    memset(best_seqs,0,(MAX_BEST+1)*sizeof(struct seq_record));

    /* re-initialize lib_buf2_list buffers */
    for (lib_bhead_p = lib_buf2_list; 
	 lib_bhead_p < lib_buf2_list+max_work_buf; ) {

      /* memset(lib_bhead_p->start-1,0,(size_t)buf_siz); */
      current_seq_p = lib_bhead_p->buf2[0].seq_b;
      /* this wipes out lib_bhead_p->buf2[0].seq_b */
      memset(lib_bhead_p->buf2,0,(size_t)(max_buf2_cnt+1)*sizeof(struct buf2_str));
      /* replace it */
      lib_bhead_p->buf2[0].seq_b = lib_bhead_p->buf2[0].seq = current_seq_p;
      lib_bhead_p->buf2[0].seq->aa1b = lib_bhead_p->start;
      lib_bhead_p++->have_results = 0;
    }

    /* re-initialize library counts */
    m_msg.ldb.length = 0l;
    m_msg.ldb.entries = m_msg.ldb.carry = 0;

    /* clean up alignment encodings */
    for (i=0; i < m_msg.nshow; i++) {
      if (bestp_arr[i]->have_ares & 0x2) {
	next_ares = bestp_arr[i]->a_res.next;
	free(bestp_arr[i]->a_res.res);
	bestp_arr[i]->a_res.res = NULL;
	bestp_arr[i]->have_ares = 0;
	while (next_ares != NULL) {
	  old_ares = next_ares;
	  next_ares = next_ares->next;
	  free(old_ares->res);
	  free(old_ares);
	}
      }
    }