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