hmmpfam.c

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

	SQD_DPRINTF1(("   ... using P7SmallViterbi(); size ~%d MB\n",
		      P7ViterbiSize(sqinfo->len, hmm->M)));
	sc = P7SmallViterbi(dsq, sqinfo->len, hmm, mx, &tr);
      }

    /* Implement do_forward; we'll override the whole_sc with a P7Forward()
     * calculation.
     * HMMER is so trace- (alignment-) dependent that this gets a bit hacky.
     * Some important implications: 
     *   1) if --do_forward is selected, the domain (Viterbi) scores do not
     *      necessarily add up to the whole sequence (Forward) score.
     *   2) The implementation of null2 for a Forward score is undefined,
     *      since the null2 correction is trace-dependent. As a total hack, 
     *      we use a null2 correction derived from the whole trace
     *      (which was the behavior of HMMER 2.1.1 and earlier, anyway).
     *      This could put the sum of domain scores and whole seq score even
     *      further in disagreement. 
     * 
     * Note that you can't move the Forward calculation into 
     * PostprocessSignificantHit(). The Forward score will exceed the
     * Viterbi score, so you can't apply thresholds until you
     * know the Forward score. Also, since PostprocessSignificantHit()
     * is wrapped by a mutex in the threaded implementation,
     * you'd destroy all useful parallelism if PostprocessSignificantHit()
     * did anything compute intensive. 
     */
    if (do_forward) {
      sc = P7Forward(dsq, sqinfo->len, hmm, NULL);
      if (do_null2) sc -= TraceScoreCorrection(hmm, tr, dsq);
    }
	
    /* Store scores/pvalue for each HMM aligned to this sequence, overall
     */
    pvalue = PValue(hmm, sc);
    evalue = thresh->Z ? (double) thresh->Z * pvalue : (double) nhmm * pvalue;
    if (sc >= thresh->globT && evalue <= thresh->globE) {
      sc = PostprocessSignificantHit(ghit, dhit, 
				     tr, hmm, dsq, sqinfo->len, 
				     sqinfo->name, NULL, NULL, /* won't need acc or desc even if we have 'em */
				     do_forward, sc,
				     do_null2,
				     thresh,
				     TRUE); /* TRUE -> hmmpfam mode */
    }
    P7FreeTrace(tr);
    FreePlan7(hmm);
    nhmm++;
  }

  FreePlan7Matrix(mx);
  free(dsq);
  *ret_nhmm = nhmm;
  return;
}


#ifdef HMMER_PVM
/*****************************************************************
 * PVM specific functions
 ****************************************************************/ 

/* Function: main_loop_pvm()
 * Date:     SRE, Fri Aug  7 13:58:34 1998 [St. Louis]
 *
 * Purpose:  Search a sequence against an HMM database;
 *           main loop for the PVM version.
 *
 * Args:     hmmfile - name of HMM file
 *           hmmfp   - open HMM file (and at start of file)
 *           seq     - sequence to search against
 *           sqinfo  - ptr to SQINFO optional info for dsq
 *           thresh  - score/evalue threshold settings
 *           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
 *           ghit    - global hits list
 *           dhit    - domain hits list
 *           nhmm    - number of HMMs searched.
 *
 * Returns:  (void)
 */
static void
main_loop_pvm(char *hmmfile, HMMFILE *hmmfp, char *seq, SQINFO *sqinfo, 
	      struct threshold_s *thresh, int do_xnu, int do_forward, int do_null2,
	      struct tophit_s *ghit, struct tophit_s *dhit, int *ret_nhmm)
{
  struct plan7_s   *hmm;        /* HMM that was searched with */
  struct p7trace_s *tr;         /* a traceback structure */
  unsigned char    *dsq;        /* digitized sequence */ 
  float  sc;			/* score of an HMM match */
  int    master_tid;		/* master's ID */
  int   *slave_tid;             /* array of slave IDs */
  int   *hmmlist;               /* array of hmm indexes being worked on by slaves */
  int    nslaves;	       	/* number of slaves in virtual machine   */
  int    nhmm;			/* number of HMMs searched */
  int    slaveidx;		/* index of a slave wanting work         */
  int    slave, msg;
  int    sent_trace;		/* TRUE if slave sent us a trace */
  char   slavename[32];		/* name of HMM that slave actually did */
  double pvalue;		/* pvalue of HMM score */
  int    arglen;

