JAL-1925 update source version in license
[jalview.git] / src / jalview / analysis / StructureFrequency.java
index 4d83019..991fe5e 100644 (file)
@@ -1,25 +1,33 @@
 /*
- * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)
- * Copyright (C) 2010 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle
+ * Jalview - A Sequence Alignment Editor and Viewer (Version 2.9.0b2)
+ * Copyright (C) 2015 The Jalview Authors
  * 
  * This file is part of Jalview.
  * 
  * Jalview is free software: you can redistribute it and/or
  * modify it under the terms of the GNU General Public License 
- * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
- * 
+ * as published by the Free Software Foundation, either version 3
+ * of the License, or (at your option) any later version.
+ *  
  * Jalview is distributed in the hope that it will be useful, but 
  * WITHOUT ANY WARRANTY; without even the implied warranty 
  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
  * PURPOSE.  See the GNU General Public License for more details.
  * 
- * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
+ * You should have received a copy of the GNU General Public License
+ * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
+ * The Jalview Authors are detailed in the 'AUTHORS' file.
  */
 package jalview.analysis;
 
-import java.util.*;
+import jalview.datamodel.AlignmentAnnotation;
+import jalview.datamodel.Annotation;
+import jalview.datamodel.SequenceFeature;
+import jalview.datamodel.SequenceI;
+import jalview.util.Format;
 
