src/jalview/analysis/PCA.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2)
   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
  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   public PCA(String[] s, boolean nucleotides, String s_m)
  78   {
  79
  80     BinarySequence[] bs = new BinarySequence[s.length];
  81     int ii = 0;
  82
  83     while ((ii < s.length) && (s[ii] != null))
  84     {
  85       bs[ii] = new BinarySequence(s[ii], nucleotides);
  86       bs[ii].encode();
  87       ii++;
  88     }
  89
  90     BinarySequence[] bs2 = new BinarySequence[s.length];
  91     ii = 0;
  92     ScoreMatrix smtrx = null;
  93     String sm=s_m;
  94     if (sm!=null)
  95     {
  96       smtrx = ResidueProperties.getScoreMatrix(sm);
  97     }
  98     if (smtrx==null)
  99     {
 100       // either we were given a non-existent score matrix or a scoremodel that isn't based on a pairwise symbol score matrix
 101       smtrx = ResidueProperties.getScoreMatrix(sm=(nucleotides ? "DNA" : "BLOSUM62"));
 102     }
 103     details.append("PCA calculation using " + sm
 104             + " sequence similarity matrix\n========\n\n");
 105     while ((ii < s.length) && (s[ii] != null))
 106     {
 107       bs2[ii] = new BinarySequence(s[ii], nucleotides);
 108       if (smtrx != null)
 109       {
 110         try
 111         {
 112           bs2[ii].matrixEncode(smtrx);
 113         } catch (InvalidSequenceTypeException x)
 114         {
 115           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
 116         }
 117       }
 118       ii++;
 119     }
 120
 121     // System.out.println("Created binary encoding");
 122     // printMemory(rt);
 123     int count = 0;
 124
 125     while ((count < bs.length) && (bs[count] != null))
 126     {
 127       count++;
 128     }
 129
 130     double[][] seqmat = new double[count][bs[0].getDBinary().length];
 131     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
 132     int i = 0;
 133
 134     while (i < count)
 135     {
 136       seqmat[i] = bs[i].getDBinary();
 137       seqmat2[i] = bs2[i].getDBinary();
 138       i++;
 139     }
 140
 141     // System.out.println("Created array");
 142     // printMemory(rt);
 143     // System.out.println(" --- Original matrix ---- ");
 144     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
 145     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
 146
 147   }
 148
 149   /**
 150    * Returns the matrix used in PCA calculation
 151    *
 152    * @return java.math.Matrix object
 153    */
 154
 155   public Matrix getM()
 156   {
 157     return m;
 158   }
 159
 160   /**
 161    * Returns Eigenvalue
 162    *
 163    * @param i
 164    *          Index of diagonal within matrix
 165    *
 166    * @return Returns value of diagonal from matrix
 167    */
 168   public double getEigenvalue(int i)
 169   {
 170     return eigenvector.d[i];
 171   }
 172
 173   /**
 174    * DOCUMENT ME!
 175    *
 176    * @param l
 177    *          DOCUMENT ME!
 178    * @param n
 179    *          DOCUMENT ME!
 180    * @param mm
 181    *          DOCUMENT ME!
 182    * @param factor
 183    *          DOCUMENT ME!
 184    *
 185    * @return DOCUMENT ME!
 186    */
 187   public float[][] getComponents(int l, int n, int mm, float factor)
 188   {
 189     float[][] out = new float[m.rows][3];
 190
 191     for (int i = 0; i < m.rows; i++)
 192     {
 193       out[i][0] = (float) component(i, l) * factor;
 194       out[i][1] = (float) component(i, n) * factor;
 195       out[i][2] = (float) component(i, mm) * factor;
 196     }
 197
 198     return out;
 199   }
 200
 201   /**
 202    * DOCUMENT ME!
 203    *
 204    * @param n
 205    *          DOCUMENT ME!
 206    *
 207    * @return DOCUMENT ME!
 208    */
 209   public double[] component(int n)
 210   {
 211     // n = index of eigenvector
 212     double[] out = new double[m.rows];
 213
 214     for (int i = 0; i < m.rows; i++)
 215     {
 216       out[i] = component(i, n);
 217     }
 218
 219     return out;
 220   }
 221
 222   /**
 223    * DOCUMENT ME!
 224    *
 225    * @param row
 226    *          DOCUMENT ME!
 227    * @param n
 228    *          DOCUMENT ME!
 229    *
 230    * @return DOCUMENT ME!
 231    */
 232   double component(int row, int n)
 233   {
 234     double out = 0.0;
 235
 236     for (int i = 0; i < symm.cols; i++)
 237     {
 238       out += (symm.value[row][i] * eigenvector.value[i][n]);
 239     }
 240
 241     return out / eigenvector.d[n];
 242   }
 243
 244   public String getDetails()
 245   {
 246     return details.toString();
 247   }
 248
 249   /**
 250    * DOCUMENT ME!
 251    */
 252   public void run()
 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     try {
 268     details.append("PCA Calculation Mode is "
 269             + (jvCalcMode ? "Jalview variant" : "Original SeqSpace") + "\n");
 270     Matrix mt = m.transpose();
 271
 272     details.append(" --- OrigT * Orig ---- \n");
 273     if (!jvCalcMode)
 274     {
 275       eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
 276     }
 277     else
 278     {
 279       eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace method
 280     }
 281
 282     eigenvector.print(ps);
 283
 284     symm = eigenvector.copy();
 285
 286     eigenvector.tred();
 287
 288     details.append(" ---Tridiag transform matrix ---\n");
 289     details.append(" --- D vector ---\n");
 290     eigenvector.printD(ps);
 291     ps.println();
 292     details.append("--- E vector ---\n");
 293     eigenvector.printE(ps);
 294     ps.println();
 295
 296     // Now produce the diagonalization matrix
 297     eigenvector.tqli();
 298     } catch (Exception q)
 299     {
 300       q.printStackTrace();
 301       details.append("\n*** Unexpected exception when performing PCA ***\n"+q.getLocalizedMessage());
 302       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 303     }
 304
 305     details.append(" --- New diagonalization matrix ---\n");
 306     eigenvector.print(ps);
 307     details.append(" --- Eigenvalues ---\n");
 308     eigenvector.printD(ps);
 309     ps.println();
 310     /*
 311      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 312      * (int ev=0;ev<symm.rows; ev++) {
 313      *
 314      * ps.print(","+component(seq, ev)); } ps.println(); }
 315      */
 316   }
 317
 318   boolean jvCalcMode = true;
 319
 320   public void setJvCalcMode(boolean calcMode)
 321   {
 322     this.jvCalcMode = calcMode;
 323   }
 324 }