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