JAL-2416 ScoreMatrix.getMatrix() returns a defensive copy to keep the
[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.math.Matrix;
24 import jalview.math.MatrixI;
25
26 import java.util.Arrays;
27
28 public class ScoreMatrix implements PairwiseScoreModelI
29 {
30   public static final short UNMAPPED = (short) -1;
31
32   private static final String BAD_ASCII_ERROR = "Unexpected character %s in getPairwiseScore";
33
34   private static final int MAX_ASCII = 127;
35
36   /*
37    * the name of the model as shown in menus
38    */
39   private String name;
40
41   /*
42    * the characters that the model provides scores for
43    */
44   private char[] symbols;
45
46   /*
47    * the score matrix; both dimensions must equal the number of symbols
48    * matrix[i][j] is the substitution score for replacing symbols[i] with symbols[j]
49    */
50   private float[][] matrix;
51
52   /*
53    * quick lookup to convert from an ascii character value to the index
54    * of the corresponding symbol in the score matrix 
55    */
56   private short[] symbolIndex;
57
58   /*
59    * true for Protein Score matrix, false for dna score matrix
60    */
61   private boolean peptide;
62
63   private boolean symmetric;
64
65   /**
66    * Constructor given a name, symbol alphabet, and matrix of scores for pairs
67    * of symbols. The matrix should be square and of the same size as the
68    * alphabet, for example 20x20 for a 20 symbol alphabet.
69    * 
70    * @param name
71    *          Unique, human readable name for the matrix
72    * @param alphabet
73    *          the symbols to which scores apply
74    * @param matrix
75    *          Pairwise scores indexed according to the symbol alphabet
76    */
77   public ScoreMatrix(String name, char[] alphabet, float[][] matrix)
78   {
79     if (alphabet.length != matrix.length)
80     {
81       throw new IllegalArgumentException(
82               "score matrix size must match alphabet size");
83     }
84     for (float[] row : matrix)
85     {
86       if (row.length != alphabet.length)
87       {
88         throw new IllegalArgumentException(
89                 "score matrix size must be square");
90       }
91     }
92
93     this.matrix = matrix;
94     this.name = name;
95     this.symbols = alphabet;
96
97     symbolIndex = buildSymbolIndex(alphabet);
98
99     /*
100      * crude heuristic for now...
101      */
102     peptide = alphabet.length >= 20;
103   }
104
105   /**
106    * Returns an array A where A[i] is the position in the alphabet array of the
107    * character whose value is i. For example if the alphabet is { 'A', 'D', 'X'
108    * } then A['D'] = A[68] = 1.
109    * <p>
110    * Unmapped characters (not in the alphabet) get an index of -1.
111    * <p>
112    * Mappings are added automatically for lower case symbols (for non case
113    * sensitive scoring), unless they are explicitly present in the alphabet (are
114    * scored separately in the score matrix).
115    * 
116    * @param alphabet
117    * @return
118    */
119   static short[] buildSymbolIndex(char[] alphabet)
120   {
121     short[] index = new short[MAX_ASCII + 1];
122     Arrays.fill(index, UNMAPPED);
123     short pos = 0;
124     for (char c : alphabet)
125     {
126       if (c <= MAX_ASCII)
127       {
128         index[c] = pos;
129       }
130
131       /*
132        * also map lower-case character (unless separately mapped)
133        */
134       if (c >= 'A' && c <= 'Z')
135       {
136         short lowerCase = (short) (c + ('a' - 'A'));
137         if (index[lowerCase] == UNMAPPED)
138         {
139           index[lowerCase] = pos;
140         }
141       }
142       pos++;
143     }
144     return index;
145   }
146
147   @Override
148   public String getName()
149   {
150     return name;
151   }
152
153   @Override
154   public boolean isDNA()
155   {
156     return !peptide;
157   }
158
159   @Override
160   public boolean isProtein()
161   {
162     return peptide;
163   }
164
165   /**
166    * Returns a copy of the score matrix as used in getPairwiseScore. If using
167    * this matrix directly, callers <em>must</em> also call
168    * <code>getMatrixIndex</code> in order to get the matrix index for each
169    * character (symbol).
170    * 
171    * @return
172    * @see #getMatrixIndex(char)
173    */
174   public float[][] getMatrix()
175   {
176     float[][] v = new float[matrix.length][matrix.length];
177     for (int i = 0; i < matrix.length; i++)
178     {
179       v[i] = Arrays.copyOf(matrix[i], matrix[i].length);
180     }
181     return v;
182   }
183
184   /**
185    * Answers the matrix index for a given character, or -1 if unmapped in the
186    * matrix. Use this method only if using <code>getMatrix</code> in order to
187    * compute scores directly (without symbol lookup) for efficiency.
188    * 
189    * @param c
190    * @return
191    * @see #getMatrix()
192    */
193   public int getMatrixIndex(char c)
194   {
195     if (c < symbolIndex.length)
196     {
197       return symbolIndex[c];
198     }
199     else
200     {
201       return UNMAPPED;
202     }
203   }
204
205   /**
206    * Returns the pairwise score for substituting c with d, or zero if c or d is
207    * an unscored or unexpected character
208    */
209   @Override
210   public float getPairwiseScore(char c, char d)
211   {
212     if (c >= symbolIndex.length)
213     {
214       System.err.println(String.format(BAD_ASCII_ERROR, c));
215       return 0;
216     }
217     if (d >= symbolIndex.length)
218     {
219       System.err.println(String.format(BAD_ASCII_ERROR, d));
220       return 0;
221     }
222
223     int cIndex = symbolIndex[c];
224     int dIndex = symbolIndex[d];
225     if (cIndex != UNMAPPED && dIndex != UNMAPPED)
226     {
227       return matrix[cIndex][dIndex];
228     }
229     return 0;
230   }
231
232   /**
233    * pretty print the matrix
234    */
235   @Override
236   public String toString()
237   {
238     return outputMatrix(false);
239   }
240
241   /**
242    * Print the score matrix, optionally formatted as html, with the alphabet
243    * symbols as column headings and at the start of each row.
244    * <p>
245    * The non-html format should give an output which can be parsed as a score
246    * matrix file
247    * 
248    * @param html
249    * @return
250    */
251   public String outputMatrix(boolean html)
252   {
253     StringBuilder sb = new StringBuilder(512);
254
255     /*
256      * heading row with alphabet
257      */
258     if (html)
259     {
260       sb.append("<table border=\"1\">");
261       sb.append(html ? "<tr><th></th>" : "");
262     }
263     else
264     {
265       sb.append("ScoreMatrix ").append(getName()).append("\n");
266       sb.append(symbols).append("\n");
267     }
268     for (char sym : symbols)
269     {
270       if (html)
271       {
272         sb.append("<th>&nbsp;").append(sym).append("&nbsp;</th>");
273       }
274       else
275       {
276         sb.append("\t").append(sym);
277       }
278     }
279     sb.append(html ? "</tr>\n" : "\n");
280
281     /*
282      * table of scores
283      */
284     for (char c1 : symbols)
285     {
286       if (html)
287       {
288         sb.append("<tr><td>");
289       }
290       sb.append(c1).append(html ? "</td>" : "");
291       for (char c2 : symbols)
292       {
293         sb.append(html ? "<td>" : "\t")
294                 .append(matrix[symbolIndex[c1]][symbolIndex[c2]])
295                 .append(html ? "</td>" : "");
296       }
297       sb.append(html ? "</tr>\n" : "\n");
298     }
299     if (html)
300     {
301       sb.append("</table>");
302     }
303     return sb.toString();
304   }
305
306   /**
307    * Answers the number of symbols coded for (also equal to the number of rows
308    * and columns of the score matrix)
309    * 
310    * @return
311    */
312   public int getSize()
313   {
314     return symbols.length;
315   }
316
317   /**
318    * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
319    * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
320    * computed using the current score matrix. For example
321    * <ul>
322    * <li>Sequences:</li>
323    * <li>FKL</li>
324    * <li>R-D</li>
325    * <li>QIA</li>
326    * <li>GWC</li>
327    * <li>Score matrix is BLOSUM62</li>
328    * <li>Gaps treated same as X (unknown)</li>
329    * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
330    * <li>product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10</li>
331    * <li>product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13</li>
332    * <li>product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26</li>
333    * <li>product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
334    * <li>and so on</li>
335    * </ul>
336    */
337   public MatrixI computePairwiseScores(String[] seqs)
338   {
339     double[][] values = new double[seqs.length][];
340     for (int row = 0; row < seqs.length; row++)
341     {
342       values[row] = new double[seqs.length];
343       for (int col = 0; col < seqs.length; col++)
344       {
345         int total = 0;
346         int width = Math.min(seqs[row].length(), seqs[col].length());
347         for (int i = 0; i < width; i++)
348         {
349           char c1 = seqs[row].charAt(i);
350           char c2 = seqs[col].charAt(i);
351           float score = getPairwiseScore(c1, c2);
352           total += score;
353         }
354         values[row][col] = total;
355       }
356     }
357     return new Matrix(values);
358   }
359
360   /**
361    * Answers a hashcode computed from the symbol alphabet and the matrix score
362    * values
363    */
364   @Override
365   public int hashCode()
366   {
367     int hs = Arrays.hashCode(symbols);
368     for (float[] row : matrix)
369     {
370       hs = hs * 31 + Arrays.hashCode(row);
371     }
372     return hs;
373   }
374
375   /**
376    * Answers true if the argument is a ScoreMatrix with the same symbol alphabet
377    * and score values, else false
378    */
379   @Override
380   public boolean equals(Object obj)
381   {
382     if (!(obj instanceof ScoreMatrix))
383     {
384       return false;
385     }
386     ScoreMatrix sm = (ScoreMatrix) obj;
387     if (Arrays.equals(symbols, sm.symbols)
388             && Arrays.deepEquals(matrix, sm.matrix))
389     {
390       return true;
391     }
392     return false;
393   }
394
395   public boolean isSymmetric()
396   {
397     return symmetric;
398   }
399
400   /**
401    * Returns the alphabet the matrix scores for, as a string of characters
402    * 
403    * @return
404    */
405   public String getSymbols()
406   {
407     return new String(symbols);
408   }
409 }