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.io.File;
  31 import java.io.FileNotFoundException;
  32 import java.io.PrintStream;
  33 import java.util.Arrays;
  34 import java.util.Random;
  35
  36 import org.testng.annotations.BeforeClass;
  37
  38 /**
  39  * Generates, and outputs in Fasta format, a random peptide or nucleotide alignment for given
  40  * sequence length and count. Will regenerate the same alignment each time if
  41  * the same random seed is used (so may be used for reproducible unit tests).
  42  * Not guaranteed to reproduce the same results between versions, as the rules
  43  * may get tweaked to produce more 'realistic' results.
  44  *
  45  * @author gmcarstairs
  46  */
  47 public class AlignmentGenerator
  48 {
  49   private static final char GAP = '-';
  50
  51   private static final char ZERO = '0';
  52
  53   private static final char[] NUCS = "GTCA".toCharArray();
  54
  55   private static final char[] PEPS = "MILVFYWHKRDEQNTCGASNP".toCharArray();
  56
  57   private static char[] BASES;
  58
  59   private Random random;
  60
  61   private PrintStream ps;
  62
  63   /**
  64    * Outputs a pseudo-randomly generated nucleotide or peptide alignment
  65    * Arguments:
  66    * <ul>
  67    * <li>n (for nucleotide) or p (for peptide)</li>
  68    * <li>length (number of bases in each sequence)</li>
  69    * <li>height (number of sequences)</li>
  70    * <li>a whole number random seed</li>
  71    * <li>percentage of gaps to include (0-100)</li>
  72    * <li>percentage chance of variation of each position (0-100)</li>
  73    * <li>(optional) path to a file to write the alignment to</li>
  74    * </ul>
  75    *
  76    *
  77    * @param args
  78    * @throws FileNotFoundException
  79    */
  80   public static void main(String[] args) throws FileNotFoundException
  81   {
  82     if (args.length != 6 && args.length != 7)
  83     {
  84       usage();
  85       return;
  86     }
  87
  88     PrintStream ps = System.out;
  89     if (args.length == 7)
  90     {
  91       ps = new PrintStream(new File(args[6]));
  92     }
  93
  94     boolean nucleotide = args[0].toLowerCase().startsWith("n");
  95     int width = Integer.parseInt(args[1]);
  96     int height = Integer.parseInt(args[2]);
  97     long randomSeed = Long.valueOf(args[3]);
  98     int gapPercentage = Integer.valueOf(args[4]);
  99     int changePercentage = Integer.valueOf(args[5]);
 100
 101     ps.println("; " + height + " sequences of " + width
 102             + " bases with " + gapPercentage + "% gaps and "
 103             + changePercentage + "% mutations (random seed = " + randomSeed
 104             + ")");
 105
 106     new AlignmentGenerator(nucleotide, ps).generate(width, height,
 107             randomSeed, gapPercentage, changePercentage);
 108
 109     if (ps != System.out)
 110     {
 111       ps.close();
 112     }
 113   }
 114
 115   /**
 116    * Prints parameter help
 117    */
 118   private static void usage()
 119   {
 120     System.out.println("Usage:");
 121     System.out.println("arg0: n (for nucleotide) or p (for peptide)");
 122     System.out.println("arg1: number of (non-gap) bases per sequence");
 123     System.out.println("arg2: number of sequences");
 124     System.out
 125             .println("arg3: an integer as random seed (same seed = same results)");
 126     System.out.println("arg4: percentage of gaps to (randomly) generate");
 127     System.out
 128             .println("arg5: percentage of 'mutations' to (randomly) generate");
 129     System.out
 130             .println("arg6: (optional) path to output file (default is sysout)");
 131     System.out.println("Example: AlignmentGenerator n 12 15 387 10 5");
 132     System.out
 133             .println("- 15 nucleotide sequences of 12 bases each, approx 10% gaps and 5% mutations, random seed = 387");
 134
 135   }
 136
 137   /**
 138    * Constructor that sets nucleotide or peptide symbol set, and also writes the
 139    * generated alignment to sysout
 140    */
 141   public AlignmentGenerator(boolean nuc)
 142   {
 143     this(nuc, System.out);
 144   }
 145
 146   /**
 147    * Constructor that sets nucleotide or peptide symbol set, and also writes the
 148    * generated alignment to the specified output stream (if not null). This can
 149    * be used to write the alignment to a file or sysout.
 150    */
 151   public AlignmentGenerator(boolean nucleotide, PrintStream printStream)
 152   {
 153     BASES = nucleotide ? NUCS : PEPS;
 154     ps = printStream;
 155   }
 156
 157   /**
 158    * Outputs an 'alignment' of given width and height, where each position is a
 159    * random choice from the symbol alphabet, or - for gap
 160    *
 161    * @param width
 162    * @param height
 163    * @param randomSeed
 164    * @param changePercentage
 165    * @param gapPercentage
 166    */
 167   public AlignmentI generate(int width, int height, long randomSeed,
 168           int gapPercentage, int changePercentage)
 169   {
 170     SequenceI[] seqs = new SequenceI[height];
 171     random = new Random(randomSeed);
 172     seqs[0] = generateSequence(1, width, gapPercentage);
 173     for (int seqno = 1; seqno < height; seqno++)
 174     {
 175       seqs[seqno] = generateAnotherSequence(seqs[0].getSequence(),
 176               seqno + 1, width, changePercentage);
 177     }
 178     AlignmentI al = new Alignment(seqs);
 179
 180     if (ps != null)
 181     {
 182       ps.println(new FastaFile().print(al.getSequencesArray(), true));
 183     }
 184
 185     return al;
 186   }
 187
 188   /**
 189    * Outputs a DNA 'sequence' of given length, with some random gaps included.
 190    *
 191    * @param seqno
 192    * @param length
 193    * @param gapPercentage
 194    */
 195   private SequenceI generateSequence(int seqno, int length,
 196           int gapPercentage)
 197   {
 198     StringBuilder seq = new StringBuilder(length);
 199
 200     /*
 201      * Loop till we've added 'length' bases (excluding gaps)
 202      */
 203     for (int count = 0; count < length;)
 204     {
 205       boolean addGap = random.nextInt(100) < gapPercentage;
 206       char c = addGap ? GAP : BASES[random.nextInt(Integer.MAX_VALUE)
 207               % BASES.length];
 208       seq.append(c);
 209       if (!addGap)
 210       {
 211         count++;
 212       }
 213     }
 214     final String seqName = "SEQ" + seqno;
 215     final String seqString = seq.toString();
 216     SequenceI sq = new Sequence(seqName, seqString);
 217     sq.createDatasetSequence();
 218     return sq;
 219   }
 220
 221   /**
 222    * Generate a sequence approximately aligned to the first one.
 223    *
 224    * @param ds
 225    * @param seqno
 226    * @param width
 227    *          number of bases
 228    * @param changePercentage
 229    * @return
 230    */
 231   private SequenceI generateAnotherSequence(char[] ds, int seqno,
 232           int width, int changePercentage)
 233   {
 234     int length = ds.length;
 235     char[] seq = new char[length];
 236     Arrays.fill(seq, ZERO);
 237     int gapsWanted = length - width;
 238     int gapsAdded = 0;
 239
 240     /*
 241      * First 'randomly' mimic gaps in model sequence.
 242      */
 243     for (int pos = 0; pos < length; pos++)
 244     {
 245       if (ds[pos] == GAP)
 246       {
 247         /*
 248          * Add a gap at the same position with changePercentage likelihood
 249          */
 250         seq[pos] = randomCharacter(GAP, changePercentage);
 251         if (seq[pos] == GAP)
 252         {
 253           gapsAdded++;
 254         }
 255       }
 256     }
 257
 258     /*
 259      * Next scatter any remaining gaps (if any) at random. This gives an even
 260      * distribution.
 261      */
 262     while (gapsAdded < gapsWanted)
 263     {
 264       boolean added = false;
 265       while (!added)
 266       {
 267         int pos = random.nextInt(length);
 268         if (seq[pos] != GAP)
 269         {
 270           seq[pos] = GAP;
 271           added = true;
 272           gapsAdded++;
 273         }
 274       }
 275     }
 276
 277     /*
 278      * Finally fill in the rest with randomly mutated bases.
 279      */
 280     for (int pos = 0; pos < length; pos++)
 281     {
 282       if (seq[pos] == ZERO)
 283       {
 284         char c = randomCharacter(ds[pos], changePercentage);
 285         seq[pos] = c;
 286       }
 287     }
 288     final String seqName = "SEQ" + seqno;
 289     final String seqString = new String(seq);
 290     SequenceI sq = new Sequence(seqName, seqString);
 291     sq.createDatasetSequence();
 292     return sq;
 293   }
 294
 295   /**
 296    * Returns a random character that is changePercentage% likely to match the
 297    * given type (as base or gap).
 298    *
 299    * @param changePercentage
 300    *
 301    * @param c
 302    * @return
 303    */
 304   private char randomCharacter(char c, int changePercentage)
 305   {
 306     final boolean mutation = random.nextInt(100) < changePercentage;
 307
 308     if (!mutation)
 309     {
 310       return c;
 311     }
 312
 313     char newchar = c;
 314     while (newchar == c)
 315     {
 316       newchar = BASES[random.nextInt(Integer.MAX_VALUE) % BASES.length];
 317     }
 318     return newchar;
 319   }
 320 }