src/jalview/util/HMMProbabilityDistributionAnalyser.java

   1 package jalview.util;
   2
   3 import jalview.datamodel.HiddenMarkovModel;
   4 import jalview.datamodel.SequenceI;
   5 import jalview.io.DataSourceType;
   6 import jalview.io.FileParse;
   7 import jalview.io.HMMFile;
   8 import jalview.io.StockholmFile;
   9 import jalview.schemes.ResidueProperties;
  10
  11 import java.io.BufferedReader;
  12 import java.io.File;
  13 import java.io.FileNotFoundException;
  14 import java.io.FileReader;
  15 import java.io.IOException;
  16 import java.io.PrintWriter;
  17 import java.util.ArrayList;
  18 import java.util.HashMap;
  19 import java.util.List;
  20 import java.util.Map;
  21 import java.util.Scanner;
  22 import java.util.Vector;
  23
  24 /**
  25  * Processes probability data. The file indexes used in this program represent
  26  * the index of the location of a family or hmm in their respective files,
  27  * starting from 0.
  28  *
  29  * @author TZVanaalten
  30  *
  31  */
  32 public class HMMProbabilityDistributionAnalyser
  33 {
  34
  35   Vector<SequenceI> sequences;
  36
  37   HiddenMarkovModel hmm;
  38
  39   // contains the raw data produced
  40   List<ArrayList<Double>> raw = new ArrayList<>();
  41
  42   // contains binned data
  43   Map<String, Double> binned = new HashMap<>();
  44
  45   // location of the family file
  46   final static String FAMILIES = "C:/Users/TZVanaalten/Pfam-A.full";
  47
  48   // location of the HMM file
  49   final static String HMMS = "H:/Desktop/PFAM/HMMs/Pfam-A.hmm";
  50
  51   // suffix for raw file
  52   final static String RAW = "/Raw.csv";
  53
  54   // suffix for binned file
  55   final static String BINNED = "/Binned.csv";
  56
  57   // normalisation scale
  58   final static double SCALE = 100000;
  59
  60   // current position in file
  61   int currentFilePosition = 0;
  62
  63   final static String NL = "\n";
  64
  65   // current directory
  66   String currentFolder;
  67
  68   /**
  69    * Sets the working directory.
  70    *
  71    * @param path
  72    */
  73   public void setFolder(String path)
  74   {
  75     currentFolder = path;
  76   }
  77
  78   /**
  79    * Moves a buffered reader to a specific location in the file, delimited by
  80    * '//'.
  81    *
  82    * @param index
  83    *          The index of the location in the file.
  84    * @param br
  85    * @throws IOException
  86    */
  87   public void moveToFile(int index, BufferedReader br) throws IOException
  88   {
  89     for (int i = 0; i < index; i++)
  90     {
  91       String line = br.readLine();
  92       while (!"//".equals(line))
  93       {
  94         line = br.readLine();
  95       }
  96     }
  97
  98   }
  99
 100   /**
 101    * Analyses a specified number of families and then saves the data. Before
 102    * analysing the data, the previous saved data will be imported and after
 103    * analysing this data is exported back into the file.
 104    *
 105    * @param increments
 106    *          The number of families to read before saving.
 107    * @throws IOException
 108    */
 109   public void run(int increments) throws IOException
 110   {
 111
 112     readPreviousData(currentFolder);
 113
 114     BufferedReader posReader = new BufferedReader(
 115             new FileReader(currentFolder + "/CurrentPosition.txt"));
 116     String line = posReader.readLine();
 117     posReader.close();
 118     currentFilePosition = Integer.parseInt(line);
 119
 120     BufferedReader inputSTO = new BufferedReader(
 121             new FileReader(FAMILIES));
 122     BufferedReader inputHMM = new BufferedReader(
 123             new FileReader(HMMS));
 124
 125     moveToFile(currentFilePosition, inputHMM);
 126     moveToFile(currentFilePosition, inputSTO);
 127
 128     int filesRead = 0;
 129     while (filesRead < increments)
 130     {
 131       FileParse parserSTO = new FileParse(inputSTO, "",
 132               DataSourceType.FILE);
 133       readStockholm(parserSTO);
 134
 135       FileParse parserHMM = new FileParse(inputHMM, "",
 136               DataSourceType.FILE);
 137       readHMM(parserHMM);
 138
 139       if (hmm.getAlphabetType().equals("amino"))
 140       {
 141         int count = countValidResidues();
 142         processData(count);
 143         filesRead++;
 144       }
 145       currentFilePosition++;
 146     }
 147
 148     PrintWriter p = new PrintWriter(
 149             new File(currentFolder + "/CurrentPosition"));
 150     p.print(currentFilePosition);
 151     p.close();
 152     exportData(currentFolder);
 153     raw.clear();
 154     binned.clear();
 155
 156   }
 157
 158   /**
 159    * Reads the previous data from both files
 160    *
 161    * @param source
 162    * @throws IOException
 163    */
 164   public void readPreviousData(String source) throws IOException
 165   {
 166     readBinned(source);
 167     readRaw(source);
 168   }
 169
 170   /**
 171    * Reads the previous data from the binned file.
 172    *
 173    * @param source
 174    * @throws IOException
 175    */
 176   public void readBinned(String source) throws IOException
 177   {
 178     BufferedReader input = new BufferedReader(
 179             new FileReader(source + BINNED));
 180     String line = input.readLine();
 181     while (!("".equals(line) || line == null))
 182     {
 183       binned = new HashMap<>();
 184       Scanner scanner = new Scanner(line);
 185       scanner.useDelimiter(",");
 186       binned.put(scanner.next(), scanner.nextDouble());
 187       scanner.close();
 188       line = input.readLine();
 189     }
 190
 191     input.close();
 192   }
 193
 194   /**
 195    * Reads the previous data from the raw file.
 196    *
 197    * @param source
 198    * @throws IOException
 199    */
 200   public void readRaw(String source) throws IOException
 201   {
 202     BufferedReader input = new BufferedReader(new FileReader(source + RAW));
 203     String line = input.readLine();
 204     if (line == null)
 205     {
 206       input.close();
 207       return;
 208     }
 209     Scanner numberScanner = new Scanner(line);
 210     numberScanner.useDelimiter(",");
 211     raw = new ArrayList<>();
 212     while (numberScanner.hasNext())
 213     {
 214       numberScanner.next();
 215       raw.add(new ArrayList<Double>());
 216     }
 217     numberScanner.close();
 218
 219     line = input.readLine();
 220     while (!("".equals(line) || line == null))
 221     {
 222       Scanner scanner = new Scanner(line);
 223       scanner.useDelimiter(",");
 224
 225       int i = 0;
 226       while (scanner.hasNext())
 227       {
 228         String value;
 229         value = scanner.next();
 230         if (!value.equals("EMPTY"))
 231         {
 232           raw.get(i).add(Double.parseDouble(value));
 233         }
 234
 235         i++;
 236       }
 237       scanner.close();
 238       line = input.readLine();
 239     }
 240
 241     input.close();
 242   }
 243
 244   /**
 245    * Counts the number of valid residues in the sequence.
 246    *
 247    * @return
 248    */
 249   public int countValidResidues()
 250   {
 251     int count = 0;
 252
 253     for (int width = 0; width < sequences.size(); width++)
 254     {
 255       for (int length = 1; length < hmm.getLength(); length++)
 256       {
 257         char symbol;
 258         int alignPos;
 259         alignPos = hmm.getNodeAlignmentColumn(length);
 260
 261         symbol = sequences.get(width).getCharAt(alignPos);
 262         if (ResidueProperties.aminoBackgroundFrequencies
 263                 .containsKey(symbol))
 264         {
 265           count++;
 266         }
 267       }
 268     }
 269
 270     return count;
 271   }
 272
 273   /**
 274    * Processes data, and stores it in both a raw and binned format.
 275    *
 276    * @param count
 277    */
 278   public void processData(int count)
 279   {
 280
 281     raw.add(new ArrayList<Double>());
 282     int rawPos = raw.size() - 1;
 283     for (int width = 0; width < sequences.size(); width++)
 284     {
 285       for (int length = 1; length < hmm.getLength(); length++)
 286       {
 287         char symbol;
 288         int alignPos;
 289         alignPos = hmm.getNodeAlignmentColumn(length);
 290
 291         symbol = sequences.get(width).getCharAt(alignPos);
 292         if (ResidueProperties.aminoBackgroundFrequencies
 293                 .containsKey(symbol))
 294         {
 295
 296           Double prob;
 297           Float bfreq;
 298           Double llr;
 299           prob = hmm.getMatchEmissionProbability(alignPos, symbol);
 300           bfreq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
 301           llr = Math.log(prob / bfreq);
 302           raw.get(rawPos).add(llr);
 303           String output;
 304           output = String.format("%.1f", llr);
 305           if ("-0.0".equals(output))
 306           {
 307             output = "0.0";
 308           }
 309           if (binned.containsKey(output))
 310           {
 311             double prev = binned.get(output);
 312             prev += (SCALE / count);
 313             binned.put(output, prev);
 314
 315           }
 316           else
 317           {
 318             binned.put(output, SCALE / count);
 319           }
 320         }
 321       }
 322     }
 323   }
 324
 325
 326   /**
 327    * Reads in the sequence data from a Stockholm file.
 328    *
 329    * @param source
 330    * @throws IOException
 331    */
 332   public void readStockholm(FileParse source) throws IOException
 333   {
 334     StockholmFile file = new StockholmFile(source);
 335     file.parse();
 336     sequences = file.getSeqs();
 337   }
 338
 339   /**
 340    * Reads in the HMM data from a HMMer file.
 341    *
 342    * @param source
 343    * @throws IOException
 344    */
 345   public void readHMM(FileParse source) throws IOException
 346   {
 347
 348     HMMFile file = new HMMFile(source);
 349     file.parse();
 350     hmm = file.getHMM();
 351
 352   }
 353
 354   /**
 355    * Exports both the binned and raw data into separate files.
 356    *
 357    * @param location
 358    * @throws FileNotFoundException
 359    */
 360   public void exportData(String location) throws FileNotFoundException
 361   {
 362     PrintWriter writerBin = new PrintWriter(new File(location + BINNED));
 363     for (Map.Entry<String, Double> entry : binned.entrySet())
 364     {
 365       writerBin.println(entry.getKey() + "," + entry.getValue());
 366     }
 367     writerBin.close();
 368
 369     PrintWriter writerRaw = new PrintWriter(new File(location + RAW));
 370
 371     StringBuilder identifier = new StringBuilder();
 372
 373     for (int i = 1; i < raw.size() + 1; i++)
 374     {
 375       identifier.append("Fam " + i + ",");
 376     }
 377
 378     writerRaw.println(identifier);
 379
 380     boolean rowIsEmpty = false;
 381     int row = 0;
 382     while (!rowIsEmpty)
 383     {
 384       rowIsEmpty = true;
 385       StringBuilder string = new StringBuilder();
 386       for (int column = 0; column < raw.size(); column++)
 387       {
 388         if (raw.get(column).size() <= row)
 389         {
 390           string.append("EMPTY,");
 391         }
 392         else
 393         {
 394           string.append(raw.get(column).get(row) + ",");
 395           rowIsEmpty = false;
 396         }
 397       }
 398       row++;
 399       writerRaw.println(string);
 400     }
 401     writerRaw.close();
 402
 403   }
 404
 405   /**
 406    * Prints the specified family on the console.
 407    *
 408    * @param index
 409    * @throws IOException
 410    */
 411   public void printFam(int index) throws IOException
 412   {
 413     BufferedReader br = new BufferedReader(new FileReader(FAMILIES));
 414
 415     moveToFile(index, br);
 416
 417     String line = br.readLine();
 418
 419     while (!"//".equals(line))
 420     {
 421       System.out.println(line);
 422       line = br.readLine();
 423     }
 424     System.out.println(line);
 425     br.close();
 426
 427   }
 428
 429   /**
 430    * Prints the specified HMM on the console.
 431    *
 432    * @param index
 433    * @throws IOException
 434    */
 435   public void printHMM(int index) throws IOException
 436   {
 437     BufferedReader br = new BufferedReader(new FileReader(HMMS));
 438
 439     moveToFile(index, br);
 440
 441     String line = br.readLine();
 442
 443     while (!"//".equals(line))
 444     {
 445       System.out.println(line);
 446       line = br.readLine();
 447     }
 448     System.out.println(line);
 449     br.close();
 450
 451   }
 452
 453   /**
 454    * Prints the specified family to a .sto file in the current directory.
 455    *
 456    * @param index
 457    * @throws IOException
 458    */
 459   public void printFamToFile(int index) throws IOException
 460   {
 461     String name;
 462
 463     BufferedReader nameFinder = new BufferedReader(
 464             new FileReader(FAMILIES));
 465
 466     moveToFile(index, nameFinder);
 467
 468     nameFinder.readLine();
 469
 470     Scanner scanner = new Scanner(nameFinder.readLine());
 471     scanner.next();
 472     scanner.next();
 473     name = scanner.next();
 474     scanner.close();
 475
 476     BufferedReader br = new BufferedReader(new FileReader(FAMILIES));
 477
 478     moveToFile(index, br);
 479
 480     String line = br.readLine();
 481     PrintWriter writer = new PrintWriter(
 482             currentFolder + "/" + name + ".sto");
 483     while (!"//".equals(line))
 484     {
 485       writer.println(line);
 486       line = br.readLine();
 487     }
 488     writer.println(line);
 489     writer.close();
 490     br.close();
 491
 492   }
 493
 494   /**
 495    * Prints the specified family to a .hmm file in the current directory.
 496    *
 497    * @param index
 498    * @throws IOException
 499    */
 500   public void printHMMToFile(int index) throws IOException
 501   {
 502
 503     String name;
 504
 505     BufferedReader nameFinder = new BufferedReader(new FileReader(HMMS));
 506
 507     moveToFile(index, nameFinder);
 508
 509     nameFinder.readLine();
 510
 511     Scanner scanner = new Scanner(nameFinder.readLine());
 512     name = scanner.next();
 513     name = scanner.next();
 514     scanner.close();
 515
 516     BufferedReader br = new BufferedReader(new FileReader(HMMS));
 517
 518     moveToFile(index, br);
 519
 520     String line = br.readLine();
 521
 522     PrintWriter writer = new PrintWriter(
 523             currentFolder + "/" + name + ".hmm");
 524     while (!"//".equals(line))
 525     {
 526       writer.println(line);
 527       line = br.readLine();
 528     }
 529     writer.println(line);
 530     writer.close();
 531     br.close();
 532
 533   }
 534
 535   /**
 536    * Clears all raw, binned and current position data in the current directory.
 537    *
 538    * @throws FileNotFoundException
 539    */
 540   public void clear() throws FileNotFoundException
 541   {
 542     PrintWriter pos = new PrintWriter(
 543             currentFolder + "/CurrentPosition.txt");
 544     pos.println("0");
 545
 546     PrintWriter raw = new PrintWriter(currentFolder + RAW);
 547
 548     PrintWriter bin = new PrintWriter(currentFolder + BINNED);
 549
 550     pos.close();
 551     bin.close();
 552     raw.close();
 553   }
 554
 555 }