hmmer-3.3/easel/esl_msaweight.c
/* collect_counts()
  * Collect a matrix of observed symbol counts in each column: ct[apos=1..alen][a=0..Kp-1].
  *   - use HMMER's fragment-definition rule, and don't count external gaps of fragments.
  *   - This matrix is O(KL) for alignment length L, alphabet size K. 
  *     It's responsible for the asymptotic memory complexity of PB weights.
  *   - "What if I've already marked fragments", you say (maybe in hmmbuild), "can I
  *     save some time here?" Not really. Even if you know a seq is a fragment, you'd
  *     need to run the lpos and rpos loops to find its start/end.
  *   - If we already know what the consensus columns are, only collect counts in them,
  *     leaving counts in nonconsensus columns zero. This is a time optimization.
  */
collect_counts()
Collect a matrix of observed symbol counts in each column: ct[apos=1..alen][a=0..Kp-1].
- use HMMER's fragment-definition rule, and don't count external gaps of fragments.
- This matrix is O(KL) for alignment length L, alphabet size K.
It's responsible for the asymptotic memory complexity of PB weights.
- "What if I've already marked fragments", you say (maybe in hmmbuild), "can I
save some time here?" Not really. Even if you know a seq is a fragment, you'd
need to run the lpos and rpos loops to find its start/end.
- If we already know what the consensus columns are, only collect counts in them,
leaving counts in nonconsensus columns zero. This is a time optimization.
collect_counts () Collect a matrix of observed symbol counts in each column : ct [ apos = 1 .. alen ] [ a = 0 .. Kp - 1 ] . - use HMMER ' s fragment - definition rule and don ' t count external gaps of fragments . - This matrix is O ( KL ) for alignment length L alphabet size K . It ' s responsible for the asymptotic memory complexity of PB weights . - " What if I ' ve already marked fragments " you say ( maybe in hmmbuild ) " can I save some time here? " Not really . Even if you know a seq is a fragment you ' d need to run the lpos and rpos loops to find its start / end . - If we already know what the consensus columns are only collect counts in them leaving counts in nonconsensus columns zero . This is a time optimization .
static int
collect_counts(const ESL_MSAWEIGHT_CFG *cfg, const ESL_MSA *msa, const int *conscols, int ncons, int **ct, ESL_MSAWEIGHT_DAT *dat)
{
  float fragthresh  = (cfg? cfg->fragthresh : eslMSAWEIGHT_FRAGTHRESH);     
  int   minspan     = (int) ceil( fragthresh * (float) msa->alen );         
  int   lpos, rpos;     
  int   idx, apos, j;        
  esl_mat_ISet(ct, msa->alen+1, msa->abc->Kp, 0);
  for (idx = 0; idx < msa->nseq; idx++)
    {
      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) { lpos = 1; rpos = msa->alen; } else if (dat) dat->all_nfrag++;   
      if (ncons) 
	{
	  for (j = 0; j < ncons && conscols[j] <= rpos; j++)
	    {
	      apos = conscols[j];
	      if (apos < lpos) continue;
	      ct[apos][msa->ax[idx][apos]]++;
	    }
	}
      else      
	{
	  for (apos = lpos; apos <= rpos; apos++)
	    ct[apos][msa->ax[idx][apos]]++;
	}
    }
  return eslOK;
}
static int collect_counts ( const ESL_MSAWEIGHT_CFG * cfg , const ESL_MSA * msa , const int * conscols , int ncons , int * * ct , ESL_MSAWEIGHT_DAT * dat ) { float fragthresh = ( cfg ? cfg -> fragthresh : eslMSAWEIGHT_FRAGTHRESH ) ; int minspan = ( int ) ceil ( fragthresh * ( float ) msa -> alen ) ; int lpos , rpos ; int idx , apos , j ; esl_mat_ISet ( ct , msa -> alen + 1 , msa -> abc -> Kp , 0 ) ; for ( idx = 0 ; idx < msa -> nseq ; idx ++ ) { 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 ) { lpos = 1 ; rpos = msa -> alen ; } else if ( dat ) dat -> all_nfrag ++ ; if ( ncons ) { for ( j = 0 ; j < ncons && conscols [ j ] <= rpos ; j ++ ) { apos = conscols [ j ] ; if ( apos < lpos ) continue ; ct [ apos ] [ msa -> ax [ idx ] [ apos ] ] ++ ; } } else { for ( apos = lpos ; apos <= rpos ; apos ++ ) ct [ apos ] [ msa -> ax [ idx ] [ apos ] ] ++ ; } } return eslOK ; }
collect_counts()
 Collect a matrix of observed symbol counts in each column: ct[apos=1..alen][a=0..Kp-1].
collect_counts () Collect a matrix of observed symbol counts in each column : ct [ apos = 1 .. alen ] [ a = 0 .. Kp - 1 ] .
// seq is fragment if (length from 1st to last aligned residue)/alen < fragthresh (i.e. span < minspan) // precalculated span length threshold using <fragthresh> // leftmost, rightmost aligned residue (1..alen) // HMMER mark_fragments() rule. Count "span" from first to last aligned residue. If alispan/alen < fragthresh, it's a fragment. // L=0 seq or alen=0? then lpos == msa->alen+1, rpos == 0 => lpos > rpos. rpos-lpos-1 <= 0 but test below still works. // full len seqs count cols 1..alen; fragments only count lpos..rpos. // if we have consensus columns already, only count symbols in those columns (faster)... // ... else, count symbols in all columns.