src/jalview/datamodel/HiddenMarkovModel.java

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