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