hmmer-3.3/src/p7_domaindef.c
/* region_trace_ensemble()
 * SRE, Fri Feb  8 11:49:44 2008 [Janelia]
 *
 * Here, we've decided that region <ireg>..<jreg> in sequence <dsq> might be
 * composed of more than one domain, and we're going to use clustering
 * of a posterior ensemble of stochastic tracebacks to sort it out.
 * 
 * Caller provides a filled Forward matrix in <fwd> for the sequence
 * region <dsq+ireg-1>, length <jreg-ireg+1>, for the model <om>
 * configured in multihit mode with its target length distribution
 * set to the total length of <dsq>: i.e., the same model
 * configuration used to score the complete sequence (if it weren't
 * multihit, we wouldn't be worried about multiple domains).
 * 
 * Caller also provides a DP matrix in <wrk> containing at least one
 * row, for use as temporary workspace. (This will typically be the
 * caller's Backwards matrix, which we haven't yet used at this point
 * in the processing pipeline.)
 * 
 * Caller provides <ddef>, which defines heuristic parameters that
 * control the clustering, and provides working space for the
 * calculation and the answers. The <ddef->sp> object must have been
 * reused (i.e., it needs to be fresh; we're going to use it here);
 * the caller needs to Reuse() it specifically, because it can't just
 * Reuse() the whole <ddef>, when it's in the process of analyzing
 * regions.
 * 
 * Upon return, <*ret_nc> contains the number of clusters that were
 * defined.
 * 
 * The caller can retrieve info on each cluster by calling
 * <p7_spensemble_GetClusterCoords(ddef->sp...)> on the
 * <P7_SPENSEMBLE> object in <ddef>.
 * 
 * Other information on what's happened in working memory:
 * 
 * <ddef->n2sc[ireg..jreg]> now contains log f'(x_i) / f(x_i) null2 scores
 *    for each residue.
 *
 * <ddef->sp> gets filled in, and upon return, it's holding the answers 
 *    (the cluster definitions). When the caller is done retrieving those
 *    answers, it needs to <esl_spensemble_Reuse()> it before calling
 *    <region_trace_ensemble()> again.
 *    
 * <ddef->tr> is used as working memory for sampled traces.
 *    
 * <wrk> has had its zero row clobbered as working space for a null2 calculation.
 */


Here, we've decided that region .. in sequence  might be
composed of more than one domain, and we're going to use clustering
of a posterior ensemble of stochastic tracebacks to sort it out.

Caller provides a filled Forward matrix in  for the sequence
region , length , for the model 
configured in multihit mode with its target length distribution
set to the total length of : i.e., the same model
configuration used to score the complete sequence (if it weren't
multihit, we wouldn't be worried about multiple domains).

Caller also provides a DP matrix in  containing at least one
row, for use as temporary workspace. (This will typically be the
caller's Backwards matrix, which we haven't yet used at this point
in the processing pipeline.)

Caller provides , which defines heuristic parameters that
control the clustering, and provides working space for the
calculation and the answers. The sp> object must have been
reused ;
the caller needs to Reuse() it specifically, because it can't just
Reuse() the whole , when it's in the process of analyzing
regions.

Upon return, <*ret_nc> contains the number of clusters that were
defined.

The caller can retrieve info on each cluster by calling
sp...)> on the
 object in .

Other information on what's happened in working memory.

n2sc[ireg..jreg]> now contains log f'(x_i) / f(x_i) null2 scores
for each residue.

sp> gets filled in, and upon return, it's holding the answers
(the cluster definitions). When the caller is done retrieving those
answers, it needs to  it before calling
 again.

tr> is used as working memory for sampled traces.

 has had its zero row clobbered as working space for a null2 calculation.
