Merge branch 'develop' into update_212_Dec_merge_with_21125_chamges
[jalview.git] / src / jalview / analysis / AAFrequency.java
index d83f138..892b1b1 100755 (executable)
@@ -50,31 +50,12 @@ import java.util.List;
  * This class is used extensively in calculating alignment colourschemes that
  * depend on the amount of conservation in each alignment column.
  * 
- * @author $author$
- * @version $Revision$
  */
 public class AAFrequency
 {
-  public static final String PROFILE = "P";
-
-  private static final String AMINO = "amino";
-
-  private static final String DNA = "DNA";
-
-  private static final String RNA = "RNA";
+  private static final double LOG2 = Math.log(2);
 
-  /*
-   * Quick look-up of String value of char 'A' to 'Z'
-   */
-  private static final String[] CHARS = new String['Z' - 'A' + 1];
-
-  static
-  {
-    for (char c = 'A'; c <= 'Z'; c++)
-    {
-      CHARS[c - 'A'] = String.valueOf(c);
-    }
-  }
+  public static final String PROFILE = "P";
 
   public static final ProfilesI calculate(List<SequenceI> list, int start,
           int end)
@@ -109,8 +90,6 @@ public class AAFrequency
     }
   }
 
-
-
   /**
    * Calculate the consensus symbol(s) for each column in the given range.
    * 
@@ -157,14 +136,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 +180,38 @@ public class AAFrequency
     // System.out.println(elapsed);
   }
 
-  public static ProfilesI calculateInformation(final HiddenMarkovModel hmm,
-          int width, int start, int end, boolean saveFullProfile,
-          boolean removeBelowBackground)
+  /**
+   * 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 removeBelowBackground
+   *          if true, symbols with a match emission probability less than
+   *          background frequency are ignored
+   * @return
+   */
+  public static ProfilesI calculateHMMProfiles(final HiddenMarkovModel hmm,
+          int width, int start, int end, boolean removeBelowBackground,
+          boolean infoLetterHeight)
   {
     ProfileI[] result = new ProfileI[width];
-    int symbolCount = hmm.getNumberOfSymbols();
-    String alph = hmm.getAlphabetType();
+    char[] symbols = hmm.getSymbols().toCharArray();
+    int symbolCount = symbols.length;
     for (int column = start; column < end; column++)
     {
       ResidueCount counts = new ResidueCount();
-      for (char symbol : hmm.getSymbols())
+      for (char symbol : symbols)
       {
-        int value = getAnalogueCount(hmm, column, removeBelowBackground,
-                alph, symbol);
+        int value = getAnalogueCount(hmm, column, symbol,
+                removeBelowBackground, infoLetterHeight);
         counts.put(symbol, value);
       }
       int maxCount = counts.getModalCount();
@@ -223,11 +219,7 @@ public class AAFrequency
       int gapCount = counts.getGapCount();
       ProfileI profile = new Profile(symbolCount, gapCount, maxCount,
               maxResidue);
-
-      if (saveFullProfile)
-      {
-        profile.setCounts(counts);
-      }
+      profile.setCounts(counts);
 
       result[column] = profile;
     }
@@ -346,71 +338,66 @@ public class AAFrequency
    * @param endCol
    *          end column (exclusive)
    * @param ignoreGaps
-   *          if true, normalise residue percentages 
+   *          if true, normalise residue percentages
    * @param showSequenceLogo
    *          if true include all information symbols, else just show modal
    *          residue
-   * @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 now = System.currentTimeMillis();
-    if (information == null || information.annotations == null
-            || information.annotations.length < endCol)
+    if (information == null || information.annotations == null)
     {
       /*
        * called with a bad alignment annotation row 
        * wait for it to be initialised properly
        */
-      return;
+      return 0;
     }
 
-    Float max = 0f;
+    float max = 0f;
+    SequenceI hmmSeq = information.sequenceRef;
 
-    for (int i = startCol; i < endCol; i++)
+    int seqLength = hmmSeq.getLength();
+    if (information.annotations.length < seqLength)
     {
-      ProfileI profile = profiles.get(i);
-      if (profile == null)
+      return 0;
+    }
+
+    HiddenMarkovModel hmm = hmmSeq.getHMM();
+
+    for (int column = startCol; column < endCol; column++)
+    {
+      if (column >= seqLength)
       {
-        /*
-         * happens if sequences calculated over were 
-         * shorter than alignment width
-         */
-        information.annotations[i] = null;
-        return;
+        // hmm consensus sequence is shorter than the alignment
+        break;
       }
-
-      HiddenMarkovModel hmm;
-      
-      SequenceI hmmSeq = information.sequenceRef;
-      
-      hmm = hmmSeq.getHMM();
       
-      Float value = getInformationContent(i, hmm);
-
-      if (value > max)
+      float value = hmm.getInformationContent(column);
+      boolean isNaN = Float.isNaN(value);
+      if (!isNaN)
       {
-        max = value;
+        max = Math.max(max, value);
       }
 
-      String description = value + " bits";
-      
-      information.annotations[i] = new Annotation(" ", description,
-              ' ', value);
+      String description = isNaN ? null
+              : String.format("%.4f bits", value);
+      information.annotations[column] = new Annotation(
+              Character.toString(
+                      Character.toUpperCase(hmmSeq.getCharAt(column))),
+              description, ' ', value);
     }
+
     information.graphMax = max;
-    // long elapsed = System.currentTimeMillis() - now;
-    // System.out.println(-elapsed);
+    return max;
   }
 
   /**
-   * Derive the gap count annotation row.
+   * Derive the occupancy count annotation
    * 
-   * @param gaprow
+   * @param occupancy
    *          the annotation row to add annotations to
    * @param profiles
    *          the source consensus data
@@ -419,11 +406,11 @@ public class AAFrequency
    * @param endCol
    *          end column (exclusive)
    */
