JAL-2416 quality calculation now only depends on loaded score matrix

diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java
index 79fa3d6..d1a3c37 100755 (executable)
--- a/src/jalview/analysis/Conservation.java
+++ b/src/jalview/analysis/Conservation.java
@@ -20,6 +20,7 @@
  */
 package jalview.analysis;
 
+import jalview.analysis.scoremodels.ScoreMatrix;
 import jalview.analysis.scoremodels.ScoreModels;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.Annotation;
@@ -163,27 +164,30 @@ public class Conservation
   }
 
   /**
-   * Translate sequence i into a numerical representation and store it in the
-   * i'th position of the seqNums array.
+   * Translate sequence i into score matrix indices and store it in the i'th
+   * position of the seqNums array.
    * 
    * @param i
+   * @param sm
    */
-  private void calcSeqNum(int i)
+  private void calcSeqNum(int i, ScoreMatrix sm)
   {
-    String sq = null; // for dumb jbuilder not-inited exception warning
-    int[] sqnum = null;
-
+    int gapIndex = sm.getMatrixIndex(' ');
     int sSize = sequences.length;
 
     if ((i > -1) && (i < sSize))
     {
-      sq = sequences[i].getSequenceAsString();
+      String sq = sequences[i].getSequenceAsString();
 
       if (seqNums.size() <= i)
       {
         seqNums.addElement(new int[sq.length() + 1]);
       }
 
+      /*
+       * the first entry in the array is the sequence's hashcode,
+       * following entries are matrix indices of sequence characters
+       */
       if (sq.hashCode() != seqNums.elementAt(i)[0])
       {
         int j;
@@ -196,14 +200,18 @@ public class Conservation
           maxLength = len;
         }
 
-        sqnum = new int[len + 1]; // better to always make a new array -
+        int[] sqnum = new int[len + 1]; // better to always make a new array -
         // sequence can change its length
         sqnum[0] = sq.hashCode();
 
         for (j = 1; j <= len; j++)
         {
-          sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
-                  .charAt(j - 1)];
+          // sqnum[j] = ResidueProperties.aaIndex[sq.charAt(j - 1)];
+          sqnum[j] = sm.getMatrixIndex(sq.charAt(j - 1));
+          if (sqnum[j] == -1)
+          {
+            sqnum[j] = gapIndex;
+          }
         }
 
         seqNums.setElementAt(sqnum, i);
@@ -533,37 +541,32 @@ public class Conservation
 
   /**
    * DOCUMENT ME!
+   * 
+   * @param sm
    */
-  private void percentIdentity2()
+  private void percentIdentity(ScoreMatrix sm)
   {
     seqNums = new Vector<int[]>();
-    // calcSeqNum(s);
     int i = 0, iSize = sequences.length;
     // Do we need to calculate this again?
     for (i = 0; i < iSize; i++)
     {
-      calcSeqNum(i);
+      calcSeqNum(i, sm);
     }
 
+    int gapIndex = sm.getMatrixIndex(' ');
+
     if ((cons2 == null) || seqNumsChanged)
     {
+      // FIXME remove magic number 24 without changing calc
+      // sm.getSize() returns 25 so doesn't quite do it...
       cons2 = new int[maxLength][24];
 
-      // Initialize the array
-      for (int j = 0; j < 24; j++)
-      {
-        for (i = 0; i < maxLength; i++)
-        {
-          cons2[i][j] = 0;
-        }
-      }
-
-      int[] sqnum;
       int j = 0;
 
       while (j < sequences.length)
       {
-        sqnum = seqNums.elementAt(j);
+        int[] sqnum = seqNums.elementAt(j);
 
         for (i = 1; i < sqnum.length; i++)
         {
@@ -572,21 +575,10 @@ public class Conservation
 
         for (i = sqnum.length - 1; i < maxLength; i++)
         {
-          cons2[i][23]++; // gap count
+          cons2[i][gapIndex]++; // gap count
         }
-
         j++;
       }
-
-      // unnecessary ?
-
-      /*
-       * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
-       * maxj = -1;
-       * 
-       * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
-       * maxi = i; maxj = j; } } }
-       */
     }
   }
 
@@ -602,14 +594,15 @@ public class Conservation
   {
     quality = new Vector<Double>();
 
-    double max = -10000;
-    float[][] BLOSUM62 = ScoreModels.getInstance().getBlosum62()
-            .getMatrix();
+    double max = -Double.MAX_VALUE;
+    ScoreMatrix blosum62 = ScoreModels.getInstance().getBlosum62();
+    float[][] blosumScores = blosum62.getMatrix();
+    int gapIndex = blosum62.getMatrixIndex(' ');
 
     // Loop over columns // JBPNote Profiling info
     // long ts = System.currentTimeMillis();
     // long te = System.currentTimeMillis();
-    percentIdentity2();
+    percentIdentity(blosum62);
 
     int size = seqNums.size();
     int[] lengths = new int[size];
@@ -622,20 +615,24 @@ public class Conservation
       lengths[l] = seqNums.elementAt(l).length - 1;
     }
 
+    // todo ? remove '*' (unused?) from score matrix and
+    // use getSize() here instead of getSize() - 1 ??
+    final int symbolCount = blosum62.getSize() - 1; // 24;
+
     for (j = startRes; j <= endRes; j++)
     {
       bigtot = 0;
 
       // First Xr = depends on column only
-      x = new double[24];
+      x = new double[symbolCount];
 
-      for (ii = 0; ii < 24; ii++)
+      for (ii = 0; ii < symbolCount; ii++)
       {
         x[ii] = 0;
 
-        for (i2 = 0; i2 < 24; i2++)
+        for (i2 = 0; i2 < symbolCount; i2++)
         {
-          x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
+          x[ii] += (((double) cons2[j][i2] * blosumScores[ii][i2]) + 4);
         }
 
         x[ii] /= size;
@@ -645,18 +642,16 @@ public class Conservation
       for (k = 0; k < size; k++)
       {
         tot = 0;
-        xx = new double[24];
-        seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : 23; // Sequence,
-                                                                      // or gap
-                                                                      // at the
-                                                                      // end
+        xx = new double[symbolCount];
+        seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : gapIndex;
+        // Sequence, or gap at the end
 
         // This is a loop over r
-        for (i = 0; i < 23; i++)
+        for (i = 0; i < symbolCount - 1; i++)
         {
           sr = 0;
 
-          sr = (double) BLOSUM62[i][seqNum] + 4;
+          sr = (double) blosumScores[i][seqNum] + 4;
 
           // Calculate X with another loop over residues
           // System.out.println("Xi " + i + " " + x[i] + " " + sr);
@@ -680,12 +675,13 @@ public class Conservation
       // Need to normalize by gaps
     }
 
-    double newmax = -10000;
+    double newmax = -Double.MAX_VALUE;
 
     for (j = startRes; j <= endRes; j++)
     {
       tmp = quality.elementAt(j).doubleValue();
-      tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+      // tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+      tmp = ((max - tmp) * (size - cons2[j][gapIndex])) / size;
 
       // System.out.println(tmp+ " " + j);
       quality.setElementAt(new Double(tmp), j);