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.addElement(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.setElementAt(sqnum, i);\r
  62        }\r
  63      } else {\r
  64        // JBPNote INFO level debug\r
  65        System.err.println("ERROR: 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         // JBPNote - have to make sure elements of the sequences vector\r
 102         //  are tested like this everywhere...\r
 103         if (sequences.elementAt(j) instanceof Sequence) {\r
 104           Sequence s = (Sequence)sequences.elementAt(j);\r
 105 \r
 106           if (s.getLength() > i) {\r
 107             String res = s.getSequence().substring(i,i+1);\r
 108 \r
 109             if (residueHash.containsKey(res)) {\r
 110               int count = ((Integer)residueHash.get(res)).intValue() ;\r
 111               count++;\r
 112               residueHash.put(res,new Integer(count));\r
 113             } else {\r
 114               residueHash.put(res,new Integer(1));\r
 115             }\r
 116           } else {\r
 117             if (residueHash.containsKey("-")) {\r
 118               int count = ((Integer)residueHash.get("-")).intValue() ;\r
 119               count++;\r
 120               residueHash.put("-",new Integer(count));\r
 121             } else {\r
 122               residueHash.put("-",new Integer(1));\r
 123             }\r
 124           }\r
 125         }\r
 126       }\r
 127 \r
 128       //What is the count threshold to count the residues in residueHash()\r
 129       int thresh = threshold*(sequences.size())/100;\r
 130 \r
 131       //loop over all the found residues\r
 132       Enumeration e = residueHash.keys();\r
 133 \r
 134       while (e.hasMoreElements()) {\r
 135 \r
 136         String res = (String)e.nextElement();\r
 137         if (((Integer)residueHash.get(res)).intValue() > thresh) {\r
 138 \r
 139           //Now loop over the properties\r
 140           Enumeration e2 = propHash.keys();\r
 141 \r
 142           while (e2.hasMoreElements()) {\r
 143             String    type = (String)e2.nextElement();\r
 144             Hashtable ht   = (Hashtable)propHash.get(type);\r
 145 \r
 146             //Have we ticked this before?\r
 147             if (! resultHash.containsKey(type)) {\r
 148               if (ht.containsKey(res)) {\r
 149                 resultHash.put(type,ht.get(res));\r
 150               } else {\r
 151                 resultHash.put(type,ht.get("-"));\r
 152               }\r
 153             } else if ( ((Integer)resultHash.get(type)).equals((Integer)ht.get(res)) == false) {\r
 154               resultHash.put(type,new Integer(-1));\r
 155             }\r
 156           }\r
 157         }\r
 158       }\r
 159 \r
 160       total.addElement(resultHash);\r
 161     }\r
 162   }\r
 163 \r
 164   public int countGaps(int j)\r
 165   {\r
 166     int count = 0;\r
 167 \r
 168     for (int i = 0; i < sequences.size();i++)\r
 169     {\r
 170       if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length())\r
 171       {  count++; continue;}\r
 172 \r
 173       char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);\r
 174       if (jalview.util.Comparison.isGap((c)))\r
 175         count++;\r
 176 \r
 177     }\r
 178     return count;\r
 179   }\r
 180   /***\r
 181    * countConsNGaps\r
 182    * returns gap count in int[0], and conserved residue count in int[1]\r
 183    */\r
 184   public int[] countConsNGaps(int j)\r
 185   {\r
 186     int count = 0;\r
 187     int cons=0;\r
 188     int nres = 0;\r
 189     int[] r = new int[2];\r
 190     char f='$';\r
 191     for (int i = 0; i < sequences.size();i++)\r
 192     {\r
 193       if( j >= ((Sequence)sequences.elementAt(i)).getSequence().length())\r
 194       {  count++;\r
 195       continue;}\r
 196 \r
 197       char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);\r
 198       if (jalview.util.Comparison.isGap((c)))\r
 199         count++;\r
 200       else {\r
 201         nres++;\r
 202         if (nres==1) {\r
 203           f = c;\r
 204           cons++;\r
 205         } else\r
 206           if (f == c) {\r
 207             cons++;\r
 208           }\r
 209       }\r
 210     }\r
 211     r[0] = (nres==cons) ? 1 : 0;\r
 212     r[1] = count;\r
 213 \r
 214 \r
 215     return r;\r
 216   }\r
 217 \r
 218   public  void  verdict(boolean consflag, float percentageGaps) {\r
 219     String consString = "";\r
 220 \r
 221     for (int i=start; i <= end; i++) {\r
 222       int[] gapcons = countConsNGaps(i);\r
 223       boolean cons = (gapcons[0]==1) ? true : false;\r
 224       int totGaps = gapcons[1];\r
 225       float pgaps = (float)totGaps*100/(float)sequences.size();\r
 226       //      System.out.println("percentage gaps = "+pgaps+"\n");\r
 227       if (percentageGaps > pgaps)\r
 228       {\r
 229         Hashtable resultHash = (Hashtable)total.elementAt(i-start);\r
 230 \r
 231         //Now find the verdict\r
 232         int         count = 0;\r
 233         Enumeration e3    = resultHash.keys();\r
 234 \r
 235         while (e3.hasMoreElements())\r
 236         {\r
 237           String type    = (String)e3.nextElement();\r
 238           Integer result = (Integer)resultHash.get(type);\r
 239 \r
 240           //Do we want to count +ve conservation or +ve and -ve cons.?\r
 241 \r
 242           if (consflag)\r
 243           {\r
 244             if (result.intValue() == 1)\r
 245               count++;\r
 246           }\r
 247           else\r
 248           {\r
 249             if (result.intValue() != -1)\r
 250               count++;\r
 251           }\r
 252         }\r
 253 \r
 254         if (count < 10)\r
 255           consString = consString + String.valueOf(count);// Conserved props!=Identity\r
 256         else\r
 257           consString = consString + ((gapcons[0]==1) ? "*" : "+");\r
 258 \r
 259       }\r
 260       else\r
 261       {\r
 262         consString = consString + "-";\r
 263       }\r
 264     }\r
 265 \r
 266     consSequence = new Sequence(name,consString,start,end);\r
 267   }\r
 268 \r
 269   public Sequence getConsSequence() {\r
 270     return consSequence;\r
 271   }\r
 272 \r
 273   // From Alignment.java in jalview118\r
 274 \r
 275   public void findQuality() {\r
 276     findQuality(0,maxLength-1);\r
 277   }\r
 278 \r
 279   int[][] cons2;\r
 280 \r
 281   private void percentIdentity2() {\r
 282     calcSeqNums(); // updates maxLength, too.\r
 283     if (cons2==null || seqNumsChanged) {\r
 284       cons2 = new int[maxLength][24];\r
 285       // Initialize the array\r
 286       for (int j=0;j<24;j++) {\r
 287         for (int i=0; i < maxLength;i++) {\r
 288           cons2[i][j] = 0;\r
 289         }\r
 290       }\r
 291 \r
 292       int sqnum[];\r
 293       int j = 0;\r
 294       while(j < sequences.size()) {\r
 295         sqnum=(int[])seqNums.elementAt(j);\r
 296         for (int i = 1; i < sqnum.length; i++) {\r
 297           cons2[i-1][sqnum[i]]++;\r
 298         }\r
 299         for (int i=sqnum.length-1; i<maxLength; i++) {\r
 300           cons2[i][23]++; // gap count\r
 301         }\r
 302         j++;\r
 303       }\r
 304 \r
 305       // unnecessary ?\r
 306       /* for (int i=start; i <= end; i++) {\r
 307            int max = -1000;\r
 308     int maxi = -1;\r
 309     int maxj = -1;\r
 310 \r
 311     for (int j=0;j<24;j++) {\r
 312         if (cons2[i][j] > max) {\r
 313         max = cons2[i][j];\r
 314         maxi = i;\r
 315         maxj = j;\r
 316       }\r
 317 \r
 318     }\r
 319   } */\r
 320     }\r
 321 \r
 322 }\r
 323 \r
 324 \r
 325 public void findQuality(int start, int end) {\r
 326     quality = new Vector();\r
 327     double max = -10000;\r
 328     String s = "";\r
 329     int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();\r
 330     //Loop over columns // JBPNote Profiling info\r
 331     //    long ts = System.currentTimeMillis();\r
 332     //long te = System.currentTimeMillis();\r
 333     percentIdentity2();\r
 334 \r
 335     int size = seqNums.size();\r
 336     int[] lengths = new int[size];\r
 337 \r
 338     for (int l = 0; l < size; l++)\r
 339       lengths[l] = ((int[]) seqNums.elementAt(l)).length-1;\r
 340 \r
 341     for (int j=start; j <= end; j++) {\r
 342       double bigtot = 0;\r
 343 \r
 344       // First Xr = depends on column only\r
 345       double x[] = new double[24];\r
 346 \r
 347       for (int ii=0; ii < 24; ii++) {\r
 348         x[ii] = 0;\r
 349         try {\r
 350           for (int i2=0; i2 < 24; i2++) {\r
 351             x[ii]  += (double)cons2[j][i2] * BLOSUM62[ii][i2]+4;\r
 352           }\r
 353         } catch (Exception e) {\r
 354           System.err.println("Exception during quality calculation.");\r
 355           e.printStackTrace();\r
 356         }\r
 357         //System.out.println("X " + ii + " " + x[ii]);\r
 358         x[ii] /= (size);\r
 359         //System.out.println("X " + ii + " " + x[ii]);\r
 360       }\r
 361       // Now calculate D for each position and sum\r
 362       for (int k=0; k < size; k++) {\r
 363         double tot = 0;\r
 364         double[] xx = new double[24];\r
 365         int seqNum =\r
 366             (j<lengths[k])\r
 367             ? ((int[]) seqNums.elementAt(k))[j+1]\r
 368             : 23; // Sequence, or gap at the end\r
 369 \r
 370         // This is a loop over r\r
 371         for (int i=0; i < 23; i++) {\r
 372           double sr = 0;\r
 373           try {\r
 374             sr = (double)BLOSUM62[i][seqNum]+4;\r
 375           } catch (Exception e) {\r
 376             System.out.println("Exception in sr: " + e);\r
 377             e.printStackTrace();\r
 378           }\r
 379           //Calculate X with another loop over residues\r
 380 \r
 381           //  System.out.println("Xi " + i + " " + x[i] + " " + sr);\r
 382           xx[i] = x[i] - sr;\r
 383 \r
 384           tot += xx[i]*xx[i];\r
 385         }\r
 386         bigtot += Math.sqrt(tot);\r
 387       }\r
 388       // This is the quality for one column\r
 389       if (max < bigtot) {max = bigtot;}\r
 390       //      bigtot  = bigtot * (size-cons2[j][23])/size;\r
 391 \r
 392       quality.addElement(new Double(bigtot));\r
 393 \r
 394       s += "-";\r
 395 \r
 396       // Need to normalize by gaps\r
 397     }\r
 398     double newmax=-10000;\r
 399     for (int j=start; j <= end; j++) {\r
 400      double tmp =  ((Double)quality.elementAt(j)).doubleValue();\r
 401      tmp = (max - tmp)*(size-cons2[j][23])/size;\r
 402      //     System.out.println(tmp+ " " + j);\r
 403      quality.setElementAt(new Double(tmp),j);\r
 404      if (tmp>newmax)\r
 405        newmax = tmp;\r
 406     }\r
 407     //    System.out.println("Quality " + s);\r
 408     qualityRange[0] = new Double(0);\r
 409     qualityRange[1] = new Double(newmax);\r
 410   }\r
 411 \r
 412 \r
 413 }\r