hmmer-3.3/easel/miniapps/esl-alimanip.c
/* find_seqs_with_given_insert
 *                   
 * Given an MSA with RF annotation, determine which sequences have inserts after column <target>
 * of at least size <min> and at most size <max>. Fill an array <useme> of size msa->nseq with 
 * TRUE seq i has the insert, FALSE if it doesn't.
 */


Given an MSA with RF annotation, determine which sequences have inserts after column 
of at least size  and at most size . Fill an array  of size msa->nseq with
TRUE seq i has the insert, FALSE if it doesn't.
Given an MSA with RF annotation determine which sequences have inserts after column of at least size and at most size . Fill an array of size msa - > nseq with TRUE seq i has the insert FALSE if it doesn ' t .
static int find_seqs_with_given_insert(ESL_MSA *msa, int *i_am_rf, char *errbuf, int target, int min, int max, int **ret_useme)
{
  int status;
  int apos, rfpos, clen;
  int **ict;
  int *useme;
  int i;
  if(! (msa->flags & eslMSA_DIGITAL)) ESL_FAIL(eslEINVAL, errbuf, "in find_seqs_with_given_insert(), msa must be digitized.");
  if(msa->rf == NULL) ESL_FAIL(eslEINVAL, errbuf, "No #=GC RF markup in alignment, it is needed for --seq-ins.");
  ESL_ALLOC(useme,sizeof(int) * (msa->nseq));
  ESL_ALLOC(ict,  sizeof(int *) * (msa->alen+2));
  for(i = 0; i <= msa->alen; i++)
    {
      ESL_ALLOC(ict[i],  sizeof(int) * (msa->nseq));
      esl_vec_ISet(ict[i], (msa->nseq), 0);
    }
  rfpos = 0;
  for(apos = 1; apos <= msa->alen; apos++)
    {
      if(i_am_rf[apos-1]) { 
	rfpos++;
      }
      else 
	for(i = 0; i < msa->nseq; i++)
	  if(! esl_abc_XIsGap(msa->abc, msa->ax[i][apos])) { 
	    ict[rfpos][i]++;
	  }	  
    }
  clen = rfpos;
  if(target > clen) ESL_FAIL(eslEINVAL, errbuf, "--seq-ins <n> enabled with <n> = %d, but non-gap RF length of alignment is only %d columns.", target, clen);
  for(i = 0; i < msa->nseq; i++) { 
    useme[i] = ((ict[target][i] >= min) && (ict[target][i] <= max)) ?  TRUE : FALSE;
  }
  *ret_useme = useme;
  for(i = 0; i <= msa->alen; i++)
    free(ict[i]);
  free(ict);
  return eslOK;
 ERROR:
  return status;
}
static int find_seqs_with_given_insert ( ESL_MSA * msa , int * i_am_rf , char * errbuf , int target , int min , int max , int * * ret_useme ) { int status ; int apos , rfpos , clen ; int * * ict ; int * useme ; int i ; if ( ! ( msa -> flags & eslMSA_DIGITAL ) ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; if ( msa -> rf == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; ESL_ALLOC ( useme , sizeof ( int ) * ( msa -> nseq ) ) ; ESL_ALLOC ( ict , sizeof ( int * ) * ( msa -> alen + 2 ) ) ; for ( i = 0 ; i <= msa -> alen ; i ++ ) { ESL_ALLOC ( ict [ i ] , sizeof ( int ) * ( msa -> nseq ) ) ; esl_vec_ISet ( ict [ i ] , ( msa -> nseq ) , 0 ) ; } rfpos = 0 ; for ( apos = 1 ; apos <= msa -> alen ; apos ++ ) { if ( i_am_rf [ apos - 1 ] ) { rfpos ++ ; } else for ( i = 0 ; i < msa -> nseq ; i ++ ) if ( ! esl_abc_XIsGap ( msa -> abc , msa -> ax [ i ] [ apos ] ) ) { ict [ rfpos ] [ i ] ++ ; } } clen = rfpos ; if ( target > clen ) ESL_FAIL ( eslEINVAL , errbuf , " " , target , clen ) ; for ( i = 0 ; i < msa -> nseq ; i ++ ) { useme [ i ] = ( ( ict [ target ] [ i ] >= min ) && ( ict [ target ] [ i ] <= max ) ) ? TRUE : FALSE ; } * ret_useme = useme ; for ( i = 0 ; i <= msa -> alen ; i ++ ) free ( ict [ i ] ) ; free ( ict ) ; return eslOK ; ERROR : return status ; }
find_seqs_with_given_insert
find_seqs_with_given_insert
/* contract check */ /* printf("ict[target:%d][i:%d]: %d, min: %d max: %d\n", target, i, ict[target][i], min, max); */