quick fix for conservation calculation when residue ambiguity codes are used.
[jalview.git] / src / jalview / analysis / Conservation.java
index 4eb1218..c216471 100755 (executable)
@@ -24,44 +24,68 @@ import jalview.datamodel.*;
 
 /**
  * Calculates conservation values for a given set of sequences
- *
+ * 
  * @author $author$
  * @version $Revision$
  */
 public class Conservation
 {
   SequenceI[] sequences;
+
   int start;
+
   int end;
+
   Vector seqNums; // vector of int vectors where first is sequence checksum
-  int maxLength = 0; //  used by quality calcs
+
+  int maxLength = 0; // used by quality calcs
+
   boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
+
   Hashtable[] total;
 
+  boolean canonicaliseAa = true; // if true then conservation calculation will
+                                  // map all symbols to canonical aa numbering
+                                  // rather than consider conservation of that
+                                  // symbol
+
   /** Stores calculated quality values */
   public Vector quality;
 
-  /** Stores maximum and minimum values of quality values  */
+  /** Stores maximum and minimum values of quality values */
   public Double[] qualityRange = new Double[2];
+
   String consString = "";
+
   Sequence consSequence;
+
   Hashtable propHash;
+
   int threshold;
+
   String name = "";
+
   int[][] cons2;
 
   /**
    * Creates a new Conservation object.
-   *
-   * @param name Name of conservation
-   * @param propHash DOCUMENT ME!
-   * @param threshold to count the residues in residueHash(). commonly used value is 3
-   * @param sequences sequences to be used in calculation
-   * @param start start residue position
-   * @param end end residue position
+   * 
+   * @param name
+   *                Name of conservation
+   * @param propHash
+   *                hash of properties for each symbol
+   * @param threshold
+   *                to count the residues in residueHash(). commonly used value
+   *                is 3
+   * @param sequences
+   *                sequences to be used in calculation
+   * @param start
+   *                start residue position
+   * @param end
+   *                end residue position
    */
   public Conservation(String name, Hashtable propHash, int threshold,
-                      Vector sequences, int start, int end)
+          Vector sequences, int start, int end)
   {
 
     this.name = name;
@@ -70,7 +94,8 @@ public class Conservation
     this.start = start;
     this.end = end;
 
-    maxLength = end - start + 1; // default width includes bounds of calculation
+    maxLength = end - start + 1; // default width includes bounds of
+                                  // calculation
 
     int s, sSize = sequences.size();
     SequenceI[] sarray = new SequenceI[sSize];
@@ -87,9 +112,10 @@ public class Conservation
   }
 
   /**
-   * DOCUMENT ME!
-   *
-   * @param i DOCUMENT ME!
+   * Translate sequence i into a numerical 
+   * representation and store it in the i'th position of the seqNums array.
+   * 
+   * @param i
    */
   private void calcSeqNum(int i)
   {
@@ -98,7 +124,7 @@ public class Conservation
 
     int sSize = sequences.length;
 
-    if ( (i > -1) && (i < sSize))
+    if ((i > -1) && (i < sSize))
     {
       sq = sequences[i].getSequenceAsString();
 
@@ -107,7 +133,7 @@ public class Conservation
         seqNums.addElement(new int[sq.length() + 1]);
       }
 
-      if (sq.hashCode() != ( (int[]) seqNums.elementAt(i))[0])
+      if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
       {
         int j;
         int len;
@@ -119,12 +145,14 @@ public class Conservation
           maxLength = len;
         }
 
-        sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
+        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] = jalview.schemes.ResidueProperties.aaIndex[sq
+                  .charAt(j - 1)];
         }
 
         seqNums.setElementAt(sqnum, i);
@@ -137,8 +165,8 @@ public class Conservation
     else
     {
       // JBPNote INFO level debug
-      System.err.println(
-          "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
+      System.err
+              .println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
     }
   }
 
