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     while (line != null && !TERMINATOR.equals(line))
 261     {
 262       HMMNode node = new HMMNode();
 263       hmm.addNode(node);
 264       Scanner scanner = new Scanner(line);
 265       String next = scanner.next();
 266
 267       /*
 268        * expect COMPO (optional) for average match emissions
 269        * or a node number followed by node's match emissions
 270        */
 271       if (COMPO.equals(next) || nodeNo > 0)
 272       {
 273         /*
 274          * parse match emissions
 275          */
 276         double[] matches = parseDoubles(scanner, numberOfSymbols);
 277         node.setMatchEmissions(matches);
 278         if (!COMPO.equals(next))
 279         {
 280           int column = parseAnnotations(scanner, node);
 281           if (column == 0)
 282           {
 283             /*
 284              * no MAP annotation provided, just number off from 0 (begin node)
 285              */
 286             column = nodeNo;
 287           }
 288           hmm.setAlignmentColumn(node, column - 1); // node 1 <==> column 0
 289         }
 290         line = input.readLine();
 291       }
 292       scanner.close();
 293
 294       /*
 295        * parse insert emissions
 296        */
 297       scanner = new Scanner(line);
 298       double[] inserts = parseDoubles(scanner, numberOfSymbols);
 299       node.setInsertEmissions(inserts);
 300       scanner.close();
 301
 302       /*
 303        * parse state transitions
 304        */
 305       line = input.readLine();
 306       scanner = new Scanner(line);
 307       double[] transitions = parseDoubles(scanner,
 308               NUMBER_OF_TRANSITIONS);
 309       node.setStateTransitions(transitions);
 310       scanner.close();
 311       line = input.readLine();
 312
 313       nodeNo++;
 314     }
 315   }
 316
 317   /**
 318    * Parses the annotations on the match emission line and add them to the node.
 319    * (See p109 of the HMMER User Guide (V3.1b2) for the specification.) Returns
 320    * the alignment column number (base 1) that the node maps to, if provided,
 321    * else zero.
 322    *
 323    * @param scanner
 324    * @param node
 325    */
 326   int parseAnnotations(Scanner scanner, HMMNode node)
 327   {
 328     /*
 329      * map from hmm node to alignment column index, if provided
 330      * HMM counts columns from 1, convert to base 0 for Jalview
 331      */
 332     int column = 0;
 333     String value;
 334     if (scanner.hasNext())
 335     {
 336       value = scanner.next();
 337       if (!"-".equals(value))
 338       {
 339         try
 340         {
 341           column = Integer.parseInt(value);
 342           node.setAlignmentColumn(column - 1);
 343         } catch (NumberFormatException e)
 344         {
 345           // ignore
 346         }
 347       }
 348     }
 349
 350     /*
 351      * hmm consensus residue if provided, else '-'
 352      */
 353     if (scanner.hasNext())
 354     {
 355       node.setConsensusResidue(scanner.next().charAt(0));
 356     }
 357
 358     /*
 359      * RF reference annotation, if provided, else '-'
 360      */
 361     if (scanner.hasNext())
 362     {
 363       node.setReferenceAnnotation(scanner.next().charAt(0));
 364     }
 365
 366     /*
 367      * 'm' for masked position, if provided, else '-'
 368      */
 369     if (scanner.hasNext())
 370     {
 371       node.setMaskValue(scanner.next().charAt(0));
 372     }
 373
 374     /*
 375      * structure consensus symbol, if provided, else '-'
 376      */
 377     if (scanner.hasNext())
 378     {
 379       node.setConsensusStructure(scanner.next().charAt(0));
 380     }
 381
 382     return column;
 383   }
 384
 385   /**
 386    * Fills an array of doubles parsed from an input line
 387    *
 388    * @param input
 389    * @param numberOfElements
 390    * @return
 391    * @throws IOException
 392    */
 393   static double[] parseDoubles(Scanner input,
 394           int numberOfElements) throws IOException
 395   {
 396     double[] values = new double[numberOfElements];
 397     for (int i = 0; i < numberOfElements; i++)
 398     {
 399       if (!input.hasNext())
 400       {
 401         throw new IOException("Incomplete data");
 402       }
 403       String next = input.next();
 404       if (next.contains("*"))
 405       {
 406         values[i] = Double.NEGATIVE_INFINITY;
 407       }
 408       else
 409       {
 410         double prob = Double.valueOf(next);
 411         prob = Math.pow(Math.E, -prob);
 412         values[i] = prob;
 413       }
 414     }
 415     return values;
 416   }
 417
 418   /**
 419    * Returns a string to be added to the StringBuilder containing the entire
 420    * output String.
 421    *
 422    * @param initialColumnSeparation
 423    *          The initial whitespace separation between the left side of the
 424    *          file and first character.
 425    * @param columnSeparation
 426    *          The separation between subsequent data entries.
 427    * @param data
 428    *          The list of data to be added to the String.
 429    * @return
 430    */
 431   String addData(int initialColumnSeparation,
 432           int columnSeparation, List<String> data)
 433   {
 434     String line = "";
 435     boolean first = true;
 436     for (String value : data)
 437     {
 438       int sep = first ? initialColumnSeparation : columnSeparation;
 439       line += String.format("%" + sep + "s", value);
 440       first = false;
 441     }
 442     return line;
 443   }
 444
 445   /**
 446    * Converts list of characters into a list of Strings.
 447    *
 448    * @param list
 449    * @return Returns the list of Strings.
 450    */
 451   List<String> charListToStringList(List<Character> list)
 452   {
 453     List<String> strList = new ArrayList<>();
 454     for (char value : list)
 455     {
 456       String strValue = Character.toString(value);
 457       strList.add(strValue);
 458     }
 459     return strList;
 460   }
 461
 462   /**
 463    * Converts an array of doubles into a list of Strings, rounded to the nearest
 464    * 5th decimal place
 465    *
 466    * @param doubles
 467    * @param noOfDecimals
 468    * @return
 469    */
 470   List<String> doublesToStringList(double[] doubles)
 471   {
 472     List<String> strList = new ArrayList<>();
 473     for (double value : doubles)
 474     {
 475       String strValue;
 476       if (value > 0)
 477       {
 478         strValue = String.format("%.5f", value);
 479       }
 480       else if (value == -0.00000d)
 481       {
 482         strValue = "0.00000";
 483       }
 484       else
 485       {
 486         strValue = "*";
 487       }
 488       strList.add(strValue);
 489     }
 490     return strList;
 491   }
 492
 493   /**
 494    * Appends model data in string format to the string builder
 495    *
 496    * @param output
 497    */
 498   void appendModelAsString(StringBuilder output)
 499   {
 500     output.append(HMM).append("  ");
 501     String charSymbols = hmm.getSymbols();
 502     for (char c : charSymbols.toCharArray())
 503     {
 504       output.append(String.format("%9s", c));
 505     }
 506     output.append(NL).append(TRANSITIONTYPELINE);
 507
 508     int length = hmm.getLength();
 509
 510     for (int nodeNo = 0; nodeNo <= length; nodeNo++)
 511     {
 512       String matchLine = String.format("%7s",
 513               nodeNo == 0 ? COMPO : Integer.toString(nodeNo));
 514
 515       double[] doubleMatches = convertToLogSpace(
 516               hmm.getNode(nodeNo).getMatchEmissions());
 517       List<String> strMatches = doublesToStringList(doubleMatches);
 518       matchLine += addData(10, 9, strMatches);
 519
 520       if (nodeNo != 0)
 521       {
 522         matchLine += SPACE + (hmm.getNodeAlignmentColumn(nodeNo) + 1);
 523         matchLine += SPACE + hmm.getConsensusResidue(nodeNo);
 524         matchLine += SPACE + hmm.getReferenceAnnotation(nodeNo);
 525         if (hmm.getFileHeader().contains("HMMER3/f"))
 526         {
 527           matchLine += SPACE + hmm.getMaskedValue(nodeNo);
 528           matchLine += SPACE + hmm.getConsensusStructure(nodeNo);
 529         }
 530       }
 531
 532       output.append(NL).append(matchLine);
 533
 534       String insertLine = "";
 535
 536       double[] doubleInserts = convertToLogSpace(
 537               hmm.getNode(nodeNo).getInsertEmissions());
 538       List<String> strInserts = doublesToStringList(doubleInserts);
 539       insertLine += addData(17, 9, strInserts);
 540
 541       output.append(NL).append(insertLine);
 542
 543       String transitionLine = "";
 544       double[] doubleTransitions = convertToLogSpace(
 545               hmm.getNode(nodeNo).getStateTransitions());
 546       List<String> strTransitions = doublesToStringList(
 547               doubleTransitions);
 548       transitionLine += addData(17, 9, strTransitions);
 549
 550       output.append(NL).append(transitionLine);
 551     }
 552   }
 553
 554   /**
 555    * Appends formatted HMM file properties to the string builder
 556    *
 557    * @param output
 558    */
 559   void appendProperties(StringBuilder output)
 560   {
 561     output.append(hmm.getFileHeader());
 562
 563     String format = "%n%-5s %1s";
 564     appendProperty(output, format, NAME);
 565     appendProperty(output, format, ACCESSION_NUMBER);
 566     appendProperty(output, format, DESCRIPTION);
 567     appendProperty(output, format, LENGTH);
 568     appendProperty(output, format, MAX_LENGTH);
 569     appendProperty(output, format, ALPHABET);
 570     appendBooleanProperty(output, format, REFERENCE_ANNOTATION);
 571     appendBooleanProperty(output, format, MASKED_VALUE);
 572     appendBooleanProperty(output, format, CONSENSUS_RESIDUE);
 573     appendBooleanProperty(output, format, CONSENSUS_STRUCTURE);
 574     appendBooleanProperty(output, format, MAP);
 575     appendProperty(output, format, DATE);
 576     appendProperty(output, format, NUMBER_OF_SEQUENCES);
 577     appendProperty(output, format, EFF_NUMBER_OF_SEQUENCES);
 578     appendProperty(output, format, CHECK_SUM);
 579     appendProperty(output, format, GATHERING_THRESHOLD);
 580     appendProperty(output, format, TRUSTED_CUTOFF);
 581     appendProperty(output, format, NOISE_CUTOFF);
 582
 583     if (hmm.getMSV() != null)
 584     {
 585       format = "%n%-19s %18s";
 586       output.append(String.format(format, "STATS LOCAL MSV", hmm.getMSV()));
 587
 588       output.append(String.format(format, "STATS LOCAL VITERBI",
 589               hmm.getViterbi()));
 590
 591       output.append(String.format(format, "STATS LOCAL FORWARD",
 592               hmm.getForward()));
 593     }
 594   }
 595
 596   /**
 597    * Appends 'yes' or 'no' for the given property, according to whether or not
 598    * it is set in the HMM
 599    *
 600    * @param output
 601    * @param format
 602    * @param propertyName
 603    */
 604   private void appendBooleanProperty(StringBuilder output, String format,
 605           String propertyName)
 606   {
 607     boolean set = hmm.getBooleanProperty(propertyName);
 608     output.append(String.format(format, propertyName,
 609             set ? HiddenMarkovModel.YES : HiddenMarkovModel.NO));
 610   }
 611
 612   /**
 613    * Appends the value of the given property to the output, if not null
 614    *
 615    * @param output
 616    * @param format
 617    * @param propertyName
 618    */
 619   private void appendProperty(StringBuilder output, String format,
 620           String propertyName)
 621   {
 622     String value = hmm.getProperty(propertyName);
 623     if (value != null)
 624     {
 625       output.append(String.format(format, propertyName, value));
 626     }
 627   }
 628
 629   @Override
 630   public String print(SequenceI[] sequences, boolean jvsuffix)
 631   {
 632     if (sequences[0].getHMM() != null)
 633     {
 634       hmm = sequences[0].getHMM();
 635     }
 636     return print();
 637   }
 638
 639   /**
 640    * Prints the .hmm file to a String.
 641    *
 642    * @return
 643    */
 644   public String print()
 645   {
 646     StringBuilder output = new StringBuilder();
 647     appendProperties(output);
 648     output.append(NL);
 649     appendModelAsString(output);
 650     output.append(NL).append(TERMINATOR).append(NL);
 651     return output.toString();
 652   }
 653
 654   /**
 655    * Converts the probabilities contained in an array into log space
 656    *
 657    * @param ds
 658    */
 659   double[] convertToLogSpace(double[] ds)
 660   {
 661     double[] converted = new double[ds.length];
 662     for (int i = 0; i < ds.length; i++)
 663     {
 664       double prob = ds[i];
 665       double logProb = -1 * Math.log(prob);
 666
 667       converted[i] = logProb;
 668     }
 669     return converted;
 670   }
 671
 672   /**
 673    * Returns the HMM sequence produced by reading a .hmm file.
 674    */
 675   @Override
 676   public SequenceI[] getSeqsAsArray()
 677   {
 678     SequenceI hmmSeq = hmm.initHMMSequence();
 679     SequenceI[] seq = new SequenceI[1];
 680     seq[0] = hmmSeq;
 681     return seq;
 682   }
 683
 684   @Override
 685   public void setNewlineString(String newLine)
 686   {
 687     NL = newLine;
 688   }
 689
 690   @Override
 691   public void setExportSettings(AlignExportSettingI exportSettings)
 692   {
 693
 694   }
 695
 696   @Override
 697   public void configureForView(AlignmentViewPanel viewpanel)
 698   {
 699
 700   }
 701
 702   @Override
 703   public boolean hasWarningMessage()
 704   {
 705     return false;
 706   }
 707
 708   @Override
 709   public String getWarningMessage()
 710   {
 711     return "warning message";
 712   }
 713
 714 }
 715