src/jalview/analysis/PCA.java

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