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.analysis.scoremodels.PairwiseDistanceModel;
  24 import jalview.analysis.scoremodels.ScoreMatrix;
  25 import jalview.analysis.scoremodels.ScoreModels;
  26 import jalview.datamodel.BinarySequence;
  27 import jalview.datamodel.BinarySequence.InvalidSequenceTypeException;
  28 import jalview.math.Matrix;
  29
  30 import java.io.PrintStream;
  31
  32 /**
  33  * Performs Principal Component Analysis on given sequences
  34  *
  35  * @author $author$
  36  * @version $Revision$
  37  */
  38 public class PCA implements Runnable
  39 {
  40   Matrix m;
  41
  42   Matrix symm;
  43
  44   Matrix m2;
  45
  46   double[] eigenvalue;
  47
  48   Matrix eigenvector;
  49
  50   StringBuffer details = new StringBuffer();
  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     ii = 0;
  94     ScoreMatrix smtrx = null;
  95     String sm = s_m;
  96     if (sm != null)
  97     {
  98       smtrx = (ScoreMatrix) ((PairwiseDistanceModel) ScoreModels
  99               .getInstance()
 100               .forName(sm)).getScoreModel();
 101     }
 102     if (smtrx == null)
 103     {
 104       // either we were given a non-existent score matrix or a scoremodel that
 105       // isn't based on a pairwise symbol score matrix
 106       smtrx = ScoreModels.getInstance().getDefaultModel(!nucleotides);
 107     }
 108     details.append("PCA calculation using " + sm
 109             + " sequence similarity matrix\n========\n\n");
 110     while ((ii < s.length) && (s[ii] != null))
 111     {
 112       bs2[ii] = new BinarySequence(s[ii], nucleotides);
 113       if (smtrx != null)
 114       {
 115         try
 116         {
 117           bs2[ii].matrixEncode(smtrx);
 118         } catch (InvalidSequenceTypeException x)
 119         {
 120           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
 121         }
 122       }
 123       ii++;
 124     }
 125
 126     // System.out.println("Created binary encoding");
 127     // printMemory(rt);
 128     int count = 0;
 129
 130     while ((count < bs.length) && (bs[count] != null))
 131     {
 132       count++;
 133     }
 134
 135     double[][] seqmat = new double[count][bs[0].getDBinary().length];
 136     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
 137     int i = 0;
 138
 139     while (i < count)
 140     {
 141       seqmat[i] = bs[i].getDBinary();
 142       seqmat2[i] = bs2[i].getDBinary();
 143       i++;
 144     }
 145
 146     // System.out.println("Created array");
 147     // printMemory(rt);
 148     // System.out.println(" --- Original matrix ---- ");
 149     m = new Matrix(seqmat);
 150     m2 = new Matrix(seqmat2);
 151
 152   }
 153
 154   /**
 155    * Returns the matrix used in PCA calculation
 156    *
 157    * @return java.math.Matrix object
 158    */
 159
 160   public Matrix getM()
 161   {
 162     return m;
 163   }
 164
 165   /**
 166    * Returns Eigenvalue
 167    *
 168    * @param i
 169    *          Index of diagonal within matrix
 170    *
 171    * @return Returns value of diagonal from matrix
 172    */
 173   public double getEigenvalue(int i)
 174   {
 175     return eigenvector.d[i];
 176   }
 177
 178   /**
 179    * DOCUMENT ME!
 180    *
 181    * @param l
 182    *          DOCUMENT ME!
 183    * @param n
 184    *          DOCUMENT ME!
 185    * @param mm
 186    *          DOCUMENT ME!
 187    * @param factor
 188    *          DOCUMENT ME!
 189    *
 190    * @return DOCUMENT ME!
 191    */
 192   public float[][] getComponents(int l, int n, int mm, float factor)
 193   {
 194     float[][] out = new float[m.rows][3];
 195
 196     for (int i = 0; i < m.rows; i++)
 197     {
 198       out[i][0] = (float) component(i, l) * factor;
 199       out[i][1] = (float) component(i, n) * factor;
 200       out[i][2] = (float) component(i, mm) * factor;
 201     }
 202
 203     return out;
 204   }
 205
 206   /**
 207    * DOCUMENT ME!
 208    *
 209    * @param n
 210    *          DOCUMENT ME!
 211    *
 212    * @return DOCUMENT ME!
 213    */
 214   public double[] component(int n)
 215   {
 216     // n = index of eigenvector
 217     double[] out = new double[m.rows];
 218
 219     for (int i = 0; i < m.rows; i++)
 220     {
 221       out[i] = component(i, n);
 222     }
 223
 224     return out;
 225   }
 226
 227   /**
 228    * DOCUMENT ME!
 229    *
 230    * @param row
 231    *          DOCUMENT ME!
 232    * @param n
 233    *          DOCUMENT ME!
 234    *
 235    * @return DOCUMENT ME!
 236    */
 237   double component(int row, int n)
 238   {
 239     double out = 0.0;
 240
 241     for (int i = 0; i < symm.cols; i++)
 242     {
 243       out += (symm.value[row][i] * eigenvector.value[i][n]);
 244     }
 245
 246     return out / eigenvector.d[n];
 247   }
 248
 249   public String getDetails()
 250   {
 251     return details.toString();
 252   }
 253
 254   /**
 255    * DOCUMENT ME!
 256    */
 257   @Override
 258   public void run()
 259   {
 260     PrintStream ps = new PrintStream(System.out)
 261     {
 262       @Override
 263       public void print(String x)
 264       {
 265         details.append(x);
 266       }
 267
 268       @Override
 269       public void println()
 270       {
 271         details.append("\n");
 272       }
 273     };
 274
 275     try
 276     {
 277       details.append("PCA Calculation Mode is "
 278               + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
 279               + "\n");
 280       Matrix mt = m.transpose();
 281
 282       details.append(" --- OrigT * Orig ---- \n");
 283       if (!jvCalcMode)
 284       {
 285         eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
 286       }
 287       else
 288       {
 289         eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace
 290                                           // method
 291       }
 292
 293       eigenvector.print(ps);
 294
 295       symm = eigenvector.copy();
 296
 297       eigenvector.tred();
 298
 299       details.append(" ---Tridiag transform matrix ---\n");
 300       details.append(" --- D vector ---\n");
 301       eigenvector.printD(ps);
 302       ps.println();
 303       details.append("--- E vector ---\n");
 304       eigenvector.printE(ps);
 305       ps.println();
 306
 307       // Now produce the diagonalization matrix
 308       eigenvector.tqli();
 309     } catch (Exception q)
 310     {
 311       q.printStackTrace();
 312       details.append("\n*** Unexpected exception when performing PCA ***\n"
 313               + q.getLocalizedMessage());
 314       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 315     }
 316
 317     details.append(" --- New diagonalization matrix ---\n");
 318     eigenvector.print(ps);
 319     details.append(" --- Eigenvalues ---\n");
 320     eigenvector.printD(ps);
 321     ps.println();
 322     /*
 323      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 324      * (int ev=0;ev<symm.rows; ev++) {
 325      *
 326      * ps.print(","+component(seq, ev)); } ps.println(); }
 327      */
 328   }
 329
 330   boolean jvCalcMode = true;
 331
 332   public void setJvCalcMode(boolean calcMode)
 333   {
 334     this.jvCalcMode = calcMode;
 335   }
 336 }