hmmer-3.3/easel/esl_msaweight.c
/* Function:  esl_msaweight_IDFilter_adv()
 * Synopsis:  Advanced version of %id filtering; digital MSAs only.
 * Incept:    SRE, Sat 30 Mar 2019
 *
 * Purpose:   Same as <esl_msaweight_IDFilter()> but customizable with
 *            optional custom parameters in <cfg>. Requires digital
 *            mode alignment.
 *            
 *            In particular, <cfg> can be used to select two other
 *            rules for sequence preference, besides the default
 *            "conscover" rule: random preference, or preferring the
 *            lower-index sequence (same as the rule used in text mode
 *            alignments).
 *            
 *            For the "conscover" rule, consensus column determination
 *            can be customized the same way as in PB weighting.
 */
esl_msaweight_IDFilter_adv()
Synopsis:  Advanced version of %id filtering; digital MSAs only.
Incept:    SRE, Sat 30 Mar 2019

Same as  but customizable with
optional custom parameters in . Requires digital
mode alignment.

In particular,  can be used to select two other
rules for sequence preference, besides the default
"conscover" rule: random preference, or preferring the
lower-index sequence (same as the rule used in text mode
alignments).

For the "conscover" rule, consensus column determination
can be customized the same way as in PB weighting.
esl_msaweight_IDFilter_adv () Synopsis : Advanced version of %id filtering ; digital MSAs only . Incept : SRE Sat 30 Mar 2019 Same as but customizable with optional custom parameters in . Requires digital mode alignment . In particular can be used to select two other rules for sequence preference besides the default " conscover " rule : random preference or preferring the lower - index sequence ( same as the rule used in text mode alignments ) . For the " conscover " rule consensus column determination can be customized the same way as in PB weighting .
int
esl_msaweight_IDFilter_adv(const ESL_MSAWEIGHT_CFG *cfg, const ESL_MSA *msa, double maxid, ESL_MSA **ret_newmsa)
{
  int     ignore_rf   = (cfg? cfg->ignore_rf  : eslMSAWEIGHT_IGNORE_RF);      
  int     allow_samp  = (cfg? cfg->allow_samp : eslMSAWEIGHT_ALLOW_SAMP);     
  int     sampthresh  = (cfg? cfg->sampthresh : eslMSAWEIGHT_SAMPTHRESH);     
  int     filterpref  = (cfg? cfg->filterpref : eslMSAWEIGHT_FILT_CONSCOVER); 
  int   **ct          = NULL;     
  int    *conscols    = NULL;     
  double *sortwgt     = NULL;     
  int    *ranked_at   = NULL;     
  int    *list        = NULL;     
  int    *useme       = NULL;     
  int     ncons       = 0;        
  int     nnew        = 0;        
  int     apos, r, i;             
  double  ident;                  
  int     status      = eslOK;
  ESL_DASSERT1(( msa->nseq >= 1 && msa->alen >= 1));
  ESL_DASSERT1(( msa->flags & eslMSA_DIGITAL ));
  ESL_ALLOC(sortwgt,   sizeof(double) * msa->nseq);
  ESL_ALLOC(ranked_at, sizeof(double) * msa->nseq);
  ESL_ALLOC(list,      sizeof(int)    * msa->nseq);
  ESL_ALLOC(useme,     sizeof(int)    * msa->nseq);
  esl_vec_ISet(useme, msa->nseq, 0); 
  if (filterpref == eslMSAWEIGHT_FILT_CONSCOVER)
    {
      ct = esl_mat_ICreate( msa->alen+1, msa->abc->Kp );      
      ESL_ALLOC(conscols,  sizeof(int)    * msa->alen);
      if      (! ignore_rf && msa->rf)                consensus_by_rf(msa, conscols, &ncons, NULL);
      else if (allow_samp  && msa->nseq > sampthresh) consensus_by_sample(cfg, msa, ct, conscols, &ncons, NULL);
      else {
	collect_counts(cfg, msa, conscols, ncons, ct, NULL);
	consensus_by_all(cfg, msa, ct, conscols, &ncons, NULL);
      }
      if (!ncons) {
	for (apos = 1; apos <= msa->alen; apos++) conscols[apos-1] = apos; 
	ncons = msa->alen; 
      }
    }
  if      (filterpref == eslMSAWEIGHT_FILT_CONSCOVER) set_preference_conscover(msa, conscols, ncons, sortwgt);
  else if (filterpref == eslMSAWEIGHT_FILT_RANDOM)    set_preference_randomly (cfg, msa->nseq,       sortwgt);
  else if (filterpref == eslMSAWEIGHT_FILT_ORIGORDER) set_preference_origorder(     msa->nseq,       sortwgt);
  else esl_fatal("bad filterpref");
  esl_quicksort(sortwgt, msa->nseq, sort_doubles_decreasing, ranked_at);
  for (r = 0; r < msa->nseq; r++)  
    {
      for (i = 0; i < nnew; i++)   
	{
	  if ((status = esl_dst_XPairId(msa->abc, msa->ax[ranked_at[r]], msa->ax[list[i]], &ident, NULL, NULL)) != eslOK) goto ERROR;
	  if (ident >= maxid) break;
	}
      if (i == nnew)  
	{
	  list[nnew++]        = ranked_at[r];  
	  useme[ranked_at[r]] = TRUE;
	}
    }
  if ((status = esl_msa_SequenceSubset(msa, useme, ret_newmsa)) != eslOK) goto ERROR;
 ERROR:
  free(useme);
  free(list);
  free(ranked_at);
  free(sortwgt);
  free(conscols);
  esl_mat_IDestroy(ct);
  return status;
}
int esl_msaweight_IDFilter_adv ( const ESL_MSAWEIGHT_CFG * cfg , const ESL_MSA * msa , double maxid , ESL_MSA * * ret_newmsa ) { int ignore_rf = ( cfg ? cfg -> ignore_rf : eslMSAWEIGHT_IGNORE_RF ) ; int allow_samp = ( cfg ? cfg -> allow_samp : eslMSAWEIGHT_ALLOW_SAMP ) ; int sampthresh = ( cfg ? cfg -> sampthresh : eslMSAWEIGHT_SAMPTHRESH ) ; int filterpref = ( cfg ? cfg -> filterpref : eslMSAWEIGHT_FILT_CONSCOVER ) ; int * * ct = NULL ; int * conscols = NULL ; double * sortwgt = NULL ; int * ranked_at = NULL ; int * list = NULL ; int * useme = NULL ; int ncons = 0 ; int nnew = 0 ; int apos , r , i ; double ident ; int status = eslOK ; ESL_DASSERT1 ( ( msa -> nseq >= 1 && msa -> alen >= 1 ) ) ; ESL_DASSERT1 ( ( msa -> flags & eslMSA_DIGITAL ) ) ; ESL_ALLOC ( sortwgt , sizeof ( double ) * msa -> nseq ) ; ESL_ALLOC ( ranked_at , sizeof ( double ) * msa -> nseq ) ; ESL_ALLOC ( list , sizeof ( int ) * msa -> nseq ) ; ESL_ALLOC ( useme , sizeof ( int ) * msa -> nseq ) ; esl_vec_ISet ( useme , msa -> nseq , 0 ) ; if ( filterpref == eslMSAWEIGHT_FILT_CONSCOVER ) { ct = esl_mat_ICreate ( msa -> alen + 1 , msa -> abc -> Kp ) ; ESL_ALLOC ( conscols , sizeof ( int ) * msa -> alen ) ; if ( ! ignore_rf && msa -> rf ) consensus_by_rf ( msa , conscols , & ncons , NULL ) ; else if ( allow_samp && msa -> nseq > sampthresh ) consensus_by_sample ( cfg , msa , ct , conscols , & ncons , NULL ) ; else { collect_counts ( cfg , msa , conscols , ncons , ct , NULL ) ; consensus_by_all ( cfg , msa , ct , conscols , & ncons , NULL ) ; } if ( ! ncons ) { for ( apos = 1 ; apos <= msa -> alen ; apos ++ ) conscols [ apos - 1 ] = apos ; ncons = msa -> alen ; } } if ( filterpref == eslMSAWEIGHT_FILT_CONSCOVER ) set_preference_conscover ( msa , conscols , ncons , sortwgt ) ; else if ( filterpref == eslMSAWEIGHT_FILT_RANDOM ) set_preference_randomly ( cfg , msa -> nseq , sortwgt ) ; else if ( filterpref == eslMSAWEIGHT_FILT_ORIGORDER ) set_preference_origorder ( msa -> nseq , sortwgt ) ; else esl_fatal ( " " ) ; esl_quicksort ( sortwgt , msa -> nseq , sort_doubles_decreasing , ranked_at ) ; for ( r = 0 ; r < msa -> nseq ; r ++ ) { for ( i = 0 ; i < nnew ; i ++ ) { if ( ( status = esl_dst_XPairId ( msa -> abc , msa -> ax [ ranked_at [ r ] ] , msa -> ax [ list [ i ] ] , & ident , NULL , NULL ) ) != eslOK ) goto ERROR ; if ( ident >= maxid ) break ; } if ( i == nnew ) { list [ nnew ++ ] = ranked_at [ r ] ; useme [ ranked_at [ r ] ] = TRUE ; } } if ( ( status = esl_msa_SequenceSubset ( msa , useme , ret_newmsa ) ) != eslOK ) goto ERROR ; ERROR : free ( useme ) ; free ( list ) ; free ( ranked_at ) ; free ( sortwgt ) ; free ( conscols ) ; esl_mat_IDestroy ( ct ) ; return status ; }
Function:  esl_msaweight_IDFilter_adv()
 Synopsis:  Advanced version of %id filtering; digital MSAs only.
