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