hmmer-3.3/easel/esl_msaweight.c
/* consensus_by_sample()
 * Use a statistical sample of seqs to define consensus columns.
 * 
 * On deep alignments, it's usually overkill to look at all sequences
 * just to determine which columns to call consensus. We can speed the
 * decision up by looking at a subsample instead.
 * 
 * We want results to be reproducible; by default, the random number
 * generator (used for taking the sample) uses a fixed seed.
 * 
 * We have to watch out for a pathological case where there could be a
 * ridiculous number of fragments piled up on just one piece of the
 * alignment. This can happen with some DNA repeat elements, for
 * example. If the sample contains too many fragments, reject this
 * approach, the caller will have to determine consensus on all
 * sequences instead.
 * 
 * In:   cfg      - optional configuration options, or NULL to use all defaults
 *       msa      - MSA to determine consensus for
 *       ct       - allocated space for [(0).1..alen][0..Kp-1] observed symbol counts; contents irrelevant
 *       conscols - allocated space for up to <alen> consensus column indices
 *       
 * Out:  ct       - [apos=1..alen][a=0..Kp-1] are counts of observed symbol <a> in column <alen> in sample. [0][] are all 0's.
 *       conscols - [0..ncons-1] list of consensus column indices 1..alen
 *       ncons    - number of consensus columns in <conscols>, or 0 if consensus determination was rejected or failed
 *       
 * Returns <eslOK> on success.
 * 
 * Returns <eslFAIL> if there were too many fragments. Now <ct> still has the counts
 * observed in the sample, but <conscols> has no valid indices and <ncons> is 0.      
 *
 * Throws <eslEMEM> on allocation failure.
 */
consensus_by_sample()
Use a statistical sample of seqs to define consensus columns.

On deep alignments, it's usually overkill to look at all sequences
just to determine which columns to call consensus. We can speed the
decision up by looking at a subsample instead.

We want results to be reproducible; by default, the random number
generator (used for taking the sample) uses a fixed seed.

We have to watch out for a pathological case where there could be a
ridiculous number of fragments piled up on just one piece of the
alignment. This can happen with some DNA repeat elements, for
example. If the sample contains too many fragments, reject this
approach, the caller will have to determine consensus on all
sequences instead.



ct       - [apos=1..alen][a=0..Kp-1] are counts of observed symbol  in column  in sample.

Returns  on success.

Returns  if there were too many fragments. Now  still has the counts
observed in the sample, but  has no valid indices and  is 0.

