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