JS optimization for too aggressive dynamic operation when entering PDB
[jalview.git] / src / jalview / analysis / scoremodels / ScoreMatrix.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.scoremodels;
22
23 import jalview.api.AlignmentViewPanel;
24 import jalview.api.analysis.PairwiseScoreModelI;
25 import jalview.api.analysis.ScoreModelI;
26 import jalview.api.analysis.SimilarityParamsI;
27 import jalview.datamodel.AlignmentView;
28 import jalview.math.Matrix;
29 import jalview.math.MatrixI;
30 import jalview.util.Comparison;
31
32 import java.util.Arrays;
33
34 /**
35  * A class that models a substitution score matrix for any given alphabet of
36  * symbols. Instances of this class are immutable and thread-safe, so the same
37  * object is returned from calls to getInstance().
38  */
39 public class ScoreMatrix extends SimilarityScoreModel
40         implements PairwiseScoreModelI
41 {
42   private static final char GAP_CHARACTER = Comparison.GAP_DASH;
43
44   /*
45    * an arbitrary score to assign for identity of an unknown symbol
46    * (this is the value on the diagonal in the * column of the NCBI matrix)
47    * (though a case could be made for using the minimum diagonal value)
48    */
49   private static final int UNKNOWN_IDENTITY_SCORE = 1;
50
51   /*
52    * Jalview 2.10.1 treated gaps as X (peptide) or N (nucleotide)
53    * for pairwise scoring; 2.10.2 uses gap score (last column) in
54    * score matrix (JAL-2397)
55    * Set this flag to true (via Groovy) for 2.10.1 behaviour
56    */
57   private static boolean scoreGapAsAny = false;
58
59   public static final short UNMAPPED = (short) -1;
60
61   private static final String BAD_ASCII_ERROR = "Unexpected character %s in getPairwiseScore";
62
63   private static final int MAX_ASCII = 127;
64
65   /*
66    * the name of the model as shown in menus
67    * each score model in use should have a unique name
68    */
69   private String name;
70
71   /*
72    * a description for the model as shown in tooltips
73    */
74   private String description;
75
76   /*
77    * the characters that the model provides scores for
78    */
79   private char[] symbols;
80
81   /*
82    * the score matrix; both dimensions must equal the number of symbols
83    * matrix[i][j] is the substitution score for replacing symbols[i] with symbols[j]
84    */
85   private float[][] matrix;
86
87   /*
88    * quick lookup to convert from an ascii character value to the index
89    * of the corresponding symbol in the score matrix 
90    */
91   private short[] symbolIndex;
92
93   /*
94    * true for Protein Score matrix, false for dna score matrix
95    */
96   private boolean peptide;
97
98   private float minValue;
99
100   private float maxValue;
101
102   /**
103    * Constructor given a name, symbol alphabet, and matrix of scores for pairs
104    * of symbols. The matrix should be square and of the same size as the
105    * alphabet, for example 20x20 for a 20 symbol alphabet.
106    * 
107    * @param theName
108    *          Unique, human readable name for the matrix
109    * @param alphabet
110    *          the symbols to which scores apply
111    * @param values
112    *          Pairwise scores indexed according to the symbol alphabet
113    */
114   public ScoreMatrix(String theName, char[] alphabet, float[][] values)
115   {
116     this(theName, null, alphabet, values);
117   }
118
119   /**
120    * Constructor given a name, description, symbol alphabet, and matrix of
121    * scores for pairs of symbols. The matrix should be square and of the same
122    * size as the alphabet, for example 20x20 for a 20 symbol alphabet.
123    * 
124    * @param theName
125    *          Unique, human readable name for the matrix
126    * @param theDescription
127    *          descriptive display name suitable for use in menus
128    * @param alphabet
129    *          the symbols to which scores apply
130    * @param values
131    *          Pairwise scores indexed according to the symbol alphabet
132    */
133   public ScoreMatrix(String theName, String theDescription, char[] alphabet,
134           float[][] values)
135   {
136     if (alphabet.length != values.length)
137     {
138       throw new IllegalArgumentException(
139               "score matrix size must match alphabet size");
140     }
141     for (float[] row : values)
142     {
143       if (row.length != alphabet.length)
144       {
145         throw new IllegalArgumentException(
146                 "score matrix size must be square");
147       }
148     }
149
150     this.matrix = values;
151     this.name = theName;
152     this.description = theDescription;
153     this.symbols = alphabet;
154
155     symbolIndex = buildSymbolIndex(alphabet);
156
157     findMinMax();
158
159     /*
160      * crude heuristic for now...
161      */
162     peptide = alphabet.length >= 20;
163   }
164
165   /**
166    * Record the minimum and maximum score values
167    */
168   protected void findMinMax()
169   {
170     float min = Float.MAX_VALUE;
171     float max = -Float.MAX_VALUE;
172     if (matrix != null)
173     {
174       for (float[] row : matrix)
175       {
176         if (row != null)
177         {
178           for (float f : row)
179           {
180             min = Math.min(min, f);
181             max = Math.max(max, f);
182           }
183         }
184       }
185     }
186     minValue = min;
187     maxValue = max;
188   }
189
190   /**
191    * Returns an array A where A[i] is the position in the alphabet array of the
192    * character whose value is i. For example if the alphabet is { 'A', 'D', 'X'
193    * } then A['D'] = A[68] = 1.
194    * <p>
195    * Unmapped characters (not in the alphabet) get an index of -1.
196    * <p>
197    * Mappings are added automatically for lower case symbols (for non case
198    * sensitive scoring), unless they are explicitly present in the alphabet (are
199    * scored separately in the score matrix).
200    * <p>
201    * the gap character (space, dash or dot) included in the alphabet (if any) is
202    * recorded in a field
203    * 
204    * @param alphabet
205    * @return
206    */
207   short[] buildSymbolIndex(char[] alphabet)
208   {
209     short[] index = new short[MAX_ASCII + 1];
210     Arrays.fill(index, UNMAPPED);
211     short pos = 0;
212     for (char c : alphabet)
213     {
214       if (c <= MAX_ASCII)
215       {
216         index[c] = pos;
217       }
218
219       /*
220        * also map lower-case character (unless separately mapped)
221        */
222       if (c >= 'A' && c <= 'Z')
223       {
224         short lowerCase = (short) (c + ('a' - 'A'));
225         if (index[lowerCase] == UNMAPPED)
226         {
227           index[lowerCase] = pos;
228         }
229       }
230       pos++;
231     }
232     return index;
233   }
234
235   @Override
236   public String getName()
237   {
238     return name;
239   }
240
241   @Override
242   public String getDescription()
243   {
244     return description;
245   }
246
247   @Override
248   public boolean isDNA()
249   {
250     return !peptide;
251   }
252
253   @Override
254   public boolean isProtein()
255   {
256     return peptide;
257   }
258
259   /**
260    * Returns a copy of the score matrix as used in getPairwiseScore. If using
261    * this matrix directly, callers <em>must</em> also call
262    * <code>getMatrixIndex</code> in order to get the matrix index for each
263    * character (symbol).
264    * 
265    * @return
266    * @see #getMatrixIndex(char)
267    */
268   public float[][] getMatrix()
269   {
270     float[][] v = new float[matrix.length][matrix.length];
271     for (int i = 0; i < matrix.length; i++)
272     {
273       v[i] = Arrays.copyOf(matrix[i], matrix[i].length);
274     }
275     return v;
276   }
277
278   /**
279    * Answers the matrix index for a given character, or -1 if unmapped in the
280    * matrix. Use this method only if using <code>getMatrix</code> in order to
281    * compute scores directly (without symbol lookup) for efficiency.
282    * 
283    * @param c
284    * @return
285    * @see #getMatrix()
286    */
287   public int getMatrixIndex(char c)
288   {
289     if (c < symbolIndex.length)
290     {
291       return symbolIndex[c];
292     }
293     else
294     {
295       return UNMAPPED;
296     }
297   }
298
299   /**
300    * Returns the pairwise score for substituting c with d. If either c or d is
301    * an unexpected character, returns 1 for identity (c == d), else the minimum
302    * score value in the matrix.
303    */
304   @Override
305   public float getPairwiseScore(char c, char d)
306   {
307     if (c >= symbolIndex.length)
308     {
309       System.err.println(String.format(BAD_ASCII_ERROR, c));
310       return 0;
311     }
312     if (d >= symbolIndex.length)
313     {
314       System.err.println(String.format(BAD_ASCII_ERROR, d));
315       return 0;
316     }
317
318     int cIndex = symbolIndex[c];
319     int dIndex = symbolIndex[d];
320     if (cIndex != UNMAPPED && dIndex != UNMAPPED)
321     {
322       return matrix[cIndex][dIndex];
323     }
324
325     /*
326      * one or both symbols not found in the matrix
327      * currently scoring as 1 (for identity) or the minimum
328      * matrix score value (otherwise)
329      * (a case could be made for using minimum row/column value instead)
330      */
331     return c == d ? UNKNOWN_IDENTITY_SCORE : getMinimumScore();
332   }
333
334   /**
335    * pretty print the matrix
336    */
337   @Override
338   public String toString()
339   {
340     return outputMatrix(false);
341   }
342
343   /**
344    * Print the score matrix, optionally formatted as html, with the alphabet
345    * symbols as column headings and at the start of each row.
346    * <p>
347    * The non-html format should give an output which can be parsed as a score
348    * matrix file
349    * 
350    * @param html
351    * @return
352    */
353   public String outputMatrix(boolean html)
354   {
355     StringBuilder sb = new StringBuilder(512);
356
357     /*
358      * heading row with alphabet
359      */
360     if (html)
361     {
362       sb.append("<table border=\"1\">");
363       sb.append(html ? "<tr><th></th>" : "");
364     }
365     else
366     {
367       sb.append("ScoreMatrix ").append(getName()).append("\n");
368     }
369     for (char sym : symbols)
370     {
371       if (html)
372       {
373         sb.append("<th>&nbsp;").append(sym).append("&nbsp;</th>");
374       }
375       else
376       {
377         sb.append("\t").append(sym);
378       }
379     }
380     sb.append(html ? "</tr>\n" : "\n");
381
382     /*
383      * table of scores
384      */
385     for (char c1 : symbols)
386     {
387       if (html)
388       {
389         sb.append("<tr><td>");
390       }
391       sb.append(c1).append(html ? "</td>" : "");
392       for (char c2 : symbols)
393       {
394         sb.append(html ? "<td>" : "\t")
395                 .append(matrix[symbolIndex[c1]][symbolIndex[c2]])
396                 .append(html ? "</td>" : "");
397       }
398       sb.append(html ? "</tr>\n" : "\n");
399     }
400     if (html)
401     {
402       sb.append("</table>");
403     }
404     return sb.toString();
405   }
406
407   /**
408    * Answers the number of symbols coded for (also equal to the number of rows
409    * and columns of the score matrix)
410    * 
411    * @return
412    */
413   public int getSize()
414   {
415     return symbols.length;
416   }
417
418   /**
419    * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
420    * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
421    * computed using the current score matrix. For example
422    * <ul>
423    * <li>Sequences:</li>
424    * <li>FKL</li>
425    * <li>R-D</li>
426    * <li>QIA</li>
427    * <li>GWC</li>
428    * <li>Score matrix is BLOSUM62</li>
429    * <li>Gaps treated same as X (unknown)</li>
430    * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
431    * <li>product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10</li>
432    * <li>product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13</li>
433    * <li>product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26</li>
434    * <li>product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
435    * <li>and so on</li>
436    * </ul>
437    * This method is thread-safe.
438    */
439   @Override
440   public MatrixI findSimilarities(AlignmentView seqstrings,
441           SimilarityParamsI options)
442   {
443     char gapChar = scoreGapAsAny ? (seqstrings.isNa() ? 'N' : 'X')
444             : GAP_CHARACTER;
445     String[] seqs = seqstrings.getSequenceStrings(gapChar);
446     return findSimilarities(seqs, options);
447   }
448
449   /**
450    * Computes pairwise similarities of a set of sequences using the given
451    * parameters
452    * 
453    * @param seqs
454    * @param params
455    * @return
456    */
457   protected MatrixI findSimilarities(String[] seqs,
458           SimilarityParamsI params)
459   {
460     double[][] values = new double[seqs.length][];
461     for (int row = 0; row < seqs.length; row++)
462     {
463       values[row] = new double[seqs.length];
464       for (int col = 0; col < seqs.length; col++)
465       {
466         double total = computeSimilarity(seqs[row], seqs[col], params);
467         values[row][col] = total;
468       }
469     }
470     return new Matrix(values);
471   }
472
473   /**
474    * Calculates the pairwise similarity of two strings using the given
475    * calculation parameters
476    * 
477    * @param seq1
478    * @param seq2
479    * @param params
480    * @return
481    */
482   protected double computeSimilarity(String seq1, String seq2,
483           SimilarityParamsI params)
484   {
485     int len1 = seq1.length();
486     int len2 = seq2.length();
487     double total = 0;
488
489     int width = Math.max(len1, len2);
490     for (int i = 0; i < width; i++)
491     {
492       if (i >= len1 || i >= len2)
493       {
494         /*
495          * off the end of one sequence; stop if we are only matching
496          * on the shorter sequence length, else treat as trailing gap
497          */
498         if (params.denominateByShortestLength())
499         {
500           break;
501         }
502       }
503
504       char c1 = i >= len1 ? GAP_CHARACTER : seq1.charAt(i);
505       char c2 = i >= len2 ? GAP_CHARACTER : seq2.charAt(i);
506       boolean gap1 = Comparison.isGap(c1);
507       boolean gap2 = Comparison.isGap(c2);
508
509       if (gap1 && gap2)
510       {
511         /*
512          * gap-gap: include if options say so, else ignore
513          */
514         if (!params.includeGappedColumns())
515         {
516           continue;
517         }
518       }
519       else if (gap1 || gap2)
520       {
521         /*
522          * gap-residue: score if options say so
523          */
524         if (!params.includeGaps())
525         {
526           continue;
527         }
528       }
529       float score = getPairwiseScore(c1, c2);
530       total += score;
531     }
532     return total;
533   }
534
535   /**
536    * Answers a hashcode computed from the symbol alphabet and the matrix score
537    * values
538    */
539   @Override
540   public int hashCode()
541   {
542     int hs = Arrays.hashCode(symbols);
543     for (float[] row : matrix)
544     {
545       hs = hs * 31 + Arrays.hashCode(row);
546     }
547     return hs;
548   }
549
550   /**
551    * Answers true if the argument is a ScoreMatrix with the same symbol alphabet
552    * and score values, else false
553    */
554   @Override
555   public boolean equals(Object obj)
556   {
557     if (!(obj instanceof ScoreMatrix))
558     {
559       return false;
560     }
561     ScoreMatrix sm = (ScoreMatrix) obj;
562     if (Arrays.equals(symbols, sm.symbols)
563             && Arrays.deepEquals(matrix, sm.matrix))
564     {
565       return true;
566     }
567     return false;
568   }
569
570   /**
571    * Returns the alphabet the matrix scores for, as a string of characters
572    * 
573    * @return
574    */
575   String getSymbols()
576   {
577     return new String(symbols);
578   }
579
580   public float getMinimumScore()
581   {
582     return minValue;
583   }
584
585   public float getMaximumScore()
586   {
587     return maxValue;
588   }
589
590   @Override
591   public ScoreModelI getInstance(AlignmentViewPanel avp)
592   {
593     return this;
594   }
595 }