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