src/jalview/analysis/PCA.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.7)
   3  * Copyright (C) 2011 J Procter, AM Waterhouse, J Engelhardt, LM Lui, 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.
  50    * By default, uses blosum62 matrix to generate sequence similarity matrices
  51    * @param s
  52    *          Set of amino acid sequences to perform PCA on
  53    */
  54   public PCA(String[] s)
  55   {
  56     this(s, false);
  57   }
  58
  59   /**
  60    * Creates a new PCA object.
  61    * By default, uses blosum62 matrix to generate sequence similarity matrices
  62    * @param s
  63    *          Set of sequences to perform PCA on
  64    *          @param nucleotides if true, uses standard DNA/RNA matrix for sequence similarity calculation.
  65    */
  66   public PCA(String[] s, boolean nucleotides)
  67   {
  68
  69     BinarySequence[] bs = new BinarySequence[s.length];
  70     int ii = 0;
  71
  72     while ((ii < s.length) && (s[ii] != null))
  73     {
  74       bs[ii] = new BinarySequence(s[ii],nucleotides);
  75       bs[ii].encode();
  76       ii++;
  77     }
  78
  79     BinarySequence[] bs2 = new BinarySequence[s.length];
  80     ii = 0;
  81
  82     String sm=nucleotides ? "DNA" : "BLOSUM62";
  83     ScoreMatrix smtrx=ResidueProperties.getScoreMatrix(sm);
  84     details.append("PCA calculation using "+sm+" sequence similarity matrix\n========\n\n");
  85
  86     while ((ii < s.length) && (s[ii] != null))
  87     {
  88       bs2[ii] = new BinarySequence(s[ii], nucleotides);
  89       if (smtrx != null)
  90       {
  91         try
  92         {
  93           bs2[ii].matrixEncode(smtrx);
  94         } catch (InvalidSequenceTypeException x)
  95         {
  96           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
  97         }
  98       }
  99       ii++;
 100     }
 101
 102     // System.out.println("Created binary encoding");
 103     // printMemory(rt);
 104     int count = 0;
 105
 106     while ((count < bs.length) && (bs[count] != null))
 107     {
 108       count++;
 109     }
 110
 111     double[][] seqmat = new double[count][bs[0].getDBinary().length];
 112     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
 113     int i = 0;
 114
 115     while (i < count)
 116     {
 117       seqmat[i] = bs[i].getDBinary();
 118       seqmat2[i] = bs2[i].getDBinary();
 119       i++;
 120     }
 121
 122     // System.out.println("Created array");
 123     // printMemory(rt);
 124     // System.out.println(" --- Original matrix ---- ");
 125     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
 126     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
 127
 128   }
 129
 130   /**
 131    * Returns the matrix used in PCA calculation
 132    *
 133    * @return java.math.Matrix object
 134    */
 135
 136   public Matrix getM()
 137   {
 138     return m;
 139   }
 140
 141   /**
 142    * Returns Eigenvalue
 143    *
 144    * @param i
 145    *          Index of diagonal within matrix
 146    *
 147    * @return Returns value of diagonal from matrix
 148    */
 149   public double getEigenvalue(int i)
 150   {
 151     return eigenvector.d[i];
 152   }
 153
 154   /**
 155    * DOCUMENT ME!
 156    *
 157    * @param l
 158    *          DOCUMENT ME!
 159    * @param n
 160    *          DOCUMENT ME!
 161    * @param mm
 162    *          DOCUMENT ME!
 163    * @param factor
 164    *          DOCUMENT ME!
 165    *
 166    * @return DOCUMENT ME!
 167    */
 168   public float[][] getComponents(int l, int n, int mm, float factor)
 169   {
 170     float[][] out = new float[m.rows][3];
 171
 172     for (int i = 0; i < m.rows; i++)
 173     {
 174       out[i][0] = (float) component(i, l) * factor;
 175       out[i][1] = (float) component(i, n) * factor;
 176       out[i][2] = (float) component(i, mm) * factor;
 177     }
 178
 179     return out;
 180   }
 181
 182   /**
 183    * DOCUMENT ME!
 184    *
 185    * @param n
 186    *          DOCUMENT ME!
 187    *
 188    * @return DOCUMENT ME!
 189    */
 190   public double[] component(int n)
 191   {
 192     // n = index of eigenvector
 193     double[] out = new double[m.rows];
 194
 195     for (int i = 0; i < m.rows; i++)
 196     {
 197       out[i] = component(i, n);
 198     }
 199
 200     return out;
 201   }
 202
 203   /**
 204    * DOCUMENT ME!
 205    *
 206    * @param row
 207    *          DOCUMENT ME!
 208    * @param n
 209    *          DOCUMENT ME!
 210    *
 211    * @return DOCUMENT ME!
 212    */
 213   double component(int row, int n)
 214   {
 215     double out = 0.0;
 216
 217     for (int i = 0; i < symm.cols; i++)
 218     {
 219       out += (symm.value[row][i] * eigenvector.value[i][n]);
 220     }
 221
 222     return out / eigenvector.d[n];
 223   }
 224
 225   public String getDetails()
 226   {
 227     return details.toString();
 228   }
 229
 230   /**
 231    * DOCUMENT ME!
 232    */
 233   public void run()
 234   {
 235     details.append("PCA Calculation Mode is "+(jvCalcMode ? "Jalview variant" : "Original SeqSpace")+"\n");
 236     Matrix mt = m.transpose();
 237
 238     details.append(" --- OrigT * Orig ---- \n");
 239     if (!jvCalcMode)
 240     {
 241       eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
 242     }
 243     else
 244     {
 245       eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace method
 246     }
 247
 248     PrintStream ps = new PrintStream(System.out)
 249     {
 250       public void print(String x)
 251       {
 252         details.append(x);
 253       }
 254
 255       public void println()
 256       {
 257         details.append("\n");
 258       }
 259     };
 260
 261     eigenvector.print(ps);
 262
 263     symm = eigenvector.copy();
 264
 265     eigenvector.tred();
 266
 267     details.append(" ---Tridiag transform matrix ---\n");
 268     details.append(" --- D vector ---\n");
 269     eigenvector.printD(ps);
 270     ps.println();
 271     details.append("--- E vector ---\n");
 272     eigenvector.printE(ps);
 273     ps.println();
 274
 275     // Now produce the diagonalization matrix
 276     eigenvector.tqli();
 277
 278     details.append(" --- New diagonalization matrix ---\n");
 279     eigenvector.print(ps);
 280     details.append(" --- Eigenvalues ---\n");
 281     eigenvector.printD(ps);
 282     ps.println();
 283     /*
 284      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 285      * (int ev=0;ev<symm.rows; ev++) {
 286      *
 287      * ps.print(","+component(seq, ev)); } ps.println(); }
 288      */
 289   }
 290
 291   boolean jvCalcMode = true;
 292
 293   public void setJvCalcMode(boolean calcMode)
 294   {
 295     this.jvCalcMode = calcMode;
 296   }
 297 }