-import jalview.datamodel.*;
+import java.util.ArrayList;
+import java.util.Hashtable;
 
 /**
  * Takes in a vector or array of sequences and column start and column end and
@@ -32,6 +40,8 @@ import jalview.datamodel.*;
  */
 public class StructureFrequency
 {
+  public static final int STRUCTURE_PROFILE_LENGTH = 74;
+
   // No need to store 1000s of strings which are not
   // visible to the user.
   public static final String MAXCOUNT = "C";
@@ -46,150 +56,26 @@ public class StructureFrequency
 
   public static final String PAIRPROFILE = "B";
 
-  public static final Hashtable[] calculate(Vector sequences, int start,
-          int end)
-  {
-    return calculate(sequences, start, end, false);
-  }
-
-  public static final Hashtable[] calculate(Vector sequences, int start,
-          int end, boolean profile)
-  {
-    SequenceI[] seqs = new SequenceI[sequences.size()];
-    int width = 0;
-    for (int i = 0; i < sequences.size(); i++)
-    {
-      seqs[i] = (SequenceI) sequences.elementAt(i);
-      if (seqs[i].getLength() > width)
-      {
-        width = seqs[i].getLength();
-      }
-    }
-
-    Hashtable[] reply = new Hashtable[width];
-
-    if (end >= width)
-    {
-      end = width;
-    }
-
-    calculate(seqs, start, end, reply, profile);
-
-    return reply;
-  }
-
-  public static final void calculate(SequenceI[] sequences, int start,
-          int end, Hashtable[] result)
-  {
-    calculate(sequences, start, end, result, false);
-  }
-
-  public static final void calculate(SequenceI[] sequences, int start,
-          int end, Hashtable[] result, boolean profile)
-  {
-    Hashtable residueHash;
-    int maxCount, nongap, i, j, v, jSize = sequences.length;
-    String maxResidue;
-    char c;
-    float percentage;
-
-    int[] values = new int[255];
-
-    char[] seq;
-
-    for (i = start; i < end; i++)
-    {
-      residueHash = new Hashtable();
-      maxCount = 0;
-      maxResidue = "";
-      nongap = 0;
-      values = new int[255];
-
-      for (j = 0; j < jSize; j++)
-      {
-        if (sequences[j] == null)
-        {
-          System.err
-                  .println("WARNING: Consensus skipping null sequence - possible race condition.");
-          continue;
-        }
-        seq = sequences[j].getSequence();
-        if (seq.length > i)
-        {
-          c = seq[i];
-
-          if (c == '.' || c == ' ')
-          {
-            c = '-';
-          }
-
-          if (c == '-')
-          {
-            values['-']++;
-            continue;
-          }
-          else if ('a' <= c && c <= 'z')
-          {
-            c -= 32; // ('a' - 'A');
-          }
-
-          nongap++;
-          values[c]++;
-
-        }
-        else
-        {
-          values['-']++;
-        }
-      }
-
-      for (v = 'A'; v < 'Z'; v++)
-      {
-        if (values[v] < 2 || values[v] < maxCount)
-        {
-          continue;
-        }
-
-        if (values[v] > maxCount)
-        {
-          maxResidue = String.valueOf((char) v);
-        }
-        else if (values[v] == maxCount)
-        {
-          maxResidue += String.valueOf((char) v);
-        }
-        maxCount = values[v];
-      }
-
-      if (maxResidue.length() == 0)
-      {
-        maxResidue = "-";
-      }
-      if (profile)
-      {
-        residueHash.put(PROFILE, new int[][]
-        { values, new int[]
-        { jSize, nongap } });
-      }
-      residueHash.put(MAXCOUNT, new Integer(maxCount));
-      residueHash.put(MAXRESIDUE, maxResidue);
-
-      percentage = ((float) maxCount * 100) / (float) jSize;
-      residueHash.put(PID_GAPS, new Float(percentage));
-
-      percentage = ((float) maxCount * 100) / (float) nongap;
-      residueHash.put(PID_NOGAPS, new Float(percentage));
-      result[i] = residueHash;
-    }
-  }
-
+  /**
+   * Returns the 3' position of a base pair
+   * 
+   * @param pairs
+   *          Secondary structure annotation
+   * @param indice
+   *          5' position of a base pair
+   * @return 3' position of a base pair
+   */
   public static int findPair(SequenceFeature[] pairs, int indice)
   {
+
     for (int i = 0; i < pairs.length; i++)
     {
       if (pairs[i].getBegin() == indice)
+
       {
+
         return pairs[i].getEnd();
+
       }
     }
     return -1;
@@ -210,17 +96,18 @@ public class StructureFrequency
           int end, Hashtable[] result, boolean profile,
           AlignmentAnnotation rnaStruc)
   {
-    // System.out.println("StructureFrequency.calculateNEW4");
+
     Hashtable residueHash;
     String maxResidue;
-    char[] seq, struc = rnaStruc.getRNAStruc().toCharArray();
+    char[] struc = rnaStruc.getRNAStruc().toCharArray();
+
     SequenceFeature[] rna = rnaStruc._rnasecstr;
     char c, s, cEnd;
-    int count, nonGap = 0, i, bpEnd = -1, j, jSize = sequences.length;
+    int count = 0, nonGap = 0, i, bpEnd = -1, j, jSize = sequences.length;
     int[] values;
     int[][] pairs;
     float percentage;
-
+    boolean wooble = true;
     for (i = start; i < end; i++) // foreach column
     {
       residueHash = new Hashtable();
@@ -228,14 +115,22 @@ public class StructureFrequency
       values = new int[255];
       pairs = new int[255][255];
       bpEnd = -1;
+      // System.out.println("s="+struc[i]);
+      if (i < struc.length)
+      {
+        s = struc[i];
 
-      s = struc[i];
+      }
+      else
+      {
+        s = '-';
+      }
       if (s == '.' || s == ' ')
       {
         s = '-';
       }
 
-      if (s != '(')
+      if (s != '(' && s != '[')
       {
         if (s == '-')
         {
@@ -244,19 +139,21 @@ public class StructureFrequency
       }
       else
       {
-        for (j = 0; j < jSize; j++) // foreach row
-        {
-          if (sequences[j] == null)
-          {
-            System.err
-                    .println("WARNING: Consensus skipping null sequence - possible race condition.");
-            continue;
-          }
-          seq = sequences[j].getSequence();
 
-          if (seq.length > i)
+        bpEnd = findPair(rna, i);
+
+        if (bpEnd > -1)
+        {
+          for (j = 0; j < jSize; j++) // foreach row
           {
-            c = seq[i];
+            if (sequences[j] == null)
+            {
+              System.err
+                      .println("WARNING: Consensus skipping null sequence - possible race condition.");
+              continue;
+            }
+            c = sequences[j].getCharAt(i);
+            // System.out.println("c="+c);
 
             // standard representation for gaps in sequence and structure
             if (c == '.' || c == ' ')
@@ -269,17 +166,26 @@ public class StructureFrequency
               values['-']++;
               continue;
             }
-            bpEnd = findPair(rna, i);
-            cEnd = seq[bpEnd];
-            if (checkBpType(c, cEnd))
+            cEnd = sequences[j].getCharAt(bpEnd);
+
+            // System.out.println("pairs ="+c+","+cEnd);
+            if (checkBpType(c, cEnd) == true)
+            {
+              values['(']++; // H means it's a helix (structured)
+              maxResidue = "(";
+              wooble = true;
+              // System.out.println("It's a pair wc");
+
+            }
+            if (checkBpType(c, cEnd) == false)
             {
-              values['H']++; // H means it's a helix (structured)
+              wooble = false;
+              values['[']++; // H means it's a helix (structured)
+              maxResidue = "[";
+
             }
-            // System.out.println("pair: "+c+","+cEnd);
             pairs[c][cEnd]++;
-            // System.out.println("pairs: "+c+","+cEnd+" - "+pairs[c][cEnd]);
 
-            maxResidue = "H";
           }
         }
         // nonGap++;
@@ -287,20 +193,24 @@ public class StructureFrequency
       // UPDATE this for new values
       if (profile)
       {
-        // System.out.println("profile");
-        residueHash.put(PROFILE, new int[][]
-        { values, new int[]
-        { jSize, (jSize - values['-']) } });
+        // TODO 1-dim array with jsize in [0], nongapped in [1]; or Pojo
+        residueHash.put(PROFILE, new int[][] { values,
+            new int[] { jSize, (jSize - values['-']) } });
 
         residueHash.put(PAIRPROFILE, pairs);
       }
-
-      count = values['H'];
-
+      if (wooble == true)
+      {
+        count = values['('];
+      }
+      if (wooble == false)
+      {
+        count = values['['];
+      }
       residueHash.put(MAXCOUNT, new Integer(count));
       residueHash.put(MAXRESIDUE, maxResidue);
 
-      percentage = ((float) count * 100) / (float) jSize;
+      percentage = ((float) count * 100) / jSize;
       residueHash.put(PID_GAPS, new Float(percentage));
 
       // percentage = ((float) count * 100) / (float) nongap;
@@ -311,7 +221,37 @@ public class StructureFrequency
       }
       if (bpEnd > 0)
       {
+        values[')'] = values['('];
+        values[']'] = values['['];
+        values['('] = 0;
+        values['['] = 0;
+        residueHash = new Hashtable();
+        if (wooble == true)
+        {
+          // System.out.println(maxResidue+","+wooble);
+          maxResidue = ")";
+        }
+        if (wooble == false)
+        {
+          // System.out.println(maxResidue+","+wooble);
+          maxResidue = "]";
+        }
+        if (profile)
+        {
+          residueHash.put(PROFILE, new int[][] { values,
+              new int[] { jSize, (jSize - values['-']) } });
+
+          residueHash.put(PAIRPROFILE, pairs);
+        }
+
+        residueHash.put(MAXCOUNT, new Integer(count));
+        residueHash.put(MAXRESIDUE, maxResidue);
+
+        percentage = ((float) count * 100) / jSize;
+        residueHash.put(PID_GAPS, new Float(percentage));
+
         result[bpEnd] = residueHash;
+
       }
     }
   }
@@ -402,21 +342,8 @@ public class StructureFrequency
   public static void completeConsensus(AlignmentAnnotation consensus,
           Hashtable[] hconsensus, int iStart, int width,
           boolean ignoreGapsInConsensusCalculation,
-          boolean includeAllConsSymbols)
-  {
-    completeConsensus(consensus, hconsensus, iStart, width,
-            ignoreGapsInConsensusCalculation, includeAllConsSymbols, null); // new
-    // char[]
-    // { 'A', 'C', 'G', 'T', 'U' });
-  }
-
-  public static void completeConsensus(AlignmentAnnotation consensus,
-          Hashtable[] hconsensus, int iStart, int width,
-          boolean ignoreGapsInConsensusCalculation,
-          boolean includeAllConsSymbols, char[] alphabet)
+          boolean includeAllConsSymbols, long nseq)
   {
-    System.out.println("StructureFrequency.completeConsensus "
-            + includeAllConsSymbols);
     float tval, value;
     if (consensus == null || consensus.annotations == null
             || consensus.annotations.length < width)
@@ -425,9 +352,23 @@ public class StructureFrequency
       // initialised properly
       return;
     }
+    String fmtstr = "%3.1f";
+    int precision = 2;
+    while (nseq > 100)
+    {
+      precision++;
+      nseq /= 10;
+    }
+    if (precision > 2)
+    {
+      fmtstr = "%" + (2 + precision) + "." + precision + "f";
+    }
+    Format fmt = new Format(fmtstr);
+
     for (int i = iStart; i < width; i++)
     {
-      if (i >= hconsensus.length)
+      Hashtable hci;
+      if (i >= hconsensus.length || ((hci = hconsensus[i]) == null))
       {
         // happens if sequences calculated over were shorter than alignment
         // width
@@ -435,78 +376,85 @@ public class StructureFrequency
         continue;
       }
       value = 0;
+      Float fv;
       if (ignoreGapsInConsensusCalculation)
       {
-        value = ((Float) hconsensus[i].get(StructureFrequency.PID_NOGAPS))
-                .floatValue();
+        fv = (Float) hci.get(StructureFrequency.PID_NOGAPS);
       }
       else
       {
-        value = ((Float) hconsensus[i].get(StructureFrequency.PID_GAPS))
-                .floatValue();
+        fv = (Float) hci.get(StructureFrequency.PID_GAPS);
       }
-
-      String maxRes = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
-              .toString();
-      String mouseOver = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
-              + " ";
+      if (fv == null)
+      {
+        consensus.annotations[i] = null;
+        // data has changed below us .. give up and
+        continue;
+      }
+      value = fv.floatValue();
+      String maxRes = hci.get(StructureFrequency.MAXRESIDUE).toString();
+      String mouseOver = hci.get(StructureFrequency.MAXRESIDUE) + " ";
       if (maxRes.length() > 1)
       {
         mouseOver = "[" + maxRes + "] ";
         maxRes = "+";
       }
-      int[][] profile = (int[][]) hconsensus[i]
-              .get(StructureFrequency.PROFILE);
-      if (profile != null && includeAllConsSymbols) // Just responsible for the
+      int[][] profile = (int[][]) hci.get(StructureFrequency.PROFILE);
+      int[][] pairs = (int[][]) hci.get(StructureFrequency.PAIRPROFILE);
+
+      if (pairs != null && includeAllConsSymbols) // Just responsible for the
       // tooltip
+      // TODO Update tooltips for Structure row
       {
-        // System.out.println("StructureFrequency.includeAllConsSymbols");
         mouseOver = "";
-        if (alphabet != null)
+
+        /*
+         * TODO It's not sure what is the purpose of the alphabet and wheter it
+         * is useful for structure?
+         * 
+         * if (alphabet != null) { for (int c = 0; c < alphabet.length; c++) {
+         * tval = ((float) profile[0][alphabet[c]]) 100f / (float)
+         * profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; mouseOver +=
+         * ((c == 0) ? "" : "; ") + alphabet[c] + " " + ((int) tval) + "%"; } }
+         * else {
+         */
+        int[][] ca = new int[625][];
+        float[] vl = new float[625];
+        int x = 0;
+        for (int c = 65; c < 90; c++)
         {
-          for (int c = 0; c < alphabet.length; c++)
+          for (int d = 65; d < 90; d++)
           {
-            tval = ((float) profile[0][alphabet[c]])
-                    * 100f
-                    / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
-                            : 0];
-            mouseOver += ((c == 0) ? "" : "; ") + alphabet[c] + " "
-                    + ((int) tval) + "%";
+            ca[x] = new int[] { c, d };
+            vl[x] = pairs[c][d];
+            x++;
           }
         }
-        else
+        jalview.util.QuickSort.sort(vl, ca);
+        int p = 0;
+
+        /*
+         * profile[1] is {total, ungappedTotal}
+         */
+        final int divisor = profile[1][ignoreGapsInConsensusCalculation ? 1
+                : 0];
+        for (int c = 624; c > 0; c--)
         {
-          // System.out.println("StructureFrequency.NOTincludeAllConsSymbols");
-          Object[] ca = new Object[profile[0].length];
-          float[] vl = new float[profile[0].length];
-          for (int c = 0; c < ca.length; c++)
-          {
-            ca[c] = new char[]
-            { (char) c };
-            vl[c] = (float) profile[0][c];
-          }
-          ;
-          jalview.util.QuickSort.sort(vl, ca);
-          for (int p = 0, c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
+          if (vl[c] > 0)
           {
-            if (((char[]) ca[c])[0] != '-')
-            {
-              tval = ((float) profile[0][((char[]) ca[c])[0]])
-                      * 100f
-                      / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
-                              : 0];
-              mouseOver += ((p == 0) ? "" : "; ") + ((char[]) ca[c])[0]
-                      + " " + ((int) tval) + "%";
-              p++;
+            tval = (vl[c] * 100f / divisor);
+            mouseOver += ((p == 0) ? "" : "; ") + (char) ca[c][0]
+                    + (char) ca[c][1] + " " + fmt.form(tval) + "%";
+            p++;
 
-            }
           }
-
         }
+
+        // }
       }
       else
       {
-        mouseOver += ((int) value + "%");
+        mouseOver += (fmt.form(value) + "%");
       }
       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
               value);
