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 symbols as column headings and at the start of each row
 235    * @param html
 236    * @return
 237    */
 238   public String outputMatrix(boolean html)
 239   {
 240     StringBuilder sb = new StringBuilder(512);
 241
 242     /*
 243      * heading row with alphabet
 244      */
 245     if (html)
 246     {
 247       sb.append("<table border=\"1\">");
 248       sb.append(html ? "<tr><th></th>" : "");
 249     }
 250     for (char sym : symbols)
 251     {
 252       if (html)
 253       {
 254         sb.append("<th>&nbsp;").append(sym).append("&nbsp;</th>");
 255       }
 256       else
 257       {
 258         sb.append("\t").append(sym);
 259       }
 260     }
 261     sb.append(html ? "</tr>\n" : "\n");
 262
 263     /*
 264      * table of scores
 265      */
 266     for (char c1 : symbols)
 267     {
 268       if (html)
 269       {
 270         sb.append("<tr><td>");
 271       }
 272       sb.append(c1).append(html ? "</td>" : "");
 273       for (char c2 : symbols)
 274       {
 275         sb.append(html ? "<td>" : "\t")
 276                 .append(matrix[symbolIndex[c1]][symbolIndex[c2]])
 277                 .append(html ? "</td>" : "");
 278       }
 279       sb.append(html ? "</tr>\n" : "\n");
 280     }
 281     if (html)
 282     {
 283       sb.append("</table>");
 284     }
 285     return sb.toString();
 286   }
 287
 288   /**
 289    * Answers the number of symbols coded for (also equal to the number of rows
 290    * and columns of the score matrix)
 291    *
 292    * @return
 293    */
 294   public int getSize()
 295   {
 296     return symbols.length;
 297   }
 298
 299   /**
 300    * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
 301    * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
 302    * computed using the current score matrix. For example
 303    * <ul>
 304    * <li>Sequences:</li>
 305    * <li>FKL</li>
 306    * <li>R-D</li>
 307    * <li>QIA</li>
 308    * <li>GWC</li>
 309    * <li>Score matrix is BLOSUM62</li>
 310    * <li>Gaps treated same as X (unknown)</li>
 311    * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
 312    * <li>product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10</li>
 313    * <li>product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13</li>
 314    * <li>product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26</li>
 315    * <li>product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
 316    * <li>and so on</li>
 317    * </ul>
 318    */
 319   public MatrixI computePairwiseScores(String[] seqs)
 320   {
 321     double[][] values = new double[seqs.length][];
 322     for (int row = 0; row < seqs.length; row++)
 323     {
 324       values[row] = new double[seqs.length];
 325       for (int col = 0; col < seqs.length; col++)
 326       {
 327         int total = 0;
 328         int width = Math.min(seqs[row].length(), seqs[col].length());
 329         for (int i = 0; i < width; i++)
 330         {
 331           char c1 = seqs[row].charAt(i);
 332           char c2 = seqs[col].charAt(i);
 333           float score = getPairwiseScore(c1, c2);
 334           total += score;
 335         }
 336         values[row][col] = total;
 337       }
 338     }
 339     return new Matrix(values);
 340   }
 341 }