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