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