hmmer-3.3/easel/miniapps/esl-alimanip.c
/* Function: msa_remove_truncated_seqs()
 * 
 * Purpose:  Remove sequences in MSA that have all gaps in the first <ntrunc> 5' leading 
 *           non-gap RF columns OR the last <ntrunc> 3' leading non-gap RF columns
 */


Remove sequences in MSA that have all gaps in the first  5' leading
non-gap RF columns OR the last  3' leading non-gap RF columns
Remove sequences in MSA that have all gaps in the first 5 ' leading non - gap RF columns OR the last 3 ' leading non - gap RF columns
static int
msa_remove_truncated_seqs(ESL_MSA *msa, char *errbuf, int ntrunc, int *i_am_rf, ESL_MSA **ret_new_msa)
{
  int  status;
  int *useme;
  int  i;
  int  leading_okay, trailing_okay;
  int  apos, rfpos_ct;
  int  nused = 0;
  ESL_MSA *new_msa;
  if(! (msa->flags & eslMSA_DIGITAL)) ESL_FAIL(eslEINVAL, errbuf, "in msa_remove_truncated_seqs(), msa must be digitized.");
  if(msa->rf == NULL) ESL_FAIL(eslEINVAL, errbuf, "No #=GC RF markup in alignment, it is needed for --detrunc.");
  if(i_am_rf == NULL) ESL_FAIL(eslEINVAL, errbuf, "internal error, msa_remove_truncated_seq() i_am_rf is NULL.");
  ESL_ALLOC(useme, sizeof(int) * msa->nseq);
  for(i = 0; i < msa->nseq; i++) { 
    leading_okay  = FALSE;
    rfpos_ct = 0; 
    apos = 1;
    while(!leading_okay && (rfpos_ct < ntrunc) && (apos <= msa->alen)) { 
      if(i_am_rf[(apos-1)]) { 
	rfpos_ct++;
	if(! esl_abc_XIsGap(msa->abc, msa->ax[i][apos])) leading_okay = TRUE;
      }
      apos++;
    }
    trailing_okay = FALSE;
    rfpos_ct = 0;
    apos = msa->alen;
    while(!trailing_okay && (rfpos_ct < ntrunc) && (apos >= 1)) { 
      if(i_am_rf[(apos-1)]) { 
	rfpos_ct++;
	if(! esl_abc_XIsGap(msa->abc, msa->ax[i][apos])) trailing_okay = TRUE;
      }
      apos--;
    }
    useme[i] = (leading_okay && trailing_okay) ? TRUE : FALSE;
    if(useme[i]) nused++;
  }
  if(nused == 0) ESL_FAIL(eslEINVAL, errbuf, "--detrunc removed ALL sequences!");
  if((status = esl_msa_SequenceSubset(msa, useme, &new_msa)) != eslOK) esl_fatal("esl_msa_SequenceSubset() had a problem.");
  free(useme);
  *ret_new_msa = new_msa;
  return eslOK;
 ERROR:
  ESL_FAIL(status, errbuf, "msa_remove_truncated_seqs(): memory allocation error.");
  return eslOK; 
}
static int msa_remove_truncated_seqs ( ESL_MSA * msa , char * errbuf , int ntrunc , int * i_am_rf , ESL_MSA * * ret_new_msa ) { int status ; int * useme ; int i ; int leading_okay , trailing_okay ; int apos , rfpos_ct ; int nused = 0 ; ESL_MSA * new_msa ; if ( ! ( msa -> flags & eslMSA_DIGITAL ) ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; if ( msa -> rf == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; if ( i_am_rf == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; ESL_ALLOC ( useme , sizeof ( int ) * msa -> nseq ) ; for ( i = 0 ; i < msa -> nseq ; i ++ ) { leading_okay = FALSE ; rfpos_ct = 0 ; apos = 1 ; while ( ! leading_okay && ( rfpos_ct < ntrunc ) && ( apos <= msa -> alen ) ) { if ( i_am_rf [ ( apos - 1 ) ] ) { rfpos_ct ++ ; if ( ! esl_abc_XIsGap ( msa -> abc , msa -> ax [ i ] [ apos ] ) ) leading_okay = TRUE ; } apos ++ ; } trailing_okay = FALSE ; rfpos_ct = 0 ; apos = msa -> alen ; while ( ! trailing_okay && ( rfpos_ct < ntrunc ) && ( apos >= 1 ) ) { if ( i_am_rf [ ( apos - 1 ) ] ) { rfpos_ct ++ ; if ( ! esl_abc_XIsGap ( msa -> abc , msa -> ax [ i ] [ apos ] ) ) trailing_okay = TRUE ; } apos -- ; } useme [ i ] = ( leading_okay && trailing_okay ) ? TRUE : FALSE ; if ( useme [ i ] ) nused ++ ; } if ( nused == 0 ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; if ( ( status = esl_msa_SequenceSubset ( msa , useme , & new_msa ) ) != eslOK ) esl_fatal ( " " ) ; free ( useme ) ; * ret_new_msa = new_msa ; return eslOK ; ERROR : ESL_FAIL ( status , errbuf , " " ) ; return eslOK ; }
Function: msa_remove_truncated_seqs()
Function : msa_remove_truncated_seqs ()
/* contract check */ /* if ALL of the first 5' <ntrunc> non-gap RF columns are gaps in this seq, we'll remove it */ /* NEVERREACHED */