hmmpfam.c

来自「这是一个基于HMM 模型的生物多序列比对算法的linux实现版本。hmmer」· C语言 代码 · 共 1,210 行 · 第 1/4 页

 *           
 *           On return, ghit and dhit contain info for all hits
 *           that satisfy the set thresholds. If an evalue
 *           cutoff is used at all, the lists will be overestimated --
 *           because the evalue will be underestimated until
 *           we know the final Z. (Thus the main program must recheck
 *           thresholds before printing any results.) If only
 *           score cutoffs are used, then the lists are correct,
 *           and may be printed exactly as they come (after 
 *           appropriate sorting, anyway). This is especially
 *           important for dynamic thresholding using Pfam
 *           score cutoffs -- the main caller cannot afford to
 *           rescan the HMM file just to get the GA/TC/NC cutoffs
 *           back out for each HMM, and neither do I want to
 *           burn the space to store them as I make a pass thru
 *           Pfam.
 *
 * Args:     hmmfile - name of HMM file
 *           hmmfp   - open HMM file (and at start of file)
 *           dsq     - digitized sequence
 *           sqinfo  - ptr to SQINFO optional info for dsq
 *           thresh  - score/evalue threshold information
 *           do_xnu     - TRUE to apply XNU filter to sequence
 *           do_forward - TRUE to use Forward() scores
 *           do_null2   - TRUE to adjust scores w/ ad hoc null2 model
 *           num_threads- number of threads, >= 1
 *           ghit    - global hits list
 *           dhit    - domain hits list
 *           ret_nhmm    - number of HMMs searched.
 *
 * Returns:  (void)
 */
static void
main_loop_threaded(char *hmmfile, HMMFILE *hmmfp, char *seq, SQINFO *sqinfo, 
		   struct threshold_s *thresh, int do_xnu, int do_forward, int do_null2,
		   int num_threads,
		   struct tophit_s *ghit, struct tophit_s *dhit, int *ret_nhmm) 
{
  unsigned char     *dsq;       /* digitized sequence      */
  int                nhmm;	/* number of HMMs searched */
  struct workpool_s *wpool;     /* pool of worker threads  */

  /* Prepare sequence.
   */
  dsq = DigitizeSequence(seq, sqinfo->len);
  if (do_xnu && Alphabet_type == hmmAMINO) XNU(dsq, sqinfo->len);

  wpool = workpool_start(hmmfile, hmmfp, dsq, sqinfo->name, sqinfo->len, 
			 do_forward, do_null2, thresh,
			 ghit, dhit, num_threads);
  workpool_stop(wpool);
  nhmm = wpool->nhmm;
  workpool_free(wpool);

  free(dsq);
  *ret_nhmm = nhmm;
  return;
}
/* Function: workpool_start()
 * Date:     SRE, Mon Sep 28 11:10:58 1998 [St. Louis]
 *
 * Purpose:  Initialize a workpool_s structure, and return it.
 *
 * Args:     hmmfile    - name of HMM file
 *           hmmfp      - open HMM file, at start
 *           dsq        - ptr to sequence to search
 *           seqname    - ptr to name of dsq
 *           L          - length of dsq
 *           do_forward - TRUE to score using Forward
 *           do_null2   - TRUE to apply null2 ad hoc correction
 *           threshold  - evalue/score threshold settings
 *           ghit       - per-seq hit list
 *           dhit       - per-domain hit list             
 *           num_threads- number of worker threads to run.
 *
 * Returns:  ptr to struct workpool_s.
 *           Caller must wait for threads to finish with workpool_stop(),
 *           then free the structure with workpool_free().
 */
