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.analysis.scoremodels.PairwiseDistanceModel;
  24 import jalview.analysis.scoremodels.ScoreMatrix;
  25 import jalview.analysis.scoremodels.ScoreModels;
  26 import jalview.math.MatrixI;
  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     scoreMatrix = null;
  81     String sm = s_m;
  82     if (sm != null)
  83     {
  84       scoreMatrix = (ScoreMatrix) ((PairwiseDistanceModel) ScoreModels
  85               .getInstance()
  86               .forName(sm)).getScoreModel();
  87     }
  88     if (scoreMatrix == null)
  89     {
  90       // either we were given a non-existent score matrix or a scoremodel that
  91       // isn't based on a pairwise symbol score matrix
  92       scoreMatrix = ScoreModels.getInstance().getDefaultModel(!nucleotides);
  93     }
  94     details.append("PCA calculation using " + sm
  95             + " sequence similarity matrix\n========\n\n");
  96   }
  97
  98   /**
  99    * Returns Eigenvalue
 100    *
 101    * @param i
 102    *          Index of diagonal within matrix
 103    *
 104    * @return Returns value of diagonal from matrix
 105    */
 106   public double getEigenvalue(int i)
 107   {
 108     return eigenvector.getD()[i];
 109   }
 110
 111   /**
 112    * DOCUMENT ME!
 113    *
 114    * @param l
 115    *          DOCUMENT ME!
 116    * @param n
 117    *          DOCUMENT ME!
 118    * @param mm
 119    *          DOCUMENT ME!
 120    * @param factor
 121    *          DOCUMENT ME!
 122    *
 123    * @return DOCUMENT ME!
 124    */
 125   public float[][] getComponents(int l, int n, int mm, float factor)
 126   {
 127     float[][] out = new float[getHeight()][3];
 128
 129     for (int i = 0; i < getHeight(); i++)
 130     {
 131       out[i][0] = (float) component(i, l) * factor;
 132       out[i][1] = (float) component(i, n) * factor;
 133       out[i][2] = (float) component(i, mm) * factor;
 134     }
 135
 136     return out;
 137   }
 138
 139   /**
 140    * DOCUMENT ME!
 141    *
 142    * @param n
 143    *          DOCUMENT ME!
 144    *
 145    * @return DOCUMENT ME!
 146    */
 147   public double[] component(int n)
 148   {
 149     // n = index of eigenvector
 150     double[] out = new double[getHeight()];
 151
 152     for (int i = 0; i < out.length; i++)
 153     {
 154       out[i] = component(i, n);
 155     }
 156
 157     return out;
 158   }
 159
 160   /**
 161    * DOCUMENT ME!
 162    *
 163    * @param row
 164    *          DOCUMENT ME!
 165    * @param n
 166    *          DOCUMENT ME!
 167    *
 168    * @return DOCUMENT ME!
 169    */
 170   double component(int row, int n)
 171   {
 172     double out = 0.0;
 173
 174     for (int i = 0; i < symm.width(); i++)
 175     {
 176       out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
 177     }
 178
 179     return out / eigenvector.getD()[n];
 180   }
 181
 182   public String getDetails()
 183   {
 184     return details.toString();
 185   }
 186
 187   /**
 188    * DOCUMENT ME!
 189    */
 190   @Override
 191   public void run()
 192   {
 193     PrintStream ps = new PrintStream(System.out)
 194     {
 195       @Override
 196       public void print(String x)
 197       {
 198         details.append(x);
 199       }
 200
 201       @Override
 202       public void println()
 203       {
 204         details.append("\n");
 205       }
 206     };
 207
 208     // long now = System.currentTimeMillis();
 209     try
 210     {
 211       details.append("PCA Calculation Mode is "
 212               + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
 213               + "\n");
 214
 215       eigenvector = scoreMatrix.computePairwiseScores(seqs);
 216
 217       details.append(" --- OrigT * Orig ---- \n");
 218       eigenvector.print(ps, "%8.2f");
 219
 220       symm = eigenvector.copy();
 221
 222       eigenvector.tred();
 223
 224       details.append(" ---Tridiag transform matrix ---\n");
 225       details.append(" --- D vector ---\n");
 226       eigenvector.printD(ps, "%15.4e");
 227       ps.println();
 228       details.append("--- E vector ---\n");
 229       eigenvector.printE(ps, "%15.4e");
 230       ps.println();
 231
 232       // Now produce the diagonalization matrix
 233       eigenvector.tqli();
 234     } catch (Exception q)
 235     {
 236       q.printStackTrace();
 237       details.append("\n*** Unexpected exception when performing PCA ***\n"
 238               + q.getLocalizedMessage());
 239       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 240     }
 241
 242     details.append(" --- New diagonalization matrix ---\n");
 243     eigenvector.print(ps, "%8.2f");
 244     details.append(" --- Eigenvalues ---\n");
 245     eigenvector.printD(ps, "%15.4e");
 246     ps.println();
 247     /*
 248      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 249      * (int ev=0;ev<symm.rows; ev++) {
 250      *
 251      * ps.print(","+component(seq, ev)); } ps.println(); }
 252      */
 253     // System.out.println(("PCA.run() took "
 254     // + (System.currentTimeMillis() - now) + "ms"));
 255   }
 256
 257   public void setJvCalcMode(boolean calcMode)
 258   {
 259     this.jvCalcMode = calcMode;
 260   }
 261
 262   /**
 263    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 264    * sequences involved in the pairwise score calculation.
 265    *
 266    * @return
 267    */
 268   public int getHeight()
 269   {
 270     // TODO can any of seqs[] be null?
 271     return seqs.length;
 272   }
 273 }