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