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.SimilarityParams;
  24 import jalview.api.analysis.DistanceScoreModelI;
  25 import jalview.api.analysis.ScoreModelI;
  26 import jalview.api.analysis.SimilarityScoreModelI;
  27 import jalview.datamodel.AlignmentView;
  28 import jalview.math.MatrixI;
  29
  30 import java.io.PrintStream;
  31
  32 /**
  33  * Performs Principal Component Analysis on given sequences
  34  */
  35 public class PCA implements Runnable
  36 {
  37   boolean jvCalcMode = true;
  38
  39   MatrixI symm;
  40
  41   double[] eigenvalue;
  42
  43   MatrixI eigenvector;
  44
  45   StringBuilder details = new StringBuilder(1024);
  46
  47   private AlignmentView seqs;
  48
  49   private ScoreModelI scoreModel;
  50
  51   public PCA(AlignmentView s, ScoreModelI sm)
  52   {
  53     this.seqs = s;
  54
  55     scoreModel = sm;
  56     details.append("PCA calculation using " + sm.getName()
  57             + " sequence similarity matrix\n========\n\n");
  58   }
  59
  60   /**
  61    * Returns Eigenvalue
  62    *
  63    * @param i
  64    *          Index of diagonal within matrix
  65    *
  66    * @return Returns value of diagonal from matrix
  67    */
  68   public double getEigenvalue(int i)
  69   {
  70     return eigenvector.getD()[i];
  71   }
  72
  73   /**
  74    * DOCUMENT ME!
  75    *
  76    * @param l
  77    *          DOCUMENT ME!
  78    * @param n
  79    *          DOCUMENT ME!
  80    * @param mm
  81    *          DOCUMENT ME!
  82    * @param factor
  83    *          DOCUMENT ME!
  84    *
  85    * @return DOCUMENT ME!
  86    */
  87   public float[][] getComponents(int l, int n, int mm, float factor)
  88   {
  89     float[][] out = new float[getHeight()][3];
  90
  91     for (int i = 0; i < getHeight(); i++)
  92     {
  93       out[i][0] = (float) component(i, l) * factor;
  94       out[i][1] = (float) component(i, n) * factor;
  95       out[i][2] = (float) component(i, mm) * factor;
  96     }
  97
  98     return out;
  99   }
 100
 101   /**
 102    * DOCUMENT ME!
 103    *
 104    * @param n
 105    *          DOCUMENT ME!
 106    *
 107    * @return DOCUMENT ME!
 108    */
 109   public double[] component(int n)
 110   {
 111     // n = index of eigenvector
 112     double[] out = new double[getHeight()];
 113
 114     for (int i = 0; i < out.length; i++)
 115     {
 116       out[i] = component(i, n);
 117     }
 118
 119     return out;
 120   }
 121
 122   /**
 123    * DOCUMENT ME!
 124    *
 125    * @param row
 126    *          DOCUMENT ME!
 127    * @param n
 128    *          DOCUMENT ME!
 129    *
 130    * @return DOCUMENT ME!
 131    */
 132   double component(int row, int n)
 133   {
 134     double out = 0.0;
 135
 136     for (int i = 0; i < symm.width(); i++)
 137     {
 138       out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
 139     }
 140
 141     return out / eigenvector.getD()[n];
 142   }
 143
 144   public String getDetails()
 145   {
 146     return details.toString();
 147   }
 148
 149   /**
 150    * DOCUMENT ME!
 151    */
 152   @Override
 153   public void run()
 154   {
 155     PrintStream ps = new PrintStream(System.out)
 156     {
 157       @Override
 158       public void print(String x)
 159       {
 160         details.append(x);
 161       }
 162
 163       @Override
 164       public void println()
 165       {
 166         details.append("\n");
 167       }
 168     };
 169
 170     // long now = System.currentTimeMillis();
 171     try
 172     {
 173       details.append("PCA Calculation Mode is "
 174               + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
 175               + "\n");
 176
 177       eigenvector = computeSimilarity(seqs);
 178
 179       details.append(" --- OrigT * Orig ---- \n");
 180       eigenvector.print(ps, "%8.2f");
 181
 182       symm = eigenvector.copy();
 183
 184       eigenvector.tred();
 185
 186       details.append(" ---Tridiag transform matrix ---\n");
 187       details.append(" --- D vector ---\n");
 188       eigenvector.printD(ps, "%15.4e");
 189       ps.println();
 190       details.append("--- E vector ---\n");
 191       eigenvector.printE(ps, "%15.4e");
 192       ps.println();
 193
 194       // Now produce the diagonalization matrix
 195       eigenvector.tqli();
 196     } catch (Exception q)
 197     {
 198       q.printStackTrace();
 199       details.append("\n*** Unexpected exception when performing PCA ***\n"
 200               + q.getLocalizedMessage());
 201       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 202     }
 203
 204     details.append(" --- New diagonalization matrix ---\n");
 205     eigenvector.print(ps, "%8.2f");
 206     details.append(" --- Eigenvalues ---\n");
 207     eigenvector.printD(ps, "%15.4e");
 208     ps.println();
 209     /*
 210      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 211      * (int ev=0;ev<symm.rows; ev++) {
 212      *
 213      * ps.print(","+component(seq, ev)); } ps.println(); }
 214      */
 215     // System.out.println(("PCA.run() took "
 216     // + (System.currentTimeMillis() - now) + "ms"));
 217   }
 218
 219   /**
 220    * Computes a pairwise similarity matrix for the given sequence regions using
 221    * the configured score model. If the score model is a similarity model, then
 222    * it computes the result directly. If it is a distance model, then use it to
 223    * compute pairwise distances, and convert these to similarity scores by
 224    * substracting from the maximum value.
 225    *
 226    * @param av
 227    * @return
 228    */
 229   MatrixI computeSimilarity(AlignmentView av)
 230   {
 231     MatrixI result = null;
 232     // TODO pass choice of params from GUI in constructo
 233     if (scoreModel instanceof SimilarityScoreModelI)
 234     {
 235       result = ((SimilarityScoreModelI) scoreModel).findSimilarities(av,
 236               SimilarityParams.SeqSpace);
 237     }
 238     else if (scoreModel instanceof DistanceScoreModelI)
 239     {
 240       result = ((DistanceScoreModelI) scoreModel).findDistances(av,
 241               SimilarityParams.SeqSpace);
 242       result.reverseRange(false);
 243     }
 244     else
 245     {
 246       System.err
 247               .println("Unexpected type of score model, cannot calculate similarity");
 248     }
 249
 250     return result;
 251   }
 252
 253   public void setJvCalcMode(boolean calcMode)
 254   {
 255     this.jvCalcMode = calcMode;
 256   }
 257
 258   /**
 259    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 260    * sequences involved in the pairwise score calculation.
 261    *
 262    * @return
 263    */
 264   public int getHeight()
 265   {
 266     // TODO can any of seqs[] be null?
 267     return seqs.getSequences().length;
 268   }
 269 }