Merge branch 'develop' into features/mchmmer
[jalview.git] / src / jalview / analysis / AAFrequency.java
index 5ea892b..85a9bb0 100755 (executable)
@@ -157,14 +157,13 @@ public class AAFrequency
       {
         if (sequences[row] == null)
         {
-          System.err
-                  .println("WARNING: Consensus skipping null sequence - possible race condition.");
+          System.err.println(
+                  "WARNING: Consensus skipping null sequence - possible race condition.");
           continue;
         }
-        char[] seq = sequences[row].getSequence();
-        if (seq.length > column)
+        if (sequences[row].getLength() > column)
         {
-          char c = seq[column];
+          char c = sequences[row].getCharAt(column);
           residueCounts.add(c);
           if (Comparison.isNucleotide(c))
           {
@@ -202,20 +201,39 @@ public class AAFrequency
     // System.out.println(elapsed);
   }
 
-  public static ProfilesI calculateInformation(final HiddenMarkovModel hmm,
+  /**
+   * Returns the full set of profiles for a hidden Markov model. The underlying
+   * data is the raw probabilities of a residue being emitted at each node,
+   * however the profiles returned by this function contain the percentage
+   * chance of a residue emission.
+   * 
+   * @param hmm
+   * @param width
+   *          The width of the Profile array (Profiles) to be returned.
+   * @param start
+   *          The alignment column on which the first profile is based.
+   * @param end
+   *          The alignment column on which the last profile is based.
+   * @param saveFullProfile
+   *          Flag for saving the counts for each profile
+   * @param removeBelowBackground
+   *          Flag for removing any characters with a match emission probability
+   *          less than its background frequency
+   * @return
+   */
+  public static ProfilesI calculateHMMProfiles(final HiddenMarkovModel hmm,
           int width, int start, int end, boolean saveFullProfile,
-          boolean removeBelowBackground)
+          boolean removeBelowBackground, boolean infoLetterHeight)
   {
     ProfileI[] result = new ProfileI[width];
     int symbolCount = hmm.getNumberOfSymbols();
-    String alph = hmm.getAlphabetType();
     for (int column = start; column < end; column++)
     {
       ResidueCount counts = new ResidueCount();
       for (char symbol : hmm.getSymbols())
       {
-        int value = getAnalogueCount(hmm, column, removeBelowBackground,
-                alph, symbol);
+        int value = getAnalogueCount(hmm, column, symbol,
+                removeBelowBackground, infoLetterHeight);
         counts.put(symbol, value);
       }
       int maxCount = counts.getModalCount();
@@ -353,10 +371,9 @@ public class AAFrequency
    * @param nseq
    *          number of sequences
    */
-  public static void completeInformation(AlignmentAnnotation information,
-          ProfilesI profiles, int startCol, int endCol,
-          boolean ignoreBelowBackground,
-          boolean showSequenceLogo, long nseq)
+  public static float completeInformation(AlignmentAnnotation information,
+          ProfilesI profiles, int startCol, int endCol, long nseq,
+          Float currentMax)
   {
     // long now = System.currentTimeMillis();
     if (information == null || information.annotations == null
@@ -366,7 +383,7 @@ public class AAFrequency
        * called with a bad alignment annotation row 
        * wait for it to be initialised properly
        */
-      return;
+      return 0;
     }
 
     Float max = 0f;
@@ -381,7 +398,7 @@ public class AAFrequency
          * shorter than alignment width
          */
         information.annotations[i] = null;
-        return;
+        return 0;
       }
 
       HiddenMarkovModel hmm;
@@ -398,13 +415,21 @@ public class AAFrequency
       }
 
       String description = value + " bits";
-      
-      information.annotations[i] = new Annotation(" ", description,
-              ' ', value);
+      information.annotations[i] = new Annotation(
+              Character.toString(Character
+                      .toUpperCase(hmm.getConsensusAtAlignColumn(i))),
+              description, ' ', value);
+    }
+    if (max > currentMax)
+    {
+      information.graphMax = max;
+      return max;
+    }
+    else
+    {
+      information.graphMax = currentMax;
+      return currentMax;
     }
-    information.graphMax = max;
-    // long elapsed = System.currentTimeMillis() - now;
-    // System.out.println(-elapsed);
   }
 
   /**
@@ -434,7 +459,7 @@ public class AAFrequency
     // always set ranges again
     gaprow.graphMax = nseq;
     gaprow.graphMin = 0;
-    double scale = 0.8/nseq;
+    double scale = 0.8 / nseq;
     for (int i = startCol; i < endCol; i++)
     {
       ProfileI profile = profiles.get(i);
@@ -452,9 +477,9 @@ public class AAFrequency
 
       String description = "" + gapped;
 
-      gaprow.annotations[i] = new Annotation("", description,
-              '\0', gapped, jalview.util.ColorUtils.bleachColour(
-                      Color.DARK_GRAY, (float) scale * gapped));
+      gaprow.annotations[i] = new Annotation("", description, '\0', gapped,
+              jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY,
+                      (float) scale * gapped));
     }
   }
 
@@ -463,7 +488,8 @@ public class AAFrequency
    * <ul>
    * <li>the full profile (percentages of all residues present), if
    * showSequenceLogo is true, or</li>
-   * <li>just the modal (most common) residue(s), if showSequenceLogo is false</li>
+   * <li>just the modal (most common) residue(s), if showSequenceLogo is
+   * false</li>
    * </ul>
    * Percentages are as a fraction of all sequence, or only ungapped sequences
    * if ignoreGaps is true.
@@ -484,8 +510,8 @@ public class AAFrequency
     String description = null;
     if (counts != null && showSequenceLogo)
     {
-      int normaliseBy = ignoreGaps ? profile.getNonGapped() : profile
-              .getHeight();
+      int normaliseBy = ignoreGaps ? profile.getNonGapped()
+              : profile.getHeight();
       description = counts.getTooltip(normaliseBy, dp);
     }
     else
@@ -542,8 +568,8 @@ public class AAFrequency
     QuickSort.sort(values, symbols);
     int nextArrayPos = 2;
     int totalPercentage = 0;
-    final int divisor = ignoreGaps ? profile.getNonGapped() : profile
-            .getHeight();
+    final int divisor = ignoreGaps ? profile.getNonGapped()
+            : profile.getHeight();
 
     /*
      * traverse the arrays in reverse order (highest counts first)
@@ -661,8 +687,8 @@ public class AAFrequency
         {
           continue;
         }
-        List<char[]> codons = MappingUtils
-                .findCodonsFor(seq, col, mappings);
+        List<char[]> codons = MappingUtils.findCodonsFor(seq, col,
+                mappings);
         for (char[] codon : codons)
         {
           int codonEncoded = CodingUtils.encodeCodon(codon);
@@ -742,10 +768,10 @@ public class AAFrequency
 
       int modalCodonEncoded = codons[codons.length - 1];
       int modalCodonCount = sortedCodonCounts[codons.length - 1];
-      String modalCodon = String.valueOf(CodingUtils
-              .decodeCodon(modalCodonEncoded));
-      if (sortedCodonCounts.length > 1
-              && sortedCodonCounts[codons.length - 2] == sortedCodonCounts[codons.length - 1])
+      String modalCodon = String
+              .valueOf(CodingUtils.decodeCodon(modalCodonEncoded));
+      if (sortedCodonCounts.length > 1 && sortedCodonCounts[codons.length
+              - 2] == sortedCodonCounts[codons.length - 1])
       {
         /*
          * two or more codons share the modal count
@@ -804,8 +830,8 @@ public class AAFrequency
           {
             if (samePercent.length() > 0)
             {
-              mouseOver.append(samePercent).append(": ")
-                      .append(lastPercent).append("% ");
+              mouseOver.append(samePercent).append(": ").append(lastPercent)
+                      .append("% ");
             }
             samePercent.setLength(0);
             samePercent.append(codon);
@@ -853,18 +879,8 @@ public class AAFrequency
     for (char symbol : hmm.getSymbols())
     {
       float freq = 0f;
-      if ("amino".equals(hmm.getAlphabetType()))
-      {
-        freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
-      }
-      if ("DNA".equals(hmm.getAlphabetType()))
-      {
-        freq = ResidueProperties.dnaBackgroundFrequencies.get(symbol);
-      }
-      if ("RNA".equals(hmm.getAlphabetType()))
-      {
-        freq = ResidueProperties.rnaBackgroundFrequencies.get(symbol);
-      }
+      freq = ResidueProperties.backgroundFrequencies
+              .get(hmm.getAlphabetType()).get(symbol);
       Double hmmProb = hmm.getMatchEmissionProbability(column, symbol);
       float prob = hmmProb.floatValue();
       informationContent += prob * (Math.log(prob / freq) / Math.log(2));
@@ -887,12 +903,11 @@ public class AAFrequency
    * @return
    */
   public static int[] extractHMMProfile(HiddenMarkovModel hmm, int column,
-          boolean removeBelowBackground)
+          boolean removeBelowBackground, boolean infoHeight)
   {
 
     if (hmm != null)
     {
-      String alph = hmm.getAlphabetType();
       int size = hmm.getNumberOfSymbols();
       char symbols[] = new char[size];
       int values[] = new int[size];
@@ -903,8 +918,8 @@ public class AAFrequency
       {
         char symbol = charList.get(i);
         symbols[i] = symbol;
-        int value = getAnalogueCount(hmm, column, removeBelowBackground,
-                alph, symbol);
+        int value = getAnalogueCount(hmm, column, symbol,
+                removeBelowBackground, infoHeight);
         values[i] = value;
         totalCount += value;
       }
@@ -944,38 +959,38 @@ public class AAFrequency
     return null;
   }
 
-  private static int getAnalogueCount(HiddenMarkovModel hmm, int column,
-          boolean removeBelowBackground, String alph, char symbol)
+  /**
+   * Converts the emission probability of a residue at a column in the alignment
+   * to a 'count' to allow for processing by the annotation renderer.
+   * 
+   * @param hmm
+   * @param column
+   * @param removeBelowBackground
+   *          When true, this method returns 0 for any symbols with a match
+   *          emission probability less than the background frequency.
+   * @param symbol
+   * @return
+   */
+  static int getAnalogueCount(HiddenMarkovModel hmm, int column,
+          char symbol, boolean removeBelowBackground, boolean infoHeight)
   {
     Double value;
 
     value = hmm.getMatchEmissionProbability(column, symbol);
     double freq;
 
-    if (AMINO.equals(alph) && removeBelowBackground)
+    freq = ResidueProperties.backgroundFrequencies
+            .get(hmm.getAlphabetType()).get(symbol);
+    if (value < freq && removeBelowBackground)
     {
-      freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
-    }
-    else if (DNA.equals(alph) && removeBelowBackground)
-    {
-      freq = ResidueProperties.dnaBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
+      return 0;
     }
-    else if (RNA.equals(alph) && removeBelowBackground)
+
+    if (infoHeight)
     {
-      freq = ResidueProperties.rnaBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
+      value = value * (Math.log(value / freq) / Math.log(2));
     }
+
     value = value * 10000;
     return Math.round(value.floatValue());
   }