src/jalview/analysis/Conservation.java

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