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