-  public static void completeGapAnnot(AlignmentAnnotation gaprow,
+  public static void completeGapAnnot(AlignmentAnnotation occupancy,
           ProfilesI profiles, int startCol, int endCol, long nseq)
   {
-    if (gaprow == null || gaprow.annotations == null
-            || gaprow.annotations.length < endCol)
+    if (occupancy == null || occupancy.annotations == null
+            || occupancy.annotations.length < endCol)
     {
       /*
        * called with a bad alignment annotation row 
@@ -432,9 +419,9 @@ public class AAFrequency
       return;
     }
     // always set ranges again
-    gaprow.graphMax = nseq;
-    gaprow.graphMin = 0;
-    double scale = 0.8/nseq;
+    occupancy.graphMax = nseq;
+    occupancy.graphMin = 0;
+    double scale = 0.8 / nseq;
     for (int i = startCol; i < endCol; i++)
     {
       ProfileI profile = profiles.get(i);
@@ -444,7 +431,7 @@ public class AAFrequency
          * happens if sequences calculated over were 
          * shorter than alignment width
          */
-        gaprow.annotations[i] = null;
+        occupancy.annotations[i] = null;
         return;
       }
 
@@ -452,9 +439,10 @@ public class AAFrequency
 
       String description = "" + gapped;
 
-      gaprow.annotations[i] = new Annotation("", description,
-              '\0', gapped, jalview.util.ColorUtils.bleachColour(
-                      Color.DARK_GRAY, (float) scale * gapped));
+      occupancy.annotations[i] = new Annotation("", description, '\0',
+              gapped,
+              jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY,
+                      (float) scale * gapped));
     }
   }
 
@@ -463,7 +451,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 +473,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
@@ -516,7 +505,7 @@ public class AAFrequency
    * contains
    * 
    * <pre>
-   *    [profileType, numberOfValues, nonGapCount, charValue1, percentage1, charValue2, percentage2, ...]
+   *    [profileType, numberOfValues, totalPercent, charValue1, percentage1, charValue2, percentage2, ...]
    * in descending order of percentage value
    * </pre>
    * 