@@ -514,43 +462,62 @@ public class StructureFrequency
   }
 
   /**
-   * get the sorted profile for the given position of the consensus
+   * get the sorted base-pair profile for the given position of the consensus
    * 
    * @param hconsensus
-   * @return
+   * @return profile of the given column
    */
   public static int[] extractProfile(Hashtable hconsensus,
-          boolean ignoreGapsInConsensusCalculation, int column)
+          boolean ignoreGapsInConsensusCalculation)
   {
-    // TODO is there a more elegant way to acces the column number?
-    /*
-     * calculate the frequence of the 16 bp variations for this column 'somehow'
-     * transfer this via getProfile and let it correctly draw
-     */
-    int[] rtnval = new int[51]; // 2*(5*5)+1
+    int[] rtnval = new int[STRUCTURE_PROFILE_LENGTH]; // 2*(5*5)+2
     int[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
     int[][] pairs = (int[][]) hconsensus
             .get(StructureFrequency.PAIRPROFILE);
 
     if (profile == null)
+    {
       return null;
-    
-    rtnval[0] = 1;
-    for (int j = 65; j <= 90; j++)
+    }
+
+    // TODO fix the object length, also do it in completeConsensus
+    // Object[] ca = new Object[625];
+    int[][] ca = new int[625][];
+    float[] vl = new float[625];
+    int x = 0;
+    for (int c = 65; c < 90; c++)
     {
-      for (int k = 65; k <= 90; k++)
+      for (int d = 65; d < 90; d++)
       {
-        if (pairs[j][k] > 0)
-        {
-          rtnval[rtnval[0]++] = j;
-          rtnval[rtnval[0]++] = k;
-          rtnval[rtnval[0]++] = (int) ((float) pairs[j][k] * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
-                  : 0]);
-        }
+        ca[x] = new int[] { c, d };
+        vl[x] = pairs[c][d];
+        x++;
+      }
+    }
+    jalview.util.QuickSort.sort(vl, ca);
+
+    int valuesCount = 0;
+    rtnval[1] = 0;
+    int offset = 2;
+    final int divisor = profile[1][ignoreGapsInConsensusCalculation ? 1 : 0];
+    for (int c = 624; c > 0; c--)
+    {
+      if (vl[c] > 0)
+      {
+        rtnval[offset++] = ca[c][0];
+        rtnval[offset++] = ca[c][1];
+        rtnval[offset] = (int) (vl[c] * 100f / divisor);
+        rtnval[1] += rtnval[offset++];
+        valuesCount++;
       }
     }
+    rtnval[0] = valuesCount;
 
-    return rtnval;
+    // insert profile type code in position 0
+    int[] result = new int[rtnval.length + 1];
+    result[0] = AlignmentAnnotation.STRUCTURE_PROFILE;
+    System.arraycopy(rtnval, 0, result, 1, rtnval.length);
+    return result;
   }
 
   public static void main(String args[])