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