src/jalview/analysis/Conservation.java

   1 /*\r
   2 * Jalview - A Sequence Alignment Editor and Viewer\r
   3 * Copyright (C) 2006 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     SequenceI [] 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     Hashtable [] total;\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 \r
  68         this.name = name;\r
  69         this.propHash = propHash;\r
  70         this.threshold = threshold;\r
  71         this.start = start;\r
  72         this.end = end;\r
  73 \r
  74 \r
  75         int s, sSize = sequences.size();\r
  76         SequenceI[] sarray = new SequenceI[sSize];\r
  77         this.sequences = sarray;\r
  78 \r
  79         for (s = 0; s < sSize; s++)\r
  80         {\r
  81           sarray[s] = (SequenceI) sequences.elementAt(s);\r
  82           if(sarray[s].getLength()>maxLength)\r
  83             maxLength = sarray[s].getLength();\r
  84         }\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         int sSize = sequences.length;\r
  99 \r
 100         if ((i > -1) && (i < sSize))\r
 101         {\r
 102             sq = sequences[i].getSequence();\r
 103 \r
 104             if (seqNums.size() <= i)\r
 105             {\r
 106                 seqNums.addElement(new int[sq.length() + 1]);\r
 107             }\r
 108 \r
 109             if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])\r
 110             {\r
 111                 int j;\r
 112                 int len;\r
 113                 seqNumsChanged = true;\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] = jalview.schemes.ResidueProperties.aaIndex[sq.charAt(j-1)];\r
 127                 }\r
 128 \r
 129 \r
 130                 seqNums.setElementAt(sqnum, i);\r
 131             }\r
 132             else\r
 133               System.out.println("NEVER THE EXCEPTION");\r
 134         }\r
 135         else\r
 136         {\r
 137             // JBPNote INFO level debug\r
 138             System.err.println(\r
 139                 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");\r
 140         }\r
 141     }\r
 142 \r
 143     /**\r
 144      * Calculates the conservation values for given set of sequences\r
 145      */\r
 146     public void calculate()\r
 147     {\r
 148       Hashtable resultHash, ht;\r
 149       int thresh, j, jSize = sequences.length;\r
 150       int[] values; // Replaces residueHash\r
 151       String type, res=null;\r
 152       char c;\r
 153       Enumeration enumeration2;\r
 154 \r
 155       total = new Hashtable[maxLength];\r
 156 \r
 157       for (int i = start; i <= end; i++)\r
 158         {\r
 159             values = new int[132];\r
 160 \r
 161             for (j = 0; j < jSize; j++)\r
 162             {\r
 163               if (sequences[j].getLength() > i)\r
 164               {\r
 165                 c = sequences[j].getCharAt(i);\r
 166 \r
 167                 // No need to check if its a '-'\r
 168                 if (c == '.' || c == ' ')\r
 169                   c = '-';\r
 170 \r
 171                 if ('a' <= c && c <= 'z')\r
 172                 {\r
 173                   c -= (32);// 32 = 'a' - 'A'\r
 174                 }\r
 175 \r
 176                 values[c]++;\r
 177               }\r
 178               else\r
 179               {\r
 180                 values['-']++;\r
 181               }\r
 182             }\r
 183 \r
 184             //What is the count threshold to count the residues in residueHash()\r
 185             thresh = (threshold * (jSize)) / 100;\r
 186 \r
 187             //loop over all the found residues\r
 188             resultHash = new Hashtable();\r
 189             for (int v = '-'; v < 'Z'; v++)\r
 190             {\r
 191 \r
 192                 if (values[v] > thresh)\r
 193                 {\r
 194                   res =  String.valueOf( (char) v);\r
 195 \r
 196                     //Now loop over the properties\r
 197                     enumeration2 = propHash.keys();\r
 198 \r
 199                     while (enumeration2.hasMoreElements())\r
 200                     {\r
 201                         type = (String) enumeration2.nextElement();\r
 202                         ht = (Hashtable) propHash.get(type);\r
 203 \r
 204                         //Have we ticked this before?\r
 205                         if (!resultHash.containsKey(type))\r
 206                         {\r
 207                             if (ht.containsKey(res))\r
 208                             {\r
 209                                 resultHash.put(type, ht.get(res));\r
 210                             }\r
 211                             else\r
 212                             {\r
 213                                 resultHash.put(type, ht.get("-"));\r
 214                             }\r
 215                         }\r
 216                         else if (((Integer) resultHash.get(type)).equals(\r
 217                                     (Integer) ht.get(res)) == false)\r
 218                         {\r
 219                             resultHash.put(type, new Integer(-1));\r
 220                         }\r
 221                     }\r
 222                 }\r
 223             }\r
 224 \r
 225             total[i] = resultHash;\r
 226         }\r
 227     }\r
 228 \r
 229 \r
 230     /***\r
 231      * countConsNGaps\r
 232      * returns gap count in int[0], and conserved residue count in int[1]\r
 233      */\r
 234     public int[] countConsNGaps(int j)\r
 235     {\r
 236         int count = 0;\r
 237         int cons = 0;\r
 238         int nres = 0;\r
 239         int[] r = new int[2];\r
 240         char f = '$';\r
 241         int i, iSize = sequences.length;\r
 242         char c;\r
 243 \r
 244         for (i = 0; i < iSize; i++)\r
 245         {\r
 246             if (j >= sequences[i].getLength())\r
 247             {\r
 248                 count++;\r
 249                 continue;\r
 250             }\r
 251 \r
 252             c = sequences[i].getCharAt(j); // gaps do not have upper/lower case\r
 253 \r
 254             if (jalview.util.Comparison.isGap((c)))\r
 255             {\r
 256                 count++;\r
 257             }\r
 258             else\r
 259             {\r
 260                 nres++;\r
 261 \r
 262                 if (nres == 1)\r
 263                 {\r
 264                     f = c;\r
 265                     cons++;\r
 266                 }\r
 267                 else if (f == c)\r
 268                 {\r
 269                     cons++;\r
 270                 }\r
 271             }\r
 272         }\r
 273 \r
 274         r[0] = (nres == cons) ? 1 : 0;\r
 275         r[1] = count;\r
 276 \r
 277         return r;\r
 278     }\r
 279 \r
 280     /**\r
 281      * Calculates the conservation sequence\r
 282      *\r
 283      * @param consflag if true, poitiveve conservation; false calculates negative conservation\r
 284      * @param percentageGaps commonly used value is 25\r
 285      */\r
 286     public void verdict(boolean consflag, float percentageGaps)\r
 287     {\r
 288         StringBuffer consString = new StringBuffer();\r
 289         String type;\r
 290         Integer result;\r
 291         int[] gapcons;\r
 292         int totGaps, count;\r
 293         float pgaps;\r
 294         Hashtable resultHash ;\r
 295         Enumeration enumeration;\r
 296 \r
 297 \r
 298         for (int i = start; i <= end; i++)\r
 299         {\r
 300             gapcons = countConsNGaps(i);\r
 301             totGaps = gapcons[1];\r
 302             pgaps = ((float) totGaps * 100) / (float) sequences.length;\r
 303 \r
 304             if (percentageGaps > pgaps)\r
 305             {\r
 306                 resultHash =  total[i - start];\r
 307 \r
 308                 //Now find the verdict\r
 309                 count = 0;\r
 310                 enumeration = resultHash.keys();\r
 311 \r
 312                 while (enumeration.hasMoreElements())\r
 313                 {\r
 314                    type = (String) enumeration.nextElement();\r
 315                    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.append(count); // Conserved props!=Identity\r
 337                 }\r
 338                 else\r
 339                 {\r
 340                     consString.append((gapcons[0] == 1) ? "*" : "+");\r
 341                 }\r
 342             }\r
 343             else\r
 344             {\r
 345                 consString.append("-");\r
 346             }\r
 347         }\r
 348 \r
 349         consSequence = new Sequence(name, consString.toString(), 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       seqNums = new Vector();\r
 374      // calcSeqNum(s);\r
 375       int i = 0, iSize = sequences.length;\r
 376     //Do we need to calculate this again?\r
 377       for (i = 0; i < iSize; i++)\r
 378       {\r
 379        calcSeqNum(i);\r
 380       }\r
 381 \r
 382 \r
 383         if ((cons2 == null) || seqNumsChanged)\r
 384         {\r
 385             cons2 = new int[maxLength][24];\r
 386 \r
 387             // Initialize the array\r
 388             for (int j = 0; j < 24; j++)\r
 389             {\r
 390                 for (i = 0; i < maxLength; i++)\r
 391                 {\r
 392                     cons2[i][j] = 0;\r
 393                 }\r
 394             }\r
 395 \r
 396             int[] sqnum;\r
 397             int j = 0;\r
 398 \r
 399             while (j < sequences.length)\r
 400             {\r
 401                 sqnum = (int[]) seqNums.elementAt(j);\r
 402 \r
 403                 for (i = 1; i < sqnum.length; i++)\r
 404                 {\r
 405                     cons2[i - 1][sqnum[i]]++;\r
 406                 }\r
 407 \r
 408                 for (i = sqnum.length - 1; i < maxLength; i++)\r
 409                 {\r
 410                     cons2[i][23]++; // gap count\r
 411                 }\r
 412 \r
 413                 j++;\r
 414             }\r
 415 \r
 416             // unnecessary ?\r
 417 \r
 418             /* for (int i=start; i <= end; i++) {\r
 419                  int max = -1000;\r
 420             int maxi = -1;\r
 421             int maxj = -1;\r
 422 \r
 423             for (int j=0;j<24;j++) {\r
 424               if (cons2[i][j] > max) {\r
 425               max = cons2[i][j];\r
 426               maxi = i;\r
 427               maxj = j;\r
 428             }\r
 429 \r
 430             }\r
 431             } */\r
 432         }\r
 433     }\r
 434 \r
 435     /**\r
 436      * Calculates the quality of the set of sequences\r
 437      *\r
 438      * @param start Start residue\r
 439      * @param end End residue\r
 440      */\r
 441     public void findQuality(int start, int end)\r
 442     {\r
 443         quality = new Vector();\r
 444 \r
 445         double max = -10000;\r
 446         int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();\r
 447 \r
 448         //Loop over columns // JBPNote Profiling info\r
 449         //long ts = System.currentTimeMillis();\r
 450         //long te = System.currentTimeMillis();\r
 451         percentIdentity2();\r
 452 \r
 453         int size = seqNums.size();\r
 454         int[] lengths = new int[size];\r
 455         double tot, bigtot, sr, tmp;\r
 456         double [] x, xx;\r
 457         int l, j, i, ii, i2, k, seqNum;\r
 458 \r
 459         for (l = 0; l < size; l++)\r
 460             lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;\r
 461 \r
 462         for (j = start; j <= end; j++)\r
 463         {\r
 464             bigtot = 0;\r
 465 \r
 466             // First Xr = depends on column only\r
 467             x = new double[24];\r
 468 \r
 469             for (ii = 0; ii < 24; ii++)\r
 470             {\r
 471                 x[ii] = 0;\r
 472 \r
 473                 for (i2 = 0; i2 < 24; i2++)\r
 474                 {\r
 475                   x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) +\r
 476                             4);\r
 477                 }\r
 478 \r
 479                 x[ii] /= size;\r
 480             }\r
 481 \r
 482             // Now calculate D for each position and sum\r
 483             for (k = 0; k < size; k++)\r
 484             {\r
 485                 tot = 0;\r
 486                 xx = new double[24];\r
 487                 seqNum = (j < lengths[k])\r
 488                     ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end\r
 489 \r
 490                 // This is a loop over r\r
 491                 for (i = 0; i < 23; i++)\r
 492                 {\r
 493                     sr = 0;\r
 494 \r
 495                     sr = (double) BLOSUM62[i][seqNum] + 4;\r
 496 \r
 497                     //Calculate X with another loop over residues\r
 498                     //  System.out.println("Xi " + i + " " + x[i] + " " + sr);\r
 499                     xx[i] = x[i] - sr;\r
 500 \r
 501                     tot += (xx[i] * xx[i]);\r
 502                 }\r
 503 \r
 504                 bigtot += Math.sqrt(tot);\r
 505             }\r
 506 \r
 507             // This is the quality for one column\r
 508             if (max < bigtot)\r
 509             {\r
 510                 max = bigtot;\r
 511             }\r
 512 \r
 513             //      bigtot  = bigtot * (size-cons2[j][23])/size;\r
 514             quality.addElement(new Double(bigtot));\r
 515 \r
 516             // Need to normalize by gaps\r
 517         }\r
 518 \r
 519         double newmax = -10000;\r
 520 \r
 521         for (j = start; j <= end; j++)\r
 522         {\r
 523             tmp = ((Double) quality.elementAt(j)).doubleValue();\r
 524             tmp = ((max - tmp) * (size - cons2[j][23])) / size;\r
 525 \r
 526             //     System.out.println(tmp+ " " + j);\r
 527             quality.setElementAt(new Double(tmp), j);\r
 528 \r
 529             if (tmp > newmax)\r
 530             {\r
 531                 newmax = tmp;\r
 532             }\r
 533         }\r
 534 \r
 535         //    System.out.println("Quality " + s);\r
 536         qualityRange[0] = new Double(0);\r
 537         qualityRange[1] = new Double(newmax);\r
 538     }\r
 539 }\r