JAL-885; Structure conservation row is just added if alignment
[jalview.git] / src / jalview / analysis / StructureFrequency.java
index 29c0a0b..9bdaa10 100644 (file)
@@ -44,300 +44,214 @@ public class StructureFrequency
 
   public static final String PROFILE = "P";
 
-  public static final Hashtable[] calculate(Vector sequences, int start,
-          int end)
-  {
-    return calculate(sequences, start, end, false);
-  }
+  public static final String PAIRPROFILE = "B";
 
-  public static final Hashtable[] calculate(Vector sequences, int start,
-          int end, boolean profile)
+/**
+ * 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)
   {
-    SequenceI[] seqs = new SequenceI[sequences.size()];
-    int width = 0;
-    for (int i = 0; i < sequences.size(); i++)
+    for (int i = 0; i < pairs.length; i++)
     {
-      seqs[i] = (SequenceI) sequences.elementAt(i);
-      if (seqs[i].getLength() > width)
+      if (pairs[i].getBegin() == indice)
       {
-        width = seqs[i].getLength();
+        return pairs[i].getEnd();
       }
     }
-
-    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);
+    return -1;
   }
 
+  /**
+   * Method to calculate a 'base pair consensus row', very similar to nucleotide
+   * consensus but takes into account a given structure
+   * 
+   * @param sequences
+   * @param start
+   * @param end
+   * @param result
+   * @param profile
+   * @param rnaStruc
+   */
   public static final void calculate(SequenceI[] sequences, int start,
-          int end, Hashtable[] result, boolean profile)
+          int end, Hashtable[] result, boolean profile,
+          AlignmentAnnotation rnaStruc)
   {
     Hashtable residueHash;
-    int maxCount, nongap, i, j, v, jSize = sequences.length;
     String maxResidue;
-    char c;
+    char[] seq, struc = rnaStruc.getRNAStruc().toCharArray();
+    SequenceFeature[] rna = rnaStruc._rnasecstr;
+    char c, s, cEnd;
+    int count, nonGap = 0, i, bpEnd = -1, j, jSize = sequences.length;
+    int[] values;
+    int[][] pairs;
     float percentage;
 
-    int[] values = new int[255];
-
-    char[] seq;
-
-    for (i = start; i < end; i++)
+    for (i = start; i < end; i++) // foreach column
     {
       residueHash = new Hashtable();
-      maxCount = 0;
-      maxResidue = "";
-      nongap = 0;
+      maxResidue = "-";
       values = new int[255];
-      
+      pairs = new int[255][255];
+      bpEnd = -1;
+
+      s = struc[i];
+      if (s == '.' || s == ' ')
+      {
+        s = '-';
+      }
 
-      for (j = 0; j < jSize; j++)
+      if (s != '(')
       {
-        if (sequences[j]==null)
+        if (s == '-')
         {
-          System.err.println("WARNING: Consensus skipping null sequence - possible race condition.");
-          continue;
+          values['-']++;
         }
-        seq = sequences[j].getSequence();
-        if (seq.length > i)
+      }
+      else
+      {
+        for (j = 0; j < jSize; j++) // foreach row
         {
-          c = seq[i];
-
-          if (c == '.' || c == ' ')
-          {
-            c = '-';
-          }
-
-          if (c == '-')
+          if (sequences[j] == null)
           {
-            values['-']++;
+            System.err
+                    .println("WARNING: Consensus skipping null sequence - possible race condition.");
             continue;
           }
-          else if ('a' <= c && c <= 'z')
-          {
-            c -= 32; // ('a' - 'A');
-          }
+          seq = sequences[j].getSequence();
 
-          nongap++;
-          values[c]++;
+          if (seq.length > i)
+          {
+            c = seq[i];
 
-        }
-        else
-        {
-          values['-']++;
-        }
-      }
+            // standard representation for gaps in sequence and structure
+            if (c == '.' || c == ' ')
+            {
+              c = '-';
+            }
 
-      for (v = 'A'; v < 'Z'; v++)
-      {
-        if (values[v] < 2 || values[v] < maxCount)
-        {
-          continue;
-        }
+            if (c == '-')
+            {
+              values['-']++;
+              continue;
+            }
+            bpEnd = findPair(rna, i);
+            cEnd = seq[bpEnd];
+            if (checkBpType(c, cEnd))
+            {
+              values['H']++; // H means it's a helix (structured)
+            }
+            pairs[c][cEnd]++;
 
-        if (values[v] > maxCount)
-        {
-          maxResidue = String.valueOf((char) v);
-        }
-        else if (values[v] == maxCount)
-        {
-          maxResidue += String.valueOf((char) v);
+            maxResidue = "H";
+          }
         }
-        maxCount = values[v];
-      }
-
-      if (maxResidue.length() == 0)
-      {
-        maxResidue = "-";
+        // nonGap++;
       }
+      // UPDATE this for new values
       if (profile)
       {
         residueHash.put(PROFILE, new int[][]
         { values, new int[]
-        { jSize, nongap } });
+        { jSize, (jSize - values['-']) } });
+
+        residueHash.put(PAIRPROFILE, pairs);
       }