@@ -166,17 +194,34 @@ public class Conservation
         {
           c = sequences[j].getCharAt(i);
 
-          // No need to check if its a '-'
-          if (c == '.' || c == ' ')
-          {
-            c = '-';
+          if (canonicaliseAa)
+          { // lookup the base aa code symbol
+            c = (char) jalview.schemes.ResidueProperties.aaIndex[sequences[j]
+                    .getCharAt(i)];
+            if (c > 20)
+            {
+              c = '-';
+            }
+            else
+            {
+              // recover canonical aa symbol
+              c = jalview.schemes.ResidueProperties.aa[c].charAt(0);
+            }
           }
-
-          if ('a' <= c && c <= 'z')
+          else
           {
-            c -= (32); // 32 = 'a' - 'A'
-          }
+            // original behaviour - operate on ascii symbols directly
+            // No need to check if its a '-'
+            if (c == '.' || c == ' ')
+            {
+              c = '-';
+            }
 
+            if (!canonicaliseAa && 'a' <= c && c <= 'z')
+            {
+              c -= (32); // 32 = 'a' - 'A'
+            }
+          }
           values[c]++;
         }
         else
@@ -185,10 +230,10 @@ public class Conservation
         }
       }
 
-      //What is the count threshold to count the residues in residueHash()
+      // What is the count threshold to count the residues in residueHash()
       thresh = (threshold * (jSize)) / 100;
 
-      //loop over all the found residues
+      // loop over all the found residues
       resultHash = new Hashtable();
       for (char v = '-'; v < 'Z'; v++)
       {
@@ -197,7 +242,7 @@ public class Conservation
         {
           res = String.valueOf(v);
 
-          //Now loop over the properties
+          // Now loop over the properties
           enumeration2 = propHash.keys();
 
           while (enumeration2.hasMoreElements())
@@ -205,7 +250,7 @@ public class Conservation
             type = (String) enumeration2.nextElement();
             ht = (Hashtable) propHash.get(type);
 
-            //Have we ticked this before?
+            // Have we ticked this before?
             if (!resultHash.containsKey(type))
             {
               if (ht.containsKey(res))
@@ -217,10 +262,10 @@ public class Conservation
                 resultHash.put(type, ht.get("-"));
               }
             }
-            else if ( ( (Integer) resultHash.get(type)).equals(
-                (Integer) ht.get(res)) == false)
+            else if (((Integer) resultHash.get(type)).equals((Integer) ht
+                    .get(res)) == false)
             {
-              resultHash.put(type, new Integer( -1));
+              resultHash.put(type, new Integer(-1));
             }
           }
         }
@@ -230,9 +275,9 @@ public class Conservation
     }
   }
 