Throws  on allocation failure.
consensus_by_sample () Use a statistical sample of seqs to define consensus columns . On deep alignments it ' s usually overkill to look at all sequences just to determine which columns to call consensus . We can speed the decision up by looking at a subsample instead . We want results to be reproducible ; by default the random number generator ( used for taking the sample ) uses a fixed seed . We have to watch out for a pathological case where there could be a ridiculous number of fragments piled up on just one piece of the alignment . This can happen with some DNA repeat elements for example . If the sample contains too many fragments reject this approach the caller will have to determine consensus on all sequences instead . ct - [ apos = 1 .. alen ] [ a = 0 .. Kp - 1 ] are counts of observed symbol in column in sample . Returns on success . Returns if there were too many fragments . Now still has the counts observed in the sample but has no valid indices and is 0 . Throws on allocation failure .
static int
consensus_by_sample(const ESL_MSAWEIGHT_CFG *cfg, const ESL_MSA *msa, int **ct, int *conscols, int *ret_ncons, ESL_MSAWEIGHT_DAT *dat)
{
  float       fragthresh  = (cfg? cfg->fragthresh : eslMSAWEIGHT_FRAGTHRESH);
  float       symfrac     = (cfg? cfg->symfrac    : eslMSAWEIGHT_SYMFRAC);
  int         nsamp       = (cfg? cfg->nsamp      : eslMSAWEIGHT_NSAMP);
  int         maxfrag     = (cfg? cfg->maxfrag    : eslMSAWEIGHT_MAXFRAG);
  ESL_RAND64 *rng         = (cfg? esl_rand64_Create(cfg->seed) : esl_rand64_Create(eslMSAWEIGHT_RNGSEED));   
  int64_t    *sampidx     = NULL;
  int         nfrag       = 0;        
  int         ncons       = 0;        
  int         tot;                    
  int         minspan;                 
  int         apos, a;
  int         lpos, rpos;
  int         i, idx;
  int         status      = eslOK;
  ESL_ALLOC(sampidx, sizeof(int64_t) * nsamp);
  esl_mat_ISet(ct, msa->alen+1, msa->abc->Kp, 0);
  if (dat) dat->seed = esl_rand64_GetSeed(rng);
  esl_rand64_Deal(rng, nsamp, (int64_t) msa->nseq, sampidx);  
  minspan = (int) ceil( fragthresh * (float) msa->alen );     
  for (i = 0; i < nsamp; i++)
    {
      idx = (int) sampidx[i];
      for (lpos = 1;         lpos <= msa->alen; lpos++) if (esl_abc_XIsResidue(msa->abc, msa->ax[idx][lpos])) break;
      for (rpos = msa->alen; rpos >= 1;         rpos--) if (esl_abc_XIsResidue(msa->abc, msa->ax[idx][rpos])) break;
      if  (rpos - lpos + 1 < minspan) nfrag++; else { lpos = 1; rpos = msa->alen; }   
      for (apos = lpos; apos <= rpos; apos++)
	ct[apos][msa->ax[idx][apos]]++;
    }
  if (dat) dat->samp_nfrag = nfrag;
  if (nfrag <= maxfrag)
    {
      for (apos = 1; apos <= msa->alen; apos++)
	{
	  for (tot = 0, a = 0; a < msa->abc->Kp-2; a++) tot += ct[apos][a];  
	  if ( ((float) ct[apos][msa->abc->K] / (float) tot) < symfrac)      
	    conscols[ncons++] = apos;
	}
      if (dat) dat->cons_by_sample = TRUE;
    }
  else
    {
      if (dat) dat->rejected_sample = TRUE;
      status = eslFAIL;
    }
 ERROR:
  free(sampidx);
  esl_rand64_Destroy(rng);
  *ret_ncons = ncons;   
  return status;
}
static int consensus_by_sample ( const ESL_MSAWEIGHT_CFG * cfg , const ESL_MSA * msa , int * * ct , int * conscols , int * ret_ncons , ESL_MSAWEIGHT_DAT * dat ) { float fragthresh = ( cfg ? cfg -> fragthresh : eslMSAWEIGHT_FRAGTHRESH ) ; float symfrac = ( cfg ? cfg -> symfrac : eslMSAWEIGHT_SYMFRAC ) ; int nsamp = ( cfg ? cfg -> nsamp : eslMSAWEIGHT_NSAMP ) ; int maxfrag = ( cfg ? cfg -> maxfrag : eslMSAWEIGHT_MAXFRAG ) ; ESL_RAND64 * rng = ( cfg ? esl_rand64_Create ( cfg -> seed ) : esl_rand64_Create ( eslMSAWEIGHT_RNGSEED ) ) ; int64_t * sampidx = NULL ; int nfrag = 0 ; int ncons = 0 ; int tot ; int minspan ; int apos , a ; int lpos , rpos ; int i , idx ; int status = eslOK ; ESL_ALLOC ( sampidx , sizeof ( int64_t ) * nsamp ) ; esl_mat_ISet ( ct , msa -> alen + 1 , msa -> abc -> Kp , 0 ) ; if ( dat ) dat -> seed = esl_rand64_GetSeed ( rng ) ; esl_rand64_Deal ( rng , nsamp , ( int64_t ) msa -> nseq , sampidx ) ; minspan = ( int ) ceil ( fragthresh * ( float ) msa -> alen ) ; for ( i = 0 ; i < nsamp ; i ++ ) { idx = ( int ) sampidx [ i ] ; for ( lpos = 1 ; lpos <= msa -> alen ; lpos ++ ) if ( esl_abc_XIsResidue ( msa -> abc , msa -> ax [ idx ] [ lpos ] ) ) break ; for ( rpos = msa -> alen ; rpos >= 1 ; rpos -- ) if ( esl_abc_XIsResidue ( msa -> abc , msa -> ax [ idx ] [ rpos ] ) ) break ; if ( rpos - lpos + 1 < minspan ) nfrag ++ ; else { lpos = 1 ; rpos = msa -> alen ; } for ( apos = lpos ; apos <= rpos ; apos ++ ) ct [ apos ] [ msa -> ax [ idx ] [ apos ] ] ++ ; } if ( dat ) dat -> samp_nfrag = nfrag ; if ( nfrag <= maxfrag ) { for ( apos = 1 ; apos <= msa -> alen ; apos ++ ) { for ( tot = 0 , a = 0 ; a < msa -> abc -> Kp - 2 ; a ++ ) tot += ct [ apos ] [ a ] ; if ( ( ( float ) ct [ apos ] [ msa -> abc -> K ] / ( float ) tot ) < symfrac ) conscols [ ncons ++ ] = apos ; } if ( dat ) dat -> cons_by_sample = TRUE ; } else { if ( dat ) dat -> rejected_sample = TRUE ; status = eslFAIL ; } ERROR : free ( sampidx ) ; esl_rand64_Destroy ( rng ) ; * ret_ncons = ncons ; return status ; }
consensus_by_sample()
 Use a statistical sample of seqs to define consensus columns.
consensus_by_sample () Use a statistical sample of seqs to define consensus columns .
// fixed seed for default reproducibility // number of fragments in sample // number of consensus columns defined // total # of residues+gaps in a column // <sampidx> is now an ordered list of <nsamp> indices in range 0..nseq-1 // define alispan as aligned length from first to last non-gap. If alispan < minspan, define seq as a fragment // i.e. <tot> symbols over esl_abc_XIsResidue() || esl_abc_XIsGap(); inclusive of degeneracies, exclusive of missing | nonresidue. // This is the rule for determining a consensus column: if the frequency of residues/(residues+gaps) >= symfrac  // will be 0 if we saw too many fragments.