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