src/jalview/analysis/PCA.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;
  22
  23 import jalview.math.Matrix;
  24 import jalview.math.MatrixI;
  25 import jalview.schemes.ResidueProperties;
  26 import jalview.schemes.ScoreMatrix;
  27
  28 import java.io.PrintStream;
  29
  30 /**
  31  * Performs Principal Component Analysis on given sequences
  32  */
  33 public class PCA implements Runnable
  34 {
  35   boolean jvCalcMode = true;
  36
  37   MatrixI symm;
  38
  39   double[] eigenvalue;
  40
  41   MatrixI eigenvector;
  42
  43   StringBuilder details = new StringBuilder(1024);
  44
  45   private String[] seqs;
  46
  47   private ScoreMatrix scoreMatrix;
  48
  49   /**
  50    * Creates a new PCA object. By default, uses blosum62 matrix to generate
  51    * sequence similarity matrices
  52    *
  53    * @param s
  54    *          Set of amino acid sequences to perform PCA on
  55    */
  56   public PCA(String[] s)
  57   {
  58     this(s, false);
  59   }
  60
  61   /**
  62    * Creates a new PCA object. By default, uses blosum62 matrix to generate
  63    * sequence similarity matrices
  64    *
  65    * @param s
  66    *          Set of sequences to perform PCA on
  67    * @param nucleotides
  68    *          if true, uses standard DNA/RNA matrix for sequence similarity
  69    *          calculation.
  70    */
  71   public PCA(String[] s, boolean nucleotides)
  72   {
  73     this(s, nucleotides, null);
  74   }
  75
  76   public PCA(String[] s, boolean nucleotides, String s_m)
  77   {
  78     this.seqs = s;
  79
  80     // BinarySequence[] bs = new BinarySequence[s.length];
  81     // int ii = 0;
  82     //
  83     // while ((ii < s.length) && (s[ii] != null))
  84     // {
  85     // bs[ii] = new BinarySequence(s[ii], nucleotides);
  86     // bs[ii].encode();
  87     // ii++;
  88     // }
  89     //
  90     // BinarySequence[] bs2 = new BinarySequence[s.length];
  91     scoreMatrix = null;
  92     String sm = s_m;
  93     if (sm != null)
  94     {
  95       scoreMatrix = ResidueProperties.getScoreMatrix(sm);
  96     }
  97     if (scoreMatrix == null)
  98     {
  99       // either we were given a non-existent score matrix or a scoremodel that
 100       // isn't based on a pairwise symbol score matrix
 101       scoreMatrix = ResidueProperties
 102               .getScoreMatrix(sm = (nucleotides ? "DNA" : "BLOSUM62"));
 103     }
 104     details.append("PCA calculation using " + sm
 105             + " sequence similarity matrix\n========\n\n");
 106     // ii = 0;
 107     // while ((ii < s.length) && (s[ii] != null))
 108     // {
 109     // bs2[ii] = new BinarySequence(s[ii], nucleotides);
 110     // if (scoreMatrix != null)
 111     // {
 112     // try
 113     // {
 114     // bs2[ii].matrixEncode(scoreMatrix);
 115     // } catch (InvalidSequenceTypeException x)
 116     // {
 117     // details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
 118     // }
 119     // }
 120     // ii++;
 121     // }
 122     //
 123     // int count = 0;
 124     // while ((count < bs.length) && (bs[count] != null))
 125     // {
 126     // count++;
 127     // }
 128     //
 129     // double[][] seqmat = new double[count][];
 130     // double[][] seqmat2 = new double[count][];
 131     //
 132     // int i = 0;
 133     // while (i < count)
 134     // {
 135     // seqmat[i] = bs[i].getDBinary();
 136     // seqmat2[i] = bs2[i].getDBinary();
 137     // i++;
 138     // }
 139     //
 140     // /*
 141     // * using a SparseMatrix to hold the encoded sequences matrix
 142     // * greatly speeds up matrix multiplication as these are mostly zero
 143     // */
 144     // m = new SparseMatrix(seqmat);
 145     // m2 = new Matrix(seqmat2);
 146
 147   }
 148
 149   /**
 150    * Returns Eigenvalue
 151    *
 152    * @param i
 153    *          Index of diagonal within matrix
 154    *
 155    * @return Returns value of diagonal from matrix
 156    */
 157   public double getEigenvalue(int i)
 158   {
 159     return eigenvector.getD()[i];
 160   }
 161
 162   /**
 163    * DOCUMENT ME!
 164    *
 165    * @param l
 166    *          DOCUMENT ME!
 167    * @param n
 168    *          DOCUMENT ME!
 169    * @param mm
 170    *          DOCUMENT ME!
 171    * @param factor
 172    *          DOCUMENT ME!
 173    *
 174    * @return DOCUMENT ME!
 175    */
 176   public float[][] getComponents(int l, int n, int mm, float factor)
 177   {
 178     float[][] out = new float[getHeight()][3];
 179
 180     for (int i = 0; i < getHeight(); i++)
 181     {
 182       out[i][0] = (float) component(i, l) * factor;
 183       out[i][1] = (float) component(i, n) * factor;
 184       out[i][2] = (float) component(i, mm) * factor;
 185     }
 186
 187     return out;
 188   }
 189
 190   /**
 191    * DOCUMENT ME!
 192    *
 193    * @param n
 194    *          DOCUMENT ME!
 195    *
 196    * @return DOCUMENT ME!
 197    */
 198   public double[] component(int n)
 199   {
 200     // n = index of eigenvector
 201     double[] out = new double[getHeight()];
 202
 203     for (int i = 0; i < out.length; i++)
 204     {
 205       out[i] = component(i, n);
 206     }
 207
 208     return out;
 209   }
 210
 211   /**
 212    * DOCUMENT ME!
 213    *
 214    * @param row
 215    *          DOCUMENT ME!
 216    * @param n
 217    *          DOCUMENT ME!
 218    *
 219    * @return DOCUMENT ME!
 220    */
 221   double component(int row, int n)
 222   {
 223     double out = 0.0;
 224
 225     for (int i = 0; i < symm.width(); i++)
 226     {
 227       out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
 228     }
 229
 230     return out / eigenvector.getD()[n];
 231   }
 232
 233   public String getDetails()
 234   {
 235     return details.toString();
 236   }
 237
 238   /**
 239    * DOCUMENT ME!
 240    */
 241   @Override
 242   public void run()
 243   {
 244     PrintStream ps = new PrintStream(System.out)
 245     {
 246       @Override
 247       public void print(String x)
 248       {
 249         details.append(x);
 250       }
 251
 252       @Override
 253       public void println()
 254       {
 255         details.append("\n");
 256       }
 257     };
 258
 259     // long now = System.currentTimeMillis();
 260     try
 261     {
 262       details.append("PCA Calculation Mode is "
 263               + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
 264               + "\n");
 265
 266       // MatrixI mt = m.transpose();
 267       // eigenvector = mt.preMultiply(jvCalcMode ? m2 : m);
 268       eigenvector = computePairwiseScores();
 269
 270       details.append(" --- OrigT * Orig ---- \n");
 271       eigenvector.print(ps, "%8.2f");
 272
 273       symm = eigenvector.copy();
 274
 275       eigenvector.tred();
 276
 277       details.append(" ---Tridiag transform matrix ---\n");
 278       details.append(" --- D vector ---\n");
 279       eigenvector.printD(ps, "%15.4e");
 280       ps.println();
 281       details.append("--- E vector ---\n");
 282       eigenvector.printE(ps, "%15.4e");
 283       ps.println();
 284
 285       // Now produce the diagonalization matrix
 286       eigenvector.tqli();
 287     } catch (Exception q)
 288     {
 289       q.printStackTrace();
 290       details.append("\n*** Unexpected exception when performing PCA ***\n"
 291               + q.getLocalizedMessage());
 292       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 293     }
 294
 295     details.append(" --- New diagonalization matrix ---\n");
 296     eigenvector.print(ps, "%8.2f");
 297     details.append(" --- Eigenvalues ---\n");
 298     eigenvector.printD(ps, "%15.4e");
 299     ps.println();
 300     /*
 301      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 302      * (int ev=0;ev<symm.rows; ev++) {
 303      *
 304      * ps.print(","+component(seq, ev)); } ps.println(); }
 305      */
 306     // System.out.println(("PCA.run() took "
 307     // + (System.currentTimeMillis() - now) + "ms"));
 308   }
 309
 310   /**
 311    * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
 312    * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
 313    * computed using the current score matrix. For example
 314    * <ul>
 315    * <li>Sequences:</li>
 316    * <li>FKL</li>
 317    * <li>RSD</li>
 318    * <li>QIA</li>
 319    * <li>GWC</li>
 320    * <li>Score matrix is BLOSUM62</li>
 321    * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
 322    * <li>product [2, 1] = R.R + S.S + D.D = 5 + 4 + 6 = 15</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.S + L.D = -3 + 0 + -3 = -7
 326    * <li>and so on</li>
 327    * </ul>
 328    */
 329   MatrixI computePairwiseScores()
 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           int score = scoreMatrix.getPairwiseScore(c1, c2);
 344           total += score;
 345         }
 346         values[row][col] = total;
 347       }
 348     }
 349     return new Matrix(values);
 350   }
 351
 352   public void setJvCalcMode(boolean calcMode)
 353   {
 354     this.jvCalcMode = calcMode;
 355   }
 356
 357   /**
 358    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 359    * sequences involved in the pairwise score calculation.
 360    *
 361    * @return
 362    */
 363   public int getHeight()
 364   {
 365     // TODO can any of seqs[] be null?
 366     return seqs.length;
 367   }
 368 }