hmmer-3.3/easel/esl_msaweight.c
/* Function:  esl_msaweight_BLOSUM()
 * Synopsis:  BLOSUM weights.
 * Incept:    SRE, Sun Nov  5 09:52:41 2006 [Janelia]
 *
 * Purpose:   Given a multiple sequence alignment <msa> and an identity
 *            threshold <maxid>, calculate sequence weights using the
 *            BLOSUM algorithm (Henikoff and Henikoff, PNAS
 *            89:10915-10919, 1992). These weights are stored
 *            internally in the <msa> object, replacing any weights
 *            that may have already been there. Weights are $\geq 0$
 *            and they sum to <msa->nseq>.
 *            
 *            The algorithm does a single linkage clustering by
 *            fractional id, defines clusters such that no two clusters
 *            have a pairwise link $\geq$ <maxid>), and assigns
 *            weights of $\frac{1}{M_i}$ to each of the $M_i$
 *            sequences in each cluster $i$. The <maxid> threshold
 *            is a fractional pairwise identity, in the range
 *            $0..1$.
 *            
 *            The <msa> may be in either digitized or text mode.
 *            Digital mode is preferred, so that the pairwise identity
 *            calculations deal with degenerate residue symbols
 *            properly.
 *
 * Returns:   <eslOK> on success, and the weights inside <msa> have been
 *            modified. 
 *            
 * Throws:    <eslEMEM> on allocation error. <eslEINVAL> if a pairwise
 *            identity calculation fails because of corrupted sequence 
 *            data. In either case, the <msa> is unmodified.
 *
 * Xref:      [Henikoff92]; squid::weight.c::BlosumWeights().
 */
esl_msaweight_BLOSUM()
Synopsis:  BLOSUM weights.

Given a multiple sequence alignment  and an identity
threshold , calculate sequence weights using the
BLOSUM algorithm (Henikoff and Henikoff, PNAS
89:10915-10919, 1992). These weights are stored
internally in the  object, replacing any weights
that may have already been there. Weights are $\geq 0$
and they sum to nseq>.



The  may be in either digitized or text mode.
Digital mode is preferred, so that the pairwise identity
calculations deal with degenerate residue symbols
properly.

 on success, and the weights inside  have been
modified.

 on allocation error.  if a pairwise
identity calculation fails because of corrupted sequence
data. In either case, the  is unmodified.


esl_msaweight_BLOSUM () Synopsis : BLOSUM weights . Given a multiple sequence alignment and an identity threshold calculate sequence weights using the BLOSUM algorithm ( Henikoff and Henikoff PNAS 89 : 10915 - 10919 1992 ) . These weights are stored internally in the object replacing any weights that may have already been there . Weights are $ \ geq 0$ and they sum to nseq > . The may be in either digitized or text mode . Digital mode is preferred so that the pairwise identity calculations deal with degenerate residue symbols properly . on success and the weights inside have been modified . on allocation error . if a pairwise identity calculation fails because of corrupted sequence data . In either case the is unmodified .
int
esl_msaweight_BLOSUM(ESL_MSA *msa, double maxid)
{
  int  *c    = NULL; 
  int  *nmem = NULL; 
  int   nc;	     
  int   i;           
  int   status;
  ESL_DASSERT1( (maxid >= 0. && maxid <= 1.) );
  ESL_DASSERT1( (msa->nseq >= 1) );
  ESL_DASSERT1( (msa->alen >= 1) );
  if (msa->nseq == 1) { msa->wgt[0] = 1.0; return eslOK; }
  if ((status = esl_msacluster_SingleLinkage(msa, maxid, &c, NULL, &nc)) != eslOK) goto ERROR;
  ESL_ALLOC(nmem, sizeof(int) * nc);
  esl_vec_ISet(nmem, nc, 0);
  for (i = 0; i < msa->nseq; i++) nmem[c[i]]++;
  for (i = 0; i < msa->nseq; i++) msa->wgt[i] = 1. / (double) nmem[c[i]];
  esl_vec_DNorm(msa->wgt, msa->nseq);
  esl_vec_DScale(msa->wgt, msa->nseq, (double) msa->nseq);	
  msa->flags |= eslMSA_HASWGTS;
  free(nmem);
  free(c);
  return eslOK;
 ERROR:
  if (c    != NULL) free(c);
  if (nmem != NULL) free(nmem);
  return status;
}
int esl_msaweight_BLOSUM ( ESL_MSA * msa , double maxid ) { int * c = NULL ; int * nmem = NULL ; int nc ; int i ; int status ; ESL_DASSERT1 ( ( maxid >= 0. && maxid <= 1. ) ) ; ESL_DASSERT1 ( ( msa -> nseq >= 1 ) ) ; ESL_DASSERT1 ( ( msa -> alen >= 1 ) ) ; if ( msa -> nseq == 1 ) { msa -> wgt [ 0 ] = 1.0 ; return eslOK ; } if ( ( status = esl_msacluster_SingleLinkage ( msa , maxid , & c , NULL , & nc ) ) != eslOK ) goto ERROR ; ESL_ALLOC ( nmem , sizeof ( int ) * nc ) ; esl_vec_ISet ( nmem , nc , 0 ) ; for ( i = 0 ; i < msa -> nseq ; i ++ ) nmem [ c [ i ] ] ++ ; for ( i = 0 ; i < msa -> nseq ; i ++ ) msa -> wgt [ i ] = 1. / ( double ) nmem [ c [ i ] ] ; esl_vec_DNorm ( msa -> wgt , msa -> nseq ) ; esl_vec_DScale ( msa -> wgt , msa -> nseq , ( double ) msa -> nseq ) ; msa -> flags |= eslMSA_HASWGTS ; free ( nmem ) ; free ( c ) ; return eslOK ; ERROR : if ( c != NULL ) free ( c ) ; if ( nmem != NULL ) free ( nmem ) ; return status ; }
Function:  esl_msaweight_BLOSUM()
 Synopsis:  BLOSUM weights.
Function : esl_msaweight_BLOSUM () Synopsis : BLOSUM weights .
/* cluster assignments for each sequence */ /* number of seqs in each cluster */ /* number of clusters  */ /* loop counter */ /* Contract checks
   */ /* Make weights normalize up to nseq, and return.
   */