src/jalview/datamodel/HiddenMarkovModel.java

   1 package jalview.datamodel;
   2
   3 import jalview.bin.Cache;
   4 import jalview.io.HMMFile;
   5 import jalview.schemes.ResidueProperties;
   6 import jalview.util.Comparison;
   7 import jalview.util.MapList;
   8
   9 import java.util.ArrayList;
  10 import java.util.Arrays;
  11 import java.util.HashMap;
  12 import java.util.List;
  13 import java.util.Map;
  14
  15 /**
  16  * Data structure which stores a hidden Markov model
  17  *
  18  * @author TZVanaalten
  19  *
  20  */
  21 public class HiddenMarkovModel
  22 {
  23   private static final char GAP_DASH = '-';
  24
  25   public final static String YES = "yes";
  26
  27   public final static String NO = "no";
  28
  29   public static final int MATCHTOMATCH = 0;
  30
  31   public static final int MATCHTOINSERT = 1;
  32
  33   public static final int MATCHTODELETE = 2;
  34
  35   public static final int INSERTTOMATCH = 3;
  36
  37   public static final int INSERTTOINSERT = 4;
  38
  39   public static final int DELETETOMATCH = 5;
  40
  41   public static final int DELETETODELETE = 6;
  42
  43   private static final double LOG2 = Math.log(2);
  44
  45   /*
  46    * properties read from HMM file header lines
  47    */
  48   private Map<String, String> fileProperties = new HashMap<>();
  49
  50   private String fileHeader;
  51
  52   /*
  53    * the symbols used in this model e.g. "ACGT"
  54    */
  55   private String alphabet;
  56
  57   /*
  58    * symbol lookup index into the alphabet for 'A' to 'Z'
  59    */
  60   private int[] symbolIndexLookup = new int['Z' - 'A' + 1];
  61
  62   /*
  63    * Nodes in the model. The begin node is at index 0, and contains
  64    * average emission probabilities for each symbol.
  65    */
  66   private List<HMMNode> nodes = new ArrayList<>();
  67
  68   /*
  69    * the aligned HMM consensus sequence extracted from the HMM profile
  70    */
  71   private SequenceI hmmSeq;
  72
  73   /*
  74    * mapping from HMM nodes to residues of the hmm consensus sequence
  75    */
  76   private Mapping mapToHmmConsensus;
  77
  78   /**
  79    * Constructor
  80    */
  81   public HiddenMarkovModel()
  82   {
  83   }
  84
  85   /**
  86    * Copy constructor given a new aligned sequence with which to associate the
  87    * HMM profile
  88    *
  89    * @param hmm
  90    * @param sq
  91    */
  92   public HiddenMarkovModel(HiddenMarkovModel hmm, SequenceI sq)
  93   {
  94     super();
  95     this.fileProperties = new HashMap<>(hmm.fileProperties);
  96     this.alphabet = hmm.alphabet;
  97     this.nodes = new ArrayList<>(hmm.nodes);
  98     this.symbolIndexLookup = hmm.symbolIndexLookup;
  99     this.fileHeader = new String(hmm.fileHeader);
 100     this.hmmSeq = sq;
 101     if (sq.getDatasetSequence() == hmm.mapToHmmConsensus.getTo())
 102     {
 103       this.mapToHmmConsensus = hmm.mapToHmmConsensus;
 104     }
 105     else
 106     {
 107       Cache.log.error("Error: HMM copied with change of mapped sequence");
 108     }
 109   }
 110
 111   /**
 112    * Returns the information content at a specified column, calculated as the
 113    * sum (over possible symbols) of the log ratio
 114    *
 115    * <pre>
 116    *  log(emission probability / background probability) / log(2)
 117    * </pre>
 118    *
 119    * @param column
 120    *          column position (base 0)
 121    * @return
 122    */
 123   public float getInformationContent(int column)
 124   {
 125     float informationContent = 0f;
 126
 127     for (char symbol : getSymbols().toCharArray())
 128     {
 129       float freq = ResidueProperties.backgroundFrequencies
 130               .get(getAlphabetType()).get(symbol);
 131       float prob = (float) getMatchEmissionProbability(column, symbol);
 132       informationContent += prob * Math.log(prob / freq);
 133     }
 134
 135     informationContent = informationContent / (float) LOG2;
 136
 137     return informationContent;
 138   }
 139
 140   /**
 141    * Gets the file header of the .hmm file this model came from
 142    *
 143    * @return
 144    */
 145   public String getFileHeader()
 146   {
 147     return fileHeader;
 148   }
 149
 150   /**
 151    * Sets the file header of this model.
 152    *
 153    * @param header
 154    */
 155   public void setFileHeader(String header)
 156   {
 157     fileHeader = header;
 158   }
 159
 160   /**
 161    * Returns the symbols used in this hidden Markov model
 162    *
 163    * @return
 164    */
 165   public String getSymbols()
 166   {
 167     return alphabet;
 168   }
 169
 170   /**
 171    * Gets the node in the hidden Markov model at the specified position.
 172    *
 173    * @param nodeIndex
 174    *          The index of the node requested. Node 0 optionally contains the
 175    *          average match emission probabilities across the entire model, and
 176    *          always contains the insert emission probabilities and state
 177    *          transition probabilities for the begin node. Node 1 contains the
 178    *          first node in the HMM that can correspond to a column in the
 179    *          alignment.
 180    * @return
 181    */
 182   public HMMNode getNode(int nodeIndex)
 183   {
 184     return nodes.get(nodeIndex);
 185   }
 186
 187   /**
 188    * Returns the name of the sequence alignment on which the HMM is based.
 189    *
 190    * @return
 191    */
 192   public String getName()
 193   {
 194     return fileProperties.get(HMMFile.NAME);
 195   }
 196
 197   /**
 198    * Answers the string value of the property (parsed from an HMM file) for the
 199    * given key, or null if the property is not present
 200    *
 201    * @param key
 202    * @return
 203    */
 204   public String getProperty(String key)
 205   {
 206     return fileProperties.get(key);
 207   }
 208
 209   /**
 210    * Answers true if the property with the given key is present with a value of
 211    * "yes" (not case-sensitive), else false
 212    *
 213    * @param key
 214    * @return
 215    */
 216   public boolean getBooleanProperty(String key)
 217   {
 218     return YES.equalsIgnoreCase(fileProperties.get(key));
 219   }
 220
 221   /**
 222    * Returns the length of the hidden Markov model. The value returned is the
 223    * LENG property if specified, else the number of nodes, excluding the begin
 224    * node (which should be the same thing).
 225    *
 226    * @return
 227    */
 228   public int getLength()
 229   {
 230     if (fileProperties.get(HMMFile.LENGTH) == null)
 231     {
 232       return nodes.size() - 1; // not counting BEGIN node
 233     }
 234     return Integer.parseInt(fileProperties.get(HMMFile.LENGTH));
 235   }
 236
 237   /**
 238    * Returns the value of mandatory property "ALPH" - "amino", "DNA", "RNA" are
 239    * the options. Other alphabets may be added.
 240    *
 241    * @return
 242    */
 243   public String getAlphabetType()
 244   {
 245     return fileProperties.get(HMMFile.ALPHABET);
 246   }
 247
 248   /**
 249    * Sets the model alphabet to the symbols in the given string (ignoring any
 250    * whitespace), and returns the number of symbols
 251    *
 252    * @param symbols
 253    */
 254   public int setAlphabet(String symbols)
 255   {
 256     String trimmed = symbols.toUpperCase().replaceAll("\\s", "");
 257     int count = trimmed.length();
 258     alphabet = trimmed;
 259     symbolIndexLookup = new int['Z' - 'A' + 1];
 260     Arrays.fill(symbolIndexLookup, -1);
 261     int ignored = 0;
 262
 263     /*
 264      * save the symbols in order, and a quick lookup of symbol position
 265      */
 266     for (short i = 0; i < count; i++)
 267     {
 268       char symbol = trimmed.charAt(i);
 269       if (symbol >= 'A' && symbol <= 'Z'
 270               && symbolIndexLookup[symbol - 'A'] == -1)
 271       {
 272         symbolIndexLookup[symbol - 'A'] = i;
 273       }
 274       else
 275       {
 276         System.err
 277                 .println(
 278                         "Unexpected or duplicated character in HMM ALPHabet: "
 279                                 + symbol);
 280         ignored++;
 281       }
 282     }
 283     return count - ignored;
 284   }
 285
 286   /**
 287    * Answers the node of the model corresponding to an aligned column position
 288    * (0...), or null if there is no such node
 289    *
 290    * @param column
 291    * @return
 292    */
 293   HMMNode getNodeForColumn(int column)
 294   {
 295     /*
 296      * if the hmm consensus is gapped at the column,
 297      * there is no corresponding node
 298      */
 299     if (Comparison.isGap(hmmSeq.getCharAt(column)))
 300     {
 301       return null;
 302     }
 303
 304     /*
 305      * find the node (if any) that is mapped to the
 306      * consensus sequence residue position at the column
 307      */
 308     int seqPos = hmmSeq.findPosition(column);
 309     int[] nodeNo = mapToHmmConsensus.getMap().locateInFrom(seqPos, seqPos);
 310     if (nodeNo != null)
 311     {
 312       return getNode(nodeNo[0]);
 313     }
 314     return null;
 315   }
 316
 317   /**
 318    * Gets the match emission probability for a given symbol at a column in the
 319    * alignment.
 320    *
 321    * @param alignColumn
 322    *          The index of the alignment column, starting at index 0. Index 0
 323    *          usually corresponds to index 1 in the HMM.
 324    * @param symbol
 325    *          The symbol for which the desired probability is being requested.
 326    * @return
 327    *
 328    */
 329   public double getMatchEmissionProbability(int alignColumn, char symbol)
 330   {
 331     HMMNode node = getNodeForColumn(alignColumn);
 332     int symbolIndex = getSymbolIndex(symbol);
 333     if (node != null && symbolIndex != -1)
 334     {
 335       return node.getMatchEmission(symbolIndex);
 336     }
 337     return 0D;
 338   }
 339
 340   /**
 341    * Gets the insert emission probability for a given symbol at a column in the
 342    * alignment.
 343    *
 344    * @param alignColumn
 345    *          The index of the alignment column, starting at index 0. Index 0
 346    *          usually corresponds to index 1 in the HMM.
 347    * @param symbol
 348    *          The symbol for which the desired probability is being requested.
 349    * @return
 350    *
 351    */
 352   public double getInsertEmissionProbability(int alignColumn, char symbol)
 353   {
 354     HMMNode node = getNodeForColumn(alignColumn);
 355     int symbolIndex = getSymbolIndex(symbol);
 356     if (node != null && symbolIndex != -1)
 357     {
 358       return node.getInsertEmission(symbolIndex);
 359     }
 360     return 0D;
 361   }
 362
 363   /**
 364    * Gets the state transition probability for a given symbol at a column in the
 365    * alignment.
 366    *
 367    * @param alignColumn
 368    *          The index of the alignment column, starting at index 0. Index 0
 369    *          usually corresponds to index 1 in the HMM.
 370    * @param symbol
 371    *          The symbol for which the desired probability is being requested.
 372    * @return
 373    *
 374    */
 375   public double getStateTransitionProbability(int alignColumn,
 376           int transition)
 377   {
 378     HMMNode node = getNodeForColumn(alignColumn);
 379     if (node != null)
 380     {
 381       return node.getStateTransition(transition);
 382     }
 383     return 0D;
 384   }
 385
 386   /**
 387    * Returns the sequence position linked to the node at the given index. This
 388    * corresponds to an aligned column position (counting from 1).
 389    *
 390    * @param nodeIndex
 391    *          The index of the node, starting from index 1. Index 0 is the begin
 392    *          node, which does not correspond to a column in the alignment.
 393    * @return
 394    */
 395   public int getNodeMapPosition(int nodeIndex)
 396   {
 397     return nodes.get(nodeIndex).getResidueNumber();
 398   }
 399
 400   /**
 401    * Returns the consensus residue at the specified node.
 402    *
 403    * @param nodeIndex
 404    *          The index of the specified node.
 405    * @return
 406    */
 407   public char getConsensusResidue(int nodeIndex)
 408   {
 409    char value = nodes.get(nodeIndex).getConsensusResidue();
 410    return value;
 411   }
 412
 413   /**
 414    * Returns the reference annotation at the specified node.
 415    *
 416    * @param nodeIndex
 417    *          The index of the specified node.
 418    * @return
 419    */
 420   public char getReferenceAnnotation(int nodeIndex)
 421   {
 422    char value = nodes.get(nodeIndex).getReferenceAnnotation();
 423    return value;
 424   }
 425
 426   /**
 427    * Returns the mask value at the specified node.
 428    *
 429    * @param nodeIndex
 430    *          The index of the specified node.
 431    * @return
 432    */
 433   public char getMaskedValue(int nodeIndex)
 434   {
 435    char value = nodes.get(nodeIndex).getMaskValue();
 436    return value;
 437   }
 438
 439   /**
 440    * Returns the consensus structure at the specified node.
 441    *
 442    * @param nodeIndex
 443    *          The index of the specified node.
 444    * @return
 445    */
 446   public char getConsensusStructure(int nodeIndex)
 447   {
 448    char value = nodes.get(nodeIndex).getConsensusStructure();
 449    return value;
 450   }
 451
 452   /**
 453    * Sets a property read from an HMM file
 454    *
 455    * @param key
 456    * @param value
 457    */
 458   public void setProperty(String key, String value)
 459   {
 460     fileProperties.put(key, value);
 461   }
 462
 463   /**
 464    * Temporary implementation, should not be used.
 465    *
 466    * @return
 467    */
 468   public String getViterbi()
 469   {
 470     String value;
 471     value = fileProperties.get(HMMFile.VITERBI);
 472     return value;
 473   }
 474
 475   /**
 476    * Temporary implementation, should not be used.
 477    *
 478    * @return
 479    */
 480   public String getMSV()
 481   {
 482     String value;
 483     value = fileProperties.get(HMMFile.MSV);
 484     return value;
 485   }
 486
 487   /**
 488    * Temporary implementation, should not be used.
 489    *
 490    * @return
 491    */
 492   public String getForward()
 493   {
 494     String value;
 495     value = fileProperties.get(HMMFile.FORWARD);
 496     return value;
 497   }
 498
 499   /**
 500    * Constructs the consensus sequence based on the most probable symbol at each
 501    * position. Gap characters are inserted for discontinuities in the node map
 502    * numbering (if provided), else an ungapped sequence is generated.
 503    * <p>
 504    * A mapping between the HMM nodes and residue positions of the sequence is
 505    * also built and saved.
 506    *
 507    * @return
 508    */
 509   void buildConsensusSequence()
 510   {
 511     List<int[]> toResidues = new ArrayList<>();
 512
 513     /*
 514      * if the HMM provided a map to sequence, use those start/end values,
 515      * else just treat it as for a contiguous sequence numbered from 1
 516      */
 517     boolean hasMap = getBooleanProperty(HMMFile.MAP);
 518     int start = hasMap ? getNode(1).getResidueNumber() : 1;
 519     int endResNo = hasMap ? getNode(nodes.size() - 1).getResidueNumber()
 520             : (start + getLength() - 1);
 521     char[] sequence = new char[endResNo - start + 1];
 522
 523     int lastResNo = start - 1;
 524     int seqOffset = 0;
 525     int gapCount = 0;
 526
 527     for (int nodeNo = 1; nodeNo < nodes.size(); nodeNo++)
 528     {
 529       HMMNode node = nodes.get(nodeNo);
 530       final int resNo = hasMap ? node.getResidueNumber() : lastResNo + 1;
 531
 532       /*
 533        * insert gaps if map numbering is not continuous
 534        */
 535       while (resNo > lastResNo + 1)
 536       {
 537         sequence[seqOffset++] = '-';
 538         lastResNo++;
 539         gapCount++;
 540       }
 541       char consensusResidue = node.getConsensusResidue();
 542       if (GAP_DASH == consensusResidue)
 543       {
 544         /*
 545          * no residue annotation in HMM - scan for the symbol
 546          * with the highest match emission probability
 547          */
 548         int symbolIndex = node.getMaxMatchEmissionIndex();
 549         consensusResidue = alphabet.charAt(symbolIndex);
 550         if (node.getMatchEmission(symbolIndex) < 0.5D)
 551         {
 552           // follow convention of lower case if match emission prob < 0.5
 553           consensusResidue = Character.toLowerCase(consensusResidue);
 554         }
 555       }
 556       sequence[seqOffset++] = consensusResidue;
 557       lastResNo = resNo;
 558     }
 559
 560     Sequence seq = new Sequence(getName(), sequence, start,
 561             lastResNo - gapCount);
 562     seq.createDatasetSequence();
 563     seq.setHMM(this);
 564     this.hmmSeq = seq;
 565
 566     /*
 567      * construct and store Mapping of nodes to residues
 568      * note as constructed this is just an identity mapping,
 569      * but it allows for greater flexibility in future
 570      */
 571     List<int[]> fromNodes = new ArrayList<>();
 572     fromNodes.add(new int[] { 1, getLength() });
 573     toResidues.add(new int[] { seq.getStart(), seq.getEnd() });
 574     MapList mapList = new MapList(fromNodes, toResidues, 1, 1);
 575     mapToHmmConsensus = new Mapping(seq.getDatasetSequence(), mapList);
 576   }
 577
 578
 579   /**
 580    * Answers the aligned consensus sequence for the profile. Note this will
 581    * return null if called before <code>setNodes</code> has been called.
 582    *
 583    * @return
 584    */
 585   public SequenceI getConsensusSequence()
 586   {
 587     return hmmSeq;
 588   }
 589
 590   /**
 591    * Answers the index position (0...) of the given symbol, or -1 if not a valid
 592    * symbol for this HMM
 593    *
 594    * @param symbol
 595    * @return
 596    */
 597   private int getSymbolIndex(char symbol)
 598   {
 599     /*
 600      * symbolIndexLookup holds the index for 'A' to 'Z'
 601      */
 602     char c = Character.toUpperCase(symbol);
 603     if ('A' <= c && c <= 'Z')
 604     {
 605       return symbolIndexLookup[c - 'A'];
 606     }
 607     return -1;
 608   }
 609
 610   /**
 611    * Sets the nodes of this HMM, and also extracts the HMM consensus sequence
 612    * and a mapping between node numbers and sequence positions
 613    *
 614    * @param nodeList
 615    */
 616   public void setNodes(List<HMMNode> nodeList)
 617   {
 618     nodes = nodeList;
 619     if (nodes.size() > 1)
 620     {
 621       buildConsensusSequence();
 622     }
 623   }
 624
 625   /**
 626    * Sets the aligned consensus sequence this HMM is the model for
 627    *
 628    * @param hmmSeq
 629    */
 630   public void setHmmSeq(SequenceI hmmSeq)
 631   {
 632     this.hmmSeq = hmmSeq;
 633   }
 634 }
 635