test HMM annotaTION

[jalview.git] / src / jalview / analysis / AAFrequency.java

diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java
index 514b300..3cc57bf 100755 (executable)
--- a/src/jalview/analysis/AAFrequency.java
+++ b/src/jalview/analysis/AAFrequency.java
@@ -24,6 +24,7 @@ import jalview.datamodel.AlignedCodonFrame;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.AlignmentI;
 import jalview.datamodel.Annotation;
+import jalview.datamodel.HiddenMarkovModel;
 import jalview.datamodel.Profile;
 import jalview.datamodel.ProfileI;
 import jalview.datamodel.Profiles;
@@ -32,6 +33,7 @@ import jalview.datamodel.ResidueCount;
 import jalview.datamodel.ResidueCount.SymbolCounts;
 import jalview.datamodel.SequenceI;
 import jalview.ext.android.SparseIntArray;
+import jalview.schemes.ResidueProperties;
 import jalview.util.Comparison;
 import jalview.util.Format;
 import jalview.util.MappingUtils;
@@ -55,6 +57,10 @@ public class AAFrequency
 {
   public static final String PROFILE = "P";
 
+  private static final String AMINO = "amino";
+
+  private static final String DNA = "DNA";
+
   /*
    * Quick look-up of String value of char 'A' to 'Z'
    */
@@ -719,4 +725,89 @@ public class AAFrequency
     }
     return scale;
   }
+
+  /**
+   * produces a HMM profile for a column in an alignment
+   * 
+   * @param aa
+   *          Alignment annotation for which the profile is being calculated
+   * @param column
+   *          column in the alignment the profile is being made for
+   * @param removeBelowBackground
+   *          boolean, indicating whether to ignore residues with probabilities
+   *          less than their background frequencies
+   * @return
+   */
+  public static int[] getHMMProfileFor(AlignmentAnnotation aa, int column,
+          boolean removeBelowBackground)
+  {
+
+    HiddenMarkovModel hmm;
+    hmm = aa.getHMM();
+    if (hmm != null)
+    {
+      String alph = hmm.getAlphabetType();
+      int size = hmm.getNumberOfSymbols();
+      char symbols[] = new char[size];
+      int values[] = new int[size];
+      List<Character> charList = hmm.getSymbols();
+      Integer totalCount = 0;
+
+      for (int i = 0; i < size; i++)
+      {
+        char symbol = charList.get(i);
+        symbols[i] = symbol;
+        Double value;
+
+        value = hmm.getMatchEmissionProbability(column, symbol);
+        double freq;
+
+        if (alph == AMINO && removeBelowBackground)
+        {
+          freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
+          if (value < freq)
+          {
+            value = 0d;
+          }
+        }
+        else if (alph == DNA && removeBelowBackground)
+        {
+          freq = ResidueProperties.nucleotideBackgroundFrequencies
+                  .get(symbol);
+          if (value < freq)
+          {
+            value = 0d;
+          }
+        }
+        value = value * 10000;
+        values[i] = value.intValue();
+        totalCount += value.intValue();
+      }
+
+      QuickSort.sort(values, symbols);
+
+      int[] profile = new int[3 + size * 2];
+
+      profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
+      profile[1] = size;
+      profile[2] = totalCount / 100;
+
+      if (totalCount != 0)
+      {
+        int arrayPos = 3;
+        for (int k = size - 1; k >= 0; k--)
+        {
+          Double percentage;
+          Integer value = values[k];
+          percentage = (value.doubleValue() / totalCount.doubleValue())
+                  * 100d;
+          profile[arrayPos] = symbols[k];
+          profile[arrayPos + 1] = percentage.intValue();
+          arrayPos += 2;
+        }
+      }
+      return profile;
+    }
+    return null;
+  }
 }