hmmer-3.3/easel/miniapps/esl-alimanip.c
/* individualize_consensus
 *                   
 * Given an MSA with a consensus structure impose it to create
 * individual secondary structures. Simple rule, for consensus
 * bp i,j if seq positions i and j are both non-gaps seq i,j are 
 * paired, if >= 1 is a gap, they're not paired.
 */


Given an MSA with a consensus structure impose it to create
individual secondary structures. Simple rule, for consensus
bp i,j if seq positions i and j are both non-gaps seq i,j are
paired, if >= 1 is a gap, they're not paired.
Given an MSA with a consensus structure impose it to create individual secondary structures . Simple rule for consensus bp i j if seq positions i and j are both non - gaps seq i j are paired if > = 1 is a gap they ' re not paired .
static int
individualize_consensus(const ESL_GETOPTS *go, char *errbuf, ESL_MSA *msa)
{
  int64_t apos;
  int     i;
  int    *cct = NULL;		   
  int    *ct  = NULL;		   
  char   *ss  = NULL;              
  char   *ss_cons_nopseudo = NULL; 
  int     status;
  if(msa->ss_cons == NULL)                                ESL_FAIL(eslEINVAL, errbuf, "--sindi requires MSA to have consensus structure annotation.\n");
  if(! (msa->flags & eslMSA_DIGITAL))                     ESL_FAIL(eslEINVAL, errbuf, "individualize_consensus() MSA is not digitized.\n");
  ESL_ALLOC(cct, sizeof(int)  * (msa->alen+1));
  ESL_ALLOC(ct,  sizeof(int)  * (msa->alen+1));
  ESL_ALLOC(ss,  sizeof(char) * (msa->alen+1));
  ESL_ALLOC(ss_cons_nopseudo, sizeof(char) * (msa->alen+1));
  esl_wuss_nopseudo(msa->ss_cons, ss_cons_nopseudo);
  if (esl_wuss2ct(ss_cons_nopseudo, msa->alen, cct) != eslOK) ESL_FAIL(eslEINVAL, errbuf, "Consensus structure string is inconsistent.");
  for (i = 0; i < msa->nseq; i++)
    {
      esl_vec_ICopy(cct, (msa->alen+1), ct);
      for (apos = 1; apos <= msa->alen; apos++)
	if (esl_abc_XIsGap(msa->abc, msa->ax[i][apos]))
	  { 
	    if (ct[apos] != 0)  ct[ct[apos]] = 0;
	    ct[apos] = 0;
	  }
      if (esl_ct2wuss(ct, msa->alen, ss) != eslOK) ESL_FAIL(eslEINVAL, errbuf, "Unexpected error converting de-knotted bp ct array to wuss notation.");
      esl_msa_AppendGR(msa, "SS", i, ss);
    }
  free(cct);
  free(ct);
  free(ss);
  free(ss_cons_nopseudo);
  return eslOK;
 ERROR:
  if (cct)               free(cct);
  if (ct)                free(ct);
  if (ss)                free(ss);
  if (ss_cons_nopseudo)  free(ss_cons_nopseudo);
  return status;
}
static int individualize_consensus ( const ESL_GETOPTS * go , char * errbuf , ESL_MSA * msa ) { int64_t apos ; int i ; int * cct = NULL ; int * ct = NULL ; char * ss = NULL ; char * ss_cons_nopseudo = NULL ; int status ; if ( msa -> ss_cons == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " \n " ) ; if ( ! ( msa -> flags & eslMSA_DIGITAL ) ) ESL_FAIL ( eslEINVAL , errbuf , " \n " ) ; ESL_ALLOC ( cct , sizeof ( int ) * ( msa -> alen + 1 ) ) ; ESL_ALLOC ( ct , sizeof ( int ) * ( msa -> alen + 1 ) ) ; ESL_ALLOC ( ss , sizeof ( char ) * ( msa -> alen + 1 ) ) ; ESL_ALLOC ( ss_cons_nopseudo , sizeof ( char ) * ( msa -> alen + 1 ) ) ; esl_wuss_nopseudo ( msa -> ss_cons , ss_cons_nopseudo ) ; if ( esl_wuss2ct ( ss_cons_nopseudo , msa -> alen , cct ) != eslOK ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; for ( i = 0 ; i < msa -> nseq ; i ++ ) { esl_vec_ICopy ( cct , ( msa -> alen + 1 ) , ct ) ; for ( apos = 1 ; apos <= msa -> alen ; apos ++ ) if ( esl_abc_XIsGap ( msa -> abc , msa -> ax [ i ] [ apos ] ) ) { if ( ct [ apos ] != 0 ) ct [ ct [ apos ] ] = 0 ; ct [ apos ] = 0 ; } if ( esl_ct2wuss ( ct , msa -> alen , ss ) != eslOK ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; esl_msa_AppendGR ( msa , " " , i , ss ) ; } free ( cct ) ; free ( ct ) ; free ( ss ) ; free ( ss_cons_nopseudo ) ; return eslOK ; ERROR : if ( cct ) free ( cct ) ; if ( ct ) free ( ct ) ; if ( ss ) free ( ss ) ; if ( ss_cons_nopseudo ) free ( ss_cons_nopseudo ) ; return status ; }
individualize_consensus
individualize_consensus
/* 0..alen-1 base pair partners array for consensus        */ /* 0..alen-1 base pair partners array for current sequence */ /* individual secondary structure we've built              */ /* no-pseudoknot version of consensus structure            */ /* go through each position of each sequence, 
     if it's a gap and it is part of a base pair, remove that base pair */ /* convert to WUSS SS string and append to MSA */