text: /* @source msbar application
  **
  ** Mutate sequence beyond all recognition
  **
  ** @author Copyright (C) Gary Williams (gwilliam@hgmp.mrc.ac.uk)
  ** @@
  **
  ** This program is free software; you can redistribute it and/or
  ** modify it under the terms of the GNU General Public License
  ** as published by the Free Software Foundation; either version 2
  ** of the License, or (at your option) any later version.
  **
  ** This program is distributed in the hope that it will be useful,
  ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  ** GNU General Public License for more details.
  **
  ** You should have received a copy of the GNU General Public License
  ** along with this program; if not, write to the Free Software
  ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  ******************************************************************************/
  
  #include "emboss.h"
  #include <ctype.h>	/* for tolower, toupper */
  
  
  
  
  static void msbar_blockmutn(AjPStr *str, AjBool isnuc,
  			    AjPStr const *blocklist,
  			    ajint min, ajint max, AjBool inframe);
  static void msbar_codonmutn(AjPStr *str, AjBool isnuc,
  			    AjPStr const *codonlist,
  			    AjBool inframe);
  static void msbar_pointmutn(AjPStr *str, AjBool isnuc,
  			    AjPStr const *pointlist);
  static void msbar_Insert(AjPStr *str, AjBool isnuc, ajint start, ajint end);
  static void msbar_Move(AjPStr *str, ajint start, ajint end, ajint destination);
  static void msbar_Duplicate(AjPStr *str, ajint start, ajint end);
  static AjBool msbar_Unlike(const AjPStr str, AjPSeqall other);
  
  
  
  
  /* @prog msbar ****************************************************************
  **
  ** Mutate sequence beyond all recognition
  **
  ******************************************************************************/
  
  int main(int argc, char **argv)
  {
  
      AjPSeqall seqall;
      AjPSeqall other;
      AjPSeqout seqout;
      AjPSeq seq = NULL;
  
      AjBool isnuc;
      AjPStr str = NULL;
  
      AjPStr *pointlist;
      AjPStr *codonlist;
      AjPStr *blocklist;
  
      ajint i;
      ajint count;
      ajint min;
      ajint max;
      AjBool inframe;
  
      /*
      ** the number of times to try to produce a sequence that is different to
      **  the 'other' sequences before we give up and output it anyway
      */
      ajint attempts = 10;
  
      /* count of tries to get a unique sequence */
      ajint try;
  
      embInit("msbar", argc, argv);
  
      seqall = ajAcdGetSeqall("sequence");
      seqout = ajAcdGetSeqoutall("outseq");
  
      pointlist = ajAcdGetList("point");
      codonlist = ajAcdGetList("codon");
      blocklist = ajAcdGetList("block");
  
      count   = ajAcdGetInt("count");
      inframe = ajAcdGetBoolean("inframe");
      min     = ajAcdGetInt("minimum");
      max     = ajAcdGetInt("maximum");
      other   = ajAcdGetSeqall("othersequence");
  
  
      str = ajStrNew();
  
      /* seed the random number generator */
      ajRandomSeed();
  
      while(ajSeqallNext(seqall, &seq))
      {
          ajSeqTrim(seq);
          
  	/* is this a protein or nucleic sequence? */
  	isnuc = ajSeqIsNuc(seq);
  
          /*
  	** try mutating until there is a result that is not the same as
  	** the 'other' sequences
  	*/
          for(try=0; try<attempts; try++)
  	{
  	    /* get a copy of the sequence string */
  	    ajStrAssignS(&str, ajSeqGetSeqS(seq));
  
  	    /* do the mutation operations */
  	    for(i=0; i<count; i++)
  	    {
  	        msbar_blockmutn(&str, isnuc, blocklist, min, max, inframe);
  	        msbar_codonmutn(&str, isnuc, codonlist, inframe);
  	        msbar_pointmutn(&str, isnuc, pointlist);
  	    }
  
              /* Is this mutated sequence different to the 'other' sequences?  */
  	    if(msbar_Unlike(str, other))
  	        break;
  	}
  
  	if(try >= attempts)
  	    ajWarn("No unique mutation found after %d attempts", attempts);
  
  	ajSeqAssignSeqS(seq, str);
  	ajSeqoutWriteSeq(seqout, seq);
      }
  
      ajSeqoutClose(seqout);
  
      ajStrDel(&str);
  
      ajSeqallDel(&seqall);
      ajSeqallDel(&other);
      ajSeqDel(&seq);
      ajSeqoutDel(&seqout);
      ajStrDelarray(&pointlist);
      ajStrDelarray(&codonlist);
      ajStrDelarray(&blocklist);
  
      embExit();
  
      return 0;
  }
  
  
  
  
  /* @funcstatic msbar_blockmutn ************************************************
  **
  ** Mutate a random block of sequence
  **
  ** @param [u] str [AjPStr*] sequence to mutate
  ** @param [r] isnuc [AjBool] TRUE if sequence is nucleic
  ** @param [r] blocklist [AjPStr const *] Types of block mutations to perform
  ** @param [r] min [ajint] minimum size of block
  ** @param [r] max [ajint] maximum size of block
  ** @param [r] inframe [AjBool] mutate blocks preserving codon frame if TRUE
  ** @@
  ******************************************************************************/
  
  static void msbar_blockmutn(AjPStr *str, AjBool isnuc, AjPStr const *blocklist,
  			    ajint min, ajint max, AjBool inframe)
  {
      ajint i = -1;
      ajint rposstart;
      ajint rposend;
      ajint rpos2;
      ajint opt;
      double xlen;
      double xdiff;
      double td;
      double tda;
      double tdb;
      ajint imax = max;
      ajint imin = min;
      
      if(inframe)
          if(imin < 3) imin = 3;
      if(imax < imin) imax = imin;
      xlen = ajStrGetLen(*str)-1;
  
      while(blocklist[++i])
      {
  	ajDebug("Next block mutation operation = %S\n", blocklist[i]);
  
  	/* None */
  	if(!ajStrCmpC(blocklist[i], "0"))
  	    return;
  
  	/* get the option value */
  	ajStrToInt(blocklist[i], &opt);
  
  	/* if 'Any' required, then choose which one */
  	if(opt == 1)
  	{
  	    td = ajRandomDouble() * 5.;
  	    opt = (ajint) td;
  	    opt += 2;
  	}
  
  	/*
  	**  get random block start and end positions in the sequence
  	**  (0 to ajStrLen - 1)
  	*/
  	if(inframe)
  	{
  	    tda = ajRandomDouble() * xlen/3.0;
  	    rposstart = (ajint) tda;
  	    xdiff = imax - imin;
  	    tdb = ajRandomDouble() * xdiff/3.0;
  	    rposend = rposstart + (imin/3) + (ajint)tdb - 1;
  	    rposstart *= 3;
  	    rposend *= 3;
  	    rposend--;
  	}
  	else
  	{
  	    tda = ajRandomDouble() * xlen;
  	    rposstart = (ajint) tda;
  	    xdiff = imax - imin;
  	    tdb = ajRandomDouble() * xdiff;
  	    rposend = rposstart + imin + (ajint) tdb - 1;
  	}
  
  	if(opt == 2)
  	{
  	    /* Insert */
  	    ajDebug("block insert from %d to %d\n", rposstart, rposend);
  	    msbar_Insert(str, isnuc, rposstart, rposend);
  	}
  
  	if(opt == 3)
  	{
  	    /* Delete */
  	    ajDebug("block deletion from %d to %d\n", rposstart,
  			   rposend);
  	    ajStrCutRange(str, rposstart, rposend);
  	}
  
  	if(opt == 4)
  	{
  	    /* Change */
  	    ajDebug("block change from %d to %d\n", rposstart, rposend);
  	    ajStrCutRange(str, rposstart, rposend);
  	    msbar_Insert(str, isnuc, rposstart, rposend);
  	}
  
  	if(opt == 5)
  	{
  	    /* Duplication */
  	    ajDebug("block duplication from %d to %d\n", rposstart,
  			   rposend);
  	    msbar_Duplicate(str, rposstart, rposend);
  	}
  
  	if(opt == 6)
  	{
  	    /* Move */
  	    if(inframe)
  	    {
  		xlen = ajStrGetLen(*str);
  		td = ajRandomDouble() * xlen/3.0;
  		rpos2 = (ajint) td;
  		rpos2 *= 3;
  	    }
  	    else
  	    {
  		xlen = ajStrGetLen(*str);
  		td = ajRandomDouble() * xlen;
  		rpos2 = (ajint) td;
  	    }
  	    ajDebug("block move from %d to %d to position %d\n",
  			   rposstart, rposend, rpos2);
  	    msbar_Move(str, rposstart, rposend, rpos2);
  	}
      }
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_codonmutn ************************************************
  **
  ** Mutate codons
  **
  ** @param [u] str [AjPStr*] Sequence to mutate
  ** @param [r] isnuc [AjBool] TRUE if sequence is nucleic
  ** @param [r] codonlist [AjPStr const *] Types of codon mutations to perform
  ** @param [r] inframe [AjBool] mutate blocks preserving codon frame if TRUE
  ** @@
  ******************************************************************************/
  
  static void msbar_codonmutn(AjPStr *str, AjBool isnuc, AjPStr const *codonlist,
  			    AjBool inframe)
  {
      ajint rpos;
      ajint rpos2;
      ajint i = -1;
      ajint opt;
      double xlen;
      double td;
      
      while(codonlist[++i])
      {
  	ajDebug("Next codon mutation operation = %S\n", codonlist[i]);
  
  	/* None */
  	if(!ajStrCmpC(codonlist[i], "0"))
  	    return;
  
  	/* get the option value */
  	ajStrToInt(codonlist[i], &opt);
  
  	/* if 'Any' required, then choose which one */
  	if(opt == 1)
  	{
  	    td = ajRandomDouble() * 5.;
  	    opt = (ajint) td;
  	    opt += 2;
  	}
  
  
  	/* get a random position in the sequence (0 to ajStrLen - 1) */
  	if(inframe)
  	{
  	    xlen = ajStrGetLen(*str);
  	    td = ajRandomDouble() * xlen/3.0;
  	    rpos = (ajint) td;
  	    rpos *= 3;
  	}
  	else
  	{
  	    xlen = ajStrGetLen(*str);
  	    td = ajRandomDouble() * xlen;
  	    rpos = (ajint) td;
  	}
  
  	if(opt == 2)
  	{
  	    /* Insert */
  	    ajDebug("codon insert at %d\n", rpos);
  	    msbar_Insert(str, isnuc, rpos, rpos+2);
  	}
  
  	if(opt == 3)
  	{
  	    /* Delete */
  	    ajDebug("codon deletion at %d\n", rpos);
  	    ajStrCutRange(str, rpos, rpos+2);
  	}
  
  	if(opt == 4)
  	{
  	    /* Change */
  	    ajDebug("codon change at %d\n", rpos);
  	    ajStrCutRange(str, rpos, rpos+2);
  	    msbar_Insert(str, isnuc, rpos, rpos+2);
  	}
  
  	if(opt == 5)
  	{
  	    /* Duplication */
  	    ajDebug("codon duplication at %d\n", rpos);
  	    msbar_Duplicate(str, rpos, rpos+2);
  
  	}
  
  	if(opt == 6)
  	{
  	    /* Move */
  	    if(inframe)
  	    {
  		xlen = ajStrGetLen(*str);
  		td = ajRandomDouble() * xlen/3.0;
  		rpos2 = (ajint) td;
  		rpos2 *= 3;
  	    }
  	    else
  	    {
  		xlen = ajStrGetLen(*str);
  		td = ajRandomDouble() * xlen;
  		rpos2 = (ajint) td;
  	    }
  	    ajDebug("codon move from %d to %d\n", rpos, rpos2);
  	    msbar_Move(str, rpos, rpos+2, rpos2);
  	}
  
  
      }
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_pointmutn ************************************************
  **
  ** Mutate random single points
  **
  ** @param [u] str [AjPStr*] sequence to mutate
  ** @param [r] isnuc [AjBool] TRUE if sequence is nucleic
  ** @param [r] pointlist [AjPStr const *] Types of point mutations to perform
  ** @@
  ******************************************************************************/
  
  static void msbar_pointmutn(AjPStr *str, AjBool isnuc, AjPStr const *pointlist)
  {
      ajint i = -1;
      ajint rpos, rpos2;
      ajint opt;
      double xlen;
      double td;
      
      while(pointlist[++i])
      {
  	ajDebug("Next point mutation operation = %S\n", pointlist[i]);
  	/* None */
  	if(!ajStrCmpC(pointlist[i], "0"))
  	    return;
  
  	/* get the option value */
  	ajStrToInt(pointlist[i], &opt);
  
  	/* if 'Any' required, then choose which one */
  	if(opt == 1)
  	{
  	    td = ajRandomDouble() * 5.;
  	    opt = (ajint) td;
  	    opt += 2;
  	}
  
  	/* get a random position in the sequence (0 to ajStrLen - 1) */
  	xlen = ajStrGetLen(*str);
  	td = ajRandomDouble() * xlen;
  	rpos = (ajint) td;
  	
  	if(opt == 2)
  	{
  	    /* Insert */
  	    ajDebug("Point insert at %d\n", rpos);
  	    msbar_Insert(str, isnuc, rpos, rpos);
  
  	}
  
  	if(opt == 3)
  	{
  	    /* Delete */
  	    ajDebug("Point deletion at %d\n", rpos);
  	    ajStrCutRange(str, rpos, rpos);
  
  	}
  
  	if(opt == 4)
  	{
  	    /* Change */
  	    ajDebug("Point change at %d\n", rpos);
  	    ajStrCutRange(str, rpos, rpos);
  	    msbar_Insert(str, isnuc, rpos, rpos);
  
  	}
  
  	if(opt == 5)
  	{
  	    /* Duplication */
  	    ajDebug("Point duplication at %d\n", rpos);
  	    msbar_Duplicate(str, rpos, rpos);
  
  	}
  
  	if(opt == 6)
  	{
  	    /* Move */
  	    xlen = ajStrGetLen(*str);
  	    td = ajRandomDouble() * xlen;
  	    rpos2 = (ajint) td;
  	    ajDebug("Point move from %d to %d\n", rpos, rpos2);
  	    msbar_Move(str, rpos, rpos, rpos2);
  	}
      }
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_Insert ***************************************************
  **
  ** Insert random sequence at a position in the main sequence
  **
  ** @param [u] str [AjPStr*] sequence to insert into
  ** @param [r] isnuc [AjBool] TRUE if sequence is nucleic
  ** @param [r] start [ajint] start position of insert
  ** @param [r] end [ajint] end of insert
  ** @@
  ******************************************************************************/
  
  static void msbar_Insert(AjPStr *str, AjBool isnuc, ajint start, ajint end)
  {
      char nuc[] = "ACGT";
      char prot[] = "ARNDCQEGHILKMFPSTWYV";
      AjPStr ins = NULL;
      ajint count;
      ajint r;
      double xlen;
      double td;
      
      count = end - start +1;
  
      while(count--)
      {
  	if(isnuc)
  	{
  	    xlen = strlen(nuc);
  	    td = ajRandomDouble() * xlen;
  	    r = (ajint) td;
  	    ajStrAppendK(&ins, nuc[r]);
  	}
  	else
  	{
  	    xlen = strlen(prot);
  	    td = ajRandomDouble() * xlen;
  	    r = (ajint) td;
  	    ajStrAppendK(&ins, prot[r]);
  	}
      }
      ajDebug("Inserting %S at %d\n", ins, start);
      ajStrInsertS(str, start, ins);
      ajStrDel(&ins);
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_Move *****************************************************
  **
  ** Move a block of sequence from one position to another
  **
  ** @param [u] str [AjPStr*] sequence to move within
  ** @param [r] start [ajint] start position of block to move
  ** @param [r] end [ajint] end position of block to move
  ** @param [r] destination [ajint] destination of move
  ** @@
  ******************************************************************************/
  
  static void msbar_Move(AjPStr *str, ajint start, ajint end, ajint destination)
  {
      AjPStr mov = NULL;
  
      ajStrAssignS(&mov, *str);
      ajStrKeepRange(&mov, start, end);
      ajStrInsertS(str, destination, mov);
      ajStrDel(&mov);
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_Duplicate ************************************************
  **
  ** Duplicate a block of sequence adjacent to the source block of sequence
  **
  ** @param [u] str [AjPStr*] sequence to duplicate within
  ** @param [r] start [ajint] start position of block to duplicate
  ** @param [r] end [ajint] end position of block to duplicate
  ** @@
  ******************************************************************************/
  
  static void msbar_Duplicate(AjPStr *str, ajint start, ajint end)
  {
      msbar_Move(str, start, end, start);
  
      return;
  }
  
  
  
  
  /* @funcstatic msbar_Unlike ************************************************
  **
  ** Check that the sequence is unlike any other input sequence
  **
  ** @param [r] str [const AjPStr] sequence to check
  ** @param [u] other [AjPSeqall] set of sequences to check against
  ** @return [AjBool] True if unlike any other sequence
  ** @@
  ******************************************************************************/
  
  static AjBool msbar_Unlike(const AjPStr str, AjPSeqall other)
  {
      AjPSeq nextother;
  
      while(ajSeqallNext(other, &nextother))
          if(ajStrMatchCaseS(str, ajSeqGetSeqS(nextother)))
              return ajFalse;
      
      return ajTrue;
  }
  