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 package jalview.analysis;\r
  20 \r
  21 import jalview.datamodel.*;\r
  22 \r
  23 import java.util.*;\r
  24 \r
  25 \r
  26 /**\r
  27  * Calculates conservation values for a given set of sequences\r
  28  *\r
  29  * @author $author$\r
  30  * @version $Revision$\r
  31  */\r
  32 public class Conservation\r
  33 {\r
  34     Vector sequences;\r
  35     int start;\r
  36     int end;\r
  37     Vector seqNums; // vector of int vectors where first is sequence checksum\r
  38     int maxLength = 0; //  used by quality calcs\r
  39     boolean seqNumsChanged = false; // updated after any change via calcSeqNum;\r
  40     Vector total = new Vector();\r
  41 \r
  42     /** Stores calculated quality values */\r
  43     public Vector quality;\r
  44 \r
  45     /** Stores maximum and minimum values of quality values  */\r
  46     public Double[] qualityRange = new Double[2];\r
  47     String consString = "";\r
  48     Sequence consSequence;\r
  49     Hashtable propHash;\r
  50     int threshold;\r
  51     String name = "";\r
  52     int[][] cons2;\r
  53 \r
  54     /**\r
  55      * Creates a new Conservation object.\r
  56      *\r
  57      * @param name Name of conservation\r
  58      * @param propHash DOCUMENT ME!\r
  59      * @param threshold to count the residues in residueHash(). commonly used value is 3\r
  60      * @param sequences sequences to be used in calculation\r
  61      * @param start start residue position\r
  62      * @param end end residue position\r
  63      */\r
  64     public Conservation(String name, Hashtable propHash, int threshold,\r
  65         Vector sequences, int start, int end)\r
  66     {\r
  67         this.name = name;\r
  68         this.propHash = propHash;\r
  69         this.threshold = threshold;\r
  70         this.sequences = sequences;\r
  71         this.start = start;\r
  72         this.end = end;\r
  73         seqNums = new Vector(sequences.size());\r
  74         calcSeqNums();\r
  75     }\r
  76 \r
  77     /**\r
  78      * DOCUMENT ME!\r
  79      */\r
  80     private void calcSeqNums()\r
  81     {\r
  82       int i=0, iSize=sequences.size();\r
  83         for (i=0; i < iSize; i++)\r
  84         {\r
  85             calcSeqNum(i);\r
  86         }\r
  87     }\r
  88 \r
  89     /**\r
  90      * DOCUMENT ME!\r
  91      *\r
  92      * @param i DOCUMENT ME!\r
  93      */\r
  94     private void calcSeqNum(int i)\r
  95     {\r
  96         String sq = null; // for dumb jbuilder not-inited exception warning\r
  97         int[] sqnum = null;\r
  98 \r
  99         if ((i > -1) && (i < sequences.size()))\r
 100         {\r
 101             sq = ((SequenceI) sequences.elementAt(i)).getSequence();\r
 102 \r
 103             if (seqNums.size() <= i)\r
 104             {\r
 105                 seqNums.addElement(new int[sq.length() + 1]);\r
 106             }\r
 107 \r
 108             if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])\r
 109             {\r
 110                 int j;\r
 111                 int len;\r
 112                 seqNumsChanged = true;\r
 113                 sq = ((SequenceI) sequences.elementAt(i)).getSequence();\r
 114                 len = sq.length();\r
 115 \r
 116                 if (maxLength < len)\r
 117                 {\r
 118                     maxLength = len;\r
 119                 }\r
 120 \r
 121                 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length\r
 122                 sqnum[0] = sq.hashCode();\r
 123 \r
 124                 for (j = 1; j <= len; j++)\r
 125                 {\r
 126                     sqnum[j] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(String.valueOf(\r
 127                                 sq.charAt(j - 1)))).intValue(); // yuk - JBPNote - case taken care of in aaHash\r
 128                 }\r
 129 \r
 130                 seqNums.setElementAt(sqnum, i);\r
 131             }\r
 132         }\r
 133         else\r
 134         {\r
 135             // JBPNote INFO level debug\r
 136             System.err.println(\r
 137                 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");\r
 138         }\r
 139     }\r
 140 \r
 141     /**\r
 142      * Calculates the conservation values for given set of sequences\r
 143      */\r
 144     public void calculate()\r
 145     {\r
 146       Hashtable resultHash, residueHash, ht;\r
 147       int count, thresh, j, jSize = sequences.size();\r
 148       String sequence, res, type;\r
 149       Enumeration enumeration, enumeration2;\r
 150 \r
 151       for (int i = start; i <= end; i++)\r
 152         {\r
 153             resultHash = new Hashtable();\r
 154             residueHash = new Hashtable();\r
 155 \r
 156             for (j = 0; j < jSize; j++)\r
 157             {\r
 158                 // JBPNote - have to make sure elements of the sequences vector\r
 159                 //  are tested like this everywhere...\r
 160                 if (sequences.elementAt(j) instanceof Sequence)\r
 161                 {\r
 162                     sequence = ((Sequence) sequences.elementAt(j)).getSequence();\r
 163 \r
 164                     if (sequence.length() > i)\r
 165                     {\r
 166                         res = String.valueOf(Character.toUpperCase(sequence.charAt(i)));\r
 167 \r
 168                         if (residueHash.containsKey(res))\r
 169                         {\r
 170                             count = ((Integer) residueHash.get(res)).intValue();\r
 171                             count++;\r
 172                             residueHash.put(res, new Integer(count));\r
 173                         }\r
 174                         else\r
 175                         {\r
 176                             residueHash.put(res, new Integer(1));\r
 177                         }\r
 178                     }\r
 179                     else\r
 180                     {\r
 181                         if (residueHash.containsKey("-"))\r
 182                         {\r
 183                             count = ((Integer) residueHash.get("-")).intValue();\r
 184                             count++;\r
 185                             residueHash.put("-", new Integer(count));\r
 186                         }\r
 187                         else\r
 188                         {\r
 189                             residueHash.put("-", new Integer(1));\r
 190                         }\r
 191                     }\r
 192                 }\r
 193             }\r
 194 \r
 195             //What is the count threshold to count the residues in residueHash()\r
 196             thresh = (threshold * (sequences.size())) / 100;\r
 197 \r
 198             //loop over all the found residues\r
 199             enumeration = residueHash.keys();\r
 200 \r
 201             while (enumeration.hasMoreElements())\r
 202             {\r
 203                 res = (String) enumeration.nextElement();\r
 204 \r
 205                 if (((Integer) residueHash.get(res)).intValue() > thresh)\r
 206                 {\r
 207                     //Now loop over the properties\r
 208                     enumeration2 = propHash.keys();\r
 209 \r
 210                     while (enumeration2.hasMoreElements())\r
 211                     {\r
 212                         type = (String) enumeration2.nextElement();\r
 213                         ht = (Hashtable) propHash.get(type);\r
 214 \r
 215                         //Have we ticked this before?\r
 216                         if (!resultHash.containsKey(type))\r
 217                         {\r
 218                             if (ht.containsKey(res))\r
 219                             {\r
 220                                 resultHash.put(type, ht.get(res));\r
 221                             }\r
 222                             else\r
 223                             {\r
 224                                 resultHash.put(type, ht.get("-"));\r
 225                             }\r
 226                         }\r
 227                         else if (((Integer) resultHash.get(type)).equals(\r
 228                                     (Integer) ht.get(res)) == false)\r
 229                         {\r
 230                             resultHash.put(type, new Integer(-1));\r
 231                         }\r
 232                     }\r
 233                 }\r
 234             }\r
 235 \r
 236             total.addElement(resultHash);\r
 237         }\r
 238     }\r
 239 \r
 240 \r
 241     /***\r
 242      * countConsNGaps\r
 243      * returns gap count in int[0], and conserved residue count in int[1]\r
 244      */\r
 245     public int[] countConsNGaps(int j)\r
 246     {\r
 247         int count = 0;\r
 248         int cons = 0;\r
 249         int nres = 0;\r
 250         int[] r = new int[2];\r
 251         char f = '$';\r
 252         int i, iSize = sequences.size();\r
 253         char c;\r
 254 \r
 255         for (i = 0; i < iSize; i++)\r
 256         {\r
 257             if (j >= ((Sequence) sequences.elementAt(i)).getLength())\r
 258             {\r
 259                 count++;\r
 260                 continue;\r
 261             }\r
 262 \r
 263             c = ((Sequence) sequences.elementAt(i)).getCharAt(j); // gaps do not have upper/lower case\r
 264 \r
 265             if (jalview.util.Comparison.isGap((c)))\r
 266             {\r
 267                 count++;\r
 268             }\r
 269             else\r
 270             {\r
 271                 nres++;\r
 272 \r
 273                 if (nres == 1)\r
 274                 {\r
 275                     f = c;\r
 276                     cons++;\r
 277                 }\r
 278                 else if (f == c)\r
 279                 {\r
 280                     cons++;\r
 281                 }\r
 282             }\r
 283         }\r
 284 \r
 285         r[0] = (nres == cons) ? 1 : 0;\r
 286         r[1] = count;\r
 287 \r
 288         return r;\r
 289     }\r
 290 \r
 291     /**\r
 292      * Calculates the conservation sequence\r
 293      *\r
 294      * @param consflag if true, poitiveve conservation; false calculates negative conservation\r
 295      * @param percentageGaps commonly used value is 25\r
 296      */\r
 297     public void verdict(boolean consflag, float percentageGaps)\r
 298     {\r
 299         StringBuffer consString = new StringBuffer();\r
 300         String type;\r
 301         Integer result;\r
 302         int[] gapcons;\r
 303         int totGaps, count;\r
 304         float pgaps;\r
 305         Hashtable resultHash ;\r
 306         Enumeration enumeration;\r
 307 \r
 308 \r
 309         for (int i = start; i <= end; i++)\r
 310         {\r
 311             gapcons = countConsNGaps(i);\r
 312             totGaps = gapcons[1];\r
 313             pgaps = ((float) totGaps * 100) / (float) sequences.size();\r
 314 \r
 315             if (percentageGaps > pgaps)\r
 316             {\r
 317                 resultHash = (Hashtable) total.elementAt(i - start);\r
 318 \r
 319                 //Now find the verdict\r
 320                 count = 0;\r
 321                 enumeration = resultHash.keys();\r
 322 \r
 323                 while (enumeration.hasMoreElements())\r
 324                 {\r
 325                    type = (String) enumeration.nextElement();\r
 326                    result = (Integer) resultHash.get(type);\r
 327 \r
 328                     //Do we want to count +ve conservation or +ve and -ve cons.?\r
 329                     if (consflag)\r
 330                     {\r
 331                         if (result.intValue() == 1)\r
 332                         {\r
 333                             count++;\r
 334                         }\r
 335                     }\r
 336                     else\r
 337                     {\r
 338                         if (result.intValue() != -1)\r
 339                         {\r
 340                             count++;\r
 341                         }\r
 342                     }\r
 343                 }\r
 344 \r
 345                 if (count < 10)\r
 346                 {\r
 347                     consString.append(count); // Conserved props!=Identity\r
 348                 }\r
 349                 else\r
 350                 {\r
 351                     consString.append((gapcons[0] == 1) ? "*" : "+");\r
 352                 }\r
 353             }\r
 354             else\r
 355             {\r
 356                 consString.append("-");\r
 357             }\r
 358         }\r
 359 \r
 360         consSequence = new Sequence(name, consString.toString(), start, end);\r
 361     }\r
 362 \r
 363     /**\r
 364      *\r
 365      *\r
 366      * @return Conservation sequence\r
 367      */\r
 368     public Sequence getConsSequence()\r
 369     {\r
 370         return consSequence;\r
 371     }\r
 372 \r
 373     // From Alignment.java in jalview118\r
 374     public void findQuality()\r
 375     {\r
 376         findQuality(0, maxLength - 1);\r
 377     }\r
 378 \r
 379     /**\r
 380      * DOCUMENT ME!\r
 381      */\r
 382     private void percentIdentity2()\r
 383     {\r
 384         calcSeqNums(); // updates maxLength, too.\r
 385 \r
 386         if ((cons2 == null) || seqNumsChanged)\r
 387         {\r
 388             cons2 = new int[maxLength][24];\r
 389 \r
 390             // Initialize the array\r
 391             for (int j = 0; j < 24; j++)\r
 392             {\r
 393                 for (int i = 0; i < maxLength; i++)\r
 394                 {\r
 395                     cons2[i][j] = 0;\r
 396                 }\r
 397             }\r
 398 \r
 399             int[] sqnum;\r
 400             int j = 0;\r
 401 \r
 402             while (j < sequences.size())\r
 403             {\r
 404                 sqnum = (int[]) seqNums.elementAt(j);\r
 405 \r
 406                 for (int i = 1; i < sqnum.length; i++)\r
 407                 {\r
 408                     cons2[i - 1][sqnum[i]]++;\r
 409                 }\r
 410 \r
 411                 for (int i = sqnum.length - 1; i < maxLength; i++)\r
 412                 {\r
 413                     cons2[i][23]++; // gap count\r
 414                 }\r
 415 \r
 416                 j++;\r
 417             }\r
 418 \r
 419             // unnecessary ?\r
 420 \r
 421             /* for (int i=start; i <= end; i++) {\r
 422                  int max = -1000;\r
 423             int maxi = -1;\r
 424             int maxj = -1;\r
 425 \r
 426             for (int j=0;j<24;j++) {\r
 427               if (cons2[i][j] > max) {\r
 428               max = cons2[i][j];\r
 429               maxi = i;\r
 430               maxj = j;\r
 431             }\r
 432 \r
 433             }\r
 434             } */\r
 435         }\r
 436     }\r
 437 \r
 438     /**\r
 439      * Calculates the quality of the set of sequences\r
 440      *\r
 441      * @param start Start residue\r
 442      * @param end End residue\r
 443      */\r
 444     public void findQuality(int start, int end)\r
 445     {\r
 446         quality = new Vector();\r
 447 \r
 448         double max = -10000;\r
 449         int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();\r
 450 \r
 451         //Loop over columns // JBPNote Profiling info\r
 452         //    long ts = System.currentTimeMillis();\r
 453         //long te = System.currentTimeMillis();\r
 454         percentIdentity2();\r
 455 \r
 456         int size = seqNums.size();\r
 457         int[] lengths = new int[size];\r
 458         double tot, bigtot, sr, tmp;\r
 459         double [] x, xx;\r
 460         int l, j, i, ii, seqNum;\r
 461 \r
 462         for (l = 0; l < size; l++)\r
 463             lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;\r
 464 \r
 465 \r
 466         for (j = start; j <= end; j++)\r
 467         {\r
 468             bigtot = 0;\r
 469 \r
 470             // First Xr = depends on column only\r
 471             x = new double[24];\r
 472 \r
 473             for (ii = 0; ii < 24; ii++)\r
 474             {\r
 475                 x[ii] = 0;\r
 476 \r
 477                 try\r
 478                 {\r
 479                     for (int i2 = 0; i2 < 24; i2++)\r
 480                     {\r
 481                         x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) +\r
 482                         4);\r
 483                     }\r
 484                 }\r
 485                 catch (Exception e)\r
 486                 {\r
 487                     System.err.println("Exception during quality calculation.");\r
 488                     e.printStackTrace();\r
 489                 }\r
 490 \r
 491                 //System.out.println("X " + ii + " " + x[ii]);\r
 492                 x[ii] /= (size);\r
 493 \r
 494                 //System.out.println("X " + ii + " " + x[ii]);\r
 495             }\r
 496 \r
 497             // Now calculate D for each position and sum\r
 498             for (int k = 0; k < size; k++)\r
 499             {\r
 500                 tot = 0;\r
 501                 xx = new double[24];\r
 502                 seqNum = (j < lengths[k])\r
 503                     ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end\r
 504 \r
 505                 // This is a loop over r\r
 506                 for (i = 0; i < 23; i++)\r
 507                 {\r
 508                     sr = 0;\r
 509 \r
 510                     try\r
 511                     {\r
 512                         sr = (double) BLOSUM62[i][seqNum] + 4;\r
 513                     }\r
 514                     catch (Exception e)\r
 515                     {\r
 516                         System.out.println("Exception in sr: " + e);\r
 517                         e.printStackTrace();\r
 518                     }\r
 519 \r
 520                     //Calculate X with another loop over residues\r
 521                     //  System.out.println("Xi " + i + " " + x[i] + " " + sr);\r
 522                     xx[i] = x[i] - sr;\r
 523 \r
 524                     tot += (xx[i] * xx[i]);\r
 525                 }\r
 526 \r
 527                 bigtot += Math.sqrt(tot);\r
 528             }\r
 529 \r
 530             // This is the quality for one column\r
 531             if (max < bigtot)\r
 532             {\r
 533                 max = bigtot;\r
 534             }\r
 535 \r
 536             //      bigtot  = bigtot * (size-cons2[j][23])/size;\r
 537             quality.addElement(new Double(bigtot));\r
 538 \r
 539 \r
 540             // Need to normalize by gaps\r
 541         }\r
 542 \r
 543         double newmax = -10000;\r
 544 \r
 545         for (j = start; j <= end; j++)\r
 546         {\r
 547             tmp = ((Double) quality.elementAt(j)).doubleValue();\r
 548             tmp = ((max - tmp) * (size - cons2[j][23])) / size;\r
 549 \r
 550             //     System.out.println(tmp+ " " + j);\r
 551             quality.setElementAt(new Double(tmp), j);\r
 552 \r
 553             if (tmp > newmax)\r
 554             {\r
 555                 newmax = tmp;\r
 556             }\r
 557         }\r
 558 \r
 559         //    System.out.println("Quality " + s);\r
 560         qualityRange[0] = new Double(0);\r
 561         qualityRange[1] = new Double(newmax);\r
 562     }\r
 563 }\r