JAL-2738 copy to spikes/mungo
[jalview.git] / src / jalview / analysis / Conservation.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
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
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.analysis;
22
23 import jalview.analysis.scoremodels.ScoreMatrix;
24 import jalview.analysis.scoremodels.ScoreModels;
25 import jalview.datamodel.AlignmentAnnotation;
26 import jalview.datamodel.Annotation;
27 import jalview.datamodel.ResidueCount;
28 import jalview.datamodel.ResidueCount.SymbolCounts;
29 import jalview.datamodel.Sequence;
30 import jalview.datamodel.SequenceI;
31 import jalview.schemes.ResidueProperties;
32 import jalview.util.Comparison;
33
34 import java.awt.Color;
35 import java.util.List;
36 import java.util.Map;
37 import java.util.Map.Entry;
38 import java.util.SortedMap;
39 import java.util.TreeMap;
40 import java.util.Vector;
41
42 /**
43  * Calculates conservation values for a given set of sequences
44  */
45 public class Conservation
46 {
47   /*
48    * need to have a minimum of 3% of sequences with a residue
49    * for it to be included in the conservation calculation
50    */
51   private static final int THRESHOLD_PERCENT = 3;
52
53   private static final int TOUPPERCASE = 'a' - 'A';
54
55   private static final int GAP_INDEX = -1;
56
57   SequenceI[] sequences;
58
59   int start;
60
61   int end;
62
63   /*
64    * a list whose i'th element is an array whose first entry is the checksum
65    * of the i'th sequence, followed by residues encoded to score matrix index
66    */
67   Vector<int[]> seqNums;
68
69   int maxLength = 0; // used by quality calcs
70
71   boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
72
73   /*
74    * a map per column with {property, conservation} where conservation value is
75    * 1 (property is conserved), 0 (absence of property is conserved) or -1
76    * (property is not conserved i.e. column has residues with and without it)
77    */
78   Map<String, Integer>[] total;
79
80   /*
81    * if true then conservation calculation will map all symbols to canonical aa
82    * numbering rather than consider conservation of that symbol
83    */
84   boolean canonicaliseAa = true;
85
86   private Vector<Double> quality;
87
88   private double qualityMinimum;
89
90   private double qualityMaximum;
91
92   private Sequence consSequence;
93
94   /*
95    * percentage of residues in a column to qualify for counting conservation
96    */
97   private int threshold;
98
99   private String name = "";
100
101   /*
102    * an array, for each column, of counts of symbols (by score matrix index)
103    */
104   private int[][] cons2;
105
106   /*
107    * gap counts for each column
108    */
109   private int[] cons2GapCounts;
110
111   private String[] consSymbs;
112
113   /**
114    * Constructor using default threshold of 3%
115    * 
116    * @param name
117    *          Name of conservation
118    * @param sequences
119    *          sequences to be used in calculation
120    * @param start
121    *          start residue position
122    * @param end
123    *          end residue position
124    */
125   public Conservation(String name, List<SequenceI> sequences, int start,
126           int end)
127   {
128     this(name, THRESHOLD_PERCENT, sequences, start, end);
129   }
130
131   /**
132    * Constructor
133    * 
134    * @param name
135    *          Name of conservation
136    * @param threshold
137    *          percentage of sequences at or below which property conservation is
138    *          ignored
139    * @param sequences
140    *          sequences to be used in calculation
141    * @param start
142    *          start column position
143    * @param end
144    *          end column position
145    */
146   public Conservation(String name, int threshold, List<SequenceI> sequences,
147           int start, int end)
148   {
149     this.name = name;
150     this.threshold = threshold;
151     this.start = start;
152     this.end = end;
153
154     maxLength = end - start + 1; // default width includes bounds of
155     // calculation
156
157     int s, sSize = sequences.size();
158     SequenceI[] sarray = new SequenceI[sSize];
159     this.sequences = sarray;
160     try
161     {
162       for (s = 0; s < sSize; s++)
163       {
164         sarray[s] = sequences.get(s);
165         if (sarray[s].getLength() > maxLength)
166         {
167           maxLength = sarray[s].getLength();
168         }
169       }
170     } catch (ArrayIndexOutOfBoundsException ex)
171     {
172       // bail - another thread has modified the sequence array, so the current
173       // calculation is probably invalid.
174       this.sequences = new SequenceI[0];
175       maxLength = 0;
176     }
177   }
178
179   /**
180    * Translate sequence i into score matrix indices and store it in the i'th
181    * position of the seqNums array.
182    * 
183    * @param i
184    * @param sm
185    */
186   private void calcSeqNum(int i, ScoreMatrix sm)
187   {
188     int sSize = sequences.length;
189
190     if ((i > -1) && (i < sSize))
191     {
192       String sq = sequences[i].getSequenceAsString();
193
194       if (seqNums.size() <= i)
195       {
196         seqNums.addElement(new int[sq.length() + 1]);
197       }
198
199       /*
200        * the first entry in the array is the sequence's hashcode,
201        * following entries are matrix indices of sequence characters
202        */
203       if (sq.hashCode() != seqNums.elementAt(i)[0])
204       {
205         int j;
206         int len;
207         seqNumsChanged = true;
208         len = sq.length();
209
210         if (maxLength < len)
211         {
212           maxLength = len;
213         }
214
215         int[] sqnum = new int[len + 1]; // better to always make a new array -
216         // sequence can change its length
217         sqnum[0] = sq.hashCode();
218
219         for (j = 1; j <= len; j++)
220         {
221           // sqnum[j] = ResidueProperties.aaIndex[sq.charAt(j - 1)];
222           char residue = sq.charAt(j - 1);
223           if (Comparison.isGap(residue))
224           {
225             sqnum[j] = GAP_INDEX;
226           }
227           else
228           {
229             sqnum[j] = sm.getMatrixIndex(residue);
230             if (sqnum[j] == -1)
231             {
232               sqnum[j] = GAP_INDEX;
233             }
234           }
235         }
236
237         seqNums.setElementAt(sqnum, i);
238       }
239       else
240       {
241         System.out.println("SEQUENCE HAS BEEN DELETED!!!");
242       }
243     }
244     else
245     {
246       // JBPNote INFO level debug
247       System.err.println(
248               "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
249     }
250   }
251
252   /**
253    * Calculates the conservation values for given set of sequences
254    */
255   public void calculate()
256   {
257     int height = sequences.length;
258
259     total = new Map[maxLength];
260
261     for (int column = start; column <= end; column++)
262     {
263       ResidueCount values = countResidues(column);
264
265       /*
266        * percentage count at or below which we ignore residues
267        */
268       int thresh = (threshold * height) / 100;
269
270       /*
271        * check observed residues in column and record whether each 
272        * physico-chemical property is conserved (+1), absence conserved (0),
273        * or not conserved (-1)
274        * Using TreeMap means properties are displayed in alphabetical order
275        */
276       SortedMap<String, Integer> resultHash = new TreeMap<String, Integer>();
277       SymbolCounts symbolCounts = values.getSymbolCounts();
278       char[] symbols = symbolCounts.symbols;
279       int[] counts = symbolCounts.values;
280       for (int j = 0; j < symbols.length; j++)
281       {
282         char c = symbols[j];
283         if (counts[j] > thresh)
284         {
285           recordConservation(resultHash, String.valueOf(c));
286         }
287       }
288       if (values.getGapCount() > thresh)
289       {
290         recordConservation(resultHash, "-");
291       }
292
293       if (total.length > 0)
294       {
295         total[column - start] = resultHash;
296       }
297     }
298   }
299
300   /**
301    * Updates the conservation results for an observed residue
302    * 
303    * @param resultMap
304    *          a map of {property, conservation} where conservation value is +1
305    *          (all residues have the property), 0 (no residue has the property)
306    *          or -1 (some do, some don't)
307    * @param res
308    */
309   protected static void recordConservation(Map<String, Integer> resultMap,
310           String res)
311   {
312     res = res.toUpperCase();
313     for (Entry<String, Map<String, Integer>> property : ResidueProperties.propHash
314             .entrySet())
315     {
316       String propertyName = property.getKey();
317       Integer residuePropertyValue = property.getValue().get(res);
318
319       if (!resultMap.containsKey(propertyName))
320       {
321         /*
322          * first time we've seen this residue - note whether it has this property
323          */
324         if (residuePropertyValue != null)
325         {
326           resultMap.put(propertyName, residuePropertyValue);
327         }
328         else
329         {
330           /*
331            * unrecognised residue - use default value for property
332            */
333           resultMap.put(propertyName, property.getValue().get("-"));
334         }
335       }
336       else
337       {
338         Integer currentResult = resultMap.get(propertyName);
339         if (currentResult.intValue() != -1
340                 && !currentResult.equals(residuePropertyValue))
341         {
342           /*
343            * property is unconserved - residues seen both with and without it
344            */
345           resultMap.put(propertyName, Integer.valueOf(-1));
346         }
347       }
348     }
349   }
350
351   /**
352    * Counts residues (upper-cased) and gaps in the given column
353    * 
354    * @param column
355    * @return
356    */
357   protected ResidueCount countResidues(int column)
358   {
359     ResidueCount values = new ResidueCount(false);
360
361     for (int row = 0; row < sequences.length; row++)
362     {
363       if (sequences[row].getLength() > column)
364       {
365         char c = sequences[row].getCharAt(column);
366         if (canonicaliseAa)
367         {
368           int index = ResidueProperties.aaIndex[c];
369           c = index > 20 ? '-' : ResidueProperties.aa[index].charAt(0);
370         }
371         else
372         {
373           c = toUpperCase(c);
374         }
375         if (Comparison.isGap(c))
376         {
377           values.addGap();
378         }
379         else
380         {
381           values.add(c);
382         }
383       }
384       else
385       {
386         values.addGap();
387       }
388     }
389     return values;
390   }
391
392   /**
393    * Counts conservation and gaps for a column of the alignment
394    * 
395    * @return { 1 if fully conserved, else 0, gap count }
396    */
397   public int[] countConservationAndGaps(int column)
398   {
399     int gapCount = 0;
400     boolean fullyConserved = true;
401     int iSize = sequences.length;
402
403     if (iSize == 0)
404     {
405       return new int[] { 0, 0 };
406     }
407
408     char lastRes = '0';
409     for (int i = 0; i < iSize; i++)
410     {
411       if (column >= sequences[i].getLength())
412       {
413         gapCount++;
414         continue;
415       }
416
417       char c = sequences[i].getCharAt(column); // gaps do not have upper/lower
418                                                // case
419
420       if (Comparison.isGap((c)))
421       {
422         gapCount++;
423       }
424       else
425       {
426         c = toUpperCase(c);
427         if (lastRes == '0')
428         {
429           lastRes = c;
430         }
431         if (c != lastRes)
432         {
433           fullyConserved = false;
434         }
435       }
436     }
437
438     int[] r = new int[] { fullyConserved ? 1 : 0, gapCount };
439     return r;
440   }
441
442   /**
443    * Returns the upper-cased character if between 'a' and 'z', else the
444    * unchanged value
445    * 
446    * @param c
447    * @return
448    */
449   char toUpperCase(char c)
450   {
451     if ('a' <= c && c <= 'z')
452     {
453       c -= TOUPPERCASE;
454     }
455     return c;
456   }
457
458   /**
459    * Calculates the conservation sequence
460    * 
461    * @param positiveOnly
462    *          if true, calculate positive conservation; else calculate both
463    *          positive and negative conservation
464    * @param maxPercentageGaps
465    *          the percentage of gaps in a column, at or above which no
466    *          conservation is asserted
467    */
468   public void verdict(boolean positiveOnly, float maxPercentageGaps)
469   {
470     // TODO call this at the end of calculate(), should not be a public method
471
472     StringBuilder consString = new StringBuilder(end);
473
474     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
475     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
476     // DOES NOT EXIST IN JALVIEW 2.1.2
477     for (int i = 0; i < start; i++)
478     {
479       consString.append('-');
480     }
481     consSymbs = new String[end - start + 1];
482     for (int i = start; i <= end; i++)
483     {
484       int[] gapcons = countConservationAndGaps(i);
485       boolean fullyConserved = gapcons[0] == 1;
486       int totGaps = gapcons[1];
487       float pgaps = (totGaps * 100f) / sequences.length;
488
489       if (maxPercentageGaps > pgaps)
490       {
491         Map<String, Integer> resultHash = total[i - start];
492         int count = 0;
493         StringBuilder positives = new StringBuilder(64);
494         StringBuilder negatives = new StringBuilder(32);
495         for (String type : resultHash.keySet())
496         {
497           int result = resultHash.get(type).intValue();
498           if (result == -1)
499           {
500             /*
501              * not conserved (present or absent)
502              */
503             continue;
504           }
505           count++;
506           if (result == 1)
507           {
508             /*
509              * positively conserved property (all residues have it)
510              */
511             positives.append(positives.length() == 0 ? "" : " ");
512             positives.append(type);
513           }
514           if (result == 0 && !positiveOnly)
515           {
516             /*
517              * absense of property is conserved (all residues lack it)
518              */
519             negatives.append(negatives.length() == 0 ? "" : " ");
520             negatives.append("!").append(type);
521           }
522         }
523         if (negatives.length() > 0)
524         {
525           positives.append(" ").append(negatives);
526         }
527         consSymbs[i - start] = positives.toString();
528
529         if (count < 10)
530         {
531           consString.append(count); // Conserved props!=Identity
532         }
533         else
534         {
535           consString.append(fullyConserved ? "*" : "+");
536         }
537       }
538       else
539       {
540         consString.append('-');
541       }
542     }
543
544     consSequence = new Sequence(name, consString.toString(), start, end);
545   }
546
547   /**
548    * 
549    * 
550    * @return Conservation sequence
551    */
552   public SequenceI getConsSequence()
553   {
554     return consSequence;
555   }
556
557   // From Alignment.java in jalview118
558   public void findQuality()
559   {
560     findQuality(0, maxLength - 1, ScoreModels.getInstance().getBlosum62());
561   }
562
563   /**
564    * DOCUMENT ME!
565    * 
566    * @param sm
567    */
568   private void percentIdentity(ScoreMatrix sm)
569   {
570     seqNums = new Vector<int[]>();
571     int i = 0, iSize = sequences.length;
572     // Do we need to calculate this again?
573     for (i = 0; i < iSize; i++)
574     {
575       calcSeqNum(i, sm);
576     }
577
578     if ((cons2 == null) || seqNumsChanged)
579     {
580       // FIXME remove magic number 24 without changing calc
581       // sm.getSize() returns 25 so doesn't quite do it...
582       cons2 = new int[maxLength][24];
583       cons2GapCounts = new int[maxLength];
584
585       int j = 0;
586
587       while (j < sequences.length)
588       {
589         int[] sqnum = seqNums.elementAt(j);
590
591         for (i = 1; i < sqnum.length; i++)
592         {
593           int index = sqnum[i];
594           if (index == GAP_INDEX)
595           {
596             cons2GapCounts[i - 1]++;
597           }
598           else
599           {
600             cons2[i - 1][index]++;
601           }
602         }
603
604         // TODO should this start from sqnum.length?
605         for (i = sqnum.length - 1; i < maxLength; i++)
606         {
607           cons2GapCounts[i]++;
608         }
609         j++;
610       }
611     }
612   }
613
614   /**
615    * Calculates the quality of the set of sequences over the given inclusive
616    * column range, using the specified substitution score matrix
617    * 
618    * @param startCol
619    * @param endCol
620    * @param scoreMatrix
621    */
622   protected void findQuality(int startCol, int endCol,
623           ScoreMatrix scoreMatrix)
624   {
625     quality = new Vector<Double>();
626
627     double max = -Double.MAX_VALUE;
628     float[][] scores = scoreMatrix.getMatrix();
629
630     percentIdentity(scoreMatrix);
631
632     int size = seqNums.size();
633     int[] lengths = new int[size];
634
635     for (int l = 0; l < size; l++)
636     {
637       lengths[l] = seqNums.elementAt(l).length - 1;
638     }
639
640     final int symbolCount = scoreMatrix.getSize();
641
642     for (int j = startCol; j <= endCol; j++)
643     {
644       double bigtot = 0;
645
646       // First Xr = depends on column only
647       double[] x = new double[symbolCount];
648
649       for (int ii = 0; ii < symbolCount; ii++)
650       {
651         x[ii] = 0;
652
653         /*
654          * todo JAL-728 currently assuming last symbol in matrix is * for gap
655          * (which we ignore as counted separately); true for BLOSUM62 but may
656          * not be once alternative matrices are supported
657          */
658         for (int i2 = 0; i2 < symbolCount - 1; i2++)
659         {
660           x[ii] += (((double) cons2[j][i2] * scores[ii][i2]) + 4D);
661         }
662         x[ii] += 4D + cons2GapCounts[j] * scoreMatrix.getMinimumScore();
663
664         x[ii] /= size;
665       }
666
667       // Now calculate D for each position and sum
668       for (int k = 0; k < size; k++)
669       {
670         double tot = 0;
671         double[] xx = new double[symbolCount];
672         // sequence character index, or implied gap if sequence too short
673         int seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1]
674                 : GAP_INDEX;
675
676         for (int i = 0; i < symbolCount - 1; i++)
677         {
678           double sr = 4D;
679           if (seqNum == GAP_INDEX)
680           {
681             sr += scoreMatrix.getMinimumScore();
682           }
683           else
684           {
685             sr += scores[i][seqNum];
686           }
687
688           xx[i] = x[i] - sr;
689
690           tot += (xx[i] * xx[i]);
691         }
692
693         bigtot += Math.sqrt(tot);
694       }
695
696       max = Math.max(max, bigtot);
697
698       quality.addElement(new Double(bigtot));
699     }
700
701     double newmax = -Double.MAX_VALUE;
702
703     for (int j = startCol; j <= endCol; j++)
704     {
705       double tmp = quality.elementAt(j).doubleValue();
706       // tmp = ((max - tmp) * (size - cons2[j][23])) / size;
707       tmp = ((max - tmp) * (size - cons2GapCounts[j])) / size;
708
709       // System.out.println(tmp+ " " + j);
710       quality.setElementAt(new Double(tmp), j);
711
712       if (tmp > newmax)
713       {
714         newmax = tmp;
715       }
716     }
717
718     qualityMinimum = 0D;
719     qualityMaximum = newmax;
720   }
721
722   /**
723    * Complete the given consensus and quuality annotation rows. Note: currently
724    * this method will enlarge the given annotation row if it is too small,
725    * otherwise will leave its length unchanged.
726    * 
727    * @param conservation
728    *          conservation annotation row
729    * @param quality2
730    *          (optional - may be null)
731    * @param istart
732    *          first column for conservation
733    * @param alWidth
734    *          extent of conservation
735    */
736   public void completeAnnotations(AlignmentAnnotation conservation,
737           AlignmentAnnotation quality2, int istart, int alWidth)
738   {
739     SequenceI cons = getConsSequence();
740
741     /*
742      * colour scale for Conservation and Quality;
743      */
744     float minR = 0.3f;
745     float minG = 0.0f;
746     float minB = 0f;
747     float maxR = 1.0f - minR;
748     float maxG = 0.9f - minG;
749     float maxB = 0f - minB;
750
751     float min = 0f;
752     float max = 11f;
753     float qmin = 0f;
754     float qmax = 0f;
755
756     if (conservation != null && conservation.annotations != null
757             && conservation.annotations.length < alWidth)
758     {
759       conservation.annotations = new Annotation[alWidth];
760     }
761
762     if (quality2 != null)
763     {
764       quality2.graphMax = (float) qualityMaximum;
765       if (quality2.annotations != null
766               && quality2.annotations.length < alWidth)
767       {
768         quality2.annotations = new Annotation[alWidth];
769       }
770       qmin = (float) qualityMinimum;
771       qmax = (float) qualityMaximum;
772     }
773
774     for (int i = istart; i < alWidth; i++)
775     {
776       float value = 0;
777
778       char c = cons.getCharAt(i);
779
780       if (Character.isDigit(c))
781       {
782         value = c - '0';
783       }
784       else if (c == '*')
785       {
786         value = 11;
787       }
788       else if (c == '+')
789       {
790         value = 10;
791       }
792
793       if (conservation != null)
794       {
795         float vprop = value - min;
796         vprop /= max;
797         int consp = i - start;
798         String conssym = (value > 0 && consp > -1
799                 && consp < consSymbs.length) ? consSymbs[consp] : "";
800         conservation.annotations[i] = new Annotation(String.valueOf(c),
801                 conssym, ' ', value, new Color(minR + (maxR * vprop),
802                         minG + (maxG * vprop), minB + (maxB * vprop)));
803       }
804
805       // Quality calc
806       if (quality2 != null)
807       {
808         value = quality.elementAt(i).floatValue();
809         float vprop = value - qmin;
810         vprop /= qmax;
811         quality2.annotations[i] = new Annotation(" ", String.valueOf(value),
812                 ' ', value, new Color(minR + (maxR * vprop),
813                         minG + (maxG * vprop), minB + (maxB * vprop)));
814       }
815     }
816   }
817
818   /**
819    * construct and call the calculation methods on a new Conservation object
820    * 
821    * @param name
822    *          - name of conservation
823    * @param seqs
824    * @param start
825    *          first column in calculation window
826    * @param end
827    *          last column in calculation window
828    * @param positiveOnly
829    *          calculate positive (true) or positive and negative (false)
830    *          conservation
831    * @param maxPercentGaps
832    *          percentage of gaps tolerated in column
833    * @param calcQuality
834    *          flag indicating if alignment quality should be calculated
835    * @return Conservation object ready for use in visualization
836    */
837   public static Conservation calculateConservation(String name,
838           List<SequenceI> seqs, int start, int end, boolean positiveOnly,
839           int maxPercentGaps, boolean calcQuality)
840   {
841     Conservation cons = new Conservation(name, seqs, start, end);
842     cons.calculate();
843     cons.verdict(positiveOnly, maxPercentGaps);
844
845     if (calcQuality)
846     {
847       cons.findQuality();
848     }
849
850     return cons;
851   }
852
853   /**
854    * Returns the computed tooltip (annotation description) for a given column.
855    * The tip is empty if the conservation score is zero, otherwise holds the
856    * conserved properties (and, optionally, properties whose absence is
857    * conserved).
858    * 
859    * @param column
860    * @return
861    */
862   String getTooltip(int column)
863   {
864     SequenceI cons = getConsSequence();
865     char val = column < cons.getLength() ? cons.getCharAt(column) : '-';
866     boolean hasConservation = val != '-' && val != '0';
867     int consp = column - start;
868     String tip = (hasConservation && consp > -1 && consp < consSymbs.length)
869             ? consSymbs[consp]
870             : "";
871     return tip;
872   }
873 }