  /* Sanity checks.
   */
  if (hmmfp->ssi == NULL)
    Die("HMM file %s needs an SSI index to use PVM. See: hmmindex.", hmmfile);

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

  /* Initialize PVM
   */
  master_tid = pvm_mytid();
#if DEBUGLEVEL >= 1  
  pvm_catchout(stderr);		/* catch output for debugging */
#endif
  SQD_DPRINTF1(("Spawning slaves...\n"));
  PVMSpawnSlaves("hmmpfam-pvm", &slave_tid, &nslaves);
  hmmlist   = MallocOrDie(sizeof(int) * nslaves);
  SQD_DPRINTF1(("Spawned a total of %d slaves...\n", nslaves));

  /* Initialize the slaves
   */
  SQD_DPRINTF1(("Broadcasting to %d slaves...\n", nslaves));
  pvm_initsend(PvmDataDefault);
  arglen = strlen(hmmfile);
  pvm_pkint(&arglen, 1, 1);
  pvm_pkstr(hmmfile);
  pvm_pkint(&(sqinfo->len), 1, 1);
  pvm_pkstr(seq);
  pvm_pkfloat(&(thresh->globT), 1, 1);
  pvm_pkdouble(&(thresh->globE), 1, 1);
  pvm_pkint(&(thresh->Z), 1, 1);
  pvm_pkint((int *)&(thresh->autocut), 1, 1);
  pvm_pkint(&do_xnu, 1, 1);
  pvm_pkint(&do_forward, 1, 1);
  pvm_pkint(&do_null2, 1, 1);
  pvm_pkint(&Alphabet_type, 1, 1);
  pvm_mcast(slave_tid, nslaves, HMMPVM_INIT);
  SQD_DPRINTF1(("Slaves should be ready...\n"));
				/* get their OK codes. */
  PVMConfirmSlaves(slave_tid, nslaves);
  SQD_DPRINTF1(("Slaves confirm that they're ok...\n"));

  /* Load the slaves.
   * For efficiency reasons, we don't want the master to
   * load HMMs from disk until she absolutely needs them.
   */
  for (nhmm = 0; nhmm < nslaves && nhmm < hmmfp->ssi->nprimary; nhmm++) {
    pvm_initsend(PvmDataDefault);
    pvm_pkint(&nhmm, 1, 1);	/* side effect: also tells him what number he is. */
    pvm_send(slave_tid[nhmm], HMMPVM_WORK);
    hmmlist[nhmm] = nhmm;
  }
  SQD_DPRINTF1(("%d slaves are loaded\n", nhmm));

  
  /* Receive/send loop
   */
  for (; nhmm < hmmfp->ssi->nprimary; nhmm++)
    {
				/* check slaves before blocking */
      PVMCheckSlaves(slave_tid, nslaves);
				/* receive output */
      SQD_DPRINTF1(("Waiting for a slave to give me output...\n"));
      pvm_recv(-1, HMMPVM_RESULTS);
      pvm_upkint(&slaveidx, 1, 1);     /* # of slave who's sending us stuff */
      pvm_upkstr(slavename);           /* name of HMM that slave did */
      pvm_upkfloat(&sc, 1, 1);         /* score   */
      pvm_upkdouble(&pvalue, 1, 1);    /* P-value */
      pvm_upkint(&sent_trace, 1, 1);   /* TRUE if trace is coming */
      tr = (sent_trace) ? PVMUnpackTrace() : NULL;
      SQD_DPRINTF1(("Slave %d finished %s for me...\n", slaveidx, slavename));

				/* send new work */
      pvm_initsend(PvmDataDefault);
      pvm_pkint(&nhmm, 1, 1);
      pvm_send(slave_tid[slaveidx], HMMPVM_WORK);
      SQD_DPRINTF1(("Assigned %d -> slave %d\n", nhmm, slaveidx));

				/* process output */
      /* 1b. Store scores/pvalue for each HMM aligned to this sequence, overall
       */
      SQD_DPRINTF1(("%15s : %2d : %f\n", slavename, slaveidx, sc));
      if (sent_trace)
	{ 
				/* now load the HMM, because the hit is significant */
	  HMMFilePositionByIndex(hmmfp, hmmlist[slaveidx]);
	  if (!HMMFileRead(hmmfp, &hmm))
	    { pvm_exit(); Die("Unexpected failure to read HMM file %s", hmmfile); }
	  if (hmm == NULL) 
	    { pvm_exit(); Die("HMM file %s may be corrupt; parse failed", hmmfile); }
	  P7Logoddsify(hmm, TRUE);
	  if (! SetAutocuts(thresh, hmm))
	    Die("HMM %s did not contain your GA, NC, or TC cutoffs", hmm->name);
	
	  sc = PostprocessSignificantHit(ghit, dhit, 
					 tr, hmm, dsq, sqinfo->len, 
					 sqinfo->name, 
					 sqinfo->flags & SQINFO_ACC ? sqinfo->acc : NULL,
					 sqinfo->flags & SQINFO_DESC ? sqinfo->desc : NULL,
					 do_forward, sc,
					 do_null2,
					 thresh,
					 TRUE); /* TRUE -> hmmpfam mode */

	  FreePlan7(hmm);
	  P7FreeTrace(tr);
	}
      hmmlist[slaveidx] = nhmm;
    }

