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