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