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