@@ -529,7 +518,6 @@ public class AAFrequency
    */
   public static int[] extractProfile(ProfileI profile, boolean ignoreGaps)
   {
-    int[] rtnval = new int[64];
     ResidueCount counts = profile.getCounts();
     if (counts == null)
     {
@@ -540,29 +528,51 @@ public class AAFrequency
     char[] symbols = symbolCounts.symbols;
     int[] values = symbolCounts.values;
     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)
      */
+    int[] result = new int[3 + 2 * symbols.length];
+    int nextArrayPos = 3;
+    int nonZeroCount = 0;
+
     for (int i = symbols.length - 1; i >= 0; i--)
     {
       int theChar = symbols[i];
       int charCount = values[i];
-
-      rtnval[nextArrayPos++] = theChar;
       final int percentage = (charCount * 100) / divisor;
-      rtnval[nextArrayPos++] = percentage;
+      if (percentage == 0)
+      {
+        /*
+         * this count (and any remaining) round down to 0% - discard
+         */
+        break;
+      }
+      nonZeroCount++;
+      result[nextArrayPos++] = theChar;
+      result[nextArrayPos++] = percentage;
       totalPercentage += percentage;
     }
-    rtnval[0] = symbols.length;
-    rtnval[1] = totalPercentage;
-    int[] result = new int[rtnval.length + 1];
+
+    /*
+     * truncate array if any zero values were discarded
+     */
+    if (nonZeroCount < symbols.length)
+    {
+      int[] tmp = new int[3 + 2 * nonZeroCount];
+      System.arraycopy(result, 0, tmp, 0, tmp.length);
+      result = tmp;
+    }
+
+    /*
+     * fill in 'header' values
+     */
     result[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
-    System.arraycopy(rtnval, 0, result, 1, rtnval.length);
+    result[1] = nonZeroCount;
+    result[2] = totalPercentage;
 
     return result;
   }
@@ -573,15 +583,15 @@ public class AAFrequency
    * contains
    * 
    * <pre>
-   *    [profileType, numberOfValues, totalCount, charValue1, percentage1, charValue2, percentage2, ...]
+   *    [profileType, numberOfValues, totalPercentage, charValue1, percentage1, charValue2, percentage2, ...]
    * in descending order of percentage value, where the character values encode codon triplets
    * </pre>
    * 
    * @param hashtable
    * @return
    */
-  public static int[] extractCdnaProfile(Hashtable hashtable,
-          boolean ignoreGaps)
+  public static int[] extractCdnaProfile(
+          Hashtable<String, Object> hashtable, boolean ignoreGaps)
   {
     // this holds #seqs, #ungapped, and then codon count, indexed by encoded
     // codon triplet
@@ -611,9 +621,16 @@ public class AAFrequency
       {
         break; // nothing else of interest here
       }
+      final int percentage = codonCount * 100 / divisor;
+      if (percentage == 0)
+      {
+        /*
+         * this (and any remaining) values rounded down to 0 - discard
+         */
+        break;
+      }
       distinctValuesCount++;
       result[j++] = codons[i];
-      final int percentage = codonCount * 100 / divisor;
       result[j++] = percentage;
       totalPercentage += percentage;
     }
@@ -637,7 +654,7 @@ public class AAFrequency
    *          the consensus data stores to be populated (one per column)
    */
   public static void calculateCdna(AlignmentI alignment,
-          Hashtable[] hconsensus)
+          Hashtable<String, Object>[] hconsensus)
   {
     final char gapCharacter = alignment.getGapCharacter();
     List<AlignedCodonFrame> mappings = alignment.getCodonFrames();
@@ -650,7 +667,7 @@ public class AAFrequency
     for (int col = 0; col < cols; col++)
     {
       // todo would prefer a Java bean for consensus data
-      Hashtable<String, int[]> columnHash = new Hashtable<>();
+      Hashtable<String, Object> columnHash = new Hashtable<>();
       // #seqs, #ungapped seqs, counts indexed by (codon encoded + 1)
       int[] codonCounts = new int[66];
       codonCounts[0] = alignment.getSequences().size();
@@ -661,8 +678,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);
@@ -670,6 +687,7 @@ public class AAFrequency
           {
             codonCounts[codonEncoded + 2]++;
             ungappedCount++;
+            break;
           }
         }
       }
