test/jalview/analysis/DnaAlignmentGenerator.java

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