src/jalview/analysis/PCA.java

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