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