hmmer-3.3/src/p7_trace.c
/* Function:  p7_trace_FauxFromMSA()
 * Synopsis:  Create array of faux tracebacks from an existing MSA.
 * Incept:    SRE, Thu May 21 08:07:25 2009 [Janelia]
 *
 * Purpose:   Given an existing <msa> and an array <matassign> that
 *            flags the alignment columns that are assigned to consensus
 *            match states (matassign[1..alen] = 1|0); create an array
 *            of faux traces <tr[0..msa->nseq-1]>. <optflags> controls 
 *            optional behavior; it can be <p7_DEFAULT> or <p7_MSA_COORDS>,
 *            as explained below.
 *            
 *            The traces are core traces: they start/end with B/E,
 *            they may use I_0,I_M, and D_1 states. Any flanking
 *            insertions (outside the first/last consensus column) are
 *            assigned to I_0 and I_M.
 *            
 *            If the input alignment contains sequence fragments,
 *            caller should first convert leading/trailing gaps to
 *            missing data symbols. This hack causes entry/exit
 *            transitions to be encoded in the trace as B->X->{MDI}k
 *            and {MDI}k->X->E, rather than B->DDDD->Mk, Mk->DDDDD->E
 *            paths involving terminal deletions, and all functions
 *            that use traces, such as <p7_trace_Count()>, (should)
 *            ignore transitions involving <p7T_X> states.
 *            
 *            By default (<optflags = p7_DEFAULT>), the <i> coordinate
 *            in the faux tracebacks is <1..L>, relative to the
 *            unaligned raw sequences in <msa>, the way most H3 traces
 *            are supposed to be. In some cases (such as model
 *            construction from an MSA) it is convenient to reference
 *            residues in the MSA cooordinate system directly; setting
 *            <optflags = p7_MSA_COORDS> makes the traces come out
 *            with <i=1..alen> coords for residues.
 *            
 *            Important: an MSA may imply DI and ID transitions that
 *            are illegal in a core model. If the only purpose of the
 *            traces is to go straight back into alignment
 *            construction through a <p7_tracealign_*> function, this
 *            is ok, because the <p7_tracealign_*> routines can handle
 *            DI and ID transitions (enabling reconstruction of almost
 *            exactly the same input alignment, modulo unaligned
 *            insertions). This is what happens for <hmmalign
 *            --mapali>, for example. However, if the caller wants to
 *            use the traces for anything else, these illegal DI and
 *            ID transitions have to be removed first, and the caller
 *            should use <p7_trace_Doctor()> to do it.
 *
 * Args:      msa       - digital alignment
 *            matassign - flag for each alignment column, whether
 *                        it is consensus or not. matassign[1..alen] = 1|0; 
 *                        matassign[0] = 0
 *            optflags  - p7_DEFAULT | p7_MSA_COORDS 
 *            tr        - RETURN: caller provides 0..nseq-1 pointer 
 *                        array for holding returned traces.
 *
 * Returns:   <eslOK> on success, and tr[0..nseq-1] now point to newly
 *            created traces; caller is responsible for freeing these.
 *
 * Throws:    <eslEMEM> on allocation error.
 *
 * Xref:      J5/17: build.c::fake_tracebacks() becomes p7_trace_FauxFromMSA();
 *                   ability to handle MSA or raw coords added.
 */
p7_trace_FauxFromMSA()
Synopsis:  Create array of faux tracebacks from an existing MSA.

Given an existing  and an array  that
flags the alignment columns that are assigned to consensus
match states (matassign[1..alen] = 1|0); create an array
of faux traces nseq-1]>.  controls
optional behavior; it can be  or ,
as explained below.

The traces are core traces: they start/end with B/E,
they may use I_0,I_M, and D_1 states. Any flanking
insertions (outside the first/last consensus column) are
assigned to I_0 and I_M.

If the input alignment contains sequence fragments,
caller should first convert leading/trailing gaps to
missing data symbols.

By default (), the  coordinate
in the faux tracebacks is <1..L>, relative to the
unaligned raw sequences in , the way most H3 traces
are supposed to be. In some cases (such as model
construction from an MSA) it is convenient to reference
residues in the MSA cooordinate system directly; setting
 makes the traces come out
with  coords for residues.

an MSA may imply DI and ID transitions that
are illegal in a core model. If the only purpose of the
traces is to go straight back into alignment
construction through a  function, this
is ok, because the  routines can handle
DI and ID transitions (enabling reconstruction of almost
exactly the same input alignment, modulo unaligned
insertions). This is what happens for , for example. However, if the caller wants to
use the traces for anything else, these illegal DI and
ID transitions have to be removed first, and the caller
should use  to do it.