Function : esl_msaweight_IDFilter_adv () Synopsis : Advanced version of %id filtering ; digital MSAs only .
// default is FALSE: use RF annotation as consensus definition, if RF is present // default is TRUE: allow subsampling speed optimization // if nseq > sampthresh, try to determine consensus on a subsample of seqs // default preference rule is "conscover" // matrix of symbol counts in each column. ct[apos=(0).1..alen][a=0..Kp-1] // list of consensus column indices [0..ncons-1] // when pair of seqs is >= maxid, retain seq w/ higher <sortwgt> // sorted seq indices 0..nseq-1 // array of seq indices in new msa  // useme[i] is TRUE if seq[i] is kept in new msa  // number of consensus column indices in <conscols> list // how many seqs have been added to <list> and <useme> so far // index over columns, ranked seq indices, seqs in <list> // pairwise fractional id (0..1) /* Allocations that we always need*/ /* initialize array */ /* Determine consensus columns with same procedure as advanced PB weights.
   *   1. If RF annotation is provided, use it.
   *         (...unless cfg->ignore_rf has been set TRUE)
   *   2. If it's a deep alignment, subsample and apply fragthresh/symfrac rules to the sample.
   *         (...unless cfg->allow_samp has been set FALSE)
   *   3. Otherwise, use standard rule on all sequences: 
   *      collect observed counts in <ct> (applying fragthresh rule) and
   *      determine consensus (applying symfrac rule)
   *   4. If all else fails, define all columns as consensus.
   *
   * We only need to do this if we're using the default "conscover" rule.
   */ /* allocations only needed for consensus determination */ // (0).1..alen; 0..Kp-1 /* Assign preferences for which sequences to retain.
   * Higher <sortwgt> means higher preference.
   */ /* Sort sequence indices by these preferences.
   * <ranked_at[]> is a ranked list of sequence indices 0..nseq-1.
   */ /* Determine which seqs will be kept, favoring highest ranked ones.
   */ // Try to include each sequence, starting with highest ranked ones // Test it against all the seqs we've already put in the list. // if the seq made it past comparisons with all <nnew> current seqs in the list... // ... add it to the growing list. /* Filter the input MSA.
   */