src/jalview/schemes/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.schemes;
  22
  23 import jalview.analysis.scoremodels.PairwiseSeqScoreModel;
  24 import jalview.math.Matrix;
  25 import jalview.math.MatrixI;
  26
  27 public class ScoreMatrix extends PairwiseSeqScoreModel
  28 {
  29   String name;
  30
  31   @Override
  32   public String getName()
  33   {
  34     return name;
  35   }
  36
  37   /**
  38    * reference to integer score matrix
  39    */
  40   int[][] matrix;
  41
  42   /**
  43    * 0 for Protein Score matrix. 1 for dna score matrix
  44    */
  45   int type;
  46
  47   /**
  48    *
  49    * @param name
  50    *          Unique, human readable name for the matrix
  51    * @param matrix
  52    *          Pairwise scores indexed according to appropriate symbol alphabet
  53    * @param type
  54    *          0 for Protein, 1 for NA
  55    */
  56   ScoreMatrix(String name, int[][] matrix, int type)
  57   {
  58     this.matrix = matrix;
  59     this.type = type;
  60     this.name = name;
  61   }
  62
  63   @Override
  64   public boolean isDNA()
  65   {
  66     return type == 1;
  67   }
  68
  69   @Override
  70   public boolean isProtein()
  71   {
  72     return type == 0;
  73   }
  74
  75   @Override
  76   public int[][] getMatrix()
  77   {
  78     return matrix;
  79   }
  80
  81   /**
  82    * Answers the score for substituting first char in A1 with first char in A2
  83    *
  84    * @param A1
  85    * @param A2
  86    * @return
  87    */
  88   public int getPairwiseScore(String A1, String A2)
  89   {
  90     return getPairwiseScore(A1.charAt(0), A2.charAt(0));
  91   }
  92
  93   @Override
  94   public int getPairwiseScore(char c, char d)
  95   {
  96     int pog = 0;
  97
  98     try
  99     {
 100       int a = (type == 0) ? ResidueProperties.aaIndex[c]
 101               : ResidueProperties.nucleotideIndex[c];
 102       int b = (type == 0) ? ResidueProperties.aaIndex[d]
 103               : ResidueProperties.nucleotideIndex[d];
 104
 105       /*
 106        * FIXME: 2.10.1 PCA treats gap as [22] or 'X', but Tree
 107        * calculation treats as [23]; which is correct?
 108        */
 109       /*
 110        * hack to convert unassigned / unknown (including gap)
 111        * to index of unknown (X for amino acids, N for nucleotide)
 112        * TODO: statically assign gap characters to this index?
 113        */
 114 //      if (type == 0)
 115 //      {
 116 //        if (a == ResidueProperties.maxProteinIndex)
 117 //        {
 118 //          a = ResidueProperties.aaIndex['X'];
 119 //        }
 120 //        if (b == ResidueProperties.maxProteinIndex)
 121 //        {
 122 //          b = ResidueProperties.aaIndex['X'];
 123 //        }
 124 //      }
 125 //      if (type != 0)
 126 //      {
 127 //        if (a == ResidueProperties.maxNucleotideIndex)
 128 //        {
 129 //          a = ResidueProperties.nucleotideIndex['N'];
 130 //        }
 131 //        if (b == ResidueProperties.maxNucleotideIndex)
 132 //        {
 133 //          b = ResidueProperties.nucleotideIndex['N'];
 134 //        }
 135 //      }
 136       pog = matrix[a][b];
 137     } catch (Exception e)
 138     {
 139       // System.out.println("Unknown residue in " + A1 + " " + A2);
 140     }
 141
 142     return pog;
 143   }
 144
 145   /**
 146    * pretty print the matrix
 147    */
 148   @Override
 149   public String toString()
 150   {
 151     return outputMatrix(false);
 152   }
 153
 154   public String outputMatrix(boolean html)
 155   {
 156     StringBuffer sb = new StringBuffer();
 157     int[] symbols = (type == 0) ? ResidueProperties.aaIndex
 158             : ResidueProperties.nucleotideIndex;
 159     int symMax = (type == 0) ? ResidueProperties.maxProteinIndex
 160             : ResidueProperties.maxNucleotideIndex;
 161     boolean header = true;
 162     if (html)
 163     {
 164       sb.append("<table border=\"1\">");
 165     }
 166     for (char sym = 'A'; sym <= 'Z'; sym++)
 167     {
 168       if (symbols[sym] >= 0 && symbols[sym] < symMax)
 169       {
 170         if (header)
 171         {
 172           sb.append(html ? "<tr><td></td>" : "");
 173           for (char sym2 = 'A'; sym2 <= 'Z'; sym2++)
 174           {
 175             if (symbols[sym2] >= 0 && symbols[sym2] < symMax)
 176             {
 177               sb.append((html ? "<td>&nbsp;" : "\t") + sym2
 178                       + (html ? "&nbsp;</td>" : ""));
 179             }
 180           }
 181           header = false;
 182           sb.append(html ? "</tr>\n" : "\n");
 183         }
 184         if (html)
 185         {
 186           sb.append("<tr>");
 187         }
 188         sb.append((html ? "<td>" : "") + sym + (html ? "</td>" : ""));
 189         for (char sym2 = 'A'; sym2 <= 'Z'; sym2++)
 190         {
 191           if (symbols[sym2] >= 0 && symbols[sym2] < symMax)
 192           {
 193             sb.append((html ? "<td>" : "\t")
 194                     + matrix[symbols[sym]][symbols[sym2]]
 195                     + (html ? "</td>" : ""));
 196           }
 197         }
 198         sb.append(html ? "</tr>\n" : "\n");
 199       }
 200     }
 201     if (html)
 202     {
 203       sb.append("</table>");
 204     }
 205     return sb.toString();
 206   }
 207
 208   /**
 209    * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
 210    * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
 211    * computed using the current score matrix. For example
 212    * <ul>
 213    * <li>Sequences:</li>
 214    * <li>FKL</li>
 215    * <li>R-D</li>
 216    * <li>QIA</li>
 217    * <li>GWC</li>
 218    * <li>Score matrix is BLOSUM62</li>
 219    * <li>Gaps treated same as X (unknown)</li>
 220    * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
 221    * <li>product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10</li>
 222    * <li>product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13</li>
 223    * <li>product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26</li>
 224    * <li>product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
 225    * <li>and so on</li>
 226    * </ul>
 227    */
 228   public MatrixI computePairwiseScores(String[] seqs)
 229   {
 230     double[][] values = new double[seqs.length][];
 231     for (int row = 0; row < seqs.length; row++)
 232     {
 233       values[row] = new double[seqs.length];
 234       for (int col = 0; col < seqs.length; col++)
 235       {
 236         int total = 0;
 237         int width = Math.min(seqs[row].length(), seqs[col].length());
 238         for (int i = 0; i < width; i++)
 239         {
 240           char c1 = seqs[row].charAt(i);
 241           char c2 = seqs[col].charAt(i);
 242           int score = getPairwiseScore(c1, c2);
 243           total += score;
 244         }
 245         values[row][col] = total;
 246       }
 247     }
 248     return new Matrix(values);
 249   }
 250 }