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   public static final int MATCHTOMATCH = 0;
 106
 107   public static final int MATCHTOINSERT = 1;
 108
 109   public static final int MATCHTODELETE = 2;
 110
 111   public static final int INSERTTOMATCH = 3;
 112
 113   public static final int INSERTTOINSERT = 4;
 114
 115   public static final int DELETETOMATCH = 5;
 116
 117   public static final int DELETETODELETE = 6;
 118
 119   public Map<Integer, Integer> getNodeLookup()
 120   {
 121     return nodeLookup;
 122   }
 123
 124   public void setNodeLookup(Map<Integer, Integer> nodeLookup)
 125   {
 126     this.nodeLookup = nodeLookup;
 127   }
 128
 129   public List<Character> getSymbols()
 130   {
 131     return symbols;
 132   }
 133
 134   public Map<String, String> getFileProperties()
 135   {
 136     return fileProperties;
 137   }
 138
 139   public HMMNode getNode(int nodeIndex)
 140   {
 141     return getNodes().get(nodeIndex);
 142   }
 143
 144   public void setSymbols(List<Character> symbolsL)
 145   {
 146     this.symbols = symbolsL;
 147   }
 148
 149   public String getName()
 150   {
 151     return fileProperties.get(NAME);
 152   }
 153   public String getAccessionNumber()
 154   {
 155     return fileProperties.get(ACCESSION_NUMBER);
 156   }
 157
 158   public void setAccessionNumber(String value)
 159   {
 160     fileProperties.put(ACCESSION_NUMBER, value);
 161   }
 162
 163   public String getDescription()
 164   {
 165     return fileProperties.get(DESCRIPTION);
 166   }
 167
 168   public void setDescription(String value)
 169   {
 170     fileProperties.put(DESCRIPTION, value);
 171   }
 172
 173   public Integer getLength()
 174   {
 175     if (fileProperties.get(LENGTH) == null)
 176     {
 177       return null;
 178     }
 179     return Integer.parseInt(fileProperties.get(LENGTH));
 180   }
 181
 182   public void setLength(int value)
 183   {
 184     fileProperties.put(LENGTH, String.valueOf(value));
 185   }
 186
 187   public Integer getMaxInstanceLength()
 188   {
 189     if (fileProperties.get(MAX_LENGTH) == null)
 190     {
 191       return null;
 192     }
 193     return Integer.parseInt(fileProperties.get(MAX_LENGTH));
 194   }
 195
 196   public void setMaxInstanceLength(int value)
 197   {
 198     fileProperties.put(MAX_LENGTH, String.valueOf(value));
 199   }
 200
 201   // gets type of symbol alphabet - "amino", "DNA", "RNA"
 202   public String getAlphabetType()
 203   {
 204     return fileProperties.get(ALPHABET);
 205   }
 206
 207   public void setAlphabetType(String value)
 208   {
 209     fileProperties.put(ALPHABET, value);
 210   }
 211
 212   // not sure whether to implement this with Date object
 213   public String getDate()
 214   {
 215     return fileProperties.get(DATE);
 216   }
 217
 218   public void setDate(String value)
 219   {
 220     fileProperties.put(DATE, value);
 221   }
 222
 223   // not sure whether to implement this
 224   public String getCommandLineLog()
 225   {
 226     return fileProperties.get(COMMAND_LOG);
 227   }
 228
 229   public void setCommandLineLog(String value)
 230   {
 231     fileProperties.put(COMMAND_LOG, value);
 232   }
 233
 234   // gets the number of sequences that the HMM was trained on
 235   public Integer getNumberOfSequences()
 236   {
 237     if (fileProperties.get(NUMBER_OF_SEQUENCES) == null)
 238     {
 239       return null;
 240     }
 241     return Integer.parseInt(fileProperties.get(NUMBER_OF_SEQUENCES));
 242   }
 243
 244   public void setNumberOfSequences(int value)
 245   {
 246     fileProperties.put(NUMBER_OF_SEQUENCES, String.valueOf(value));
 247   }
 248
 249   // gets the effective number determined during sequence weighting
 250   public Double getEffectiveNumberOfSequences()
 251   {
 252     if (fileProperties.get(LENGTH) == null)
 253     {
 254       return null;
 255     }
 256     return Double.parseDouble(fileProperties.get(EFF_NUMBER_OF_SEQUENCES));
 257   }
 258
 259   public void setEffectiveNumberOfSequences(double value)
 260   {
 261     fileProperties.put(EFF_NUMBER_OF_SEQUENCES, String.valueOf(value));
 262   }
 263
 264   public Long getCheckSum()
 265   {
 266     if (fileProperties.get(LENGTH) == null)
 267     {
 268       return null;
 269     }
 270     return Long.parseLong(fileProperties.get(CHECK_SUM));
 271   }
 272
 273   public void setCheckSum(long value)
 274   {
 275     fileProperties.put(CHECK_SUM, String.valueOf(value));
 276   }
 277
 278   public List<HMMNode> getNodes()
 279   {
 280     return nodes;
 281   }
 282
 283   public void setNodes(List<HMMNode> nodes)
 284   {
 285     this.nodes = nodes;
 286   }
 287
 288   /**
 289    * get match emission probability for a given symbol at a column in the
 290    * alignment
 291    *
 292    * @param alignColumn
 293    * @param symbol
 294    * @return
 295    *
 296    */
 297   public Double getMatchEmissionProbability(int alignColumn, char symbol)
 298   {
 299     int symbolIndex;
 300     int nodeIndex;
 301     Double probability;
 302     if (!symbolIndexLookup.containsKey(symbol))
 303     {
 304       return 0d;
 305     }
 306     symbolIndex = symbolIndexLookup.get(symbol);
 307     if (nodeLookup.containsKey(alignColumn + 1))
 308     {
 309       nodeIndex = nodeLookup.get(alignColumn + 1);
 310       probability = getNode(nodeIndex).getMatchEmissions().get(symbolIndex);
 311       return probability;
 312     }
 313     else
 314     {
 315       return 0d;
 316     }
 317
 318   }
 319
 320   /**
 321    * get insert emission probability for a given symbol at a column in the
 322    * alignment
 323    *
 324    * @param alignColumn
 325    * @param symbol
 326    * @return
 327    */
 328   public Double getInsertEmissionProbability(int alignColumn, char symbol)
 329   {
 330     int symbolIndex;
 331     int nodeIndex;
 332     Double probability;
 333     if (!symbolIndexLookup.containsKey(symbol))
 334     {
 335       return 0d;
 336     }
 337     symbolIndex = symbolIndexLookup.get(symbol);
 338     if (nodeLookup.containsKey(alignColumn + 1))
 339     {
 340       nodeIndex = nodeLookup.get(alignColumn + 1);
 341       probability = getNode(nodeIndex).getInsertEmissions()
 342               .get(symbolIndex);
 343       return probability;
 344     }
 345     else
 346     {
 347       return 0d;
 348     }
 349
 350   }
 351
 352   /**
 353    * get state transition probability for a given transition type at a column in
 354    * the alignment
 355    *
 356    * @param alignColumn
 357    * @param transition
 358    * @return
 359    */
 360   public Double getStateTransitionProbability(int alignColumn,
 361           int transition)
 362   {
 363     int transitionIndex;
 364     int nodeIndex;
 365     Double probability;
 366     if (nodeLookup.containsKey(alignColumn + 1))
 367     {
 368       nodeIndex = nodeLookup.get(alignColumn + 1);
 369       probability = getNode(nodeIndex).getStateTransitions()
 370               .get(transition);
 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    * find the index of the node in a hidden Markov model based on the column in
 707    * the alignment
 708    *
 709    * @param alignmentColumn
 710    */
 711
 712   public Integer findNodeIndex(int alignmentColumn)
 713   {
 714     Integer index;
 715     index = nodeLookup.get(alignmentColumn);
 716     return index;
 717   }
 718
 719   public static String findStringFromBoolean(boolean value)
 720   {
 721     if (value)
 722     {
 723       return YES;
 724     }
 725     else
 726     {
 727       return NO;
 728     }
 729   }
 730
 731   /**
 732    * @return
 733    */
 734   public AlignmentAnnotation createAnnotation(int length)
 735   {
 736     Annotation[] annotations = new Annotation[length];
 737     float max = 0f;
 738     for (int alignPos = 0; alignPos < length; alignPos++)
 739     {
 740       Float content = getInformationContent(alignPos);
 741       if (content > max)
 742       {
 743         max = content;
 744       }
 745
 746       Character cons;
 747       cons = getConsensusAtAlignColumn(alignPos);
 748       cons = Character.toUpperCase(cons);
 749
 750       String description = String.format("%.3f", content);
 751       description += " bits";
 752       annotations[alignPos] = new Annotation(cons.toString(), description,
 753               ' ',
 754               content);
 755
 756     }
 757     AlignmentAnnotation annotation = new AlignmentAnnotation(
 758             "Information Content",
 759             "The information content of each column, measured in bits",
 760             annotations,
 761             0f, max, AlignmentAnnotation.BAR_GRAPH);
 762     return annotation;
 763   }
 764
 765   public float getInformationContent(int column)
 766   {
 767     float informationContent = 0f;
 768
 769     for (char symbol : symbols)
 770     {
 771       float freq = 0f;
 772       if (symbols.size() == 20)
 773       {
 774         freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
 775       }
 776       if (symbols.size() == 4)
 777       {
 778         freq = ResidueProperties.nucleotideBackgroundFrequencies
 779                 .get(symbol);
 780       }
 781       Double hmmProb = getMatchEmissionProbability(column, symbol);
 782       float prob = hmmProb.floatValue();
 783       informationContent += prob * (Math.log(prob / freq) / Math.log(2));
 784
 785     }
 786
 787     return informationContent;
 788   }
 789
 790 }
 791