src/jalview/io/HMMFile.java

   1 package jalview.io;
   2
   3 import jalview.api.AlignExportSettingI;
   4 import jalview.api.AlignmentViewPanel;
   5 import jalview.datamodel.HMMNode;
   6 import jalview.datamodel.HiddenMarkovModel;
   7 import jalview.datamodel.SequenceI;
   8
   9 import java.io.BufferedReader;
  10 import java.io.IOException;
  11 import java.util.ArrayList;
  12 import java.util.List;
  13 import java.util.Scanner;
  14
  15
  16 /**
  17  * Adds capability to read in and write out HMMER3 files. .
  18  *
  19  *
  20  * @author TZVanaalten
  21  *
  22  */
  23 public class HMMFile extends AlignFile
  24         implements AlignmentFileReaderI, AlignmentFileWriterI
  25 {
  26   private static final String TERMINATOR = "//";
  27
  28   /*
  29    * keys to data in HMM file, used to store as properties of the HiddenMarkovModel
  30    */
  31   public static final String HMM = "HMM";
  32
  33   public static final String NAME = "NAME";
  34
  35   public static final String ACCESSION_NUMBER = "ACC";
  36
  37   public static final String DESCRIPTION = "DESC";
  38
  39   public static final String LENGTH = "LENG";
  40
  41   public static final String MAX_LENGTH = "MAXL";
  42
  43   public static final String ALPHABET = "ALPH";
  44
  45   public static final String DATE = "DATE";
  46
  47   public static final String COMMAND_LOG = "COM";
  48
  49   public static final String NUMBER_OF_SEQUENCES = "NSEQ";
  50
  51   public static final String EFF_NUMBER_OF_SEQUENCES = "EFFN";
  52
  53   public static final String CHECK_SUM = "CKSUM";
  54
  55   public static final String STATISTICS = "STATS";
  56
  57   public static final String COMPO = "COMPO";
  58
  59   public static final String GATHERING_THRESHOLD = "GA";
  60
  61   public static final String TRUSTED_CUTOFF = "TC";
  62
  63   public static final String NOISE_CUTOFF = "NC";
  64
  65   public static final String VITERBI = "VITERBI";
  66
  67   public static final String MSV = "MSV";
  68
  69   public static final String FORWARD = "FORWARD";
  70
  71   public static final String MAP = "MAP";
  72
  73   public static final String REFERENCE_ANNOTATION = "RF";
  74
  75   public static final String CONSENSUS_RESIDUE = "CONS";
  76
  77   public static final String CONSENSUS_STRUCTURE = "CS";
  78
  79   public static final String MASKED_VALUE = "MM";
  80
  81   private static final String ALPH_AMINO = "amino";
  82
  83   private static final String ALPH_DNA = "DNA";
  84
  85   private static final String ALPH_RNA = "RNA";
  86
  87   private static final String ALPHABET_AMINO = "ACDEFGHIKLMNPQRSTVWY";
  88
  89   private static final String ALPHABET_DNA = "ACGT";
  90
  91   private static final String ALPHABET_RNA = "ACGU";
  92
  93   private static final int NUMBER_OF_TRANSITIONS = 7;
  94
  95   private static final String SPACE = " ";
  96
  97   /*
  98    * optional guide line added to an output HMMER file, purely for readability
  99    */
 100   private static final String TRANSITIONTYPELINE = "            m->m     m->i     m->d     i->m     i->i     d->m     d->d";
 101
 102   private static String NL = System.lineSeparator();
 103
 104   private HiddenMarkovModel hmm;
 105
 106   // number of symbols in the alphabet used in the hidden Markov model
 107   private int numberOfSymbols;
 108
 109   /**
 110    * Constructor that parses immediately
 111    *
 112    * @param inFile
 113    * @param type
 114    * @throws IOException
 115    */
 116   public HMMFile(String inFile, DataSourceType type) throws IOException
 117   {
 118     super(inFile, type);
 119   }
 120
 121   /**
 122    * Constructor that parses immediately
 123    *
 124    * @param source
 125    * @throws IOException
 126    */
 127   public HMMFile(FileParse source) throws IOException
 128   {
 129     super(source);
 130   }
 131
 132   /**
 133    * Default constructor
 134    */
 135   public HMMFile()
 136   {
 137   }
 138
 139   /**
 140    * Constructor for HMMFile used for exporting
 141    *
 142    * @param hmm
 143    */
 144   public HMMFile(HiddenMarkovModel markov)
 145   {
 146     hmm = markov;
 147   }
 148
 149   /**
 150    * Returns the HMM produced by parsing a HMMER3 file
 151    *
 152    * @return
 153    */
 154   public HiddenMarkovModel getHMM()
 155   {
 156     return hmm;
 157   }
 158
 159   /**
 160    * Gets the name of the hidden Markov model
 161    *
 162    * @return
 163    */
 164   public String getName()
 165   {
 166     return hmm.getName();
 167   }
 168
 169   /**
 170    * Reads the data from HMM file into the HMM model
 171    */
 172   @Override
 173   public void parse()
 174   {
 175     try
 176     {
 177       hmm = new HiddenMarkovModel();
 178       parseHeaderLines(dataIn);
 179       parseModel(dataIn);
 180     } catch (Exception e)
 181     {
 182       e.printStackTrace();
 183     }
 184   }
 185
 186   /**
 187    * Reads the header properties from a HMMER3 file and saves them in the
 188    * HiddeMarkovModel. This method exits after reading the next line after the
 189    * HMM line.
 190    *
 191    * @param input
 192    * @throws IOException
 193    */
 194   void parseHeaderLines(BufferedReader input) throws IOException
 195   {
 196     boolean readingHeaders = true;
 197     hmm.setFileHeader(input.readLine());
 198     String line = input.readLine();
 199     while (readingHeaders && line != null)
 200     {
 201       Scanner parser = new Scanner(line);
 202       String next = parser.next();
 203       if (ALPHABET.equals(next))
 204       {
 205         String alphabetType = parser.next();
 206         hmm.setProperty(ALPHABET, alphabetType);
 207         String alphabet = ALPH_DNA.equalsIgnoreCase(alphabetType)
 208                 ? ALPHABET_DNA
 209                 : (ALPH_RNA.equalsIgnoreCase(alphabetType) ? ALPHABET_RNA
 210                         : ALPHABET_AMINO);
 211         numberOfSymbols = hmm.setAlphabet(alphabet);
 212       }
 213       else if (HMM.equals(next))
 214       {
 215         readingHeaders = false;
 216         String symbols = line.substring(line.indexOf(HMM) + HMM.length());
 217         numberOfSymbols = hmm.setAlphabet(symbols);
 218       }
 219       else if (STATISTICS.equals(next))
 220       {
 221         parser.next();
 222         String key;
 223         String value;
 224         key = parser.next();
 225         value = parser.next() + SPACE + SPACE + parser.next();
 226         hmm.setProperty(key, value);
 227       }
 228       else
 229       {
 230         String key = next;
 231         String value = parser.next();
 232         while (parser.hasNext())
 233         {
 234           value = value + SPACE + parser.next();
 235         }
 236         hmm.setProperty(key, value);
 237       }
 238       parser.close();
 239       line = input.readLine();
 240     }
 241   }
 242
 243   /**
 244    * Parses the model data from the HMMER3 file. The input buffer should be
 245    * positioned at the (optional) COMPO line if there is one, else at the insert
 246    * emissions line for the BEGIN node of the model.
 247    *
 248    * @param input
 249    * @throws IOException
 250    */
 251   void parseModel(BufferedReader input) throws IOException
 252   {
 253     /*
 254      * specification says there must always be an HMM header (already read)
 255      * and one more header (guide headings) which is skipped here
 256      */
 257     int nodeNo = 0;
 258     String line = input.readLine();
 259     List<HMMNode> nodes = new ArrayList<>();
 260
 261     while (line != null && !TERMINATOR.equals(line))
 262     {
 263       HMMNode node = new HMMNode();
 264       nodes.add(node);
 265       Scanner scanner = new Scanner(line);
 266       String next = scanner.next();
 267
 268       /*
 269        * expect COMPO (optional) for average match emissions
 270        * or a node number followed by node's match emissions
 271        */
 272       if (COMPO.equals(next) || nodeNo > 0)
 273       {
 274         /*
 275          * parse match emissions
 276          */
 277         double[] matches = parseDoubles(scanner, numberOfSymbols);
 278         node.setMatchEmissions(matches);
 279         if (!COMPO.equals(next))
 280         {
 281           int resNo = parseAnnotations(scanner, node);
 282           if (resNo == 0)
 283           {
 284             /*
 285              * no MAP annotation provided, just number off from 0 (begin node)
 286              */
 287             resNo = nodeNo;
 288           }
 289           node.setResidueNumber(resNo);
 290         }
 291         line = input.readLine();
 292       }
 293       scanner.close();
 294
 295       /*
 296        * parse insert emissions
 297        */
 298       scanner = new Scanner(line);
 299       double[] inserts = parseDoubles(scanner, numberOfSymbols);
 300       node.setInsertEmissions(inserts);
 301       scanner.close();
 302
 303       /*
 304        * parse state transitions
 305        */
 306       line = input.readLine();
 307       scanner = new Scanner(line);
 308       double[] transitions = parseDoubles(scanner,
 309               NUMBER_OF_TRANSITIONS);
 310       node.setStateTransitions(transitions);
 311       scanner.close();
 312       line = input.readLine();
 313
 314       nodeNo++;
 315     }
 316
 317     hmm.setNodes(nodes);
 318   }
 319
 320   /**
 321    * Parses the annotations on the match emission line and add them to the node.
 322    * (See p109 of the HMMER User Guide (V3.1b2) for the specification.) Returns
 323    * the residue position that the node maps to, if provided, else zero.
 324    *
 325    * @param scanner
 326    * @param node
 327    */
 328   int parseAnnotations(Scanner scanner, HMMNode node)
 329   {
 330     int mapTo = 0;
 331
 332     /*
 333      * map from hmm node to sequence position, if provided
 334      */
 335     if (scanner.hasNext())
 336     {
 337       String value = scanner.next();
 338       if (!"-".equals(value))
 339       {
 340         try
 341         {
 342           mapTo = Integer.parseInt(value);
 343           node.setResidueNumber(mapTo);
 344         } catch (NumberFormatException e)
 345         {
 346           // ignore
 347         }
 348       }
 349     }
 350
 351     /*
 352      * hmm consensus residue if provided, else '-'
 353      */
 354     if (scanner.hasNext())
 355     {
 356       node.setConsensusResidue(scanner.next().charAt(0));
 357     }
 358
 359     /*
 360      * RF reference annotation, if provided, else '-'
 361      */
 362     if (scanner.hasNext())
 363     {
 364       node.setReferenceAnnotation(scanner.next().charAt(0));
 365     }
 366
 367     /*
 368      * 'm' for masked position, if provided, else '-'
 369      */
 370     if (scanner.hasNext())
 371     {
 372       node.setMaskValue(scanner.next().charAt(0));
 373     }
 374
 375     /*
 376      * structure consensus symbol, if provided, else '-'
 377      */
 378     if (scanner.hasNext())
 379     {
 380       node.setConsensusStructure(scanner.next().charAt(0));
 381     }
 382
 383     return mapTo;
 384   }
 385
 386   /**
 387    * Fills an array of doubles parsed from an input line
 388    *
 389    * @param input
 390    * @param numberOfElements
 391    * @return
 392    * @throws IOException
 393    */
 394   static double[] parseDoubles(Scanner input,
 395           int numberOfElements) throws IOException
 396   {
 397     double[] values = new double[numberOfElements];
 398     for (int i = 0; i < numberOfElements; i++)
 399     {
 400       if (!input.hasNext())
 401       {
 402         throw new IOException("Incomplete data");
 403       }
 404       String next = input.next();
 405       if (next.contains("*"))
 406       {
 407         values[i] = Double.NEGATIVE_INFINITY;
 408       }
 409       else
 410       {
 411         double prob = Double.valueOf(next);
 412         prob = Math.pow(Math.E, -prob);
 413         values[i] = prob;
 414       }
 415     }
 416     return values;
 417   }
 418
 419   /**
 420    * Returns a string to be added to the StringBuilder containing the entire
 421    * output String.
 422    *
 423    * @param initialColumnSeparation
 424    *          The initial whitespace separation between the left side of the
 425    *          file and first character.
 426    * @param columnSeparation
 427    *          The separation between subsequent data entries.
 428    * @param data
 429    *          The list of data to be added to the String.
 430    * @return
 431    */
 432   String addData(int initialColumnSeparation,
 433           int columnSeparation, List<String> data)
 434   {
 435     String line = "";
 436     boolean first = true;
 437     for (String value : data)
 438     {
 439       int sep = first ? initialColumnSeparation : columnSeparation;
 440       line += String.format("%" + sep + "s", value);
 441       first = false;
 442     }
 443     return line;
 444   }
 445
 446   /**
 447    * Converts list of characters into a list of Strings.
 448    *
 449    * @param list
 450    * @return Returns the list of Strings.
 451    */
 452   List<String> charListToStringList(List<Character> list)
 453   {
 454     List<String> strList = new ArrayList<>();
 455     for (char value : list)
 456     {
 457       String strValue = Character.toString(value);
 458       strList.add(strValue);
 459     }
 460     return strList;
 461   }
 462
 463   /**
 464    * Converts an array of doubles into a list of Strings, rounded to the nearest
 465    * 5th decimal place
 466    *
 467    * @param doubles
 468    * @param noOfDecimals
 469    * @return
 470    */
 471   List<String> doublesToStringList(double[] doubles)
 472   {
 473     List<String> strList = new ArrayList<>();
 474     for (double value : doubles)
 475     {
 476       String strValue;
 477       if (value > 0)
 478       {
 479         strValue = String.format("%.5f", value);
 480       }
 481       else if (value == -0.00000d)
 482       {
 483         strValue = "0.00000";
 484       }
 485       else
 486       {
 487         strValue = "*";
 488       }
 489       strList.add(strValue);
 490     }
 491     return strList;
 492   }
 493
 494   /**
 495    * Appends model data in string format to the string builder
 496    *
 497    * @param output
 498    */
 499   void appendModelAsString(StringBuilder output)
 500   {
 501     output.append(HMM).append("  ");
 502     String charSymbols = hmm.getSymbols();
 503     for (char c : charSymbols.toCharArray())
 504     {
 505       output.append(String.format("%9s", c));
 506     }
 507     output.append(NL).append(TRANSITIONTYPELINE);
 508
 509     int length = hmm.getLength();
 510
 511     for (int nodeNo = 0; nodeNo <= length; nodeNo++)
 512     {
 513       String matchLine = String.format("%7s",
 514               nodeNo == 0 ? COMPO : Integer.toString(nodeNo));
 515
 516       double[] doubleMatches = convertToLogSpace(
 517               hmm.getNode(nodeNo).getMatchEmissions());
 518       List<String> strMatches = doublesToStringList(doubleMatches);
 519       matchLine += addData(10, 9, strMatches);
 520
 521       if (nodeNo != 0)
 522       {
 523         matchLine += SPACE + (hmm.getNodeMapPosition(nodeNo));
 524         matchLine += SPACE + hmm.getConsensusResidue(nodeNo);
 525         matchLine += SPACE + hmm.getReferenceAnnotation(nodeNo);
 526         if (hmm.getFileHeader().contains("HMMER3/f"))
 527         {
 528           matchLine += SPACE + hmm.getMaskedValue(nodeNo);
 529           matchLine += SPACE + hmm.getConsensusStructure(nodeNo);
 530         }
 531       }
 532
 533       output.append(NL).append(matchLine);
 534
 535       String insertLine = "";
 536
 537       double[] doubleInserts = convertToLogSpace(
 538               hmm.getNode(nodeNo).getInsertEmissions());
 539       List<String> strInserts = doublesToStringList(doubleInserts);
 540       insertLine += addData(17, 9, strInserts);
 541
 542       output.append(NL).append(insertLine);
 543
 544       String transitionLine = "";
 545       double[] doubleTransitions = convertToLogSpace(
 546               hmm.getNode(nodeNo).getStateTransitions());
 547       List<String> strTransitions = doublesToStringList(
 548               doubleTransitions);
 549       transitionLine += addData(17, 9, strTransitions);
 550
 551       output.append(NL).append(transitionLine);
 552     }
 553   }
 554
 555   /**
 556    * Appends formatted HMM file properties to the string builder
 557    *
 558    * @param output
 559    */
 560   void appendProperties(StringBuilder output)
 561   {
 562     output.append(hmm.getFileHeader());
 563
 564     String format = "%n%-5s %1s";
 565     appendProperty(output, format, NAME);
 566     appendProperty(output, format, ACCESSION_NUMBER);
 567     appendProperty(output, format, DESCRIPTION);
 568     appendProperty(output, format, LENGTH);
 569     appendProperty(output, format, MAX_LENGTH);
 570     appendProperty(output, format, ALPHABET);
 571     appendBooleanProperty(output, format, REFERENCE_ANNOTATION);
 572     appendBooleanProperty(output, format, MASKED_VALUE);
 573     appendBooleanProperty(output, format, CONSENSUS_RESIDUE);
 574     appendBooleanProperty(output, format, CONSENSUS_STRUCTURE);
 575     appendBooleanProperty(output, format, MAP);
 576     appendProperty(output, format, DATE);
 577     appendProperty(output, format, NUMBER_OF_SEQUENCES);
 578     appendProperty(output, format, EFF_NUMBER_OF_SEQUENCES);
 579     appendProperty(output, format, CHECK_SUM);
 580     appendProperty(output, format, GATHERING_THRESHOLD);
 581     appendProperty(output, format, TRUSTED_CUTOFF);
 582     appendProperty(output, format, NOISE_CUTOFF);
 583
 584     if (hmm.getMSV() != null)
 585     {
 586       format = "%n%-19s %18s";
 587       output.append(String.format(format, "STATS LOCAL MSV", hmm.getMSV()));
 588
 589       output.append(String.format(format, "STATS LOCAL VITERBI",
 590               hmm.getViterbi()));
 591
 592       output.append(String.format(format, "STATS LOCAL FORWARD",
 593               hmm.getForward()));
 594     }
 595   }
 596
 597   /**
 598    * Appends 'yes' or 'no' for the given property, according to whether or not
 599    * it is set in the HMM
 600    *
 601    * @param output
 602    * @param format
 603    * @param propertyName
 604    */
 605   private void appendBooleanProperty(StringBuilder output, String format,
 606           String propertyName)
 607   {
 608     boolean set = hmm.getBooleanProperty(propertyName);
 609     output.append(String.format(format, propertyName,
 610             set ? HiddenMarkovModel.YES : HiddenMarkovModel.NO));
 611   }
 612
 613   /**
 614    * Appends the value of the given property to the output, if not null
 615    *
 616    * @param output
 617    * @param format
 618    * @param propertyName
 619    */
 620   private void appendProperty(StringBuilder output, String format,
 621           String propertyName)
 622   {
 623     String value = hmm.getProperty(propertyName);
 624     if (value != null)
 625     {
 626       output.append(String.format(format, propertyName, value));
 627     }
 628   }
 629
 630   @Override
 631   public String print(SequenceI[] sequences, boolean jvsuffix)
 632   {
 633     if (sequences[0].getHMM() != null)
 634     {
 635       hmm = sequences[0].getHMM();
 636     }
 637     return print();
 638   }
 639
 640   /**
 641    * Prints the .hmm file to a String.
 642    *
 643    * @return
 644    */
 645   public String print()
 646   {
 647     StringBuilder output = new StringBuilder();
 648     appendProperties(output);
 649     output.append(NL);
 650     appendModelAsString(output);
 651     output.append(NL).append(TERMINATOR).append(NL);
 652     return output.toString();
 653   }
 654
 655   /**
 656    * Converts the probabilities contained in an array into log space
 657    *
 658    * @param ds
 659    */
 660   double[] convertToLogSpace(double[] ds)
 661   {
 662     double[] converted = new double[ds.length];
 663     for (int i = 0; i < ds.length; i++)
 664     {
 665       double prob = ds[i];
 666       double logProb = -1 * Math.log(prob);
 667
 668       converted[i] = logProb;
 669     }
 670     return converted;
 671   }
 672
 673   /**
 674    * Returns the HMM sequence produced by reading a .hmm file.
 675    */
 676   @Override
 677   public SequenceI[] getSeqsAsArray()
 678   {
 679     SequenceI hmmSeq = hmm.getConsensusSequence();
 680     SequenceI[] seq = new SequenceI[1];
 681     seq[0] = hmmSeq;
 682     return seq;
 683   }
 684
 685   @Override
 686   public void setNewlineString(String newLine)
 687   {
 688     NL = newLine;
 689   }
 690
 691   @Override
 692   public void setExportSettings(AlignExportSettingI exportSettings)
 693   {
 694
 695   }
 696
 697   @Override
 698   public void configureForView(AlignmentViewPanel viewpanel)
 699   {
 700
 701   }
 702
 703   @Override
 704   public boolean hasWarningMessage()
 705   {
 706     return false;
 707   }
 708
 709   @Override
 710   public String getWarningMessage()
 711   {
 712     return "warning message";
 713   }
 714
 715 }
 716