-  /***
-   * countConsNGaps
-   * returns gap count in int[0], and conserved residue count in int[1]
+  /*****************************************************************************
+   * count conservation  for the j'th column of the alignment
+   * @return { gap count, conserved residue count}
    */
   public int[] countConsNGaps(int j)
   {
@@ -254,7 +299,7 @@ public class Conservation
 
       c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
 
-      if (jalview.util.Comparison.isGap( (c)))
+      if (jalview.util.Comparison.isGap((c)))
       {
         count++;
       }
@@ -282,9 +327,12 @@ public class Conservation
 
   /**
    * Calculates the conservation sequence
-   *
-   * @param consflag if true, poitiveve conservation; false calculates negative conservation
-   * @param percentageGaps commonly used value is 25
+   * 
+   * @param consflag
+   *                if true, poitiveve conservation; false calculates negative
+   *                conservation
+   * @param percentageGaps
+   *                commonly used value is 25
    */
   public void verdict(boolean consflag, float percentageGaps)
   {
@@ -297,9 +345,9 @@ public class Conservation
     Hashtable resultHash;
     Enumeration enumeration;
 
-    //NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
-    //EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
-    //DOES NOT EXIST IN JALVIEW 2.1.2
+    // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
+    // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
+    // DOES NOT EXIST IN JALVIEW 2.1.2
     for (int i = 0; i < start; i++)
     {
       consString.append('-');
@@ -309,13 +357,13 @@ public class Conservation
     {
       gapcons = countConsNGaps(i);
       totGaps = gapcons[1];
-      pgaps = ( (float) totGaps * 100) / (float) sequences.length;
+      pgaps = ((float) totGaps * 100) / (float) sequences.length;
 
       if (percentageGaps > pgaps)
       {
         resultHash = total[i - start];
 
-        //Now find the verdict
+        // Now find the verdict
         count = 0;
         enumeration = resultHash.keys();
 
@@ -324,7 +372,7 @@ public class Conservation
           type = (String) enumeration.nextElement();
           result = (Integer) resultHash.get(type);
 
-          //Do we want to count +ve conservation or +ve and -ve cons.?
+          // Do we want to count +ve conservation or +ve and -ve cons.?
           if (consflag)
           {
             if (result.intValue() == 1)
@@ -347,7 +395,7 @@ public class Conservation
         }
         else
         {
-          consString.append( (gapcons[0] == 1) ? "*" : "+");
+          consString.append((gapcons[0] == 1) ? "*" : "+");
         }
       }
       else
@@ -360,8 +408,8 @@ public class Conservation
   }
 
   /**
-   *
-   *
+   * 
+   * 
    * @return Conservation sequence
    */
   public Sequence getConsSequence()
@@ -383,13 +431,13 @@ public class Conservation
     seqNums = new Vector();
     // calcSeqNum(s);
     int i = 0, iSize = sequences.length;
-    //Do we need to calculate this again?
+    // Do we need to calculate this again?
     for (i = 0; i < iSize; i++)
     {
       calcSeqNum(i);
     }
 
-    if ( (cons2 == null) || seqNumsChanged)
+    if ((cons2 == null) || seqNumsChanged)
     {
       cons2 = new int[maxLength][24];
 
@@ -424,28 +472,24 @@ public class Conservation
 
       // 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;
-                   }
-
-                   }
-                   } */
+      /*
+       * 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; }
+       *  } }
+       */
     }
   }
 
   /**
    * Calculates the quality of the set of sequences
-   *
-   * @param start Start residue
-   * @param end End residue
+   * 
+   * @param start
+   *                Start residue
+   * @param end
+   *                End residue
    */
   public void findQuality(int start, int end)
   {
@@ -454,9 +498,9 @@ public class Conservation
     double max = -10000;
     int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
 
-    //Loop over columns // JBPNote Profiling info
-    //long ts = System.currentTimeMillis();
-    //long te = System.currentTimeMillis();
+    // Loop over columns // JBPNote Profiling info
+    // long ts = System.currentTimeMillis();
+    // long te = System.currentTimeMillis();
     percentIdentity2();
 
     int size = seqNums.size();
@@ -467,7 +511,7 @@ public class Conservation
 
     for (l = 0; l < size; l++)
     {
-      lengths[l] = ( (int[]) seqNums.elementAt(l)).length - 1;
+      lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
     }
 
     for (j = start; j <= end; j++)
@@ -483,8 +527,7 @@ public class Conservation
 
         for (i2 = 0; i2 < 24; i2++)
         {
-          x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) +
-                    4);
+          x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
         }
 
         x[ii] /= size;
@@ -495,8 +538,8 @@ public class Conservation
       {
         tot = 0;
         xx = new double[24];
-        seqNum = (j < lengths[k])
-            ? ( (int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
+        seqNum = (j < lengths[k]) ? ((int[]) seqNums.elementAt(k))[j + 1]
+                : 23; // Sequence, or gap at the end
 
         // This is a loop over r
         for (i = 0; i < 23; i++)
@@ -505,8 +548,8 @@ public class Conservation
 
           sr = (double) BLOSUM62[i][seqNum] + 4;
 
-          //Calculate X with another loop over residues
-          //  System.out.println("Xi " + i + " " + x[i] + " " + sr);
+          // Calculate X with another loop over residues
+          // System.out.println("Xi " + i + " " + x[i] + " " + sr);
           xx[i] = x[i] - sr;
 
           tot += (xx[i] * xx[i]);
@@ -521,7 +564,7 @@ public class Conservation
         max = bigtot;
       }
 
-      //      bigtot  = bigtot * (size-cons2[j][23])/size;
+      // bigtot = bigtot * (size-cons2[j][23])/size;
       quality.addElement(new Double(bigtot));
 
       // Need to normalize by gaps
@@ -531,10 +574,10 @@ public class Conservation
 
     for (j = start; j <= end; j++)
     {
-      tmp = ( (Double) quality.elementAt(j)).doubleValue();
-      tmp = ( (max - tmp) * (size - cons2[j][23])) / size;
+      tmp = ((Double) quality.elementAt(j)).doubleValue();
+      tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 
-      //     System.out.println(tmp+ " " + j);
+      // System.out.println(tmp+ " " + j);
       quality.setElementAt(new Double(tmp), j);
 
       if (tmp > newmax)
@@ -543,7 +586,7 @@ public class Conservation
       }
     }
 
-    //    System.out.println("Quality " + s);
+    // System.out.println("Quality " + s);
     qualityRange[0] = new Double(0);
     qualityRange[1] = new Double(newmax);
   }