Here we ' ve decided that region .. in sequence might be composed of more than one domain and we ' re going to use clustering of a posterior ensemble of stochastic tracebacks to sort it out . Caller provides a filled Forward matrix in for the sequence region length for the model configured in multihit mode with its target length distribution set to the total length of : i . e . the same model configuration used to score the complete sequence ( if it weren ' t multihit we wouldn ' t be worried about multiple domains ) . Caller also provides a DP matrix in containing at least one row for use as temporary workspace . ( This will typically be the caller ' s Backwards matrix which we haven ' t yet used at this point in the processing pipeline . ) Caller provides which defines heuristic parameters that control the clustering and provides working space for the calculation and the answers . The sp > object must have been reused ; the caller needs to Reuse () it specifically because it can ' t just Reuse () the whole when it ' s in the process of analyzing regions . Upon return < * ret_nc > contains the number of clusters that were defined . The caller can retrieve info on each cluster by calling sp ... ) > on the object in . Other information on what ' s happened in working memory . n2sc [ ireg .. jreg ] > now contains log f ' ( x_i ) / f ( x_i ) null2 scores for each residue . sp > gets filled in and upon return it ' s holding the answers ( the cluster definitions ) . When the caller is done retrieving those answers it needs to it before calling again . tr > is used as working memory for sampled traces . has had its zero row clobbered as working space for a null2 calculation .
static int
region_trace_ensemble(P7_DOMAINDEF *ddef, const P7_OPROFILE *om, const ESL_DSQ *dsq, int ireg, int jreg, 
		      const P7_OMX *fwd, P7_OMX *wrk, int *ret_nc)
{
  int    Lr  = jreg-ireg+1;
  int    t, d, d2;
  int    nov, n;
  int    nc;
  int    pos;
  float  null2[p7_MAXCODE];
  esl_vec_FSet(ddef->n2sc+ireg, Lr, 0.0); 
  if (ddef->do_reseeding) 
    esl_randomness_Init(ddef->r, esl_randomness_GetSeed(ddef->r));
  for (t = 0; t < ddef->nsamples; t++)
    {
      p7_StochasticTrace(ddef->r, dsq+ireg-1, Lr, om, fwd, ddef->tr);
      p7_trace_Index(ddef->tr);
      pos = 1;
      for (d = 0; d < ddef->tr->ndom; d++)
	{
	  p7_spensemble_Add(ddef->sp, t, ddef->tr->sqfrom[d]+ireg-1, ddef->tr->sqto[d]+ireg-1, ddef->tr->hmmfrom[d], ddef->tr->hmmto[d]);
	  p7_Null2_ByTrace(om, ddef->tr, ddef->tr->tfrom[d], ddef->tr->tto[d], wrk, null2);
	  for (; pos <= ddef->tr->sqfrom[d]; pos++) ddef->n2sc[ireg+pos-1] += 1.0;
	  for (; pos <= ddef->tr->sqto[d];   pos++) ddef->n2sc[ireg+pos-1] += null2[dsq[ireg+pos-1]];
	}
      for (; pos <= Lr; pos++)  ddef->n2sc[ireg+pos-1] += 1.0;
      p7_trace_Reuse(ddef->tr);        
    }
  for (pos = ireg; pos <= jreg; pos++)
    ddef->n2sc[pos] = logf(ddef->n2sc[pos] / (float) ddef->nsamples);
  p7_spensemble_Cluster(ddef->sp, ddef->min_overlap, ddef->of_smaller, ddef->max_diagdiff, ddef->min_posterior, ddef->min_endpointp, &nc);
  for (d = 0; d < nc; d++) 
    ddef->sp->assignment[d] = 0; 
  for (d = 0; d < nc; d++)
    {
      for (d2 = d+1; d2 < nc; d2++)
	{
	  nov = ESL_MIN(ddef->sp->sigc[d].j, ddef->sp->sigc[d2].j) - ESL_MAX(ddef->sp->sigc[d].i, ddef->sp->sigc[d2].i) + 1;
	  if (nov == 0) break;
	  n   = ESL_MIN(ddef->sp->sigc[d].j - ddef->sp->sigc[d].i + 1,  ddef->sp->sigc[d2].j - ddef->sp->sigc[d2].i + 1);
	  if ((float) nov / (float) n >= 0.8) 
	    {
	      if (ddef->sp->sigc[d].prob > ddef->sp->sigc[d2].prob) ddef->sp->assignment[d2] = 1;
	      else                                                  ddef->sp->assignment[d]  = 1;
	    }
	}
    }
  d = 0;
  for (d2 = 0; d2 < nc; d2++)
    {
      if (ddef->sp->assignment[d2]) continue; 
      if (d != d2) memcpy(ddef->sp->sigc + d, ddef->sp->sigc + d2, sizeof(struct p7_spcoord_s));
      d++;
    }
  ddef->sp->nc = d;
  *ret_nc = d;
  return eslOK;
}
static int region_trace_ensemble ( P7_DOMAINDEF * ddef , const P7_OPROFILE * om , const ESL_DSQ * dsq , int ireg , int jreg , const P7_OMX * fwd , P7_OMX * wrk , int * ret_nc ) { int Lr = jreg - ireg + 1 ; int t , d , d2 ; int nov , n ; int nc ; int pos ; float null2 [ p7_MAXCODE ] ; esl_vec_FSet ( ddef -> n2sc + ireg , Lr , 0.0 ) ; if ( ddef -> do_reseeding ) esl_randomness_Init ( ddef -> r , esl_randomness_GetSeed ( ddef -> r ) ) ; for ( t = 0 ; t < ddef -> nsamples ; t ++ ) { p7_StochasticTrace ( ddef -> r , dsq + ireg - 1 , Lr , om , fwd , ddef -> tr ) ; p7_trace_Index ( ddef -> tr ) ; pos = 1 ; for ( d = 0 ; d < ddef -> tr -> ndom ; d ++ ) { p7_spensemble_Add ( ddef -> sp , t , ddef -> tr -> sqfrom [ d ] + ireg - 1 , ddef -> tr -> sqto [ d ] + ireg - 1 , ddef -> tr -> hmmfrom [ d ] , ddef -> tr -> hmmto [ d ] ) ; p7_Null2_ByTrace ( om , ddef -> tr , ddef -> tr -> tfrom [ d ] , ddef -> tr -> tto [ d ] , wrk , null2 ) ; for ( ; pos <= ddef -> tr -> sqfrom [ d ] ; pos ++ ) ddef -> n2sc [ ireg + pos - 1 ] += 1.0 ; for ( ; pos <= ddef -> tr -> sqto [ d ] ; pos ++ ) ddef -> n2sc [ ireg + pos - 1 ] += null2 [ dsq [ ireg + pos - 1 ] ] ; } for ( ; pos <= Lr ; pos ++ ) ddef -> n2sc [ ireg + pos - 1 ] += 1.0 ; p7_trace_Reuse ( ddef -> tr ) ; } for ( pos = ireg ; pos <= jreg ; pos ++ ) ddef -> n2sc [ pos ] = logf ( ddef -> n2sc [ pos ] / ( float ) ddef -> nsamples ) ; p7_spensemble_Cluster ( ddef -> sp , ddef -> min_overlap , ddef -> of_smaller , ddef -> max_diagdiff , ddef -> min_posterior , ddef -> min_endpointp , & nc ) ; for ( d = 0 ; d < nc ; d ++ ) ddef -> sp -> assignment [ d ] = 0 ; for ( d = 0 ; d < nc ; d ++ ) { for ( d2 = d + 1 ; d2 < nc ; d2 ++ ) { nov = ESL_MIN ( ddef -> sp -> sigc [ d ] . j , ddef -> sp -> sigc [ d2 ] . j ) - ESL_MAX ( ddef -> sp -> sigc [ d ] . i , ddef -> sp -> sigc [ d2 ] . i ) + 1 ; if ( nov == 0 ) break ; n = ESL_MIN ( ddef -> sp -> sigc [ d ] . j - ddef -> sp -> sigc [ d ] . i + 1 , ddef -> sp -> sigc [ d2 ] . j - ddef -> sp -> sigc [ d2 ] . i + 1 ) ; if ( ( float ) nov / ( float ) n >= 0.8 ) { if ( ddef -> sp -> sigc [ d ] . prob > ddef -> sp -> sigc [ d2 ] . prob ) ddef -> sp -> assignment [ d2 ] = 1 ; else ddef -> sp -> assignment [ d ] = 1 ; } } } d = 0 ; for ( d2 = 0 ; d2 < nc ; d2 ++ ) { if ( ddef -> sp -> assignment [ d2 ] ) continue ; if ( d != d2 ) memcpy ( ddef -> sp -> sigc + d , ddef -> sp -> sigc + d2 , sizeof ( struct p7_spcoord_s ) ) ; d ++ ; } ddef -> sp -> nc = d ; * ret_nc = d ; return eslOK ; }
region_trace_ensemble()
 SRE, Fri Feb  8 11:49:44 2008 [Janelia]
region_trace_ensemble () SRE Fri Feb 8 11 : 49 : 44 2008 [ Janelia ]
/* zero the null2 scores in region */ /* By default, we make results reproducible by forcing a reset of
   * the RNG to its originally seeded state.
   */ /* Collect an ensemble of sampled traces; calculate null2 odds ratios from these */ /* residues outside domains get bumped +1: because f'(x) = f(x), so f'(x)/f(x) = 1 in these segments */ /* Residues inside domains get bumped by their null2 ratio */ /* the remaining residues in the region outside any domains get +1 */ /* Convert the accumulated n2sc[] ratios in this region to log odds null2 scores on each residue. */ /* Cluster the ensemble of traces to break region into envelopes. */ /* A little hacky now. Remove "dominated" domains relative to seq coords. */ /* overload <assignment> to flag that a domain is dominated */ /* who dominates who? (by post prob) */ /* overlap */ /* shrink the sigc list, removing dominated domains */ /* skip domain d2, it's dominated. */