static struct workpool_s *
workpool_start(char *hmmfile, HMMFILE *hmmfp, unsigned char *dsq, char *seqname, int L,
	       int do_forward, int do_null2, struct threshold_s *thresh,
	       struct tophit_s *ghit, struct tophit_s *dhit, 
	       int num_threads)
{
  struct workpool_s *wpool;
  pthread_attr_t    attr;
  int i;
  int rtn;

  wpool             = MallocOrDie(sizeof(struct workpool_s));
  wpool->thread     = MallocOrDie(num_threads * sizeof(pthread_t));
  wpool->hmmfile    = hmmfile;
  wpool->dsq        = dsq;
  wpool->L          = L;
  wpool->seqname    = seqname;
  wpool->do_forward = do_forward;
  wpool->do_null2   = do_null2;
  wpool->thresh     = thresh;

  wpool->hmmfp      = hmmfp;
  wpool->nhmm       = 0;
  if ((rtn = pthread_mutex_init(&(wpool->input_lock), NULL)) != 0)
    Die("pthread_mutex_init FAILED; %s\n", strerror(rtn));

  wpool->ghit       = ghit;
  wpool->dhit       = dhit;
  if ((rtn = pthread_mutex_init(&(wpool->output_lock), NULL)) != 0)
    Die("pthread_mutex_init FAILED; %s\n", strerror(rtn));

  wpool->num_threads= num_threads;

  /* Create slave threads. See comments in hmmcalibrate.c at 
   * this step regarding concurrency and system scope.
   */
  pthread_attr_init(&attr);
#ifndef __sgi
#ifdef HAVE_PTHREAD_ATTR_SETSCOPE
  pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
#endif
#endif
#ifdef HAVE_PTHREAD_SETCONCURRENCY
  pthread_setconcurrency(num_threads+1);
#endif
  for (i = 0; i < num_threads; i++)
    if ((rtn = pthread_create(&(wpool->thread[i]), &attr,
			      worker_thread , (void *) wpool)) != 0)
      Die("Failed to create thread %d; return code %d\n", i, rtn);

  pthread_attr_destroy(&attr);
  return wpool;
}

/* Function: workpool_stop()
 * Date:     SRE, Thu Jul 16 11:20:16 1998 [St. Louis]
 *
 * Purpose:  Waits for threads in a workpool to finish.
 *
 * Args:     wpool -- ptr to the workpool structure
 *
 * Returns:  (void)
 */
static void
workpool_stop(struct workpool_s *wpool)
{
  int i;
				/* wait for threads to stop */
  for (i = 0; i < wpool->num_threads; i++)
    if (pthread_join(wpool->thread[i],NULL) != 0)
      Die("pthread_join failed");
  return;
}

/* Function: workpool_free()
 * Date:     SRE, Thu Jul 16 11:26:27 1998 [St. Louis]
 *
 * Purpose:  Free a workpool_s structure, after the threads
 *           have finished.
 *
 * Args:     wpool -- ptr to the workpool.
 *
 * Returns:  (void)
 */
static void
workpool_free(struct workpool_s *wpool)
{
  free(wpool->thread);
  free(wpool);
  return;
}


/* Function: worker_thread()
 * Date:     SRE, Mon Sep 28 10:48:29 1998 [St. Louis]
 *
 * Purpose:  The procedure executed by the worker threads.
 *
 * Args:     ptr  - (void *) that is recast to a pointer to
 *                  the workpool.
 *
 * Returns:  (void *)
 */
