MaxResidue must be at least 2
[jalview.git] / src / jalview / analysis / Conservation.java
index 0fe4f5e..d1b508b 100755 (executable)
@@ -1,6 +1,6 @@
 /*\r
 * Jalview - A Sequence Alignment Editor and Viewer\r
-* Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
+* Copyright (C) 2006 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
 *\r
 * This program is free software; you can redistribute it and/or\r
 * modify it under the terms of the GNU General Public License\r
@@ -31,13 +31,13 @@ import java.util.*;
  */\r
 public class Conservation\r
 {\r
-    Vector sequences;\r
+    SequenceI [] sequences;\r
     int start;\r
     int end;\r
     Vector seqNums; // vector of int vectors where first is sequence checksum\r
     int maxLength = 0; //  used by quality calcs\r
     boolean seqNumsChanged = false; // updated after any change via calcSeqNum;\r
-    Vector total = new Vector();\r
+    Hashtable [] total;\r
 \r
     /** Stores calculated quality values */\r
     public Vector quality;\r
@@ -64,27 +64,27 @@ public class Conservation
     public Conservation(String name, Hashtable propHash, int threshold,\r
         Vector sequences, int start, int end)\r
     {\r
+\r
         this.name = name;\r
         this.propHash = propHash;\r
         this.threshold = threshold;\r
-        this.sequences = sequences;\r
         this.start = start;\r
         this.end = end;\r
-        seqNums = new Vector(sequences.size());\r
-        calcSeqNums();\r
-    }\r
 \r
-    /**\r
-     * DOCUMENT ME!\r
-     */\r
-    private void calcSeqNums()\r
-    {\r
-        for (int i = 0; i < sequences.size(); i++)\r
+\r
+        int s, sSize = sequences.size();\r
+        SequenceI[] sarray = new SequenceI[sSize];\r
+        this.sequences = sarray;\r
+\r
+        for (s = 0; s < sSize; s++)\r
         {\r
-            calcSeqNum(i);\r
+          sarray[s] = (SequenceI) sequences.elementAt(s);\r
+          if(sarray[s].getLength()>maxLength)\r
+            maxLength = sarray[s].getLength();\r
         }\r
     }\r
 \r
+\r
     /**\r
      * DOCUMENT ME!\r
      *\r
@@ -95,9 +95,11 @@ public class Conservation
         String sq = null; // for dumb jbuilder not-inited exception warning\r
         int[] sqnum = null;\r
 \r
-        if ((i > -1) && (i < sequences.size()))\r
+        int sSize = sequences.length;\r
+\r
+        if ((i > -1) && (i < sSize))\r
         {\r
-            sq = ((SequenceI) sequences.elementAt(i)).getSequence();\r
+            sq = sequences[i].getSequence();\r
 \r
             if (seqNums.size() <= i)\r
             {\r
@@ -109,7 +111,6 @@ public class Conservation
                 int j;\r
                 int len;\r
                 seqNumsChanged = true;\r
-                sq = ((SequenceI) sequences.elementAt(i)).getSequence();\r
                 len = sq.length();\r
 \r
                 if (maxLength < len)\r
@@ -122,12 +123,14 @@ public class Conservation
 \r
                 for (j = 1; j <= len; j++)\r
                 {\r
-                    sqnum[j] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(new String(\r
-                                sq.substring(j - 1, j)))).intValue(); // yuk\r
+                    sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq.charAt(j-1)];\r
                 }\r
 \r
+\r
                 seqNums.setElementAt(sqnum, i);\r
             }\r
+            else\r
+              System.out.println("NEVER THE EXCEPTION");\r
         }\r
         else\r
         {\r
@@ -142,72 +145,61 @@ public class Conservation
      */\r
     public void calculate()\r
     {\r
-        for (int i = start; i <= end; i++)\r
-        {\r
-            Hashtable resultHash = null;\r
-            Hashtable residueHash = null;\r
+      Hashtable resultHash, ht;\r
+      int thresh, j, jSize = sequences.length;\r
+      int[] values; // Replaces residueHash\r
+      String type, res=null;\r
+      char c;\r
+      Enumeration enumeration2;\r
 \r
-            resultHash = new Hashtable();\r
-            residueHash = new Hashtable();\r
+      total = new Hashtable[maxLength];\r
+\r
+      for (int i = start; i <= end; i++)\r
+        {\r
+            values = new int[132];\r
 \r
-            for (int j = 0; j < sequences.size(); j++)\r
+            for (j = 0; j < jSize; j++)\r
             {\r
-                // JBPNote - have to make sure elements of the sequences vector\r
-                //  are tested like this everywhere...\r
-                if (sequences.elementAt(j) instanceof Sequence)\r
-                {\r
-                    Sequence s = (Sequence) sequences.elementAt(j);\r
+              if (sequences[j].getLength() > i)\r
+              {\r
+                c = sequences[j].getCharAt(i);\r
 \r
-                    if (s.getLength() > i)\r
-                    {\r
-                        String res = s.getSequence().substring(i, i + 1);\r
+                // No need to check if its a '-'\r
+                if (c == '.' || c == ' ')\r
+                  c = '-';\r
 \r
-                        if (residueHash.containsKey(res))\r
-                        {\r
-                            int count = ((Integer) residueHash.get(res)).intValue();\r
-                            count++;\r
-                            residueHash.put(res, new Integer(count));\r
-                        }\r
-                        else\r
-                        {\r
-                            residueHash.put(res, new Integer(1));\r
-                        }\r
-                    }\r
-                    else\r
-                    {\r
-                        if (residueHash.containsKey("-"))\r
-                        {\r
-                            int count = ((Integer) residueHash.get("-")).intValue();\r
-                            count++;\r
-                            residueHash.put("-", new Integer(count));\r
-                        }\r
-                        else\r
-                        {\r
-                            residueHash.put("-", new Integer(1));\r
-                        }\r
-                    }\r
+                if ('a' <= c && c <= 'z')\r
+                {\r
+                  c -= (32);// 32 = 'a' - 'A'\r
                 }\r
+\r
+                values[c]++;\r
+              }\r
+              else\r
+              {\r
+                values['-']++;\r
+              }\r
             }\r
 \r
             //What is the count threshold to count the residues in residueHash()\r
-            int thresh = (threshold * (sequences.size())) / 100;\r
+            thresh = (threshold * (jSize)) / 100;\r
 \r
             //loop over all the found residues\r
-            Enumeration e = residueHash.keys();\r
-\r
-            while (e.hasMoreElements())\r
+            resultHash = new Hashtable();\r
+            for (int v = '-'; v < 'Z'; v++)\r
             {\r
-                String res = (String) e.nextElement();\r
 \r
-                if (((Integer) residueHash.get(res)).intValue() > thresh)\r
+                if (values[v] > thresh)\r
                 {\r
+                  res =  String.valueOf( (char) v);\r
+\r
                     //Now loop over the properties\r
-                    Enumeration e2 = propHash.keys();\r
+                    enumeration2 = propHash.keys();\r
 \r
-                    while (e2.hasMoreElements())\r
+                    while (enumeration2.hasMoreElements())\r
                     {\r
-                        String type = (String) e2.nextElement();\r
-                        Hashtable ht = (Hashtable) propHash.get(type);\r
+                        type = (String) enumeration2.nextElement();\r
+                        ht = (Hashtable) propHash.get(type);\r
 \r
                         //Have we ticked this before?\r
                         if (!resultHash.containsKey(type))\r
@@ -230,7 +222,7 @@ public class Conservation
                 }\r
             }\r
 \r
-            total.addElement(resultHash);\r
+            total[i] = resultHash;\r
         }\r
     }\r
 \r
@@ -246,17 +238,18 @@ public class Conservation
         int nres = 0;\r
         int[] r = new int[2];\r
         char f = '$';\r
+        int i, iSize = sequences.length;\r
+        char c;\r
 \r
-        for (int i = 0; i < sequences.size(); i++)\r
+        for (i = 0; i < iSize; i++)\r
         {\r
-            if (j >= ((Sequence) sequences.elementAt(i)).getSequence().length())\r
+            if (j >= sequences[i].getLength())\r
             {\r
                 count++;\r
-\r
                 continue;\r
             }\r
 \r
-            char c = ((Sequence) sequences.elementAt(i)).getSequence().charAt(j);\r
+            c = sequences[i].getCharAt(j); // gaps do not have upper/lower case\r
 \r
             if (jalview.util.Comparison.isGap((c)))\r
             {\r
@@ -292,27 +285,34 @@ public class Conservation
      */\r
     public void verdict(boolean consflag, float percentageGaps)\r
     {\r
-        String consString = "";\r
+        StringBuffer consString = new StringBuffer();\r
+        String type;\r
+        Integer result;\r
+        int[] gapcons;\r
+        int totGaps, count;\r
+        float pgaps;\r
+        Hashtable resultHash ;\r
+        Enumeration enumeration;\r
+\r
 \r
         for (int i = start; i <= end; i++)\r
         {\r
-            int[] gapcons = countConsNGaps(i);\r
-            int totGaps = gapcons[1];\r
-            float pgaps = ((float) totGaps * 100) / (float) sequences.size();\r
+            gapcons = countConsNGaps(i);\r
+            totGaps = gapcons[1];\r
+            pgaps = ((float) totGaps * 100) / (float) sequences.length;\r
 \r
-            //      System.out.println("percentage gaps = "+pgaps+"\n");\r
             if (percentageGaps > pgaps)\r
             {\r
-                Hashtable resultHash = (Hashtable) total.elementAt(i - start);\r
+                resultHash =  total[i - start];\r
 \r
                 //Now find the verdict\r
-                int count = 0;\r
-                Enumeration e3 = resultHash.keys();\r
+                count = 0;\r
+                enumeration = resultHash.keys();\r
 \r
-                while (e3.hasMoreElements())\r
+                while (enumeration.hasMoreElements())\r
                 {\r
-                    String type = (String) e3.nextElement();\r
-                    Integer result = (Integer) resultHash.get(type);\r
+                   type = (String) enumeration.nextElement();\r
+                   result = (Integer) resultHash.get(type);\r
 \r
                     //Do we want to count +ve conservation or +ve and -ve cons.?\r
                     if (consflag)\r
@@ -333,20 +333,20 @@ public class Conservation
 \r
                 if (count < 10)\r
                 {\r
-                    consString = consString + String.valueOf(count); // Conserved props!=Identity\r
+                    consString.append(count); // Conserved props!=Identity\r
                 }\r
                 else\r
                 {\r
-                    consString = consString + ((gapcons[0] == 1) ? "*" : "+");\r
+                    consString.append((gapcons[0] == 1) ? "*" : "+");\r
                 }\r
             }\r
             else\r
             {\r
-                consString = consString + "-";\r
+                consString.append("-");\r
             }\r
         }\r
 \r
-        consSequence = new Sequence(name, consString, start, end);\r
+        consSequence = new Sequence(name, consString.toString(), start, end);\r
     }\r
 \r
     /**\r
@@ -370,7 +370,15 @@ public class Conservation
      */\r
     private void percentIdentity2()\r
     {\r
-        calcSeqNums(); // updates maxLength, too.\r
+      seqNums = new Vector();\r
+     // calcSeqNum(s);\r
+      int i = 0, iSize = sequences.length;\r
+    //Do we need to calculate this again?\r
+      for (i = 0; i < iSize; i++)\r
+      {\r
+       calcSeqNum(i);\r
+      }\r
+\r
 \r
         if ((cons2 == null) || seqNumsChanged)\r
         {\r
@@ -379,7 +387,7 @@ public class Conservation
             // Initialize the array\r
             for (int j = 0; j < 24; j++)\r
             {\r
-                for (int i = 0; i < maxLength; i++)\r
+                for (i = 0; i < maxLength; i++)\r
                 {\r
                     cons2[i][j] = 0;\r
                 }\r
@@ -388,16 +396,16 @@ public class Conservation
             int[] sqnum;\r
             int j = 0;\r
 \r
-            while (j < sequences.size())\r
+            while (j < sequences.length)\r
             {\r
                 sqnum = (int[]) seqNums.elementAt(j);\r
 \r
-                for (int i = 1; i < sqnum.length; i++)\r
+                for (i = 1; i < sqnum.length; i++)\r
                 {\r
                     cons2[i - 1][sqnum[i]]++;\r
                 }\r
 \r
-                for (int i = sqnum.length - 1; i < maxLength; i++)\r
+                for (i = sqnum.length - 1; i < maxLength; i++)\r
                 {\r
                     cons2[i][23]++; // gap count\r
                 }\r
@@ -435,73 +443,56 @@ public class Conservation
         quality = new Vector();\r
 \r
         double max = -10000;\r
-        String s = "";\r
         int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();\r
 \r
         //Loop over columns // JBPNote Profiling info\r
-        //    long ts = System.currentTimeMillis();\r
+        //long ts = System.currentTimeMillis();\r
         //long te = System.currentTimeMillis();\r
         percentIdentity2();\r
 \r
         int size = seqNums.size();\r
         int[] lengths = new int[size];\r
+        double tot, bigtot, sr, tmp;\r
+        double [] x, xx;\r
+        int l, j, i, ii, i2, k, seqNum;\r
 \r
-        for (int l = 0; l < size; l++)\r
+        for (l = 0; l < size; l++)\r
             lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;\r
 \r
-        for (int j = start; j <= end; j++)\r
+        for (j = start; j <= end; j++)\r
         {\r
-            double bigtot = 0;\r
+            bigtot = 0;\r
 \r
             // First Xr = depends on column only\r
-            double[] x = new double[24];\r
+            x = new double[24];\r
 \r
-            for (int ii = 0; ii < 24; ii++)\r
+            for (ii = 0; ii < 24; ii++)\r
             {\r
                 x[ii] = 0;\r
 \r
-                try\r
-                {\r
-                    for (int i2 = 0; i2 < 24; i2++)\r
-                    {\r
-                        x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) +\r
-                        4);\r
-                    }\r
-                }\r
-                catch (Exception e)\r
+                for (i2 = 0; i2 < 24; i2++)\r
                 {\r
-                    System.err.println("Exception during quality calculation.");\r
-                    e.printStackTrace();\r
+                  x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) +\r
+                            4);\r
                 }\r
 \r
-                //System.out.println("X " + ii + " " + x[ii]);\r
-                x[ii] /= (size);\r
-\r
-                //System.out.println("X " + ii + " " + x[ii]);\r
+                x[ii] /= size;\r
             }\r
 \r
             // Now calculate D for each position and sum\r
-            for (int k = 0; k < size; k++)\r
+            for (k = 0; k < size; k++)\r
             {\r
-                double tot = 0;\r
-                double[] xx = new double[24];\r
-                int seqNum = (j < lengths[k])\r
+                tot = 0;\r
+                xx = new double[24];\r
+                seqNum = (j < lengths[k])\r
                     ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end\r
 \r
                 // This is a loop over r\r
-                for (int i = 0; i < 23; i++)\r
+                for (i = 0; i < 23; i++)\r
                 {\r
-                    double sr = 0;\r
+                    sr = 0;\r
 \r
-                    try\r
-                    {\r
-                        sr = (double) BLOSUM62[i][seqNum] + 4;\r
-                    }\r
-                    catch (Exception e)\r
-                    {\r
-                        System.out.println("Exception in sr: " + e);\r
-                        e.printStackTrace();\r
-                    }\r
+                    sr = (double) BLOSUM62[i][seqNum] + 4;\r
 \r
                     //Calculate X with another loop over residues\r
                     //  System.out.println("Xi " + i + " " + x[i] + " " + sr);\r
@@ -522,16 +513,14 @@ public class Conservation
             //      bigtot  = bigtot * (size-cons2[j][23])/size;\r
             quality.addElement(new Double(bigtot));\r
 \r
-            s += "-";\r
-\r
             // Need to normalize by gaps\r
         }\r
 \r
         double newmax = -10000;\r
 \r
-        for (int j = start; j <= end; j++)\r
+        for (j = start; j <= end; j++)\r
         {\r
-            double tmp = ((Double) quality.elementAt(j)).doubleValue();\r
+            tmp = ((Double) quality.elementAt(j)).doubleValue();\r
             tmp = ((max - tmp) * (size - cons2[j][23])) / size;\r
 \r
             //     System.out.println(tmp+ " " + j);\r