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.math.*;
  24
  25 /**
  26  * Performs Principal Component Analysis on given sequences
  27  *
  28  * @author $author$
  29  * @version $Revision$
  30  */
  31 public class PCA implements Runnable
  32 {
  33   Matrix m;
  34
  35   Matrix symm;
  36
  37   Matrix m2;
  38
  39   double[] eigenvalue;
  40
  41   Matrix eigenvector;
  42
  43   StringBuffer details = new StringBuffer();
  44
  45   /**
  46    * Creates a new PCA object.
  47    *
  48    * @param s
  49    *          Set of sequences to perform PCA on
  50    */
  51   public PCA(String[] s)
  52   {
  53
  54     BinarySequence[] bs = new BinarySequence[s.length];
  55     int ii = 0;
  56
  57     while ((ii < s.length) && (s[ii] != null))
  58     {
  59       bs[ii] = new BinarySequence(s[ii]);
  60       bs[ii].encode();
  61       ii++;
  62     }
  63
  64     BinarySequence[] bs2 = new BinarySequence[s.length];
  65     ii = 0;
  66
  67     while ((ii < s.length) && (s[ii] != null))
  68     {
  69       bs2[ii] = new BinarySequence(s[ii]);
  70       bs2[ii].blosumEncode();
  71       ii++;
  72     }
  73
  74     // System.out.println("Created binary encoding");
  75     // printMemory(rt);
  76     int count = 0;
  77
  78     while ((count < bs.length) && (bs[count] != null))
  79     {
  80       count++;
  81     }
  82
  83     double[][] seqmat = new double[count][bs[0].getDBinary().length];
  84     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
  85     int i = 0;
  86
  87     while (i < count)
  88     {
  89       seqmat[i] = bs[i].getDBinary();
  90       seqmat2[i] = bs2[i].getDBinary();
  91       i++;
  92     }
  93
  94     // System.out.println("Created array");
  95     // printMemory(rt);
  96     // System.out.println(" --- Original matrix ---- ");
  97     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
  98     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
  99
 100   }
 101
 102   /**
 103    * Returns the matrix used in PCA calculation
 104    *
 105    * @return java.math.Matrix object
 106    */
 107
 108   public Matrix getM()
 109   {
 110     return m;
 111   }
 112
 113   /**
 114    * Returns Eigenvalue
 115    *
 116    * @param i
 117    *          Index of diagonal within matrix
 118    *
 119    * @return Returns value of diagonal from matrix
 120    */
 121   public double getEigenvalue(int i)
 122   {
 123     return eigenvector.d[i];
 124   }
 125
 126   /**
 127    * DOCUMENT ME!
 128    *
 129    * @param l
 130    *          DOCUMENT ME!
 131    * @param n
 132    *          DOCUMENT ME!
 133    * @param mm
 134    *          DOCUMENT ME!
 135    * @param factor
 136    *          DOCUMENT ME!
 137    *
 138    * @return DOCUMENT ME!
 139    */
 140   public float[][] getComponents(int l, int n, int mm, float factor)
 141   {
 142     float[][] out = new float[m.rows][3];
 143
 144     for (int i = 0; i < m.rows; i++)
 145     {
 146       out[i][0] = (float) component(i, l) * factor;
 147       out[i][1] = (float) component(i, n) * factor;
 148       out[i][2] = (float) component(i, mm) * factor;
 149     }
 150
 151     return out;
 152   }
 153
 154   /**
 155    * DOCUMENT ME!
 156    *
 157    * @param n
 158    *          DOCUMENT ME!
 159    *
 160    * @return DOCUMENT ME!
 161    */
 162   public double[] component(int n)
 163   {
 164     // n = index of eigenvector
 165     double[] out = new double[m.rows];
 166
 167     for (int i = 0; i < m.rows; i++)
 168     {
 169       out[i] = component(i, n);
 170     }
 171
 172     return out;
 173   }
 174
 175   /**
 176    * DOCUMENT ME!
 177    *
 178    * @param row
 179    *          DOCUMENT ME!
 180    * @param n
 181    *          DOCUMENT ME!
 182    *
 183    * @return DOCUMENT ME!
 184    */
 185   double component(int row, int n)
 186   {
 187     double out = 0.0;
 188
 189     for (int i = 0; i < symm.cols; i++)
 190     {
 191       out += (symm.value[row][i] * eigenvector.value[i][n]);
 192     }
 193
 194     return out / eigenvector.d[n];
 195   }
 196
 197   public String getDetails()
 198   {
 199     return details.toString();
 200   }
 201
 202   /**
 203    * DOCUMENT ME!
 204    */
 205   public void run()
 206   {
 207     Matrix mt = m.transpose();
 208
 209     details.append(" --- OrigT * Orig ---- \n");
 210     // eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
 211     eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace method
 212
 213     PrintStream ps = new PrintStream(System.out)
 214     {
 215       public void print(String x)
 216       {
 217         details.append(x);
 218       }
 219
 220       public void println()
 221       {
 222         details.append("\n");
 223       }
 224     };
 225
 226     eigenvector.print(ps);
 227
 228     symm = eigenvector.copy();
 229
 230     eigenvector.tred();
 231
 232     details.append(" ---Tridiag transform matrix ---\n");
 233     details.append(" --- D vector ---\n");
 234     eigenvector.printD(ps);
 235     ps.println();
 236     details.append("--- E vector ---\n");
 237     eigenvector.printE(ps);
 238     ps.println();
 239
 240     // Now produce the diagonalization matrix
 241     eigenvector.tqli();
 242
 243     details.append(" --- New diagonalization matrix ---\n");
 244     eigenvector.print(ps);
 245     details.append(" --- Eigenvalues ---\n");
 246     eigenvector.printD(ps);
 247     ps.println();
 248     /*
 249      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 250      * (int ev=0;ev<symm.rows; ev++) {
 251      *
 252      * ps.print(","+component(seq, ev)); } ps.println(); }
 253      */
 254   }
 255 }