src/jalview/io/HMMFile.java

   1 package jalview.io;
   2
   3 import jalview.datamodel.HMMNode;
   4 import jalview.datamodel.HiddenMarkovModel;
   5 import jalview.datamodel.SequenceI;
   6
   7 import java.io.BufferedReader;
   8 import java.io.FileNotFoundException;
   9 import java.io.IOException;
  10 import java.io.PrintWriter;
  11 import java.io.UnsupportedEncodingException;
  12 import java.util.ArrayList;
  13 import java.util.List;
  14 import java.util.Scanner;
  15
  16
  17 /**
  18  * Adds capability to read in and write out HMMER3 files. Currently only supports HMMER3/f.
  19  *
  20  *
  21  * @author TZVanaalten
  22  *
  23  */
  24 public class HMMFile extends AlignFile
  25         implements AlignmentFileReaderI, AlignmentFileWriterI
  26 {
  27   // HMM to store file data
  28   private HiddenMarkovModel hmm = new HiddenMarkovModel();
  29
  30   // number of possible transitions
  31   private final int NUMBER_OF_TRANSITIONS = 7;
  32
  33   private final String NEW_LINE = "\n";
  34
  35   //number of symbols in the alphabet used in the hidden Markov model
  36   int numberOfSymbols;
  37
  38   private final String SPACE = " ";
  39
  40   private final String COMPO = "COMPO";
  41
  42   private final String EMPTY = "";
  43
  44   //This is a line that needs to be added to each HMMER£ file. It is purely for readability.
  45   private static final String TRANSITIONTYPELINE = "m->m     m->i     m->d     i->m     i->i     d->m     d->d";
  46
  47   /**
  48    * Constructor for HMMFile, parses immediately
  49    *
  50    * @param source
  51    * @throws IOException
  52    */
  53   public HMMFile(FileParse source) throws IOException
  54   {
  55     super(false, source);
  56     parse();
  57   }
  58
  59   /**
  60    * Default constructor, do not use!
  61    */
  62   public HMMFile()
  63   {
  64
  65   }
  66
  67   /**
  68    * Constructor for HMMFile used for exporting.
  69    *
  70    * @param hmm
  71    * @param exportImmediately
  72    */
  73   public HMMFile(HiddenMarkovModel markov)
  74   {
  75     hmm = markov;
  76   }
  77
  78   /**
  79    * Returns the HMM produced by reading in a HMMER3 file.
  80    *
  81    * @return
  82    */
  83   public HiddenMarkovModel getHMM()
  84   {
  85     return hmm;
  86   }
  87
  88   /**
  89    * Sets the HMM used in this file.
  90    *
  91    * @param model
  92    */
  93   public void setHMM(HiddenMarkovModel model)
  94   {
  95     this.hmm = model;
  96   }
  97
  98   /**
  99    * Gets the name of the hidden Markov model.
 100    *
 101    * @return
 102    */
 103   public String getName()
 104   {
 105     return hmm.getName();
 106   }
 107
 108   /**
 109    * Reads the data from HMM file into the HMM field on this object.
 110    *
 111    * @throws IOException
 112    */
 113   @Override
 114   public void parse() throws IOException
 115   {
 116     parseFileProperties(dataIn);
 117     parseModel(dataIn);
 118   }
 119
 120   /**
 121    * Reads the data from HMM file into the HMM field on this object.
 122    *
 123    * @throws IOException
 124    */
 125
 126   public void parse(BufferedReader br) throws IOException
 127   {
 128     parseFileProperties(br);
 129     parseModel(br);
 130   }
 131
 132
 133
 134   /**
 135    * Imports the file properties from a HMMER3 file.
 136    *
 137    * @param input
 138    *          The buffered reader used to read in the file.
 139    * @throws IOException
 140    */
 141   void parseFileProperties(BufferedReader input) throws IOException
 142   {
 143     boolean readingFile = true;
 144     hmm.setFileHeader(input.readLine());
 145     String line = input.readLine();
 146     while (readingFile)
 147     {
 148       if (line != null)
 149       {
 150         Scanner parser = new Scanner(line);
 151         String next = parser.next();
 152         if ("HMM".equals(next)) // indicates start of HMM data (end of file
 153                               // properties)
 154         {
 155           readingFile = false;
 156           hmm.fillSymbols(parser);
 157           numberOfSymbols = hmm.getNumberOfSymbols();
 158         }
 159         else if ("STATS".equals(next))
 160         {
 161           parser.next();
 162           String key;
 163           String value;
 164           key = parser.next();
 165           value = parser.next() + SPACE + SPACE + parser.next();
 166           hmm.addFileProperty(key, value);
 167         }
 168         else
 169         {
 170           String key = next;
 171           String value = parser.next();
 172           while (parser.hasNext())
 173           {
 174             value = value + SPACE + parser.next();
 175           }
 176           hmm.addFileProperty(key, value);
 177         }
 178         parser.close();
 179       }
 180       line = input.readLine();
 181       if (line == null)
 182       {
 183         readingFile = false;
 184       }
 185     }
 186
 187   }
 188
 189   /**
 190    * Parses the model data from the HMMER3 file
 191    *
 192    * @param input
 193    *          The buffered reader used to read the file.
 194    * @throws IOException
 195    */
 196   void parseModel(BufferedReader input) throws IOException
 197   {
 198     String line = input.readLine();
 199     int node = 0;
 200     while (!"//".equals(line))
 201     {
 202       hmm.getNodes().add(new HMMNode());
 203       String next;
 204       Scanner matchReader = new Scanner(line);
 205       next = matchReader.next();
 206       if (next.equals(COMPO) || node > 0)
 207       {
 208         // stores match emission line in list
 209         List<Double> matches = new ArrayList<>();
 210         matches = fillList(matchReader, numberOfSymbols);
 211         hmm.getNodes().get(node).setMatchEmissions(matches);
 212         if (node > 0)
 213         {
 214           parseAnnotations(matchReader, node);
 215         }
 216       }
 217       matchReader.close();
 218       // stores insert emission line in list
 219       line = input.readLine();
 220       Scanner insertReader = new Scanner(line);
 221       List<Double> inserts = new ArrayList<>();
 222       inserts = fillList(insertReader, numberOfSymbols);
 223       hmm.getNodes().get(node).setInsertEmissions(inserts);
 224       insertReader.close();
 225
 226       // stores state transition line in list
 227       line = input.readLine();
 228       Scanner transitionReader = new Scanner(line);
 229       List<Double> transitions = new ArrayList<>();
 230       transitions = fillList(transitionReader, NUMBER_OF_TRANSITIONS);
 231       hmm.getNodes().get(node).setStateTransitions(transitions);
 232       transitionReader.close();
 233       line = input.readLine();
 234       node++;
 235     }
 236
 237   }
 238
 239   /**
 240    * Parses the annotations on the match emission line.
 241    *
 242    * @param scanner
 243    *          The scanner which is processing match emission line.
 244    * @param index
 245    *          The index of node which is being scanned.
 246    */
 247   void parseAnnotations(Scanner scanner, int index)
 248   {
 249     if (hmm.mapIsActive())
 250     {
 251       int column;
 252       column = scanner.nextInt();
 253       hmm.getNodes().get(index).setAlignmentColumn(column);
 254       hmm.getNodeLookup().put(column, index);
 255     }
 256     else
 257     {
 258       scanner.next();
 259     }
 260
 261     if (scanner.hasNext())
 262     {
 263       char consensusR;
 264       consensusR = charValue(scanner.next());
 265       hmm.getNodes().get(index).setConsensusResidue(consensusR);
 266     }
 267
 268     if (scanner.hasNext())
 269     {
 270       char reference;
 271       reference = charValue(scanner.next());
 272       hmm.getNodes().get(index).setReferenceAnnotation(reference);
 273     }
 274
 275     if (scanner.hasNext())
 276     {
 277       char value;
 278       value = charValue(scanner.next());
 279       hmm.getNodes().get(index).setMaskValue(value);
 280     }
 281     if (scanner.hasNext())
 282     {
 283       char consensusS;
 284       consensusS = charValue(scanner.next());
 285       hmm.getNodes().get(index).setConsensusStructure(consensusS);
 286     }
 287   }
 288
 289
 290
 291   /**
 292    * Fills a list of doubles based on an input line.
 293    *
 294    * @param input
 295    *          The scanner for the line containing the data to be transferred to
 296    *          the list.
 297    * @param numberOfElements
 298    *          The number of elements in the list to be filled.
 299    * @return filled list Returns the list of doubles.
 300    */
 301   static List<Double> fillList(Scanner input,
 302           int numberOfElements)
 303   {
 304     List<Double> list = new ArrayList<>();
 305     for (int i = 0; i < numberOfElements; i++)
 306     {
 307
 308       String next = input.next();
 309       if (next.contains("*")) // state transitions to or from delete states
 310                               // occasionally have values of -infinity. These
 311                               // values are represented by an * in the .hmm
 312                               // file.
 313       {
 314         list.add(Double.NEGATIVE_INFINITY);
 315       }
 316       else
 317       {
 318         double prob = Double.valueOf(next);
 319         prob = Math.pow(Math.E, -prob);
 320         list.add(prob);
 321       }
 322     }
 323     return list;
 324   }
 325
 326
 327   /**
 328    * Writes a HMM to a file/
 329    *
 330    * @param exportLocation
 331    *          Filename, URL or Pasted String to write to.
 332    * @throws FileNotFoundException
 333    * @throws UnsupportedEncodingException
 334    *
 335    **/
 336
 337   public void exportFile(String exportLocation) throws IOException
 338   {
 339     PrintWriter writer = new PrintWriter(exportLocation);
 340     appendFileProperties(writer);
 341     appendModel(writer);
 342     writer.println("//");
 343
 344     writer.close();
 345
 346   }
 347
 348   /**
 349    * Returns a string to be added to the StringBuilder containing the entire
 350    * output String.
 351    *
 352    * @param initialColumnSeparation
 353    *          The initial whitespace separation between the left side of the
 354    *          file and first character.
 355    * @param columnSeparation
 356    *          The separation between subsequent data entries.
 357    * @param data
 358    *          The list fo data to be added to the String.
 359    * @return
 360    */
 361   String addData(int initialColumnSeparation,
 362           int columnSeparation, List<String> data)
 363   {
 364     String line = EMPTY;
 365     int index = 0;
 366     for (String value : data)
 367     {
 368       if (index == 0)
 369       {
 370         line += String.format("%" + initialColumnSeparation + "s", value);
 371       }
 372       else
 373       {
 374         line += String.format("%" + columnSeparation + "s", value);
 375       }
 376       index++;
 377     }
 378     return line;
 379   }
 380
 381   /**
 382    * Converts list of characters into a list of Strings.
 383    *
 384    * @param list
 385    * @return Returns the list of Strings.
 386    */
 387   List<String> charListToStringList(List<Character> list)
 388   {
 389     List<String> strList = new ArrayList<>();
 390     for (char value : list)
 391     {
 392       String strValue = Character.toString(value);
 393       strList.add(strValue);
 394     }
 395     return strList;
 396   }
 397
 398   /**
 399    * Converts a list of doubles into a list of Strings, rounded to the nearest
 400    * 5th decimal place.
 401    *
 402    * @param list
 403    * @param noOfDecimals
 404    * @return
 405    */
 406   List<String> doubleListToStringList(List<Double> list)
 407   {
 408     List<String> strList = new ArrayList<>();
 409     for (double value : list)
 410     {
 411       String strValue;
 412       if (value > 0)
 413       {
 414         strValue = String.format("%.5f", value);
 415
 416       }
 417       else if (value == -0.00000d)
 418       {
 419         strValue = "0.00000";
 420       }
 421       else
 422       {
 423         strValue = "*";
 424       }
 425
 426       strList.add(strValue);
 427     }
 428     return strList;
 429   }
 430
 431   /**
 432    * Converts a primitive array of Strings to a list of Strings.
 433    *
 434    * @param array
 435    * @return
 436    */
 437   List<String> stringArrayToStringList(String[] array)
 438   {
 439     List<String> list = new ArrayList<>();
 440     for (String value : array)
 441     {
 442       list.add(value);
 443     }
 444
 445     return list;
 446   }
 447
 448   /**
 449    * Appends the hidden Markov model data to the StringBuilder containing the
 450    * output
 451    *
 452    * @param file
 453    *          The StringBuilder containing the output.
 454    */
 455   void appendModel(PrintWriter writer)
 456   {
 457     String symbolLine = "HMM";
 458     List<Character> charSymbols = hmm.getSymbols();
 459     List<String> strSymbols;
 460     strSymbols = charListToStringList(charSymbols);
 461     symbolLine += addData(11, 9, strSymbols);
 462     writer.println(symbolLine);
 463     writer.println(TRANSITIONTYPELINE);
 464
 465     int length = hmm.getLength();
 466
 467     for (int node = 0; node <= length; node++)
 468     {
 469       String matchLine;
 470       if (node == 0)
 471       {
 472         matchLine = String.format("%7s", "COMPO");
 473       }
 474       else
 475       {
 476         matchLine = String.format("%7s", node);
 477       }
 478
 479       List<String> strMatches;
 480       List<Double> doubleMatches;
 481       doubleMatches = convertListToLogSpace(
 482               hmm.getNode(node).getMatchEmissions());
 483       strMatches = doubleListToStringList(doubleMatches);
 484       matchLine += addData(10, 9, strMatches);
 485
 486
 487       if (node != 0)
 488       {
 489         matchLine += SPACE + (hmm.getNodeAlignmentColumn(node) + 1);
 490         matchLine += SPACE + hmm.getConsensusResidue(node);
 491         matchLine += SPACE + hmm.getReferenceAnnotation(node);
 492         if (hmm.getFileHeader().contains("HMMER3/f"))
 493         {
 494           matchLine += SPACE + hmm.getMaskedValue(node);
 495           matchLine += SPACE + hmm.getConsensusStructure(node);
 496         }
 497
 498       }
 499
 500       writer.println(matchLine);
 501
 502       String insertLine = EMPTY;
 503       List<String> strInserts;
 504       List<Double> doubleInserts;
 505       doubleInserts = convertListToLogSpace(
 506               hmm.getNode(node).getInsertEmissions());
 507       strInserts = doubleListToStringList(doubleInserts);
 508       insertLine += addData(17, 9, strInserts);
 509
 510       writer.println(insertLine);
 511
 512       String transitionLine = EMPTY;
 513       List<String> strTransitions;
 514       List<Double> doubleTransitions;
 515       doubleTransitions = convertListToLogSpace(
 516               hmm.getNode(node).getStateTransitions());
 517       strTransitions = doubleListToStringList(doubleTransitions);
 518       transitionLine += addData(17, 9, strTransitions);
 519
 520       writer.println(transitionLine);
 521     }
 522   }
 523
 524   /**
 525    * Appends the hidden Markov model file properties to the StringBuilder
 526    * containing the output
 527    *
 528    * @param file
 529    *          The StringBuilder containing the output.
 530    */
 531   void appendFileProperties(PrintWriter writer)
 532   {
 533     String line;
 534
 535     writer.println(hmm.getFileHeader());
 536
 537     line = String.format("%-5s %1s", "NAME", hmm.getName());
 538     writer.println((line));
 539
 540     if (hmm.getAccessionNumber() != null)
 541     {
 542     line = String.format("%-5s %1s", "ACC", hmm.getAccessionNumber());
 543       writer.println((line));
 544     }
 545
 546     if (hmm.getDescription() != null)
 547     {
 548     line = String.format("%-5s %1s", "DESC", hmm.getDescription());
 549       writer.println((line));
 550     }
 551     line = String.format("%-5s %1s", "LENG", hmm.getLength());
 552     writer.println((line));
 553
 554     if (hmm.getMaxInstanceLength() != null)
 555     {
 556     line = String.format("%-5s %1s", "MAXL", hmm.getMaxInstanceLength());
 557       writer.println((line));
 558     }
 559     line = String.format("%-5s %1s", "ALPH", hmm.getAlphabetType());
 560     writer.println((line));
 561
 562     boolean status;
 563     String statusStr;
 564
 565     status = hmm.referenceAnnotationIsActive();
 566     statusStr = HiddenMarkovModel.findStringFromBoolean(status);
 567     line = String.format("%-5s %1s", "RF",
 568             statusStr);
 569     writer.println((line));
 570
 571     status = hmm.maskValueIsActive();
 572     statusStr = HiddenMarkovModel.findStringFromBoolean(status);
 573     line = String.format("%-5s %1s", "MM",
 574             statusStr);
 575     writer.println((line));
 576
 577     status = hmm.consensusResidueIsActive();
 578     statusStr = HiddenMarkovModel.findStringFromBoolean(status);
 579     line = String.format("%-5s %1s", "CONS",
 580             statusStr);
 581     writer.println((line));
 582
 583     status = hmm.consensusStructureIsActive();
 584     statusStr = HiddenMarkovModel.findStringFromBoolean(status);
 585     line = String.format("%-5s %1s", "CS",
 586             statusStr);
 587     writer.println((line));
 588
 589     status = hmm.mapIsActive();
 590     statusStr = HiddenMarkovModel.findStringFromBoolean(status);
 591     line = String.format("%-5s %1s", "MAP",
 592             statusStr);
 593     writer.println((line));
 594
 595
 596     if (hmm.getDate() != null)
 597     {
 598     line = String.format("%-5s %1s", "DATE", hmm.getDate());
 599       writer.println((line));
 600     }
 601     if (hmm.getNumberOfSequences() != null)
 602     {
 603     line = String.format("%-5s %1s", "NSEQ", hmm.getNumberOfSequences());
 604       writer.println((line));
 605     }
 606     if (hmm.getEffectiveNumberOfSequences() != null)
 607     {
 608     line = String.format("%-5s %1s", "EFFN",
 609             hmm.getEffectiveNumberOfSequences());
 610       writer.println((line));
 611     }
 612     if (hmm.getCheckSum() != null)
 613     {
 614     line = String.format("%-5s %1s", "CKSUM", hmm.getCheckSum());
 615       writer.println((line));
 616     }
 617     if (hmm.getGatheringThreshold() != null)
 618     {
 619     line = String.format("%-5s %1s", "GA", hmm.getGatheringThreshold());
 620       writer.println((line));
 621     }
 622
 623     if (hmm.getTrustedCutoff() != null)
 624     {
 625     line = String.format("%-5s %1s", "TC", hmm.getTrustedCutoff());
 626       writer.println((line));
 627     }
 628     if (hmm.getNoiseCutoff() != null)
 629     {
 630     line = String.format("%-5s %1s", "NC", hmm.getNoiseCutoff());
 631       writer.println((line));
 632     }
 633     if (hmm.getMSV() != null)
 634     {
 635       line = String.format("%-19s %18s", "STATS LOCAL MSV", hmm.getMSV());
 636       writer.println((line));
 637
 638       line = String.format("%-19s %18s", "STATS LOCAL VITERBI",
 639               hmm.getViterbi());
 640       writer.println((line));
 641
 642       line = String.format("%-19s %18s", "STATS LOCAL FORWARD",
 643               hmm.getForward());
 644       writer.println((line));
 645     }
 646   }
 647
 648
 649   /**
 650    * Returns the char value of a single lettered String.
 651    *
 652    * @param string
 653    * @return
 654    */
 655   char charValue(String string)
 656   {
 657     char character;
 658     character = string.charAt(0);
 659     return character;
 660
 661   }
 662
 663   @Override
 664   public String print(SequenceI[] seqs, boolean jvsuffix)
 665   {
 666
 667     return null;
 668   }
 669
 670   /**
 671    * Converts the probabilities contained in a list into log space.
 672    *
 673    * @param list
 674    */
 675   List<Double> convertListToLogSpace(List<Double> list)
 676   {
 677
 678     List<Double> convertedList = new ArrayList<>();
 679     for (int i = 0; i < list.size(); i++)
 680     {
 681       double prob = list.get(i);
 682       double logProb = -1 * Math.log(prob);
 683
 684       convertedList.add(logProb);
 685     }
 686     return convertedList;
 687
 688
 689   }
 690 }
 691