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