src/jalview/analysis/Conservation.java

   1 /* Jalview - a java multiple alignment editor\r
   2  * Copyright (C) 1998  Michele Clamp\r
   3  *\r
   4  * This program is free software; you can redistribute it and/or\r
   5  * modify it under the terms of the GNU General Public License\r
   6  * as published by the Free Software Foundation; either version 2\r
   7  * of the License, or (at your option) any later version.\r
   8  *\r
   9  * This program is distributed in the hope that it will be useful,\r
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  12  * GNU General Public License for more details.\r
  13  *\r
  14  * You should have received a copy of the GNU General Public License\r
  15  * along with this program; if not, write to the Free Software\r
  16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.\r
  17  */\r
  18 package jalview.analysis;\r
  19 \r
  20 \r
  21 import java.util.*;\r
  22 import jalview.gui.*;\r
  23 import jalview.datamodel.*;\r
  24 \r
  25 \r
  26 public class Conservation {\r
  27   Vector sequences;\r
  28   int    start;\r
  29   int    end;\r
  30   Vector seqNums; // vector of int vectors where first is sequence checksum\r
  31   int maxLength=0; //  used by quality calcs\r
  32   boolean seqNumsChanged = false; // updated after any change via calcSeqNum;\r
  33   private void calcSeqNums() {\r
  34     for (int i=0; i<sequences.size(); i++) {\r
  35       calcSeqNum(i);\r
  36     }\r
  37   }\r
  38   private void calcSeqNum(int i) {\r
  39     String sq=null; // for dumb jbuilder not-inited exception warning\r
  40     int[] sqnum = null;\r
  41     if (i>-1 && i<sequences.size()) {\r
  42       sq = ((SequenceI)sequences.elementAt(i)).getSequence();\r
  43       if (seqNums.size()<=i) {\r
  44         seqNums.add(new int[sq.length()+1]);\r
  45       }\r
  46       if (sq.hashCode()!=((int[])seqNums.elementAt(i))[0])\r
  47        {\r
  48          int j, len;\r
  49          seqNumsChanged = true;\r
  50          sq = ((SequenceI) sequences.elementAt(i)).getSequence();\r
  51          len = sq.length();\r
  52          if (maxLength < len)\r
  53            maxLength = len;\r
  54          sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length\r
  55          sqnum[0] = sq.hashCode();\r
  56          for (j = 1; j <= len; j++)\r
  57          {\r
  58            sqnum[j] = ( (Integer) jalview.schemes.ResidueProperties.aaHash.get(new\r
  59                String(sq.substring(j - 1, j)))).intValue(); // yuk\r
  60          }\r
  61          seqNums.set(i, sqnum);\r
  62        }\r
  63      } else {\r
  64        // JBPNote INFO level debug\r
  65        System.out.println("calcSeqNum called with out of range sequence index for Alignment\n");\r
  66      }\r
  67    }\r
  68    Vector total = new Vector();\r
  69    public Vector quality;\r
  70    public Double[] qualityRange = new Double[2];\r
  71    String consString = "";\r
  72 \r
  73    Sequence consSequence;\r
  74    Hashtable        propHash;\r
  75    int              threshold;\r
  76 \r
  77   String name = "";\r
  78 \r
  79   public Conservation(String name,Hashtable propHash, int threshold, Vector sequences, int start, int end) {\r
  80     this.name      = name;\r
  81     this.propHash  = propHash;\r
  82     this.threshold = threshold;\r
  83     this.sequences = sequences;\r
  84     this.start     = start;\r
  85     this.end       = end;\r
  86     seqNums = new Vector(sequences.size());\r
  87     calcSeqNums();\r
  88   }\r
  89 \r
  90 \r
  91   public void  calculate() {\r
  92 \r
  93     for (int i = start;i <= end; i++) {\r
  94       Hashtable resultHash  = null;\r
  95       Hashtable residueHash = null;\r
  96 \r
  97       resultHash  = new Hashtable();\r
  98       residueHash = new Hashtable();\r
  99 \r
 100       for (int j=0; j < sequences.size(); j++) {\r
 101 \r
 102         if (sequences.elementAt(j) instanceof Sequence) {\r
 103           Sequence s = (Sequence)sequences.elementAt(j);\r
 104 \r
 105           if (s.getLength() > i) {\r
 106             String res = s.getSequence().substring(i,i+1);\r
 107 \r
 108             if (residueHash.containsKey(res)) {\r
 109               int count = ((Integer)residueHash.get(res)).intValue() ;\r
 110               count++;\r
 111               residueHash.put(res,new Integer(count));\r
 112             } else {\r
 113               residueHash.put(res,new Integer(1));\r
 114             }\r
 115           } else {\r
 116             if (residueHash.containsKey("-")) {\r
 117               int count = ((Integer)residueHash.get("-")).intValue() ;\r
 118               count++;\r
 119               residueHash.put("-",new Integer(count));\r
 120             } else {\r
 121               residueHash.put("-",new Integer(1));\r
 122             }\r
 123           }\r
 124         }\r
 125       }\r
 126 \r
 127       //What is the count threshold to count the residues in residueHash()\r
 128       int thresh = threshold*(sequences.size())/100;\r
 129 \r
 130       //loop over all the found residues\r
 131       Enumeration e = residueHash.keys();\r
 132 \r
 133       while (e.hasMoreElements()) {\r
 134 \r
 135         String res = (String)e.nextElement();\r
 136         if (((Integer)residueHash.get(res)).intValue() > thresh) {\r
 137 \r
 138           //Now loop over the properties\r
 139           Enumeration e2 = propHash.keys();\r
 140 \r
 141           while (e2.hasMoreElements()) {\r
 142             String    type = (String)e2.nextElement();\r
 143             Hashtable ht   = (Hashtable)propHash.get(type);\r
 144 \r
 145             //Have we ticked this before?\r
 146             if (! resultHash.containsKey(type)) {\r
 147               if (ht.containsKey(res)) {\r
 148                 resultHash.put(type,ht.get(res));\r
 149               } else {\r
 150                 resultHash.put(type,ht.get("-"));\r
 151               }\r
 152             } else if ( ((Integer)resultHash.get(type)).equals((Integer)ht.get(res)) == false) {\r
 153               resultHash.put(type,new Integer(-1));\r
 154             }\r
 155           }\r
 156         }\r
 157       }\r
 158       total.addElement(resultHash);\r
 159     }\r
 160   }\r
 161 \r
 162   public int countGaps(int j)\r
 163   {\r
 164     int count = 0;\r
 165 \r
 166     for (int i = 0; i < sequences.size();i++)\r
 167     {\r
 168       if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length())\r
 169       {  count++; continue;}\r
 170 \r
 171       char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);\r
 172       if (jalview.util.Comparison.isGap((c)))\r
 173         count++;\r
 174 \r
 175     }\r
 176     return count;\r
 177   }\r
 178 \r
 179   public  void  verdict(boolean consflag, float percentageGaps) {\r
 180     String consString = "";\r
 181 \r
 182     for (int i=start; i <= end; i++) {\r
 183       int totGaps = countGaps(i);\r
 184       float pgaps = (float)totGaps*100/(float)sequences.size();\r
 185 \r
 186       if (percentageGaps > pgaps)\r
 187       {\r
 188         Hashtable resultHash = (Hashtable)total.elementAt(i-start);\r
 189 \r
 190         //Now find the verdict\r
 191         int         count = 0;\r
 192         Enumeration e3    = resultHash.keys();\r
 193 \r
 194         while (e3.hasMoreElements())\r
 195         {\r
 196           String type    = (String)e3.nextElement();\r
 197           Integer result = (Integer)resultHash.get(type);\r
 198 \r
 199           //Do we want to count +ve conservation or +ve and -ve cons.?\r
 200 \r
 201           if (consflag)\r
 202           {\r
 203             if (result.intValue() == 1)\r
 204               count++;\r
 205           }\r
 206           else\r
 207           {\r
 208             if (result.intValue() != -1)\r
 209               count++;\r
 210           }\r
 211         }\r
 212 \r
 213         if (count < 10)\r
 214           consString = consString + String.valueOf(count);\r
 215         else\r
 216           consString = consString + "*";\r
 217 \r
 218       }\r
 219       else\r
 220       {\r
 221         consString = consString + "-";\r
 222       }\r
 223     }\r
 224 \r
 225     consSequence = new Sequence(name,consString,start,end);\r
 226   }\r
 227 \r
 228   public Sequence getConsSequence() {\r
 229     return consSequence;\r
 230   }\r
 231 \r
 232   // From Alignment.java in jalview118\r
 233 \r
 234   public void findQuality() {\r
 235     findQuality(0,maxLength-1);\r
 236   }\r
 237 \r
 238   int[][] cons2;\r
 239 \r
 240   private void percentIdentity2() {\r
 241     calcSeqNums(); // updates maxLength, too.\r
 242     if (cons2==null || seqNumsChanged) {\r
 243       cons2 = new int[maxLength][24];\r
 244       // Initialize the array\r
 245       for (int j=0;j<24;j++) {\r
 246         for (int i=0; i < maxLength;i++) {\r
 247           cons2[i][j] = 0;\r
 248         }\r
 249       }\r
 250 \r
 251       int sqnum[];\r
 252       int j = 0;\r
 253       while(j < sequences.size()) {\r
 254         sqnum=(int[])seqNums.elementAt(j);\r
 255         for (int i = 1; i < sqnum.length; i++) {\r
 256           cons2[i-1][sqnum[i]]++;\r
 257         }\r
 258         j++;\r
 259       }\r
 260 \r
 261       // unnecessary ?\r
 262       /* for (int i=start; i <= end; i++) {\r
 263            int max = -1000;\r
 264     int maxi = -1;\r
 265     int maxj = -1;\r
 266 \r
 267     for (int j=0;j<24;j++) {\r
 268         if (cons2[i][j] > max) {\r
 269         max = cons2[i][j];\r
 270         maxi = i;\r
 271         maxj = j;\r
 272       }\r
 273 \r
 274     }\r
 275   } */\r
 276     }\r
 277 \r
 278 }\r
 279 \r
 280 \r
 281 public void findQuality(int start, int end) {\r
 282     quality = new Vector();\r
 283     double max = -10000;\r
 284     String s = "";\r
 285     int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();\r
 286     //Loop over columns // JBPNote Profiling info\r
 287     //    long ts = System.currentTimeMillis();\r
 288     //long te = System.currentTimeMillis();\r
 289     percentIdentity2();\r
 290 \r
 291     int size = seqNums.size();\r
 292     int[] lengths = new int[size];\r
 293 \r
 294     for (int l = 0; l < size; l++)\r
 295       lengths[l] = ((int[]) seqNums.get(l)).length;\r
 296 \r
 297     for (int j=start; j <= end; j++) {\r
 298       double bigtot = 0;\r
 299 \r
 300       // First Xr = depends on column only\r
 301       double x[] = new double[24];\r
 302 \r
 303       for (int ii=0; ii < 24; ii++) {\r
 304         x[ii] = 0;\r
 305         try {\r
 306           for (int i2=0; i2 < 24; i2++) {\r
 307             x[ii]  += (double)cons2[j][i2] * BLOSUM62[ii][i2]+4;\r
 308           }\r
 309         } catch (Exception e) {\r
 310           System.out.println("Exception : "  + e);\r
 311         }\r
 312         //System.out.println("X " + ii + " " + x[ii]);\r
 313         x[ii] /= (size);\r
 314         //System.out.println("X " + ii + " " + x[ii]);\r
 315       }\r
 316       // Now calculate D for each position and sum\r
 317       for (int k=0; k < size; k++) {\r
 318         double tot = 0;\r
 319         double[] xx = new double[24];\r
 320         int seqNum =\r
 321             ((j+1)<lengths[size])\r
 322             ? ((int[]) seqNums.get(k))[j+1]\r
 323             : 23; // Sequence, or gap at the end\r
 324 \r
 325         // This is a loop over r\r
 326         for (int i=0; i < 23; i++) {\r
 327           double sr = 0;\r
 328           try {\r
 329             sr = (double)BLOSUM62[i][seqNum]+4;\r
 330           } catch (Exception e) {\r
 331             System.out.println("Exception in sr " + e);\r
 332           }\r
 333           //Calculate X with another loop over residues\r
 334 \r
 335           //  System.out.println("Xi " + i + " " + x[i] + " " + sr);\r
 336           xx[i] = x[i] - sr;\r
 337 \r
 338           tot += xx[i]*xx[i];\r
 339         }\r
 340         bigtot += Math.sqrt(tot);\r
 341       }\r
 342       // This is the quality for one column\r
 343       if (max < bigtot) {max = bigtot;}\r
 344       //      bigtot  = bigtot * (size-cons2[j][23])/size;\r
 345 \r
 346       quality.addElement(new Double(bigtot));\r
 347 \r
 348       s += "-";\r
 349 \r
 350       // Need to normalize by gaps\r
 351     }\r
 352     double newmax=-10000;\r
 353     for (int j=start; j <= end; j++) {\r
 354      double tmp =  ((Double)quality.elementAt(j)).doubleValue();\r
 355      tmp = (max - tmp)*(size-cons2[j][23])/size;\r
 356      //     System.out.println(tmp+ " " + j);\r
 357      quality.setElementAt(new Double(tmp),j);\r
 358      if (tmp>newmax)\r
 359        newmax = tmp;\r
 360     }\r
 361     //    System.out.println("Quality " + s);\r
 362     qualityRange[0] = new Double(0);\r
 363     qualityRange[1] = new Double(newmax);\r
 364   }\r
 365 \r
 366 \r
 367 }\r