test/jalview/datamodel/HiddenMarkovModelTest.java

   1 package jalview.datamodel;
   2
   3 import static org.testng.Assert.assertEquals;
   4
   5 import jalview.io.DataSourceType;
   6 import jalview.io.FileParse;
   7 import jalview.io.HMMFile;
   8 import jalview.schemes.ResidueProperties;
   9
  10 import java.io.IOException;
  11 import java.net.MalformedURLException;
  12 import java.util.Map;
  13
  14 import org.testng.annotations.BeforeClass;
  15 import org.testng.annotations.Test;
  16
  17 public class HiddenMarkovModelTest {
  18
  19   HiddenMarkovModel hmm;
  20
  21   HiddenMarkovModel alignmentHmm;
  22
  23   @BeforeClass(alwaysRun = true)
  24   public void setUp() throws MalformedURLException, IOException
  25   {
  26     /*
  27      * load hmm model of a Kinase domain to a HiddenMarkovModel
  28      */
  29     HMMFile file = new HMMFile(new FileParse(
  30             "test/jalview/io/test_PKinase_hmm.txt", DataSourceType.FILE));
  31     hmm = file.getHMM();
  32
  33     // used to check if consensus sequence is automatically aligned with alignment
  34     HMMFile alignmentTest = new HMMFile(
  35             new FileParse("test/jalview/io/HMMAlignmentTestHMM.hmm",
  36                     DataSourceType.FILE));
  37     alignmentHmm = alignmentTest.getHMM();
  38   }
  39
  40   @Test(groups = "Functional")
  41   public void testGetMatchEmissionProbabilities()
  42           throws MalformedURLException, IOException
  43   {
  44     /*
  45      * raw value in file is 3.67403
  46      * saved as probability e^-X = 0.05259
  47      */
  48     double mep = hmm.getMatchEmissionProbability(0, 'R');
  49     assertEquals(mep, 0.02537400637, 0.0001d);
  50     assertEquals(mep, Math.pow(Math.E, -3.67403), 0.0001d);
  51
  52     mep = hmm.getMatchEmissionProbability(19, 'W');
  53     assertEquals(mep, 0.00588228492, 0.0001d);
  54     assertEquals(mep, Math.pow(Math.E, -5.13581), 0.0001d);
  55
  56     // column 160 is a gapped region of the model
  57     mep = hmm.getMatchEmissionProbability(160, 'G');
  58     assertEquals(mep, 0D, 0.0001d);
  59
  60     mep = hmm.getMatchEmissionProbability(475, 'A');
  61     assertEquals(mep, 0.04995163708, 0.0001d);
  62     assertEquals(mep, Math.pow(Math.E, -2.99670), 0.0001d);
  63   }
  64
  65   @Test(groups = "Functional")
  66   public void testGetInsertEmissionProbabilities()
  67   {
  68     double iep = hmm.getInsertEmissionProbability(2, 'A');
  69     assertEquals(iep, Math.pow(Math.E, -2.68618), 0.0001d);
  70     // symbol is not case-sensitive
  71     assertEquals(iep, hmm.getInsertEmissionProbability(2, 'a'));
  72
  73     iep = hmm.getInsertEmissionProbability(5, 'T');
  74     assertEquals(iep, Math.pow(Math.E, -2.77519), 0.0001d);
  75
  76     // column 161 is gapped in the hmm
  77     iep = hmm.getInsertEmissionProbability(161, 'K');
  78     assertEquals(iep, 0D, 0.0001d);
  79
  80     iep = hmm.getInsertEmissionProbability(472, 'L');
  81     assertEquals(iep, Math.pow(Math.E, -2.69355), 0.0001d);
  82   }
  83
  84   @Test(groups = "Functional")
  85   public void testGetStateTransitionProbabilities()
  86   {
  87     // * in model file treated as negative infinity
  88     double stp = hmm.getStateTransitionProbability(475,
  89             HiddenMarkovModel.MATCHTODELETE);
  90     assertEquals(stp, Double.NEGATIVE_INFINITY);
  91
  92     // file value is 5.01631, saved as e^-5.01631
  93     stp = hmm.getStateTransitionProbability(8,
  94             HiddenMarkovModel.MATCHTOINSERT);
  95     assertEquals(stp, Math.pow(Math.E, -5.01631), 0.0001D);
  96
  97     stp = hmm.getStateTransitionProbability(36,
  98             HiddenMarkovModel.MATCHTODELETE);
  99     assertEquals(stp, Math.pow(Math.E, -5.73865), 0.0001D);
 100
 101     stp = hmm.getStateTransitionProbability(22,
 102             HiddenMarkovModel.INSERTTOMATCH);
 103     assertEquals(stp, Math.pow(Math.E, -0.61958), 0.0001D);
 104
 105     stp = hmm.getStateTransitionProbability(80,
 106             HiddenMarkovModel.INSERTTOINSERT);
 107     assertEquals(stp, Math.pow(Math.E, -0.77255), 0.0001D);
 108
 109     stp = hmm.getStateTransitionProbability(475,
 110             HiddenMarkovModel.DELETETOMATCH);
 111     assertEquals(stp, 1D, 0.0001D);
 112
 113     stp = hmm.getStateTransitionProbability(218,
 114             HiddenMarkovModel.DELETETODELETE);
 115     assertEquals(stp, Math.pow(Math.E, -0.95510), 0.0001D);
 116   }
 117
 118   @Test(groups = "Functional")
 119   public void testGetConsensusSequence()
 120   {
 121     SequenceI seq = hmm.getConsensusSequence();
 122     String subStr = seq.getSequenceAsString().substring(0, 10);
 123     assertEquals(subStr, "yelleklGsG");
 124     subStr = seq.getSequenceAsString().substring(150, 161);
 125     assertEquals(subStr, "-dllk------");
 126
 127     // test to see if consensus sequence maps to alignment correctly
 128     // see HMMAlignmentTest.sto for corresponding alignment file
 129     SequenceI seq2 = alignmentHmm.getConsensusSequence();
 130     assertEquals(seq2.getCharAt(0), '-');
 131     assertEquals(seq2.getCharAt(7), '-');
 132     assertEquals(seq2.getCharAt(8), 's');
 133
 134   }
 135
 136   /**
 137    * A rather pointless test that reproduces the code implemented and asserts
 138    * the result is the same...
 139    */
 140   @Test(groups = "Functional")
 141   public void testGetInformationContent()
 142   {
 143     /*
 144      * information measure is sum over all symbols of
 145      * emissionProb * log(emissionProb / background) / log(2)
 146      */
 147     Map<Character, Float> uniprotFreqs = ResidueProperties.backgroundFrequencies
 148             .get("amino");
 149     int col = 4;
 150     float expected = 0f;
 151     for (char aa : hmm.getSymbols().toCharArray())
 152     {
 153       double mep = hmm.getMatchEmissionProbability(col, aa);
 154       float background = uniprotFreqs.get(aa);
 155       expected += mep * Math.log(mep / background);
 156     }
 157     expected /= Math.log(2D);
 158
 159     float actual = hmm.getInformationContent(col);
 160     assertEquals(actual, expected, 0.0001d);
 161   }
 162 }