formatting
[jalview.git] / src / jalview / analysis / Conservation.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.7)
3  * Copyright (C) 2011 J Procter, AM Waterhouse, J Engelhardt, LM Lui, G Barton, M Clamp, S Searle
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
10  * 
11  * Jalview is distributed in the hope that it will be useful, but 
12  * WITHOUT ANY WARRANTY; without even the implied warranty 
13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14  * PURPOSE.  See the GNU General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 package jalview.analysis;
19
20 import java.awt.Color;
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
35   int start;
36
37   int end;
38
39   Vector seqNums; // vector of int vectors where first is sequence checksum
40
41   int maxLength = 0; // used by quality calcs
42
43   boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
44
45   Hashtable[] total;
46
47   boolean canonicaliseAa = true; // if true then conservation calculation will
48
49   // map all symbols to canonical aa numbering
50   // rather than consider conservation of that
51   // symbol
52
53   /** Stores calculated quality values */
54   public Vector quality;
55
56   /** Stores maximum and minimum values of quality values */
57   public Double[] qualityRange = new Double[2];
58
59   String consString = "";
60
61   Sequence consSequence;
62
63   Hashtable propHash;
64
65   int threshold;
66
67   String name = "";
68
69   int[][] cons2;
70
71   /**
72    * Creates a new Conservation object.
73    * 
74    * @param name
75    *          Name of conservation
76    * @param propHash
77    *          hash of properties for each symbol
78    * @param threshold
79    *          to count the residues in residueHash(). commonly used value is 3
80    * @param sequences
81    *          sequences to be used in calculation
82    * @param start
83    *          start residue position
84    * @param end
85    *          end residue position
86    */
87   public Conservation(String name, Hashtable propHash, int threshold,
88           List<SequenceI> sequences, int start, int end)
89   {
90     this.name = name;
91     this.propHash = propHash;
92     this.threshold = threshold;
93     this.start = start;
94     this.end = end;
95
96     maxLength = end - start + 1; // default width includes bounds of
97     // calculation
98
99     int s, sSize = sequences.size();
100     SequenceI[] sarray = new SequenceI[sSize];
101     this.sequences = sarray;
102     try {
103     for (s = 0; s < sSize; s++)
104     {
105       sarray[s] = (SequenceI) sequences.get(s);
106       if (sarray[s].getLength() > maxLength)
107       {
108         maxLength = sarray[s].getLength();
109       }
110     }
111     } catch (ArrayIndexOutOfBoundsException ex)
112     {
113       // bail - another thread has modified the sequence array, so the current calculation is probably invalid. 
114       this.sequences=new SequenceI[0];
115       maxLength=0;
116     }
117   }
118
119   /**
120    * Translate sequence i into a numerical representation and store it in the
121    * i'th position of the seqNums array.
122    * 
123    * @param i
124    */
125   private void calcSeqNum(int i)
126   {
127     String sq = null; // for dumb jbuilder not-inited exception warning
128     int[] sqnum = null;
129
130     int sSize = sequences.length;
131
132     if ((i > -1) && (i < sSize))
133     {
134       sq = sequences[i].getSequenceAsString();
135
136       if (seqNums.size() <= i)
137       {
138         seqNums.addElement(new int[sq.length() + 1]);
139       }
140
141       if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
142       {
143         int j;
144         int len;
145         seqNumsChanged = true;
146         len = sq.length();
147
148         if (maxLength < len)
149         {
150           maxLength = len;
151         }
152
153         sqnum = new int[len + 1]; // better to always make a new array -
154         // sequence can change its length
155         sqnum[0] = sq.hashCode();
156
157         for (j = 1; j <= len; j++)
158         {
159           sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
160                   .charAt(j - 1)];
161         }
162
163         seqNums.setElementAt(sqnum, i);
164       }
165       else
166       {
167         System.out.println("SEQUENCE HAS BEEN DELETED!!!");
168       }
169     }
170     else
171     {
172       // JBPNote INFO level debug
173       System.err
174               .println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
175     }
176   }
177
178   /**
179    * Calculates the conservation values for given set of sequences
180    */
181   public void calculate()
182   {
183     Hashtable resultHash, ht;
184     int thresh, j, jSize = sequences.length;
185     int[] values; // Replaces residueHash
186     String type, res = null;
187     char c;
188     Enumeration enumeration2;
189
190     total = new Hashtable[maxLength];
191
192     for (int i = start; i <= end; i++)
193     {
194       values = new int[255];
195
196       for (j = 0; j < jSize; j++)
197       {
198         if (sequences[j].getLength() > i)
199         {
200           c = sequences[j].getCharAt(i);
201
202           if (canonicaliseAa)
203           { // lookup the base aa code symbol
204             c = (char) jalview.schemes.ResidueProperties.aaIndex[sequences[j]
205                     .getCharAt(i)];
206             if (c > 20)
207             {
208               c = '-';
209             }
210             else
211             {
212               // recover canonical aa symbol
213               c = jalview.schemes.ResidueProperties.aa[c].charAt(0);
214             }
215           }
216           else
217           {
218             // original behaviour - operate on ascii symbols directly
219             // No need to check if its a '-'
220             if (c == '.' || c == ' ')
221             {
222               c = '-';
223             }
224
225             if (!canonicaliseAa && 'a' <= c && c <= 'z')
226             {
227               c -= (32); // 32 = 'a' - 'A'
228             }
229           }
230           values[c]++;
231         }
232         else
233         {
234           values['-']++;
235         }
236       }
237
238       // What is the count threshold to count the residues in residueHash()
239       thresh = (threshold * (jSize)) / 100;
240
241       // loop over all the found residues
242       resultHash = new Hashtable();
243       for (char v = '-'; v < 'Z'; v++)
244       {
245
246         if (values[v] > thresh)
247         {
248           res = String.valueOf(v);
249
250           // Now loop over the properties
251           enumeration2 = propHash.keys();
252
253           while (enumeration2.hasMoreElements())
254           {
255             type = (String) enumeration2.nextElement();
256             ht = (Hashtable) propHash.get(type);
257
258             // Have we ticked this before?
259             if (!resultHash.containsKey(type))
260             {
261               if (ht.containsKey(res))
262               {
263                 resultHash.put(type, ht.get(res));
264               }
265               else
266               {
267                 resultHash.put(type, ht.get("-"));
268               }
269             }
270             else if (((Integer) resultHash.get(type)).equals((Integer) ht
271                     .get(res)) == false)
272             {
273               resultHash.put(type, new Integer(-1));
274             }
275           }
276         }
277       }
278
279       if (total.length>0) {
280         total[i - start] = resultHash;
281       }
282     }
283   }
284
285   /*****************************************************************************
286    * count conservation for the j'th column of the alignment
287    * 
288    * @return { gap count, conserved residue count}
289    */
290   public int[] countConsNGaps(int j)
291   {
292     int count = 0;
293     int cons = 0;
294     int nres = 0;
295     int[] r = new int[2];
296     char f = '$';
297     int i, iSize = sequences.length;
298     char c;
299
300     for (i = 0; i < iSize; i++)
301     {
302       if (j >= sequences[i].getLength())
303       {
304         count++;
305         continue;
306       }
307
308       c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
309
310       if (jalview.util.Comparison.isGap((c)))
311       {
312         count++;
313       }
314       else
315       {
316         nres++;
317
318         if (nres == 1)
319         {
320           f = c;
321           cons++;
322         }
323         else if (f == c)
324         {
325           cons++;
326         }
327       }
328     }
329
330     r[0] = (nres == cons) ? 1 : 0;
331     r[1] = count;
332
333     return r;
334   }
335
336   /**
337    * Calculates the conservation sequence
338    * 
339    * @param consflag
340    *          if true, poitiveve conservation; false calculates negative
341    *          conservation
342    * @param percentageGaps
343    *          commonly used value is 25
344    */
345   public void verdict(boolean consflag, float percentageGaps)
346   {
347     StringBuffer consString = new StringBuffer();
348     String type;
349     Integer result;
350     int[] gapcons;
351     int totGaps, count;
352     float pgaps;
353     Hashtable resultHash;
354     Enumeration enumeration;
355
356     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
357     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
358     // DOES NOT EXIST IN JALVIEW 2.1.2
359     for (int i = 0; i < start; i++)
360     {
361       consString.append('-');
362     }
363
364     for (int i = start; i <= end; i++)
365     {
366       gapcons = countConsNGaps(i);
367       totGaps = gapcons[1];
368       pgaps = ((float) totGaps * 100) / (float) sequences.length;
369
370       if (percentageGaps > pgaps)
371       {
372         resultHash = total[i - start];
373
374         // Now find the verdict
375         count = 0;
376         enumeration = resultHash.keys();
377
378         while (enumeration.hasMoreElements())
379         {
380           type = (String) enumeration.nextElement();
381           result = (Integer) resultHash.get(type);
382
383           // Do we want to count +ve conservation or +ve and -ve cons.?
384           if (consflag)
385           {
386             if (result.intValue() == 1)
387             {
388               count++;
389             }
390           }
391           else
392           {
393             if (result.intValue() != -1)
394             {
395               count++;
396             }
397           }
398         }
399
400         if (count < 10)
401         {
402           consString.append(count); // Conserved props!=Identity
403         }
404         else
405         {
406           consString.append((gapcons[0] == 1) ? "*" : "+");
407         }
408       }
409       else
410       {
411         consString.append('-');
412       }
413     }
414
415     consSequence = new Sequence(name, consString.toString(), start, end);
416   }
417
418   /**
419    * 
420    * 
421    * @return Conservation sequence
422    */
423   public Sequence getConsSequence()
424   {
425     return consSequence;
426   }
427
428   // From Alignment.java in jalview118
429   public void findQuality()
430   {
431     findQuality(0, maxLength - 1);
432   }
433
434   /**
435    * DOCUMENT ME!
436    */
437   private void percentIdentity2()
438   {
439     seqNums = new Vector();
440     // calcSeqNum(s);
441     int i = 0, iSize = sequences.length;
442     // Do we need to calculate this again?
443     for (i = 0; i < iSize; i++)
444     {
445       calcSeqNum(i);
446     }
447
448     if ((cons2 == null) || seqNumsChanged)
449     {
450       cons2 = new int[maxLength][24];
451
452       // Initialize the array
453       for (int j = 0; j < 24; j++)
454       {
455         for (i = 0; i < maxLength; i++)
456         {
457           cons2[i][j] = 0;
458         }
459       }
460
461       int[] sqnum;
462       int j = 0;
463
464       while (j < sequences.length)
465       {
466         sqnum = (int[]) seqNums.elementAt(j);
467
468         for (i = 1; i < sqnum.length; i++)
469         {
470           cons2[i - 1][sqnum[i]]++;
471         }
472
473         for (i = sqnum.length - 1; i < maxLength; i++)
474         {
475           cons2[i][23]++; // gap count
476         }
477
478         j++;
479       }
480
481       // unnecessary ?
482
483       /*
484        * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
485        * maxj = -1;
486        * 
487        * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
488        * maxi = i; maxj = j; } } }
489        */
490     }
491   }
492
493   /**
494    * Calculates the quality of the set of sequences
495    * 
496    * @param start
497    *          Start residue
498    * @param end
499    *          End residue
500    */
501   public void findQuality(int start, int end)
502   {
503     quality = new Vector();
504
505     double max = -10000;
506     int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
507
508     // Loop over columns // JBPNote Profiling info
509     // long ts = System.currentTimeMillis();
510     // long te = System.currentTimeMillis();
511     percentIdentity2();
512
513     int size = seqNums.size();
514     int[] lengths = new int[size];
515     double tot, bigtot, sr, tmp;
516     double[] x, xx;
517     int l, j, i, ii, i2, k, seqNum;
518
519     for (l = 0; l < size; l++)
520     {
521       lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
522     }
523
524     for (j = start; j <= end; j++)
525     {
526       bigtot = 0;
527
528       // First Xr = depends on column only
529       x = new double[24];
530
531       for (ii = 0; ii < 24; ii++)
532       {
533         x[ii] = 0;
534
535         for (i2 = 0; i2 < 24; i2++)
536         {
537           x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
538         }
539
540         x[ii] /= size;
541       }
542
543       // Now calculate D for each position and sum
544       for (k = 0; k < size; k++)
545       {
546         tot = 0;
547         xx = new double[24];
548         seqNum = (j < lengths[k]) ? ((int[]) seqNums.elementAt(k))[j + 1]
549                 : 23; // Sequence, or gap at the end
550
551         // This is a loop over r
552         for (i = 0; i < 23; i++)
553         {
554           sr = 0;
555
556           sr = (double) BLOSUM62[i][seqNum] + 4;
557
558           // Calculate X with another loop over residues
559           // System.out.println("Xi " + i + " " + x[i] + " " + sr);
560           xx[i] = x[i] - sr;
561
562           tot += (xx[i] * xx[i]);
563         }
564
565         bigtot += Math.sqrt(tot);
566       }
567
568       // This is the quality for one column
569       if (max < bigtot)
570       {
571         max = bigtot;
572       }
573
574       // bigtot = bigtot * (size-cons2[j][23])/size;
575       quality.addElement(new Double(bigtot));
576
577       // Need to normalize by gaps
578     }
579
580     double newmax = -10000;
581
582     for (j = start; j <= end; j++)
583     {
584       tmp = ((Double) quality.elementAt(j)).doubleValue();
585       tmp = ((max - tmp) * (size - cons2[j][23])) / size;
586
587       // System.out.println(tmp+ " " + j);
588       quality.setElementAt(new Double(tmp), j);
589
590       if (tmp > newmax)
591       {
592         newmax = tmp;
593       }
594     }
595
596     // System.out.println("Quality " + s);
597     qualityRange[0] = new Double(0);
598     qualityRange[1] = new Double(newmax);
599   }
600
601   /**
602    * complete the given consensus and quuality annotation rows. Note: currently
603    * this method will enlarge the given annotation row if it is too small,
604    * otherwise will leave its length unchanged.
605    * 
606    * @param conservation
607    *          conservation annotation row
608    * @param quality2
609    *          (optional - may be null)
610    * @param istart
611    *          first column for conservation
612    * @param alWidth
613    *          extent of conservation
614    */
615   public void completeAnnotations(AlignmentAnnotation conservation,
616           AlignmentAnnotation quality2, int istart, int alWidth)
617   {
618     char[] sequence = getConsSequence().getSequence();
619     float minR;
620     float minG;
621     float minB;
622     float maxR;
623     float maxG;
624     float maxB;
625     minR = 0.3f;
626     minG = 0.0f;
627     minB = 0f;
628     maxR = 1.0f - minR;
629     maxG = 0.9f - minG;
630     maxB = 0f - minB; // scalable range for colouring both Conservation and
631     // Quality
632
633     float min = 0f;
634     float max = 11f;
635     float qmin = 0f;
636     float qmax = 0f;
637
638     char c;
639
640     if (conservation.annotations != null
641             && conservation.annotations.length < alWidth)
642     {
643       conservation.annotations = new Annotation[alWidth];
644     }
645
646     if (quality2 != null)
647     {
648       quality2.graphMax = qualityRange[1].floatValue();
649       if (quality2.annotations != null
650               && quality2.annotations.length < alWidth)
651       {
652         quality2.annotations = new Annotation[alWidth];
653       }
654       qmin = qualityRange[0].floatValue();
655       qmax = qualityRange[1].floatValue();
656     }
657
658     for (int i = 0; i < alWidth; i++)
659     {
660       float value = 0;
661
662       c = sequence[i];
663
664       if (Character.isDigit(c))
665       {
666         value = (int) (c - '0');
667       }
668       else if (c == '*')
669       {
670         value = 11;
671       }
672       else if (c == '+')
673       {
674         value = 10;
675       }
676
677       float vprop = value - min;
678       vprop /= max;
679       conservation.annotations[i] = new Annotation(String.valueOf(c),
680               String.valueOf(value), ' ', value, new Color(minR
681                       + (maxR * vprop), minG + (maxG * vprop), minB
682                       + (maxB * vprop)));
683
684       // Quality calc
685       if (quality2 != null)
686       {
687         value = ((Double) quality.elementAt(i)).floatValue();
688         vprop = value - qmin;
689         vprop /= qmax;
690         quality2.annotations[i] = new Annotation(" ",
691                 String.valueOf(value), ' ', value, new Color(minR
692                         + (maxR * vprop), minG + (maxG * vprop), minB
693                         + (maxB * vprop)));
694       }
695     }
696   }
697
698   /**
699    * construct and call the calculation methods on a new Conservation object
700    * @param name - name of conservation
701    * @param consHash - hash table of properties for each amino acid (normally ResidueProperties.propHash)
702    * @param threshold - minimum number of conserved residues needed to indicate conservation (typically 3)
703    * @param seqs
704    * @param start first column in calculation window
705    * @param end last column in calculation window
706    * @param posOrNeg positive (true) or negative (false) conservation 
707    * @param consPercGaps percentage of gaps tolerated in column
708    * @param calcQuality flag indicating if alignment quality should be calculated  
709    * @return Conservation object ready for use in visualization
710    */
711   public static Conservation calculateConservation(String name,
712           Hashtable consHash, int threshold, List<SequenceI> seqs, int start, int end, boolean posOrNeg, int consPercGaps, boolean calcQuality)
713   {
714     Conservation cons = new Conservation(name, consHash, threshold, seqs, start,end);
715     return calculateConservation(cons, posOrNeg, consPercGaps, calcQuality);
716   }
717   /**
718   * @param b positive (true) or negative (false) conservation 
719   * @param consPercGaps percentage of gaps tolerated in column
720   * @param calcQuality flag indicating if alignment quality should be calculated  
721   * @return Conservation object ready for use in visualization
722   */
723  public static Conservation calculateConservation(Conservation cons, boolean b, int consPercGaps, boolean calcQuality)
724  {
725    cons.calculate();
726     cons.verdict(b, consPercGaps); 
727
728     if (calcQuality)
729     {
730       cons.findQuality();
731     }
732
733     return cons;
734   }
735 }