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