src/jalview/datamodel/HiddenMarkovModel.java

   1 package jalview.datamodel;
   2
   3 import jalview.schemes.ResidueProperties;
   4
   5 import java.util.ArrayList;
   6 import java.util.HashMap;
   7 import java.util.List;
   8 import java.util.Map;
   9 import java.util.Scanner;
  10
  11 /**
  12  * Data structure which stores a hidden Markov model. Currently contains file properties as well, not sure whether these should be transferred to the HMMFile class
  13  *
  14  * @author TZVanaalten
  15  *
  16  */
  17 public class HiddenMarkovModel
  18 {
  19   // Stores file properties. Do not directly access this field as it contains
  20   // only string value - use the getter methods. For example, to find the length
  21   // of theHMM, use getModelLength()to return an int value
  22   Map<String, String> fileProperties = new HashMap<>();
  23
  24   //contains all of the symbols used in this model. The index of each symbol represents its lookup value
  25   List<Character> symbols = new ArrayList<>();
  26
  27   // contains information for each node in the model. The begin node is at index
  28   // 0. Node 0 contains average emission probabilities for each symbol
  29   List<HMMNode> nodes = new ArrayList<>();
  30
  31   // contains the HMM node for each alignment column
  32   Map<Integer, Integer> nodeLookup = new HashMap<>();
  33
  34   //contains the symbol index for each symbol
  35   Map<Character, Integer> symbolIndexLookup = new HashMap<>();
  36
  37   Map<Character, Double> backgroundFrequencies = new HashMap();
  38
  39
  40   final static String YES = "yes";
  41
  42   final static String NO = "no";
  43
  44   int numberOfSymbols;
  45
  46   //keys for file properties hashmap
  47   private final String NAME = "NAME";
  48
  49   private final String ACCESSION_NUMBER = "ACC";
  50
  51   private final String DESCRIPTION = "DESC";
  52
  53   private final String LENGTH = "LENG";
  54
  55   private final String MAX_LENGTH = "MAXL";
  56
  57   private final String ALPHABET = "ALPH";
  58
  59   private final String DATE = "DATE";
  60
  61   private final String COMMAND_LOG = "COM";
  62
  63   private final String NUMBER_OF_SEQUENCES = "NSEQ";
  64
  65   private final String EFF_NUMBER_OF_SEQUENCES = "EFFN";
  66
  67   private final String CHECK_SUM = "CKSUM";
  68
  69   private final String GATHERING_THRESHOLDS = "GA";
  70
  71   private final String TRUSTED_CUTOFFS = "TC";
  72
  73   private final String NOISE_CUTOFFS = "NC";
  74
  75   private final String STATISTICS = "STATS";
  76
  77   private final String COMPO = "COMPO";
  78
  79   private final String GATHERING_THRESHOLD = "GA";
  80
  81   private final String TRUSTED_CUTOFF = "TC";
  82
  83   private final String NOISE_CUTOFF = "NC";
  84
  85   private final String VITERBI = "VITERBI";
  86
  87   private final String MSV = "MSV";
  88
  89   private final String FORWARD = "FORWARD";
  90
  91   private final String MAP = "MAP";
  92
  93   private final String REFERENCE_ANNOTATION = "RF";
  94
  95   private final String CONSENSUS_RESIDUE = "CONS";
  96
  97   private final String CONSENSUS_STRUCTURE = "CS";
  98
  99   private final String MASKED_VALUE = "MM";
 100
 101   final static String[] TRANSITION_TYPES = new String[] { "m->m", "m->i",
 102       "m->d", "i->m", "i->i", "d->m", "d->d" };
 103
 104   public String getTransitionType(int index)
 105   {
 106     return TRANSITION_TYPES[index];
 107   }
 108
 109   public Map<Integer, Integer> getNodeLookup()
 110   {
 111     return nodeLookup;
 112   }
 113
 114   public void setNodeLookup(Map<Integer, Integer> nodeLookup)
 115   {
 116     this.nodeLookup = nodeLookup;
 117   }
 118
 119   public String[] getTransitionTypes()
 120   {
 121     return TRANSITION_TYPES;
 122   }
 123
 124   public List<Character> getSymbols()
 125   {
 126     return symbols;
 127   }
 128
 129   public Map<String, String> getFileProperties()
 130   {
 131     return fileProperties;
 132   }
 133
 134   public HMMNode getNode(int nodeIndex)
 135   {
 136     return getNodes().get(nodeIndex);
 137   }
 138
 139   public void setSymbols(List<Character> symbolsL)
 140   {
 141     this.symbols = symbolsL;
 142   }
 143
 144   public String getName()
 145   {
 146     return fileProperties.get(NAME);
 147   }
 148   public String getAccessionNumber()
 149   {
 150     return fileProperties.get(ACCESSION_NUMBER);
 151   }
 152
 153   public void setAccessionNumber(String value)
 154   {
 155     fileProperties.put(ACCESSION_NUMBER, value);
 156   }
 157
 158   public String getDescription()
 159   {
 160     return fileProperties.get(DESCRIPTION);
 161   }
 162
 163   public void setDescription(String value)
 164   {
 165     fileProperties.put(DESCRIPTION, value);
 166   }
 167
 168   public Integer getLength()
 169   {
 170     if (fileProperties.get(LENGTH) == null)
 171     {
 172       return null;
 173     }
 174     return Integer.parseInt(fileProperties.get(LENGTH));
 175   }
 176
 177   public void setLength(int value)
 178   {
 179     fileProperties.put(LENGTH, String.valueOf(value));
 180   }
 181
 182   public Integer getMaxInstanceLength()
 183   {
 184     if (fileProperties.get(MAX_LENGTH) == null)
 185     {
 186       return null;
 187     }
 188     return Integer.parseInt(fileProperties.get(MAX_LENGTH));
 189   }
 190
 191   public void setMaxInstanceLength(int value)
 192   {
 193     fileProperties.put(MAX_LENGTH, String.valueOf(value));
 194   }
 195
 196   // gets type of symbol alphabet - "amino", "DNA", "RNA"
 197   public String getAlphabetType()
 198   {
 199     return fileProperties.get(ALPHABET);
 200   }
 201
 202   public void setAlphabetType(String value)
 203   {
 204     fileProperties.put(ALPHABET, value);
 205   }
 206
 207   // not sure whether to implement this with Date object
 208   public String getDate()
 209   {
 210     return fileProperties.get(DATE);
 211   }
 212
 213   public void setDate(String value)
 214   {
 215     fileProperties.put(DATE, value);
 216   }
 217
 218   // not sure whether to implement this
 219   public String getCommandLineLog()
 220   {
 221     return fileProperties.get(COMMAND_LOG);
 222   }
 223
 224   public void setCommandLineLog(String value)
 225   {
 226     fileProperties.put(COMMAND_LOG, value);
 227   }
 228
 229   // gets the number of sequences that the HMM was trained on
 230   public Integer getNumberOfSequences()
 231   {
 232     if (fileProperties.get(NUMBER_OF_SEQUENCES) == null)
 233     {
 234       return null;
 235     }
 236     return Integer.parseInt(fileProperties.get(NUMBER_OF_SEQUENCES));
 237   }
 238
 239   public void setNumberOfSequences(int value)
 240   {
 241     fileProperties.put(NUMBER_OF_SEQUENCES, String.valueOf(value));
 242   }
 243
 244   // gets the effective number determined during sequence weighting
 245   public Double getEffectiveNumberOfSequences()
 246   {
 247     if (fileProperties.get(LENGTH) == null)
 248     {
 249       return null;
 250     }
 251     return Double.parseDouble(fileProperties.get(EFF_NUMBER_OF_SEQUENCES));
 252   }
 253
 254   public void setEffectiveNumberOfSequences(double value)
 255   {
 256     fileProperties.put(EFF_NUMBER_OF_SEQUENCES, String.valueOf(value));
 257   }
 258
 259   public Long getCheckSum()
 260   {
 261     if (fileProperties.get(LENGTH) == null)
 262     {
 263       return null;
 264     }
 265     return Long.parseLong(fileProperties.get(CHECK_SUM));
 266   }
 267
 268   public void setCheckSum(long value)
 269   {
 270     fileProperties.put(CHECK_SUM, String.valueOf(value));
 271   }
 272
 273   public List<HMMNode> getNodes()
 274   {
 275     return nodes;
 276   }
 277
 278   public void setNodes(List<HMMNode> nodes)
 279   {
 280     this.nodes = nodes;
 281   }
 282
 283   /**
 284    * get match emission probability for a given symbol at a column in the
 285    * alignment
 286    *
 287    * @param alignColumn
 288    * @param symbol
 289    * @return
 290    *
 291    */
 292   public Double getMatchEmissionProbability(int alignColumn, char symbol)
 293   {
 294     int symbolIndex;
 295     int nodeIndex;
 296     Double probability;
 297     if (!symbolIndexLookup.containsKey(symbol))
 298     {
 299       return 0d;
 300     }
 301     symbolIndex = symbolIndexLookup.get(symbol);
 302     if (nodeLookup.containsKey(alignColumn + 1))
 303     {
 304       nodeIndex = nodeLookup.get(alignColumn + 1);
 305       probability = getNode(nodeIndex).getMatchEmissions().get(symbolIndex);
 306       return probability;
 307     }
 308     else
 309     {
 310       return 0d;
 311     }
 312
 313   }
 314
 315   /**
 316    * get insert emission probability for a given symbol at a column in the
 317    * alignment
 318    *
 319    * @param alignColumn
 320    * @param symbol
 321    * @return
 322    */
 323   public Double getInsertEmissionProbability(int alignColumn, char symbol)
 324   {
 325     int symbolIndex;
 326     int nodeIndex;
 327     Double probability;
 328     if (!symbolIndexLookup.containsKey(symbol))
 329     {
 330       return 0d;
 331     }
 332     symbolIndex = symbolIndexLookup.get(symbol);
 333     if (nodeLookup.containsKey(alignColumn + 1))
 334     {
 335       nodeIndex = nodeLookup.get(alignColumn + 1);
 336       probability = getNode(nodeIndex).getInsertEmissions()
 337               .get(symbolIndex);
 338       return probability;
 339     }
 340     else
 341     {
 342       return 0d;
 343     }
 344
 345   }
 346
 347   /**
 348    * get state transition probability for a given transition type at a column in
 349    * the alignment
 350    *
 351    * @param alignColumn
 352    * @param transition
 353    * @return
 354    */
 355   public Double getStateTransitionProbability(int alignColumn,
 356           String transition)
 357   {
 358     int transitionIndex;
 359     int nodeIndex;
 360     Double probability;
 361     transitionIndex = getTransitionType(transition);
 362     if (nodeLookup.containsKey(alignColumn + 1))
 363     {
 364       nodeIndex = nodeLookup.get(alignColumn + 1);
 365       probability = getNode(nodeIndex).getStateTransitions()
 366             .get(transitionIndex);
 367       return probability;
 368     }
 369     else
 370     {
 371       return 0d;
 372     }
 373
 374   }
 375
 376   public Integer getNodeAlignmentColumn(int nodeIndex)
 377   {
 378     Integer value = nodes.get(nodeIndex).getAlignmentColumn();
 379    return value;
 380   }
 381
 382   public char getConsensusResidue(int nodeIndex)
 383   {
 384    char value = nodes.get(nodeIndex).getConsensusResidue();
 385    return value;
 386   }
 387
 388   public char getConsensus(int columnIndex)
 389   {
 390     char value;
 391     Integer index = findNodeIndex(columnIndex + 1);
 392     if (index == null)
 393     {
 394       return '-';
 395     }
 396     value = getNodes().get(index).getConsensusResidue();
 397     return value;
 398   }
 399
 400   public char getReferenceAnnotation(int nodeIndex)
 401   {
 402    char value = nodes.get(nodeIndex).getReferenceAnnotation();
 403    return value;
 404   }
 405
 406   public char getMaskedValue(int nodeIndex)
 407   {
 408    char value = nodes.get(nodeIndex).getMaskValue();
 409    return value;
 410   }
 411
 412   public char getConsensusStructure(int nodeIndex)
 413   {
 414    char value = nodes.get(nodeIndex).getConsensusStructure();
 415    return value;
 416   }
 417
 418   /**
 419    * returns the average match emission for a given symbol
 420    * @param symbolIndex
 421    * index of symbol
 422    * @return
 423    * average negative log propbability of a match emission of the given symbol
 424    */
 425   public double getAverageMatchEmission(int symbolIndex)
 426   {
 427     double value = nodes.get(0).getMatchEmissions().get(symbolIndex);
 428     return value;
 429   }
 430
 431   public int getNumberOfSymbols()
 432   {
 433     return numberOfSymbols;
 434   }
 435
 436   public void setNumberOfSymbols(int numberOfSymbols)
 437   {
 438     this.numberOfSymbols = numberOfSymbols;
 439   }
 440
 441
 442
 443   /**
 444    * fills symbol array and also finds numberOfSymbols
 445    *
 446    * @param parser
 447    *          scanner scanning symbol line in file
 448    */
 449   public void fillSymbols(Scanner parser)
 450   {
 451     int i = 0;
 452     while (parser.hasNext())
 453     {
 454       String strSymbol = parser.next();
 455       char[] symbol = strSymbol.toCharArray();
 456       symbols.add(symbol[0]);
 457       symbolIndexLookup.put(symbol[0], i);
 458       i++;
 459     }
 460     numberOfSymbols = symbols.size();
 461   }
 462
 463   /**
 464    * adds file property
 465    *
 466    * @param key
 467    * @param value
 468    */
 469   public void addFileProperty(String key, String value)
 470   {
 471     fileProperties.put(key, value);
 472   }
 473
 474   public boolean referenceAnnotationIsActive()
 475   {
 476     String status;
 477     status = fileProperties.get(REFERENCE_ANNOTATION);
 478     if (status == null)
 479     {
 480       return false;
 481     }
 482     switch (status)
 483     {
 484     case YES:
 485       return true;
 486     case NO:
 487       return false;
 488     default:
 489       return false;
 490     }
 491
 492   }
 493
 494   public boolean maskValueIsActive()
 495   {
 496     String status;
 497     status = fileProperties.get(MASKED_VALUE);
 498     if (status == null)
 499     {
 500       return false;
 501     }
 502     switch (status)
 503     {
 504     case YES:
 505       return true;
 506     case NO:
 507       return false;
 508     default:
 509       return false;
 510     }
 511
 512   }
 513
 514   public boolean consensusResidueIsActive()
 515   {
 516     String status;
 517     status = fileProperties.get(CONSENSUS_RESIDUE);
 518     if (status == null)
 519     {
 520       return false;
 521     }
 522     switch (status)
 523     {
 524     case YES:
 525       return true;
 526     case NO:
 527       return false;
 528     default:
 529       return false;
 530     }
 531
 532   }
 533
 534   public boolean consensusStructureIsActive()
 535   {
 536     String status;
 537     status = fileProperties.get(CONSENSUS_STRUCTURE);
 538     if (status == null)
 539     {
 540       return false;
 541     }
 542     switch (status)
 543     {
 544     case YES:
 545       return true;
 546     case NO:
 547       return false;
 548     default:
 549       return false;
 550     }
 551
 552   }
 553
 554   public boolean mapIsActive()
 555   {
 556     String status;
 557     status = fileProperties.get(MAP);
 558     if (status == null)
 559     {
 560       return false;
 561     }
 562     switch (status)
 563     {
 564     case YES:
 565       return true;
 566     case NO:
 567       return false;
 568     default:
 569       return false;
 570     }
 571
 572   }
 573
 574   public void setAlignmentColumn(int nodeIndex, int column)
 575   {
 576     nodes.get(nodeIndex).setAlignmentColumn(column);
 577   }
 578
 579   public void setReferenceAnnotation(int nodeIndex, char value)
 580   {
 581     nodes.get(nodeIndex).setReferenceAnnotation(value);
 582   }
 583
 584   public void setConsensusResidue(int nodeIndex, char value)
 585   {
 586     nodes.get(nodeIndex).setConsensusResidue(value);
 587   }
 588
 589   public void setConsensusStructure(int nodeIndex, char value)
 590   {
 591     nodes.get(nodeIndex).setConsensusStructure(value);
 592   }
 593
 594   public void setMaskValue(int nodeIndex, char value)
 595   {
 596     nodes.get(nodeIndex).setMaskValue(value);
 597   }
 598
 599   public String getGatheringThreshold()
 600   {
 601     String value;
 602     value = fileProperties.get("GA");
 603     return value;
 604   }
 605
 606   public String getNoiseCutoff()
 607   {
 608     String value;
 609     value = fileProperties.get("NC");
 610     return value;
 611   }
 612
 613   public String getTrustedCutoff()
 614   {
 615     String value;
 616     value = fileProperties.get("TC");
 617     return value;
 618   }
 619
 620   public String getViterbi()
 621   {
 622     String value;
 623     value = fileProperties.get(VITERBI);
 624     return value;
 625   }
 626
 627   public String getMSV()
 628   {
 629     String value;
 630     value = fileProperties.get(MSV);
 631     return value;
 632   }
 633
 634   public String getForward()
 635   {
 636     String value;
 637     value = fileProperties.get(FORWARD);
 638     return value;
 639   }
 640
 641   public void setMAPStatus(boolean status)
 642   {
 643     if (status == true)
 644     {
 645       fileProperties.put(MAP, YES);
 646     }
 647     else
 648     {
 649       fileProperties.put(MAP, NO);
 650     }
 651   }
 652
 653   public void setReferenceAnnotationStatus(boolean status)
 654   {
 655     if (status == true)
 656     {
 657       fileProperties.put(REFERENCE_ANNOTATION, YES);
 658     }
 659     else
 660     {
 661       fileProperties.put(REFERENCE_ANNOTATION, NO);
 662     }
 663   }
 664
 665   public void setMaskedValueStatus(boolean status)
 666   {
 667     if (status == true)
 668     {
 669       fileProperties.put(MASKED_VALUE, YES);
 670     }
 671     else
 672     {
 673       fileProperties.put(MASKED_VALUE, NO);
 674     }
 675   }
 676
 677   public void setConsensusResidueStatus(boolean status)
 678   {
 679     if (status == true)
 680     {
 681       fileProperties.put(CONSENSUS_RESIDUE, YES);
 682     }
 683     else
 684     {
 685       fileProperties.put(CONSENSUS_RESIDUE, NO);
 686     }
 687   }
 688
 689   public void setConsensusStructureStatus(boolean status)
 690   {
 691     if (status == true)
 692     {
 693       fileProperties.put(CONSENSUS_STRUCTURE, YES);
 694     }
 695     else
 696     {
 697       fileProperties.put(CONSENSUS_STRUCTURE, NO);
 698     }
 699   }
 700
 701   /**
 702    *
 703    * @param transition
 704    *          type of transition occuring
 705    * @return index value representing position along stateTransition array.
 706    */
 707   public Integer getTransitionType(String transition)
 708   {
 709     Integer index;
 710     switch (transition)
 711     {
 712     case "mm":
 713       index = 0;
 714       break;
 715     case "mi":
 716       index = 1;
 717       break;
 718     case "md":
 719       index = 2;
 720       break;
 721     case "im":
 722       index = 3;
 723       break;
 724     case "ii":
 725       index = 4;
 726       break;
 727     case "dm":
 728       index = 5;
 729       break;
 730     case "dd":
 731       index = 6;
 732       break;
 733     default:
 734       index = null;
 735     }
 736     return index;
 737   }
 738
 739   /**
 740    * find the index of the node in a hidden Markov model based on the column in
 741    * the alignment
 742    *
 743    * @param alignmentColumn
 744    */
 745
 746   public Integer findNodeIndex(int alignmentColumn)
 747   {
 748     Integer index;
 749     index = nodeLookup.get(alignmentColumn);
 750     return index;
 751   }
 752
 753   public static String findStringFromBoolean(boolean value)
 754   {
 755     if (value)
 756     {
 757       return YES;
 758     }
 759     else
 760     {
 761       return NO;
 762     }
 763   }
 764
 765   /**
 766    * creates the HMM annotation
 767    *
 768    * @return
 769    */
 770   public AlignmentAnnotation createAnnotation(int length)
 771   {
 772     Annotation[] annotations = new Annotation[length];
 773     float max = 0f;
 774     for (int i = 0; i < length; i++)
 775     {
 776       Float content = getInformationContent(i);
 777       if (content > max)
 778       {
 779         max = content;
 780       }
 781       String description = String.format("%.3f", content);
 782       description += " bits";
 783       annotations[i] = new Annotation(null, description, ' ', content);
 784
 785     }
 786     AlignmentAnnotation annotation = new AlignmentAnnotation(
 787             "Information Content",
 788             "The information content of each column, measured in bits",
 789             annotations,
 790             0f, max, AlignmentAnnotation.BAR_GRAPH);
 791     return annotation;
 792   }
 793
 794   public float getInformationContent(int column)
 795   {
 796     float informationContent = 0f;
 797
 798     for (char symbol : symbols)
 799     {
 800       float freq = 0f;
 801       if (symbols.size() == 20)
 802       {
 803         freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
 804       }
 805       if (symbols.size() == 4)
 806       {
 807         freq = ResidueProperties.nucleotideBackgroundFrequencies
 808                 .get(symbol);
 809       }
 810       Double hmmProb = getMatchEmissionProbability(column, symbol);
 811       float prob = hmmProb.floatValue();
 812       informationContent += prob * (Math.log(prob / freq) / Math.log(2));
 813
 814     }
 815
 816     return informationContent;
 817   }
 818 }
 819