src/jalview/analysis/PCA.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8)
   3  * Copyright (C) 2012 J Procter, AM Waterhouse, LM Lui, J Engelhardt, G Barton, M Clamp, S Searle
   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 of the License, or (at your option) any later version.
  10  *
  11  * Jalview is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14  * PURPOSE.  See the GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  17  */
  18 package jalview.analysis;
  19
  20 import java.io.*;
  21
  22 import jalview.datamodel.*;
  23 import jalview.datamodel.BinarySequence.InvalidSequenceTypeException;
  24 import jalview.math.*;
  25 import jalview.schemes.ResidueProperties;
  26 import jalview.schemes.ScoreMatrix;
  27
  28 /**
  29  * Performs Principal Component Analysis on given sequences
  30  *
  31  * @author $author$
  32  * @version $Revision$
  33  */
  34 public class PCA implements Runnable
  35 {
  36   Matrix m;
  37
  38   Matrix symm;
  39
  40   Matrix m2;
  41
  42   double[] eigenvalue;
  43
  44   Matrix eigenvector;
  45
  46   StringBuffer details = new StringBuffer();
  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
  73     BinarySequence[] bs = new BinarySequence[s.length];
  74     int ii = 0;
  75
  76     while ((ii < s.length) && (s[ii] != null))
  77     {
  78       bs[ii] = new BinarySequence(s[ii], nucleotides);
  79       bs[ii].encode();
  80       ii++;
  81     }
  82
  83     BinarySequence[] bs2 = new BinarySequence[s.length];
  84     ii = 0;
  85
  86     String sm = nucleotides ? "DNA" : "BLOSUM62";
  87     ScoreMatrix smtrx = ResidueProperties.getScoreMatrix(sm);
  88     details.append("PCA calculation using " + sm
  89             + " sequence similarity matrix\n========\n\n");
  90     while ((ii < s.length) && (s[ii] != null))
  91     {
  92       bs2[ii] = new BinarySequence(s[ii], nucleotides);
  93       if (smtrx != null)
  94       {
  95         try
  96         {
  97           bs2[ii].matrixEncode(smtrx);
  98         } catch (InvalidSequenceTypeException x)
  99         {
 100           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
 101         }
 102       }
 103       ii++;
 104     }
 105
 106     // System.out.println("Created binary encoding");
 107     // printMemory(rt);
 108     int count = 0;
 109
 110     while ((count < bs.length) && (bs[count] != null))
 111     {
 112       count++;
 113     }
 114
 115     double[][] seqmat = new double[count][bs[0].getDBinary().length];
 116     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
 117     int i = 0;
 118
 119     while (i < count)
 120     {
 121       seqmat[i] = bs[i].getDBinary();
 122       seqmat2[i] = bs2[i].getDBinary();
 123       i++;
 124     }
 125
 126     // System.out.println("Created array");
 127     // printMemory(rt);
 128     // System.out.println(" --- Original matrix ---- ");
 129     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
 130     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
 131
 132   }
 133
 134   /**
 135    * Returns the matrix used in PCA calculation
 136    *
 137    * @return java.math.Matrix object
 138    */
 139
 140   public Matrix getM()
 141   {
 142     return m;
 143   }
 144
 145   /**
 146    * Returns Eigenvalue
 147    *
 148    * @param i
 149    *          Index of diagonal within matrix
 150    *
 151    * @return Returns value of diagonal from matrix
 152    */
 153   public double getEigenvalue(int i)
 154   {
 155     return eigenvector.d[i];
 156   }
 157
 158   /**
 159    * DOCUMENT ME!
 160    *
 161    * @param l
 162    *          DOCUMENT ME!
 163    * @param n
 164    *          DOCUMENT ME!
 165    * @param mm
 166    *          DOCUMENT ME!
 167    * @param factor
 168    *          DOCUMENT ME!
 169    *
 170    * @return DOCUMENT ME!
 171    */
 172   public float[][] getComponents(int l, int n, int mm, float factor)
 173   {
 174     float[][] out = new float[m.rows][3];
 175
 176     for (int i = 0; i < m.rows; i++)
 177     {
 178       out[i][0] = (float) component(i, l) * factor;
 179       out[i][1] = (float) component(i, n) * factor;
 180       out[i][2] = (float) component(i, mm) * factor;
 181     }
 182
 183     return out;
 184   }
 185
 186   /**
 187    * DOCUMENT ME!
 188    *
 189    * @param n
 190    *          DOCUMENT ME!
 191    *
 192    * @return DOCUMENT ME!
 193    */
 194   public double[] component(int n)
 195   {
 196     // n = index of eigenvector
 197     double[] out = new double[m.rows];
 198
 199     for (int i = 0; i < m.rows; i++)
 200     {
 201       out[i] = component(i, n);
 202     }
 203
 204     return out;
 205   }
 206
 207   /**
 208    * DOCUMENT ME!
 209    *
 210    * @param row
 211    *          DOCUMENT ME!
 212    * @param n
 213    *          DOCUMENT ME!
 214    *
 215    * @return DOCUMENT ME!
 216    */
 217   double component(int row, int n)
 218   {
 219     double out = 0.0;
 220
 221     for (int i = 0; i < symm.cols; i++)
 222     {
 223       out += (symm.value[row][i] * eigenvector.value[i][n]);
 224     }
 225
 226     return out / eigenvector.d[n];
 227   }
 228
 229   public String getDetails()
 230   {
 231     return details.toString();
 232   }
 233
 234   /**
 235    * DOCUMENT ME!
 236    */
 237   public void run()
 238   {
 239     details.append("PCA Calculation Mode is "
 240             + (jvCalcMode ? "Jalview variant" : "Original SeqSpace") + "\n");
 241     Matrix mt = m.transpose();
 242
 243     details.append(" --- OrigT * Orig ---- \n");
 244     if (!jvCalcMode)
 245     {
 246       eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
 247     }
 248     else
 249     {
 250       eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace method
 251     }
 252
 253     PrintStream ps = new PrintStream(System.out)
 254     {
 255       public void print(String x)
 256       {
 257         details.append(x);
 258       }
 259
 260       public void println()
 261       {
 262         details.append("\n");
 263       }
 264     };
 265
 266     eigenvector.print(ps);
 267
 268     symm = eigenvector.copy();
 269
 270     eigenvector.tred();
 271
 272     details.append(" ---Tridiag transform matrix ---\n");
 273     details.append(" --- D vector ---\n");
 274     eigenvector.printD(ps);
 275     ps.println();
 276     details.append("--- E vector ---\n");
 277     eigenvector.printE(ps);
 278     ps.println();
 279
 280     // Now produce the diagonalization matrix
 281     eigenvector.tqli();
 282
 283     details.append(" --- New diagonalization matrix ---\n");
 284     eigenvector.print(ps);
 285     details.append(" --- Eigenvalues ---\n");
 286     eigenvector.printD(ps);
 287     ps.println();
 288     /*
 289      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 290      * (int ev=0;ev<symm.rows; ev++) {
 291      *
 292      * ps.print(","+component(seq, ev)); } ps.println(); }
 293      */
 294   }
 295
 296   boolean jvCalcMode = true;
 297
 298   public void setJvCalcMode(boolean calcMode)
 299   {
 300     this.jvCalcMode = calcMode;
 301   }
 302 }