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