-      residueHash.put(MAXCOUNT, new Integer(maxCount));
+
+      count = values['H'];
+
+      residueHash.put(MAXCOUNT, new Integer(count));
       residueHash.put(MAXRESIDUE, maxResidue);
 
-      percentage = ((float) maxCount * 100) / (float) jSize;
+      percentage = ((float) count * 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;
+      // percentage = ((float) count * 100) / (float) nongap;
+      // residueHash.put(PID_NOGAPS, new Float(percentage));
+      if (result[i] == null)
+      {
+        result[i] = residueHash;
+      }
+      if (bpEnd > 0)
+      {
+        result[bpEnd] = residueHash;
+      }
     }
   }
-  
-  /**
-   * Method to calculate a 'base pair consensus row', very similar 
-   * to nucleotide consensus but takes into account a given structure
-   * @param sequences
-   * @param start
-   * @param end
-   * @param result
-   * @param profile
-   * @param rnaStruc
-   */
-  public static final void calculate(SequenceI[] sequences, int start,
-                 int end, Hashtable[] result, boolean profile, AlignmentAnnotation rnaStruc){
-
-         Hashtable residueHash;
-         String maxResidue;
-         char[] seq, struc=rnaStruc.getRNAStruc().toCharArray();
-         SequenceFeature[] rna =rnaStruc._rnasecstr;
-         char c,s,cEnd;
-         int count,nonGap=0,i,bpEnd=-1,j,jSize = sequences.length;
-         int[] values = new int[255];
-         float percentage;
-         
-         for (i = start; i < end; i++) //foreach column
-         {
-                 residueHash = new Hashtable();
-                 maxResidue="";
-                 values = new int[255];
-                 bpEnd=-1;
-
-                 s = struc[i];
-                 if (s == '.' || s == ' ')
-                 {
-                         s = '-';
-                 }
-
-                 if(s == '-'){
-                         values['-']++;
-                 }
-                 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)
-                                 {
-                                         c = seq[i];
-
-                                         //standard representation for gaps in sequence and structure
-                                         if (c == '.' || c == ' ')
-                                         {
-                                                 c = '-';
-                                         }
-
-                                         if (c == '-')
-                                         {
-                                                 values['-']++;
-                                                 continue;
-                                         }
-                                         if(s == '('){
-                                                 bpEnd=rna[(rna.length-1-nonGap)].getEnd();              
-                                                 cEnd=seq[bpEnd];
-                                                 if(checkBpType(c,cEnd)){
-                                                         values['H']++; //H means it's a helix (structured)
-                                                 }
-                                         }
-                                 }
-                         }
-                         nonGap++;
-                 }
-                 /*UPDATE this for new values
-             if (profile)
-             {
-               residueHash.put(PROFILE, new int[][]
-               { values, new int[]
-               { jSize, nongap } });
-             }
-                  */
-
-                 count=values['H'];
-                 
-                 if (count == 0)
-                 {
-                         maxResidue = "-";
-                 }else{
-                         maxResidue="H";
-                 }
-                 
-                 residueHash.put(MAXCOUNT, new Integer(count));
-             residueHash.put(MAXRESIDUE, maxResidue);
-             
-                 percentage = ((float) count * 100) / (float) jSize;
-                 residueHash.put(PID_GAPS, new Float(percentage));
-
-                 //percentage = ((float) count * 100) / (float) nongap;
-                 //residueHash.put(PID_NOGAPS, new Float(percentage));
-                 if(result[i]==null){
-                         result[i] = residueHash;
-                 }
-                 if(bpEnd>0){
-                         result[bpEnd]=residueHash;
-                 }
-
-         }
-  }
 
-        
   /**
-   * Method to check if a base-pair is a canonical or a wobble bp 
-   * @param up 5' base
-   * @param down 3' base
+   * Method to check if a base-pair is a canonical or a wobble bp
+   * 
+   * @param up
+   *          5' base
+   * @param down
+   *          3' base
    * @return True if it is a canonical/wobble bp
    */
