caf7bafc197178741eb39cb1dc6a2b07c04f6c89
[jalview.git] / 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 java.io.*;
24
25 import jalview.datamodel.*;
26 import jalview.datamodel.BinarySequence.InvalidSequenceTypeException;
27 import jalview.math.*;
28 import jalview.schemes.ResidueProperties;
29 import jalview.schemes.ScoreMatrix;
30
31 /**
32  * Performs Principal Component Analysis on given sequences
33  * 
34  * @author $author$
35  * @version $Revision$
36  */
37 public class PCA implements Runnable
38 {
39   Matrix m;
40
41   Matrix symm;
42
43   Matrix m2;
44
45   double[] eigenvalue;
46
47   Matrix eigenvector;
48
49   StringBuffer details = new StringBuffer();
50
51   /**
52    * Creates a new PCA object. By default, uses blosum62 matrix to generate
53    * sequence similarity matrices
54    * 
55    * @param s
56    *          Set of amino acid sequences to perform PCA on
57    */
58   public PCA(String[] s)
59   {
60     this(s, false);
61   }
62
63   /**
64    * Creates a new PCA object. By default, uses blosum62 matrix to generate
65    * sequence similarity matrices
66    * 
67    * @param s
68    *          Set of sequences to perform PCA on
69    * @param nucleotides
70    *          if true, uses standard DNA/RNA matrix for sequence similarity
71    *          calculation.
72    */
73   public PCA(String[] s, boolean nucleotides)
74   {
75     this(s, nucleotides, null);
76   }
77
78   public PCA(String[] s, boolean nucleotides, String s_m)
79   {
80
81     BinarySequence[] bs = new BinarySequence[s.length];
82     int ii = 0;
83
84     while ((ii < s.length) && (s[ii] != null))
85     {
86       bs[ii] = new BinarySequence(s[ii], nucleotides);
87       bs[ii].encode();
88       ii++;
89     }
90
91     BinarySequence[] bs2 = new BinarySequence[s.length];
92     ii = 0;
93     ScoreMatrix smtrx = null;
94     String sm = s_m;
95     if (sm != null)
96     {
97       smtrx = ResidueProperties.getScoreMatrix(sm);
98     }
99     if (smtrx == null)
100     {
101       // either we were given a non-existent score matrix or a scoremodel that
102       // isn't based on a pairwise symbol score matrix
103       smtrx = ResidueProperties.getScoreMatrix(sm = (nucleotides ? "DNA"
104               : "BLOSUM62"));
105     }
106     details.append("PCA calculation using " + sm
107             + " sequence similarity matrix\n========\n\n");
108     while ((ii < s.length) && (s[ii] != null))
109     {
110       bs2[ii] = new BinarySequence(s[ii], nucleotides);
111       if (smtrx != null)
112       {
113         try
114         {
115           bs2[ii].matrixEncode(smtrx);
116         } catch (InvalidSequenceTypeException x)
117         {
118           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
119         }
120       }
121       ii++;
122     }
123
124     // System.out.println("Created binary encoding");
125     // printMemory(rt);
126     int count = 0;
127
128     while ((count < bs.length) && (bs[count] != null))
129     {
130       count++;
131     }
132
133     double[][] seqmat = new double[count][bs[0].getDBinary().length];
134     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
135     int i = 0;
136
137     while (i < count)
138     {
139       seqmat[i] = bs[i].getDBinary();
140       seqmat2[i] = bs2[i].getDBinary();
141       i++;
142     }
143
144     // System.out.println("Created array");
145     // printMemory(rt);
146     // System.out.println(" --- Original matrix ---- ");
147     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
148     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
149
150   }
151
152   /**
153    * Returns the matrix used in PCA calculation
154    * 
155    * @return java.math.Matrix object
156    */
157
158   public Matrix getM()
159   {
160     return m;
161   }
162
163   /**
164    * Returns Eigenvalue
165    * 
166    * @param i
167    *          Index of diagonal within matrix
168    * 
169    * @return Returns value of diagonal from matrix
170    */
171   public double getEigenvalue(int i)
172   {
173     return eigenvector.d[i];
174   }
175
176   /**
177    * DOCUMENT ME!
178    * 
179    * @param l
180    *          DOCUMENT ME!
181    * @param n
182    *          DOCUMENT ME!
183    * @param mm
184    *          DOCUMENT ME!
185    * @param factor
186    *          DOCUMENT ME!
187    * 
188    * @return DOCUMENT ME!
189    */
190   public float[][] getComponents(int l, int n, int mm, float factor)
191   {
192     float[][] out = new float[m.rows][3];
193
194     for (int i = 0; i < m.rows; i++)
195     {
196       out[i][0] = (float) component(i, l) * factor;
197       out[i][1] = (float) component(i, n) * factor;
198       out[i][2] = (float) component(i, mm) * factor;
199     }
200
201     return out;
202   }
203
204   /**
205    * DOCUMENT ME!
206    * 
207    * @param n
208    *          DOCUMENT ME!
209    * 
210    * @return DOCUMENT ME!
211    */
212   public double[] component(int n)
213   {
214     // n = index of eigenvector
215     double[] out = new double[m.rows];
216
217     for (int i = 0; i < m.rows; i++)
218     {
219       out[i] = component(i, n);
220     }
221
222     return out;
223   }
224
225   /**
226    * DOCUMENT ME!
227    * 
228    * @param row
229    *          DOCUMENT ME!
230    * @param n
231    *          DOCUMENT ME!
232    * 
233    * @return DOCUMENT ME!
234    */
235   double component(int row, int n)
236   {
237     double out = 0.0;
238
239     for (int i = 0; i < symm.cols; i++)
240     {
241       out += (symm.value[row][i] * eigenvector.value[i][n]);
242     }
243
244     return out / eigenvector.d[n];
245   }
246
247   public String getDetails()
248   {
249     return details.toString();
250   }
251
252   /**
253    * DOCUMENT ME!
254    */
255   public void run()
256   {
257     PrintStream ps = new PrintStream(System.out)
258     {
259       public void print(String x)
260       {
261         details.append(x);
262       }
263
264       public void println()
265       {
266         details.append("\n");
267       }
268     };
269
270     try
271     {
272       details.append("PCA Calculation Mode is "
273               + (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
274               + "\n");
275       Matrix mt = m.transpose();
276
277       details.append(" --- OrigT * Orig ---- \n");
278       if (!jvCalcMode)
279       {
280         eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
281       }
282       else
283       {
284         eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace
285                                           // method
286       }
287
288       eigenvector.print(ps);
289
290       symm = eigenvector.copy();
291
292       eigenvector.tred();
293
294       details.append(" ---Tridiag transform matrix ---\n");
295       details.append(" --- D vector ---\n");
296       eigenvector.printD(ps);
297       ps.println();
298       details.append("--- E vector ---\n");
299       eigenvector.printE(ps);
300       ps.println();
301
302       // Now produce the diagonalization matrix
303       eigenvector.tqli();
304     } catch (Exception q)
305     {
306       q.printStackTrace();
307       details.append("\n*** Unexpected exception when performing PCA ***\n"
308               + q.getLocalizedMessage());
309       details.append("*** Matrices below may not be fully diagonalised. ***\n");
310     }
311
312     details.append(" --- New diagonalization matrix ---\n");
313     eigenvector.print(ps);
314     details.append(" --- Eigenvalues ---\n");
315     eigenvector.printD(ps);
316     ps.println();
317     /*
318      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
319      * (int ev=0;ev<symm.rows; ev++) {
320      * 
321      * ps.print(","+component(seq, ev)); } ps.println(); }
322      */
323   }
324
325   boolean jvCalcMode = true;
326
327   public void setJvCalcMode(boolean calcMode)
328   {
329     this.jvCalcMode = calcMode;
330   }
331 }