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