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