-  public static boolean checkBpType(char up, char down){
-         if(up>'Z'){up-=32;}
-         if(down>'Z'){down-=32;}
-         
-         switch (up){
-               case 'A': 
-                       switch (down){
-                               case 'T': return true;
-                               case 'U': return true;
-                       }
-               break;
-               case 'C': 
-                       switch (down){
-                               case 'G': return true;
-                       }
-               break;
-               case 'T': 
-                       switch (down){
-                               case 'A': return true;
-                               case 'G': return true;
-                               }
-               break;
-               case 'G': 
-                       switch (down){
-                               case 'C': return true;
-                               case 'T': return true;
-                               case 'U': return true;
-                       }
-               break;
-               case 'U': 
-                       switch (down){
-                               case 'A': return true;
-                               case 'G': return true;
-                       }
-               break;
-         }       
-         return false;
+  public static boolean checkBpType(char up, char down)
+  {
+    if (up > 'Z')
+    {
+      up -= 32;
+    }
+    if (down > 'Z')
+    {
+      down -= 32;
+    }
+
+    switch (up)
+    {
+    case 'A':
+      switch (down)
+      {
+      case 'T':
+        return true;
+      case 'U':
+        return true;
+      }
+      break;
+    case 'C':
+      switch (down)
+      {
+      case 'G':
+        return true;
+      }
+      break;
+    case 'T':
+      switch (down)
+      {
+      case 'A':
+        return true;
+      case 'G':
+        return true;
+      }
+      break;
+    case 'G':
+      switch (down)
+      {
+      case 'C':
+        return true;
+      case 'T':
+        return true;
+      case 'U':
+        return true;
+      }
+      break;
+    case 'U':
+      switch (down)
+      {
+      case 'A':
+        return true;
+      case 'G':
+        return true;
+      }
+      break;
+    }
+    return false;
   }
-  
+
   /**
    * Compute all or part of the annotation row from the given consensus
    * hashtable
@@ -357,7 +271,7 @@ public class StructureFrequency
   {
     completeConsensus(consensus, hconsensus, iStart, width,
             ignoreGapsInConsensusCalculation, includeAllConsSymbols, null); // new
-                                                                            // char[]
+    // char[]
     // { 'A', 'C', 'G', 'T', 'U' });
   }
 
@@ -394,16 +308,20 @@ public class StructureFrequency
         value = ((Float) hconsensus[i].get(StructureFrequency.PID_GAPS))
                 .floatValue();
       }
-      
-      String maxRes = hconsensus[i].get(StructureFrequency.MAXRESIDUE).toString();
-      String mouseOver = hconsensus[i].get(StructureFrequency.MAXRESIDUE) + " ";
+
+      String maxRes = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
+              .toString();
+      String mouseOver = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
+              + " ";
       if (maxRes.length() > 1)
       {
         mouseOver = "[" + maxRes + "] ";
         maxRes = "+";
       }
-      int[][] profile = (int[][]) hconsensus[i].get(StructureFrequency.PROFILE);
-      if (profile != null && includeAllConsSymbols)
+      int[][] profile = (int[][]) hconsensus[i]
+              .get(StructureFrequency.PROFILE);
+      if (profile != null && includeAllConsSymbols) // Just responsible for the
+      // tooltip
       {
         mouseOver = "";
         if (alphabet != null)
@@ -457,56 +375,61 @@ 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)
+          boolean ignoreGapsInConsensusCalculation, int column)
   {
-    int[] rtnval = new int[64];
+    // 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[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
+    int[][] pairs = (int[][]) hconsensus
+            .get(StructureFrequency.PAIRPROFILE);
+
     if (profile == null)
       return null;
-    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);
+    
     rtnval[0] = 1;
-    for (int c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
+    for (int j = 65; j <= 90; j++)
     {
-      if (((char[]) ca[c])[0] != '-')
+      for (int k = 65; k <= 90; k++)
       {
-        rtnval[rtnval[0]++] = ((char[]) ca[c])[0];
-        rtnval[rtnval[0]++] = (int) (((float) profile[0][((char[]) ca[c])[0]]) * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
-                : 0]);
+        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]);
+        }
       }
     }
+
     return rtnval;
   }
 
-  enum base {A,T,g,C};
-  
-  public static void main(String args[]){
-         //Short test to see if checkBpType works
-         ArrayList<String> test = new ArrayList<String>();
-         test.add("A");
-         test.add("c");
-         test.add("g");
-         test.add("T");
-         test.add("U");
-         for (String i : test) {
-                 for (String j : test) {
-                         System.out.println(i+"-"+j+": "+StructureFrequency.checkBpType(i.charAt(0),j.charAt(0)));
-                 }
-         }
+  public static void main(String args[])
+  {
+    // Short test to see if checkBpType works
+    ArrayList<String> test = new ArrayList<String>();
+    test.add("A");
+    test.add("c");
+    test.add("g");
+    test.add("T");
+    test.add("U");
+    for (String i : test)
+    {
+      for (String j : test)
+      {
+        System.out.println(i + "-" + j + ": "
+                + StructureFrequency.checkBpType(i.charAt(0), j.charAt(0)));
+      }
+    }
   }
 }