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   @BeforeClass(alwaysRun = true)
  22   public void setUp() throws MalformedURLException, IOException
  23   {
  24     /*
  25      * load hmm model of a Kinase domain to a HiddenMarkovModel
  26      */
  27     HMMFile file = new HMMFile(new FileParse(
  28             "test/jalview/io/test_PKinase_hmm.txt", DataSourceType.FILE));
  29     hmm = file.getHMM();
  30   }
  31
  32   @Test(groups = "Functional")
  33   public void testGetMatchEmissionProbabilities()
  34           throws MalformedURLException, IOException
  35   {
  36     /*
  37      * raw value in file is 3.67403
  38      * saved as probability e^-X = 0.05259
  39      */
  40     double mep = hmm.getMatchEmissionProbability(0, 'R');
  41     assertEquals(mep, 0.02537400637, 0.0001d);
  42     assertEquals(mep, Math.pow(Math.E, -3.67403), 0.0001d);
  43
  44     mep = hmm.getMatchEmissionProbability(19, 'W');
  45     assertEquals(mep, 0.00588228492, 0.0001d);
  46     assertEquals(mep, Math.pow(Math.E, -5.13581), 0.0001d);
  47
  48     // column 160 is a gapped region of the model
  49     mep = hmm.getMatchEmissionProbability(160, 'G');
  50     assertEquals(mep, 0D, 0.0001d);
  51
  52     mep = hmm.getMatchEmissionProbability(475, 'A');
  53     assertEquals(mep, 0.04995163708, 0.0001d);
  54     assertEquals(mep, Math.pow(Math.E, -2.99670), 0.0001d);
  55   }
  56
  57   @Test(groups = "Functional")
  58   public void testGetInsertEmissionProbabilities()
  59   {
  60     double iep = hmm.getInsertEmissionProbability(2, 'A');
  61     assertEquals(iep, Math.pow(Math.E, -2.68618), 0.0001d);
  62
  63     iep = hmm.getInsertEmissionProbability(5, 'T');
  64     assertEquals(iep, Math.pow(Math.E, -2.77519), 0.0001d);
  65
  66     // column 161 is gapped in the hmm
  67     iep = hmm.getInsertEmissionProbability(161, 'K');
  68     assertEquals(iep, 0D, 0.0001d);
  69
  70     iep = hmm.getInsertEmissionProbability(472, 'L');
  71     assertEquals(iep, Math.pow(Math.E, -2.69355), 0.0001d);
  72   }
  73
  74   @Test(groups = "Functional")
  75   public void testGetStateTransitionProbabilities()
  76   {
  77     // * in model file treated as negative infinity
  78     double stp = hmm.getStateTransitionProbability(475,
  79             HiddenMarkovModel.MATCHTODELETE);
  80     assertEquals(stp, Double.NEGATIVE_INFINITY);
  81
  82     // file value is 5.01631, saved as e^-5.01631
  83     stp = hmm.getStateTransitionProbability(8,
  84             HiddenMarkovModel.MATCHTOINSERT);
  85     assertEquals(stp, Math.pow(Math.E, -5.01631), 0.0001D);
  86
  87     stp = hmm.getStateTransitionProbability(36,
  88             HiddenMarkovModel.MATCHTODELETE);
  89     assertEquals(stp, Math.pow(Math.E, -5.73865), 0.0001D);
  90
  91     stp = hmm.getStateTransitionProbability(22,
  92             HiddenMarkovModel.INSERTTOMATCH);
  93     assertEquals(stp, Math.pow(Math.E, -0.61958), 0.0001D);
  94
  95     stp = hmm.getStateTransitionProbability(80,
  96             HiddenMarkovModel.INSERTTOINSERT);
  97     assertEquals(stp, Math.pow(Math.E, -0.77255), 0.0001D);
  98
  99     stp = hmm.getStateTransitionProbability(475,
 100             HiddenMarkovModel.DELETETOMATCH);
 101     assertEquals(stp, 1D, 0.0001D);
 102
 103     stp = hmm.getStateTransitionProbability(218,
 104             HiddenMarkovModel.DELETETODELETE);
 105     assertEquals(stp, Math.pow(Math.E, -0.95510), 0.0001D);
 106   }
 107
 108   @Test(groups = "Functional")
 109   public void testGetConsensusAtAlignColumn()
 110   {
 111     assertEquals(hmm.getConsensusAtAlignColumn(10), 's');
 112     assertEquals(hmm.getConsensusAtAlignColumn(50), 'k');
 113     hmm.setConsensusResidueStatus(false);
 114     assertEquals(hmm.getConsensusAtAlignColumn(100), 'l');
 115     assertEquals(hmm.getConsensusAtAlignColumn(400), 'k');
 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
 128   /**
 129    * A rather pointless test that reproduces the code implemented and asserts
 130    * the result is the same...
 131    */
 132   @Test(groups = "Functional")
 133   public void testGetInformationContent()
 134   {
 135     /*
 136      * information measure is sum over all symbols of
 137      * emissionProb * log(emissionProb / background) / log(2)
 138      */
 139     Map<Character, Float> uniprotFreqs = ResidueProperties.backgroundFrequencies
 140             .get("amino");
 141     int col = 4;
 142     float expected = 0f;
 143     for (char aa : hmm.getSymbols())
 144     {
 145       double mep = hmm.getMatchEmissionProbability(col, aa);
 146       float background = uniprotFreqs.get(aa);
 147       expected += mep * Math.log(mep / background);
 148     }
 149     expected /= Math.log(2D);
 150
 151     float actual = hmm.getInformationContent(col);
 152     assertEquals(actual, expected, 0.0001d);
 153   }
 154 }