  /* Collect the output. all n slaves are still working, so wait for them.
   */
  for (slave = 0; slave < nslaves && slave < nhmm; slave++)
    {
				/* don't check slaves (they're exiting normally);
				   window of vulnerability here to slave crashes */
				/* receive output */
      pvm_recv(-1, HMMPVM_RESULTS);
      pvm_upkint(&slaveidx, 1, 1);     /* slave who's sending us stuff */
      pvm_upkstr(slavename);
      pvm_upkfloat(&sc, 1, 1);         /* one score */
      pvm_upkdouble(&pvalue, 1, 1);	   /* P-value */
      pvm_upkint(&sent_trace, 1, 1);   /* TRUE if trace is coming */
      tr = (sent_trace) ? PVMUnpackTrace() : NULL;

				/* process output */
      SQD_DPRINTF1(("%15s : %2d : %f\n", slavename, slaveidx, sc));
      if (sent_trace)
	{ 
	  /* now load the HMM, because the hit is significant */
	  HMMFilePositionByIndex(hmmfp, hmmlist[slaveidx]);
	  if (!HMMFileRead(hmmfp, &hmm))
	    { pvm_exit(); Die("Unexpected failure to read HMM file %s", hmmfile);}
	  if (hmm == NULL) 
	    { pvm_exit(); Die("HMM file %s may be corrupt; parse failed", hmmfile); }
	  P7Logoddsify(hmm, TRUE);
	  if (! SetAutocuts(thresh, hmm))
	    Die("HMM %s did not contain your GA, NC, or TC cutoffs", hmm->name);
	  
	  sc = PostprocessSignificantHit(ghit, dhit, 
					 tr, hmm, dsq, sqinfo->len, 
					 sqinfo->name, NULL, NULL, /* won't need acc or desc even if we have 'em */
					 do_forward, sc,
					 do_null2,
					 thresh,
					 TRUE); /* TRUE -> hmmpfam mode */

	  FreePlan7(hmm);
	  P7FreeTrace(tr);
	}
    }

  /* Shut down all the slaves. (xref bug #15/STL5 p.66).
   */
  for (slave = 0; slave < nslaves; slave++)
    {
      pvm_initsend(PvmDataDefault);
      msg = -1;			/* the special "shut down now" flag */
      pvm_pkint(&msg, 1, 1);
      pvm_send(slave_tid[slave], HMMPVM_WORK);
    }

  /* Cleanup; quit the PVM; and return
   */
  free(slave_tid);
  free(hmmlist);
  free(dsq);
  pvm_exit();
  *ret_nhmm = nhmm;
  return;
}
#else /* HMMER_PVM off; insert dummy stub function: */
static void
main_loop_pvm(char *hmmfile, HMMFILE *hmmfp, char *seq, SQINFO *sqinfo, 
	      struct threshold_s *thresh, int do_xnu, int do_forward, int do_null2,
	      struct tophit_s *ghit, struct tophit_s *dhit, int *ret_nhmm)
{
  Die("no PVM support");
}
#endif /*HMMER_PVM*/


#ifdef HMMER_THREADS
/*****************************************************************
 * POSIX threads implementation.
 * 
 * API:
 *      workpool_start()   (makes a workpool_s structure. Starts calculations.)
 *      workpool_stop()    (waits for threads to finish.)
 *      workpool_free()    (destroys the structure)
 *      
 * Threads:
 *      worker_thread()    (the actual parallelized worker thread).
 *****************************************************************/
/* Function: main_loop_threaded()
 * Date:     SRE, Wed Feb 27 11:10:13 2002 [St. Louis]
 *
 * Purpose:  Search a sequence against an HMM database;
 *           main loop for the threaded version.