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