@@ -695,7 +713,8 @@ public class AAFrequency
    */
   public static void completeCdnaConsensus(
           AlignmentAnnotation consensusAnnotation,
-          Hashtable[] consensusData, boolean showProfileLogo, int nseqs)
+          Hashtable<String, Object>[] consensusData,
+          boolean showProfileLogo, int nseqs)
   {
     if (consensusAnnotation == null
             || consensusAnnotation.annotations == null
@@ -710,7 +729,7 @@ public class AAFrequency
     consensusAnnotation.scaleColLabel = true;
     for (int col = 0; col < consensusData.length; col++)
     {
-      Hashtable hci = consensusData[col];
+      Hashtable<String, Object> hci = consensusData[col];
       if (hci == null)
       {
         // gapped protein column?
@@ -742,10 +761,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 +823,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);
@@ -839,136 +858,116 @@ public class AAFrequency
   }
 
   /**
-   * Returns the information content at a specified column.
+   * Returns the sorted HMM profile for the given column of the alignment. The
+   * returned array contains
    * 
-   * @param column
-   *          Index of the column, starting from 0.
-   * @return
-   */
-  public static float getInformationContent(int column,
-          HiddenMarkovModel hmm)
-  {
-    float informationContent = 0f;
-
-    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);
-      }
-      Double hmmProb = hmm.getMatchEmissionProbability(column, symbol);
-      float prob = hmmProb.floatValue();
-      informationContent += prob * (Math.log(prob / freq) / Math.log(2));
-
-    }
-
-    return informationContent;
-  }
-
-  /**
-   * Produces a HMM profile for a column in an alignment
+   * <pre>
+   *    [profileType=0, numberOfValues, 100, charValue1, percentage1, charValue2, percentage2, ...]
+   * in descending order of percentage value
+   * </pre>
    * 
-   * @param aa
-   *          Alignment annotation for which the profile is being calculated.
+   * @param hmm
    * @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.
+   *          if true, ignores residues with probability less than their
+   *          background frequency
+   * @param infoHeight
+   *          if true, uses the log ratio 'information' measure to scale the
+   *          value
    * @return
    */
   public static int[] extractHMMProfile(HiddenMarkovModel hmm, int column,
-          boolean removeBelowBackground)
+          boolean removeBelowBackground, boolean infoHeight)
   {
+    if (hmm == null)
+    {
+      return null;
+    }
+    String alphabet = hmm.getSymbols();
+    int size = alphabet.length();
+    char symbols[] = new char[size];
+    int values[] = new int[size];
+    int totalCount = 0;
 
-    if (hmm != null)
+    for (int i = 0; i < size; i++)
     {
-      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;
-        int value = getAnalogueCount(hmm, column, removeBelowBackground,
-                alph, symbol);
-        values[i] = value;
-        totalCount += value;
-      }
+      char symbol = alphabet.charAt(i);
+      symbols[i] = symbol;
+      int value = getAnalogueCount(hmm, column, symbol,
+              removeBelowBackground, infoHeight);
+      values[i] = value;
+      totalCount += value;
+    }
 
-      QuickSort.sort(values, symbols);
+    /*
+     * sort symbols by increasing emission probability
+     */
+    QuickSort.sort(values, symbols);
 
-      int[] profile = new int[3 + size * 2];
+    int[] profile = new int[3 + size * 2];
 
-      profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
-      profile[1] = size;
-      profile[2] = totalCount / 100;
+    profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
+    profile[1] = size;
+    profile[2] = 100;
 
-      if (totalCount != 0)
+    /*
+     * order symbol/count profile by decreasing emission probability
+     */
+    if (totalCount != 0)
+    {
+      int arrayPos = 3;
+      for (int k = size - 1; k >= 0; k--)
       {
-        int arrayPos = 3;
-        for (int k = size - 1; k >= 0; k--)
+        Float percentage;
+        int value = values[k];
+        if (removeBelowBackground)
+        {
+          percentage = ((float) value) / totalCount * 100f;
+        }
+        else
         {
-          Double percentage;
-          Integer value = values[k];
-          percentage = (value.doubleValue() / totalCount.doubleValue())
-                  * 100d;
-          profile[arrayPos] = symbols[k];
-          profile[arrayPos + 1] = percentage.intValue();
-          arrayPos += 2;
+          percentage = value / 100f;
         }
+        int intPercent = Math.round(percentage);
+        profile[arrayPos] = symbols[k];
+        profile[arrayPos + 1] = intPercent;
+        arrayPos += 2;
       }
-      return profile;
     }
-    return null;
+    return profile;
   }
 
-  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', suitable for rendering as an annotation value
+   * 
+   * @param hmm
+   * @param column
+   * @param symbol
+   * @param removeBelowBackground
+   *          if true, returns 0 for any symbol with a match emission
+   *          probability less than the background frequency
+   * @infoHeight if true, uses the log ratio 'information content' to scale the
+   *             value
+   * @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)
+    double value = hmm.getMatchEmissionProbability(column, symbol);
+    double freq = ResidueProperties.backgroundFrequencies
+            .get(hmm.getAlphabetType()).get(symbol);
+    if (value < freq && removeBelowBackground)
     {
-      freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
+      return 0;
     }
-    else if (DNA.equals(alph) && removeBelowBackground)
-    {
-      freq = ResidueProperties.dnaBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
-    }
-    else if (RNA.equals(alph) && removeBelowBackground)
+
+    if (infoHeight)
     {
-      freq = ResidueProperties.rnaBackgroundFrequencies.get(symbol);
-      if (value < freq)
-      {
-        value = 0d;
-      }
+      value = value * (Math.log(value / freq) / LOG2);
     }
-    value = value * 10000;
-    return value.intValue();
+
+    value = value * 10000d;
+    return Math.round((float) value);
   }
 }