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