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