void *
worker_thread(void *ptr)
{
  struct workpool_s *wpool;     /* our working threads structure   */
  struct plan7_s    *hmm;       /* an HMM to search with           */
  struct p7trace_s  *tr;        /* traceback from an alignment     */
  float  sc;			/* score of an alignment           */
  int    rtn;			/* a return code from pthreads lib */
  double pvalue;		/* P-value of score                */
  double evalue;		/* E-value of a score              */
  struct threshold_s thresh;	/* a local copy of thresholds      */
  struct dpmatrix_s *mx;        /* growable DP matrix              */

  wpool = (struct workpool_s *) ptr;
  /* Because we might dynamically change the thresholds using
   * Pfam GA/NC/TC cutoffs, we make a local copy of the threshold
   * structure in this thread.
   */ 
  thresh.globT   = wpool->thresh->globT;
  thresh.globE   = wpool->thresh->globE;
  thresh.domT    = wpool->thresh->domT;
  thresh.domE    = wpool->thresh->domE;
  thresh.autocut = wpool->thresh->autocut;
  thresh.Z       = wpool->thresh->Z;
  
  /* 
   * We'll create for at least N=300xM=300, and thus consume at least 1 MB,
   * regardless of RAMLIMIT -- this helps us avoid reallocating some weird
   * asymmetric matrices.
   * 
   * We're growable in both M and N, because inside of P7SmallViterbi,
   * we're going to be calling P7Viterbi on subpieces that vary in size,
   * and for different models.
   */
  mx = CreatePlan7Matrix(300, 300, 25, 25);

  for (;;) {

    /* 1. acquire lock on HMM input, and get
     *    the next HMM to work on.
     */
				/* acquire a lock */
    if ((rtn = pthread_mutex_lock(&(wpool->input_lock))) != 0)
      Die("pthread_mutex_lock failure: %s\n", strerror(rtn));
    
    if (! HMMFileRead(wpool->hmmfp, &hmm)) 
      {	/* we're done. release lock, exit thread */
	if ((rtn = pthread_mutex_unlock(&(wpool->input_lock))) != 0)
	  Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));
	FreePlan7Matrix(mx);
	pthread_exit(NULL);
      }
    wpool->nhmm++;
    SQD_DPRINTF1(("a thread is working on %s\n", hmm->name));
				/* release the lock */
    if ((rtn = pthread_mutex_unlock(&(wpool->input_lock))) != 0)
      Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));

    if (hmm == NULL)
      Die("HMM file %s may be corrupt or in incorrect format; parse failed", wpool->hmmfile);
    P7Logoddsify(hmm, !(wpool->do_forward));

    if (!SetAutocuts(&thresh, hmm)) 
      Die("HMM %s did not have the right GA, NC, or TC cutoffs", hmm->name);

    /* 2. We have an HMM in score form.
     *    Score the sequence.
     */
    if (P7ViterbiSpaceOK(wpool->L, hmm->M, mx))
      sc = P7Viterbi(wpool->dsq, wpool->L, hmm, mx, &tr);
    else
      sc = P7SmallViterbi(wpool->dsq, wpool->L, hmm, mx, &tr);

    
    /* The Forward score override (see comments in serial vers)
     */
    if (wpool->do_forward) {
      sc  = P7Forward(wpool->dsq, wpool->L, hmm, NULL);
      if (wpool->do_null2)   sc -= TraceScoreCorrection(hmm, tr, wpool->dsq);
    }
    
    /* 3. Save the output in tophits structures, after acquiring a lock
     */
    if ((rtn = pthread_mutex_lock(&(wpool->output_lock))) != 0)
      Die("pthread_mutex_lock failure: %s\n", strerror(rtn));
    SQD_DPRINTF1(("model %s scores %f\n", hmm->name, sc));

    pvalue = PValue(hmm, sc);
    evalue = thresh.Z ? (double) thresh.Z * pvalue : (double) wpool->nhmm * pvalue;
    if (sc >= thresh.globT && evalue <= thresh.globE) 
      { 
	sc = PostprocessSignificantHit(wpool->ghit, wpool->dhit, 
				       tr, hmm, wpool->dsq, wpool->L, 
				       wpool->seqname, 
				       NULL, NULL, /* won't need seq's acc or desc */
				       wpool->do_forward, sc,
				       wpool->do_null2,
				       &thresh,
				       TRUE); /* TRUE -> hmmpfam mode */
      }
    if ((rtn = pthread_mutex_unlock(&(wpool->output_lock))) != 0)
      Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));

    P7FreeTrace(tr);
    FreePlan7(hmm);

  } /* end 'infinite' loop over HMMs in this thread */
}
#else /* HMMER_THREADS off; put in a dummy stub: */
static void
main_loop_threaded(char *hmmfile, HMMFILE *hmmfp, char *seq, SQINFO *sqinfo, 
		   struct threshold_s *thresh, int do_xnu, int do_forward, int do_null2,
		   int num_threads,
		   struct tophit_s *ghit, struct tophit_s *dhit, int *ret_nhmm) 
{
  Die("no threads support");
}
#endif /* HMMER_THREADS */