src/jalview/analysis/PaSiMap.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.ScoreModelI;
  24 import jalview.api.analysis.SimilarityParamsI;
  25 import jalview.bin.Console;
  26 import jalview.datamodel.AlignmentView;
  27 import jalview.datamodel.Point;
  28 import jalview.datamodel.SequenceI;
  29 import jalview.gui.PairwiseAlignPanel;
  30 import jalview.gui.PaSiMapPanel;
  31 import jalview.math.Matrix;
  32 import jalview.math.MatrixI;
  33 import jalview.viewmodel.AlignmentViewport;
  34
  35
  36 import java.io.PrintStream;
  37 import java.util.Hashtable;
  38 import java.util.Enumeration;
  39
  40 /**
  41  * Performs Principal Component Analysis on given sequences
  42  * @AUTHOR MorellThomas
  43  */
  44 public class PaSiMap implements Runnable
  45 {
  46   /*
  47    * inputs
  48    */
  49   final private AlignmentViewport seqs;
  50
  51   final private ScoreModelI scoreModel;
  52
  53   final private SimilarityParamsI similarityParams;
  54
  55   final private byte dim = 3;
  56
  57   /*
  58    * outputs
  59    */
  60   private MatrixI pairwiseScores;
  61
  62   private MatrixI tridiagonal;
  63
  64   private MatrixI eigenMatrix;
  65
  66   /**
  67    * Constructor given the sequences to compute for, the similarity model to
  68    * use, and a set of parameters for sequence comparison
  69    *
  70    * @param sequences
  71    * @param sm
  72    * @param options
  73    */
  74   //public PaSiMap(AlignmentView sequences, ScoreModelI sm,
  75   //&! viewport or panel?
  76   public PaSiMap(AlignmentViewport sequences, ScoreModelI sm,
  77           SimilarityParamsI options)
  78   {
  79     this.seqs = sequences;
  80     this.scoreModel = sm;
  81     this.similarityParams = options;
  82   }
  83
  84   /**
  85    * Returns Eigenvalue
  86    *
  87    * @param i
  88    *          Index of diagonal within matrix
  89    *
  90    * @return Returns value of diagonal from matrix
  91    */
  92   public double getEigenvalue(int i)
  93   {
  94     return eigenMatrix.getD()[i];
  95   }
  96
  97   /**
  98    * DOCUMENT ME!
  99    *
 100    * @param l
 101    *          DOCUMENT ME!
 102    * @param n
 103    *          DOCUMENT ME!
 104    * @param mm
 105    *          DOCUMENT ME!
 106    * @param factor
 107    *          DOCUMENT ME!
 108    *
 109    * @return DOCUMENT ME!
 110    */
 111   public Point[] getComponents(int l, int n, int mm, float factor)
 112   {
 113     Point[] out = new Point[getHeight()];
 114
 115     for (int i = 0; i < out.length; i++)
 116     {
 117       float x = (float) component(i, l) * factor;
 118       float y = (float) component(i, n) * factor;
 119       float z = (float) component(i, mm) * factor;
 120       out[i] = new Point(x, y, z);
 121     }
 122     //&!
 123     System.out.println("Points:");
 124     for (Point point : out)
 125     {
 126       System.out.println(point.toString());
 127     }
 128
 129     return out;
 130   }
 131
 132   /**
 133    * DOCUMENT ME!
 134    *
 135    * @param n
 136    *          DOCUMENT ME!
 137    *
 138    * @return DOCUMENT ME!
 139    */
 140   public double[] component(int n)
 141   {
 142     // n = index of eigenvector
 143     double[] out = new double[getHeight()];
 144
 145     for (int i = 0; i < out.length; i++)
 146     {
 147       out[i] = component(i, n);
 148     }
 149
 150     return out;
 151   }
 152
 153   /**
 154    * DOCUMENT ME!
 155    *
 156    * @param row
 157    *          DOCUMENT ME!
 158    * @param n
 159    *          DOCUMENT ME!
 160    *
 161    * @return DOCUMENT ME!
 162    */
 163   double component(int row, int n)
 164   {
 165     /*
 166     double out = 0.0;
 167
 168     for (int i = 0; i < pairwiseScores.width(); i++)
 169     {
 170       double pairwiseScore = pairwiseScores.getValue(row, i);
 171       double eigenScore = eigenMatrix.getValue(i, n);
 172       out += (pairwiseScores.getValue(row, i) * eigenMatrix.getValue(i, n));
 173     }
 174
 175     return out / eigenMatrix.getD()[n];
 176     */
 177     System.out.println(String.format("row %d, col %d", row, n));
 178     return eigenMatrix.getValue(row, n);
 179   }
 180
 181   /**
 182    * Answers a formatted text report of the PaSiMap calculation results (matrices
 183    * and eigenvalues) suitable for display
 184    *
 185    * @return
 186    */
 187   public String getDetails()
 188   {
 189     StringBuilder sb = new StringBuilder(1024);
 190     sb.append("PaSiMap calculation using ").append(scoreModel.getName())
 191             .append(" sequence similarity matrix\n========\n\n");
 192     PrintStream ps = wrapOutputBuffer(sb);
 193
 194     /*
 195      * pairwise similarity scores
 196      */
 197     sb.append(" --- OrigT * Orig ---- \n");
 198     pairwiseScores.print(ps, "%8.2f");
 199
 200     /*
 201      * tridiagonal matrix, with D and E vectors
 202      */
 203     /*
 204     sb.append(" ---Tridiag transform matrix ---\n");
 205     sb.append(" --- D vector ---\n");
 206     tridiagonal.printD(ps, "%15.4e");
 207     ps.println();
 208     sb.append("--- E vector ---\n");
 209     tridiagonal.printE(ps, "%15.4e");
 210     ps.println();
 211     */
 212
 213     /*
 214      * eigenvalues matrix, with D vector
 215      */
 216     sb.append(" --- New diagonalization matrix ---\n");
 217     eigenMatrix.print(ps, "%8.2f");
 218     sb.append(" --- Eigenvalues ---\n");
 219     eigenMatrix.printD(ps, "%15.4e");
 220     ps.println();
 221
 222     return sb.toString();
 223   }
 224
 225   /**
 226    * Performs the PaSiMap calculation
 227    *
 228    * creates a new gui/PairwiseAlignPanel with the input sequences (<++>/AlignmentViewport)
 229    * uses analysis/AlignSeq to creatue the pairwise alignments and calculate the AlignmentScores (float for each pair)
 230    * gets all float[][] scores from the gui/PairwiseAlignPanel
 231    * checks the connections for each sequence with <++>/AlignmentViewport seqs.calculateConnectivity(float[][] scores, int dim) (from analysis/Connectivity) -- throws an Exception if insufficient
 232    * creates a math/MatrixI pairwiseScores of the float[][] scores
 233    * copys the scores and fills the diagonal to create a symmetric matrix using math/Matrix.fillDiagonal()
 234    * performs the analysis/ccAnalysis with the symmetric matrix
 235    * gets the eigenmatrix and the eigenvalues using math/Matrix.tqli()
 236    */
 237   @Override
 238   public void run()
 239   {
 240     try
 241     {
 242       // run needleman regardless if aligned or not
 243       // gui.PairwiseAlignPanel <++>
 244       PairwiseAlignPanel alignment = new PairwiseAlignPanel(seqs);
 245       float[][] scores = alignment.getAlignmentScores();        //bigger index first -- eg scores[14][13]
 246
 247       Hashtable<SequenceI, Integer> connectivity = seqs.calculateConnectivity(scores, dim);
 248
 249       pairwiseScores = new Matrix(scores);
 250       pairwiseScores.fillDiagonal();
 251
 252       ccAnalysis cc = new ccAnalysis(pairwiseScores, dim);
 253       pairwiseScores = cc.run();
 254       tridiagonal = pairwiseScores.copy();
 255       //tridiagonal.tred();
 256
 257       /** perform the eigendecomposition for the plot */
 258       eigenMatrix = tridiagonal.copy();
 259       eigenMatrix.setD(pairwiseScores.getD());
 260       //eigenMatrix.tqli();
 261       System.out.println(getDetails());
 262
 263
 264     } catch (Exception q)
 265     {
 266       Console.error("Error computing PaSiMap:  " + q.getMessage());
 267       q.printStackTrace();
 268     }
 269   }
 270
 271   /**
 272    * Returns a PrintStream that wraps (appends its output to) the given
 273    * StringBuilder
 274    *
 275    * @param sb
 276    * @return
 277    */
 278   protected PrintStream wrapOutputBuffer(StringBuilder sb)
 279   {
 280     PrintStream ps = new PrintStream(System.out)
 281     {
 282       @Override
 283       public void print(String x)
 284       {
 285         sb.append(x);
 286       }
 287
 288       @Override
 289       public void println()
 290       {
 291         sb.append("\n");
 292       }
 293     };
 294     return ps;
 295   }
 296
 297   /**
 298    * Answers the N dimensions of the NxM PaSiMap matrix. This is the number of
 299    * sequences involved in the pairwise score calculation.
 300    *
 301    * @return
 302    */
 303   public int getHeight()
 304   {
 305     // TODO can any of seqs[] be null?
 306     return pairwiseScores.height();// seqs.getSequences().length;
 307   }
 308
 309   /**
 310    * Answers the M dimensions of the NxM PaSiMap matrix. This is the number of
 311    * sequences involved in the pairwise score calculation.
 312    *
 313    * @return
 314    */
 315   public int getWidth()
 316   {
 317     // TODO can any of seqs[] be null?
 318     return pairwiseScores.width();// seqs.getSequences().length;
 319   }
 320
 321   /**
 322    * Answers the sequence pairwise similarity scores which were the first step
 323    * of the PaSiMap calculation
 324    *
 325    * @return
 326    */
 327   public MatrixI getPairwiseScores()
 328   {
 329     return pairwiseScores;
 330   }
 331
 332   public void setPairwiseScores(MatrixI m)
 333   {
 334     pairwiseScores = m;
 335   }
 336
 337   public MatrixI getEigenmatrix()
 338   {
 339     return eigenMatrix;
 340   }
 341
 342   public void setEigenmatrix(MatrixI m)
 343   {
 344     eigenMatrix = m;
 345   }
 346
 347   public MatrixI getTridiagonal()
 348   {
 349     return tridiagonal;
 350   }
 351
 352   public void setTridiagonal(MatrixI tridiagonal)
 353   {
 354     this.tridiagonal = tridiagonal;
 355   }
 356 }