msa       - digital alignment
matassign - flag for each alignment column, whether
it is consensus or not. matassign[1..alen] = 1|0;
matassign[0] = 0
optflags  - p7_DEFAULT | p7_MSA_COORDS
tr        - RETURN: caller provides 0..nseq-1 pointer
array for holding returned traces.

 on success, and tr[0..nseq-1] now point to newly
created traces; caller is responsible for freeing these.

 on allocation error.

J5/17: build.c::fake_tracebacks() becomes p7_trace_FauxFromMSA();
ability to handle MSA or raw coords added.
p7_trace_FauxFromMSA () Synopsis : Create array of faux tracebacks from an existing MSA . Given an existing and an array that flags the alignment columns that are assigned to consensus match states ( matassign [ 1 .. alen ] = 1|0 ) ; create an array of faux traces nseq - 1 ] > . controls optional behavior ; it can be or as explained below . The traces are core traces : they start / end with B / E they may use I_0 I_M and D_1 states . Any flanking insertions ( outside the first / last consensus column ) are assigned to I_0 and I_M . If the input alignment contains sequence fragments caller should first convert leading / trailing gaps to missing data symbols . By default () the coordinate in the faux tracebacks is <1 .. L > relative to the unaligned raw sequences in the way most H3 traces are supposed to be . In some cases ( such as model construction from an MSA ) it is convenient to reference residues in the MSA cooordinate system directly ; setting makes the traces come out with coords for residues . an MSA may imply DI and ID transitions that are illegal in a core model . If the only purpose of the traces is to go straight back into alignment construction through a function this is ok because the routines can handle DI and ID transitions ( enabling reconstruction of almost exactly the same input alignment modulo unaligned insertions ) . This is what happens for for example . However if the caller wants to use the traces for anything else these illegal DI and ID transitions have to be removed first and the caller should use to do it . msa - digital alignment matassign - flag for each alignment column whether it is consensus or not . matassign [ 1 .. alen ] = 1|0 ; matassign [ 0 ] = 0 optflags - p7_DEFAULT | p7_MSA_COORDS tr - RETURN : caller provides 0 .. nseq - 1 pointer array for holding returned traces . on success and tr [ 0 .. nseq - 1 ] now point to newly created traces ; caller is responsible for freeing these . on allocation error . J5 / 17 : build . c :: fake_tracebacks () becomes p7_trace_FauxFromMSA () ; ability to handle MSA or raw coords added .
int
p7_trace_FauxFromMSA(ESL_MSA *msa, int *matassign, int optflags, P7_TRACE **tr)
{		      
  int  idx;			
  int  k;                       
  int  apos;                    
  int  rpos;			
  int  showpos;			
  int  status = eslOK;
  for (idx = 0; idx < msa->nseq; idx++) tr[idx] = NULL;
  for (idx = 0; idx < msa->nseq; idx++)
    {
      if ((tr[idx] = p7_trace_Create())                      == NULL) goto ERROR; 
      if ((status  = p7_trace_Append(tr[idx], p7T_B, 0, 0)) != eslOK) goto ERROR;
      for (k = 0, rpos = 1, apos = 1; apos <= msa->alen; apos++)
	{
	  showpos = (optflags & p7_MSA_COORDS) ? apos : rpos;
	  if (matassign[apos]) 
	    {			
	      k++;
	      if (esl_abc_XIsResidue(msa->abc, msa->ax[idx][apos])) 
		status = p7_trace_Append(tr[idx], p7T_M, k, showpos);
	      else if (esl_abc_XIsGap    (msa->abc, msa->ax[idx][apos])) 
		status = p7_trace_Append(tr[idx], p7T_D, k, 0);          
	      else if (esl_abc_XIsNonresidue(msa->abc, msa->ax[idx][apos]))
		status = p7_trace_Append(tr[idx], p7T_M, k, showpos); 
	      else if (esl_abc_XIsMissing(msa->abc, msa->ax[idx][apos]))
		{
		  if (tr[idx]->st[tr[idx]->N-1] != p7T_X)
		    status = p7_trace_Append(tr[idx], p7T_X, k, 0); 
		}
	      else ESL_XEXCEPTION(eslEINCONCEIVABLE, "can't happen");
	    }
	  else
	    { 			
	      if (esl_abc_XIsResidue(msa->abc, msa->ax[idx][apos]))
		status = p7_trace_Append(tr[idx], p7T_I, k, showpos);
	      else if (esl_abc_XIsNonresidue(msa->abc, msa->ax[idx][apos]))
		status = p7_trace_Append(tr[idx], p7T_I, k, showpos); 
	      else if (esl_abc_XIsMissing(msa->abc, msa->ax[idx][apos]))
		{ 
		  if (tr[idx]->st[tr[idx]->N-1] != p7T_X)
		    status = p7_trace_Append(tr[idx], p7T_X, k, 0);
		}
	      else if (! esl_abc_XIsGap(msa->abc, msa->ax[idx][apos]))
		ESL_XEXCEPTION(eslEINCONCEIVABLE, "can't happen");
	    }
	  if (esl_abc_XIsResidue(msa->abc, msa->ax[idx][apos])) rpos++; 
	  if (status != eslOK) goto ERROR;
	}
      if ((status = p7_trace_Append(tr[idx], p7T_E, 0, 0)) != eslOK) goto ERROR;
      tr[idx]->M = k;
      tr[idx]->L = msa->alen;
    }
  return eslOK;
 ERROR:
  for (idx = 0; idx < msa->nseq; idx++) { p7_trace_Destroy(tr[idx]); tr[idx] = NULL; }
  return status; 
}
int p7_trace_FauxFromMSA ( ESL_MSA * msa , int * matassign , int optflags , P7_TRACE * * tr ) { int idx ; int k ; int apos ; int rpos ; int showpos ; int status = eslOK ; for ( idx = 0 ; idx < msa -> nseq ; idx ++ ) tr [ idx ] = NULL ; for ( idx = 0 ; idx < msa -> nseq ; idx ++ ) { if ( ( tr [ idx ] = p7_trace_Create ( ) ) == NULL ) goto ERROR ; if ( ( status = p7_trace_Append ( tr [ idx ] , p7T_B , 0 , 0 ) ) != eslOK ) goto ERROR ; for ( k = 0 , rpos = 1 , apos = 1 ; apos <= msa -> alen ; apos ++ ) { showpos = ( optflags & p7_MSA_COORDS ) ? apos : rpos ; if ( matassign [ apos ] ) { k ++ ; if ( esl_abc_XIsResidue ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) status = p7_trace_Append ( tr [ idx ] , p7T_M , k , showpos ) ; else if ( esl_abc_XIsGap ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) status = p7_trace_Append ( tr [ idx ] , p7T_D , k , 0 ) ; else if ( esl_abc_XIsNonresidue ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) status = p7_trace_Append ( tr [ idx ] , p7T_M , k , showpos ) ; else if ( esl_abc_XIsMissing ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) { if ( tr [ idx ] -> st [ tr [ idx ] -> N - 1 ] != p7T_X ) status = p7_trace_Append ( tr [ idx ] , p7T_X , k , 0 ) ; } else ESL_XEXCEPTION ( eslEINCONCEIVABLE , " " ) ; } else { if ( esl_abc_XIsResidue ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) status = p7_trace_Append ( tr [ idx ] , p7T_I , k , showpos ) ; else if ( esl_abc_XIsNonresidue ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) status = p7_trace_Append ( tr [ idx ] , p7T_I , k , showpos ) ; else if ( esl_abc_XIsMissing ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) { if ( tr [ idx ] -> st [ tr [ idx ] -> N - 1 ] != p7T_X ) status = p7_trace_Append ( tr [ idx ] , p7T_X , k , 0 ) ; } else if ( ! esl_abc_XIsGap ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) ESL_XEXCEPTION ( eslEINCONCEIVABLE , " " ) ; } if ( esl_abc_XIsResidue ( msa -> abc , msa -> ax [ idx ] [ apos ] ) ) rpos ++ ; if ( status != eslOK ) goto ERROR ; } if ( ( status = p7_trace_Append ( tr [ idx ] , p7T_E , 0 , 0 ) ) != eslOK ) goto ERROR ; tr [ idx ] -> M = k ; tr [ idx ] -> L = msa -> alen ; } return eslOK ; ERROR : for ( idx = 0 ; idx < msa -> nseq ; idx ++ ) { p7_trace_Destroy ( tr [ idx ] ) ; tr [ idx ] = NULL ; } return status ; }
Function:  p7_trace_FauxFromMSA()
 Synopsis:  Create array of faux tracebacks from an existing MSA.
Function : p7_trace_FauxFromMSA () Synopsis : Create array of faux tracebacks from an existing MSA .
/* counter over seqs in MSA */ /* position in HMM                 */ /* position in alignment columns 1..alen */ /* position in unaligned sequence residues 1..L */ /* coord to actually record: apos or rpos */ /* match or delete */ /* treat * as a residue! */ /* allow only one X in a row */ /* insert or nothing */ /* treat * as a residue! */ /* k == M by construction; set tr->L = msa->alen since coords are w.r.t. ax */