hmmer-3.3/easel/miniapps/esl-alimanip.c
/* trim_msa
 *                   
 * Given an MSA and unaligned 'trimmed' versions (subsequences) of all seqs in that MSA, 
 * replace all chars that have been trimmed away (not in subsequences) with gaps in the MSA.
 * 
 * We remove all GR and GC markup from the msa, except for possibly #=GC RF if <do_keeprf> 
 * is TRUE, else we remove that too.
 */


Given an MSA and unaligned 'trimmed' versions (subsequences) of all seqs in that MSA,
replace all chars that have been trimmed away (not in subsequences) with gaps in the MSA.

We remove all GR and GC markup from the msa, except for possibly #=GC RF if 
is TRUE, else we remove that too.
Given an MSA and unaligned ' trimmed ' versions ( subsequences ) of all seqs in that MSA replace all chars that have been trimmed away ( not in subsequences ) with gaps in the MSA . We remove all GR and GC markup from the msa except for possibly # = GC RF if is TRUE else we remove that too .
static int trim_msa(ESL_MSA *msa, ESL_SQ **sq, int do_keeprf, char *errbuf)
{
  int status;
  int i, r;
  int apos, uapos;
  int astart,  aend;
  int uastart, uaend;
  char *offset;
  char *aseq;
  char *uaseq;
  char *uasubseq;
  int *a2ua_map;
  int *ua2a_map;
  int ualen;
  if(! (msa->flags & eslMSA_DIGITAL))
    ESL_FAIL(eslEINVAL, errbuf, "in trim_msa(), msa must be digitized.");
  ESL_ALLOC(aseq,  sizeof(char) * (msa->alen+1));
  for(i = 0; i < msa->nseq; i++)
    {
      if (sq[i]->dsq == NULL) ESL_FAIL(eslEINVAL, errbuf, "in trim_msa(), sq's must be digitized.");
      if (sq[i]->n   == 0)    ESL_FAIL(eslEINVAL, errbuf, "in trim_msa(), sq[%d] is zero-length\n", i);
      ESL_ALLOC(a2ua_map, sizeof(int) * (msa->alen+1));
      esl_vec_ISet(a2ua_map, (msa->alen+1), -1);
      uapos = apos = 1;
      while(apos <= msa->alen)
	{
	  while(apos <= msa->alen && esl_abc_XIsGap(msa->abc, msa->ax[i][apos])) apos++;
	  if(apos <= msa->alen) a2ua_map[apos] = uapos++;
	  apos++;
	}
      ualen = uapos;
      ESL_ALLOC(ua2a_map, sizeof(int) * (ualen+1));
      ua2a_map[0] = -1;
      for(apos = 1; apos <= msa->alen; apos++)
	if(a2ua_map[apos] != -1)
	  ua2a_map[a2ua_map[apos]] = apos;
      ESL_ALLOC(uasubseq, sizeof(char) * (sq[i]->n+1));
      esl_abc_Textize(msa->abc, sq[i]->dsq, sq[i]->n, uasubseq);
      esl_abc_Textize(msa->abc, msa->ax[i], msa->alen, aseq);
      esl_strdup(aseq, -1, &(uaseq));
      esl_strdealign(uaseq, uaseq, "-_.~", NULL);
      offset = strstr(uaseq, uasubseq); 
      if(offset == NULL) ESL_FAIL(eslEINVAL, errbuf, "in trim_msa(), sq[%d] is not a subseq of msa seq %d\n", i, i);
      uastart = offset  - uaseq + 1;
      uaend   = uastart + strlen(uasubseq) - 1;
      astart  = ua2a_map[uastart];
      aend    = ua2a_map[uaend];
      free(ua2a_map);
      free(a2ua_map);
      for(apos = 1;        apos <  astart;    apos++) msa->ax[i][apos] = msa->abc->K; 
      for(apos = aend + 1; apos <= msa->alen; apos++) msa->ax[i][apos] = msa->abc->K; 
      free(uaseq);
      free(uasubseq);
    }
  if(! do_keeprf && msa->rf != NULL) { free(msa->rf); msa->rf = NULL; }
  if(msa->ss_cons != NULL) { free(msa->ss_cons); msa->ss_cons = NULL; }
  if(msa->sa_cons != NULL) { free(msa->sa_cons); msa->sa_cons = NULL; }
  if(msa->pp_cons != NULL) { free(msa->pp_cons); msa->pp_cons = NULL; }
  if(msa->ss != NULL) { 
    for(i = 0; i < msa->nseq; i++) if(msa->ss[i] != NULL) { free(msa->ss[i]); }
    free(msa->ss); 
    msa->ss = NULL;
  }
  if(msa->sa != NULL) { 
    for(i = 0; i < msa->nseq; i++) if(msa->sa[i] != NULL) { free(msa->sa[i]); }
    free(msa->sa); 
    msa->sa = NULL;
  }
  if(msa->pp != NULL) { 
    for(i = 0; i < msa->nseq; i++) if(msa->pp[i] != NULL) { free(msa->pp[i]); }
    free(msa->pp); 
    msa->pp = NULL;
  }
  if(msa->ngr > 0) { 
    for(r = 0; r < msa->ngr; r++) { 
      for(i = 0; i < msa->nseq; i++) if(msa->gr[r][i] != NULL) { free(msa->gr[r][i]); }
      free(msa->gr[r]);
    }
    if(msa->gr_idx != NULL) { esl_keyhash_Destroy(msa->gr_idx); msa->gr_idx = NULL; }
    msa->gr_idx = NULL;
    free(msa->gr);
    msa->gr = NULL;
    msa->ngr = 0;
  }
  free(aseq);
  return eslOK;
 ERROR:
  return status;
}
static int trim_msa ( ESL_MSA * msa , ESL_SQ * * sq , int do_keeprf , char * errbuf ) { int status ; int i , r ; int apos , uapos ; int astart , aend ; int uastart , uaend ; char * offset ; char * aseq ; char * uaseq ; char * uasubseq ; int * a2ua_map ; int * ua2a_map ; int ualen ; if ( ! ( msa -> flags & eslMSA_DIGITAL ) ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; ESL_ALLOC ( aseq , sizeof ( char ) * ( msa -> alen + 1 ) ) ; for ( i = 0 ; i < msa -> nseq ; i ++ ) { if ( sq [ i ] -> dsq == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " " ) ; if ( sq [ i ] -> n == 0 ) ESL_FAIL ( eslEINVAL , errbuf , " \n " , i ) ; ESL_ALLOC ( a2ua_map , sizeof ( int ) * ( msa -> alen + 1 ) ) ; esl_vec_ISet ( a2ua_map , ( msa -> alen + 1 ) , -1 ) ; uapos = apos = 1 ; while ( apos <= msa -> alen ) { while ( apos <= msa -> alen && esl_abc_XIsGap ( msa -> abc , msa -> ax [ i ] [ apos ] ) ) apos ++ ; if ( apos <= msa -> alen ) a2ua_map [ apos ] = uapos ++ ; apos ++ ; } ualen = uapos ; ESL_ALLOC ( ua2a_map , sizeof ( int ) * ( ualen + 1 ) ) ; ua2a_map [ 0 ] = -1 ; for ( apos = 1 ; apos <= msa -> alen ; apos ++ ) if ( a2ua_map [ apos ] != -1 ) ua2a_map [ a2ua_map [ apos ] ] = apos ; ESL_ALLOC ( uasubseq , sizeof ( char ) * ( sq [ i ] -> n + 1 ) ) ; esl_abc_Textize ( msa -> abc , sq [ i ] -> dsq , sq [ i ] -> n , uasubseq ) ; esl_abc_Textize ( msa -> abc , msa -> ax [ i ] , msa -> alen , aseq ) ; esl_strdup ( aseq , -1 , & ( uaseq ) ) ; esl_strdealign ( uaseq , uaseq , " " , NULL ) ; offset = strstr ( uaseq , uasubseq ) ; if ( offset == NULL ) ESL_FAIL ( eslEINVAL , errbuf , " \n " , i , i ) ; uastart = offset - uaseq + 1 ; uaend = uastart + strlen ( uasubseq ) - 1 ; astart = ua2a_map [ uastart ] ; aend = ua2a_map [ uaend ] ; free ( ua2a_map ) ; free ( a2ua_map ) ; for ( apos = 1 ; apos < astart ; apos ++ ) msa -> ax [ i ] [ apos ] = msa -> abc -> K ; for ( apos = aend + 1 ; apos <= msa -> alen ; apos ++ ) msa -> ax [ i ] [ apos ] = msa -> abc -> K ; free ( uaseq ) ; free ( uasubseq ) ; } if ( ! do_keeprf && msa -> rf != NULL ) { free ( msa -> rf ) ; msa -> rf = NULL ; } if ( msa -> ss_cons != NULL ) { free ( msa -> ss_cons ) ; msa -> ss_cons = NULL ; } if ( msa -> sa_cons != NULL ) { free ( msa -> sa_cons ) ; msa -> sa_cons = NULL ; } if ( msa -> pp_cons != NULL ) { free ( msa -> pp_cons ) ; msa -> pp_cons = NULL ; } if ( msa -> ss != NULL ) { for ( i = 0 ; i < msa -> nseq ; i ++ ) if ( msa -> ss [ i ] != NULL ) { free ( msa -> ss [ i ] ) ; } free ( msa -> ss ) ; msa -> ss = NULL ; } if ( msa -> sa != NULL ) { for ( i = 0 ; i < msa -> nseq ; i ++ ) if ( msa -> sa [ i ] != NULL ) { free ( msa -> sa [ i ] ) ; } free ( msa -> sa ) ; msa -> sa = NULL ; } if ( msa -> pp != NULL ) { for ( i = 0 ; i < msa -> nseq ; i ++ ) if ( msa -> pp [ i ] != NULL ) { free ( msa -> pp [ i ] ) ; } free ( msa -> pp ) ; msa -> pp = NULL ; } if ( msa -> ngr > 0 ) { for ( r = 0 ; r < msa -> ngr ; r ++ ) { for ( i = 0 ; i < msa -> nseq ; i ++ ) if ( msa -> gr [ r ] [ i ] != NULL ) { free ( msa -> gr [ r ] [ i ] ) ; } free ( msa -> gr [ r ] ) ; } if ( msa -> gr_idx != NULL ) { esl_keyhash_Destroy ( msa -> gr_idx ) ; msa -> gr_idx = NULL ; } msa -> gr_idx = NULL ; free ( msa -> gr ) ; msa -> gr = NULL ; msa -> ngr = 0 ; } free ( aseq ) ; return eslOK ; ERROR : return status ; }
trim_msa
trim_msa
/* we'll replace the first occurence of uasubseq in uaseq */ /* make it a gap */ /* make it a gap */ /* Free all per-column annotation that might now be invalid */ /* Free all per-residue annotation (alternatively, we could just add gaps to the gaps we've created by 
   * trimming, but this would cause problems for SS annotation, for example. */