test/jalview/analysis/DnaAlignmentGenerator.java

   1 package jalview.analysis;
   2
   3 import jalview.datamodel.Alignment;
   4 import jalview.datamodel.AlignmentI;
   5 import jalview.datamodel.Sequence;
   6 import jalview.datamodel.SequenceI;
   7 import jalview.io.FastaFile;
   8
   9 import java.util.Arrays;
  10 import java.util.Random;
  11
  12 /**
  13  * Generates, and outputs in Fasta format, a random DNA alignment for given
  14  * sequence length and count. Will regenerate the same alignment each time if
  15  * the same random seed is used (so may be used for reproducible unit tests).
  16  * Not guaranteed to reproduce the same results between versions, as the rules
  17  * may get tweaked to produce more 'realistic' results.
  18  *
  19  * Arguments:
  20  * <ul>
  21  * <li>length (number of bases in each sequence)</li>
  22  * <li>height (number of sequences)</li>
  23  * <li>a whole number random seed</li>
  24  * <li>percentage of gaps to include (0-100)</li>
  25  * <li>percentage chance of variation of each position (0-100)</li>
  26  * </ul>
  27  *
  28  * @author gmcarstairs
  29  *
  30  */
  31 public class DnaAlignmentGenerator
  32 {
  33   private static final char GAP = '-';
  34
  35   private static final char ZERO = '0';
  36
  37   private static final char[] BASES = new char[] { 'G', 'T', 'C', 'A' };
  38
  39   private Random random;
  40
  41   /**
  42    * Outputs a DNA 'alignment' where each position is a random choice from
  43    * 'GTCA-'.
  44    *
  45    * @param args
  46    */
  47   public static void main(String[] args)
  48   {
  49     if (args.length != 5)
  50     {
  51       usage();
  52       return;
  53     }
  54     int width = Integer.parseInt(args[0]);
  55     int height = Integer.parseInt(args[1]);
  56     long randomSeed = Long.valueOf(args[2]);
  57     int gapPercentage = Integer.valueOf(args[3]);
  58     int changePercentage = Integer.valueOf(args[4]);
  59     AlignmentI al = new DnaAlignmentGenerator().generate(width, height,
  60             randomSeed, gapPercentage, changePercentage);
  61
  62     System.out.println("; " + height + " sequences of " + width
  63             + " bases with " + gapPercentage + "% gaps and "
  64             + changePercentage + "% mutations (random seed = " + randomSeed
  65             + ")");
  66     System.out.println(new FastaFile().print(al.getSequencesArray()));
  67   }
  68
  69   /**
  70    * Print parameter help.
  71    */
  72   private static void usage()
  73   {
  74     System.out.println("Usage:");
  75     System.out.println("arg0: number of (non-gap) bases per sequence");
  76     System.out.println("arg1: number sequences");
  77     System.out
  78             .println("arg2: an integer as random seed (same seed = same results)");
  79     System.out.println("arg3: percentage of gaps to (randomly) generate");
  80     System.out
  81             .println("arg4: percentage of 'mutations' to (randomly) generate");
  82     System.out.println("Example: DnaAlignmentGenerator 12 15 387 10 5");
  83     System.out
  84             .println("- 15 sequences of 12 bases each, approx 10% gaps and 5% mutations, random seed = 387");
  85
  86   }
  87
  88   /**
  89    * Default constructor
  90    */
  91   public DnaAlignmentGenerator()
  92   {
  93
  94   }
  95
  96   /**
  97    * Outputs a DNA 'alignment' of given width and height, where each position is
  98    * a random choice from 'GTCA-'.
  99    *
 100    * @param width
 101    * @param height
 102    * @param randomSeed
 103    * @param changePercentage
 104    * @param gapPercentage
 105    */
 106   public AlignmentI generate(int width, int height, long randomSeed,
 107           int gapPercentage, int changePercentage)
 108   {
 109     SequenceI[] seqs = new SequenceI[height];
 110     random = new Random(randomSeed);
 111     seqs[0] = generateSequence(1, width, gapPercentage);
 112     for (int seqno = 1; seqno < height; seqno++)
 113     {
 114       seqs[seqno] = generateAnotherSequence(seqs[0].getSequence(),
 115               seqno + 1, width, changePercentage);
 116     }
 117     AlignmentI al = new Alignment(seqs);
 118     return al;
 119   }
 120
 121   /**
 122    * Outputs a DNA 'sequence' of given length, with some random gaps included.
 123    *
 124    * @param seqno
 125    * @param length
 126    * @param gapPercentage
 127    */
 128   private SequenceI generateSequence(int seqno, int length,
 129           int gapPercentage)
 130   {
 131     StringBuilder seq = new StringBuilder(length);
 132
 133     /*
 134      * Loop till we've added 'length' bases (excluding gaps)
 135      */
 136     for (int count = 0; count < length;)
 137     {
 138       boolean addGap = random.nextInt(100) < gapPercentage;
 139       char c = addGap ? GAP : BASES[random.nextInt(Integer.MAX_VALUE) % 4];
 140       seq.append(c);
 141       if (!addGap)
 142       {
 143         count++;
 144       }
 145     }
 146     final String seqName = "SEQ" + seqno;
 147     final String seqString = seq.toString();
 148     SequenceI sq = new Sequence(seqName, seqString);
 149     sq.createDatasetSequence();
 150     return sq;
 151   }
 152
 153   /**
 154    * Generate a sequence approximately aligned to the first one.
 155    *
 156    * @param ds
 157    * @param seqno
 158    * @param width
 159    *          number of bases
 160    * @param changePercentage
 161    * @return
 162    */
 163   private SequenceI generateAnotherSequence(char[] ds, int seqno,
 164           int width, int changePercentage)
 165   {
 166     int length = ds.length;
 167     char[] seq = new char[length];
 168     Arrays.fill(seq, ZERO);
 169     int gapsWanted = length - width;
 170     int gapsAdded = 0;
 171
 172     /*
 173      * First 'randomly' mimic gaps in model sequence.
 174      */
 175     for (int pos = 0; pos < length; pos++)
 176     {
 177       if (ds[pos] == GAP)
 178       {
 179         /*
 180          * Add a gap at the same position with changePercentage likelihood
 181          */
 182         seq[pos] = randomCharacter(GAP, changePercentage);
 183         if (seq[pos] == GAP)
 184         {
 185           gapsAdded++;
 186         }
 187       }
 188     }
 189
 190     /*
 191      * Next scatter any remaining gaps (if any) at random. This gives an even
 192      * distribution.
 193      */
 194     while (gapsAdded < gapsWanted)
 195     {
 196       boolean added = false;
 197       while (!added)
 198       {
 199         int pos = random.nextInt(length);
 200         if (seq[pos] != GAP)
 201         {
 202           seq[pos] = GAP;
 203           added = true;
 204           gapsAdded++;
 205         }
 206       }
 207     }
 208
 209     /*
 210      * Finally fill in the rest with randomly mutated bases.
 211      */
 212     for (int pos = 0; pos < length; pos++)
 213     {
 214       if (seq[pos] == ZERO)
 215       {
 216         char c = randomCharacter(ds[pos], changePercentage);
 217         seq[pos] = c;
 218       }
 219     }
 220     final String seqName = "SEQ" + seqno;
 221     final String seqString = new String(seq);
 222     SequenceI sq = new Sequence(seqName, seqString);
 223     sq.createDatasetSequence();
 224     return sq;
 225   }
 226
 227   /**
 228    * Returns a random character that is changePercentage% likely to match the
 229    * given type (as base or gap).
 230    *
 231    * @param changePercentage
 232    *
 233    * @param c
 234    * @return
 235    */
 236   private char randomCharacter(char c, int changePercentage)
 237   {
 238     final boolean mutation = random.nextInt(100) < changePercentage;
 239
 240     if (!mutation)
 241     {
 242       return c;
 243     }
 244
 245     char newchar = c;
 246     while (newchar == c)
 247     {
 248       newchar = BASES[random.nextInt(Integer.MAX_VALUE) % 4];
 249     }
 250     return newchar;
 251   }
 252 }