/*
* Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
* Copyright (C) $$Year-Rel$$ The Jalview Authors
*
* This file is part of Jalview.
*
* Jalview is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation, either version 3
* of the License, or (at your option) any later version.
*
* Jalview is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR
* PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Jalview. If not, see .
* The Jalview Authors are detailed in the 'AUTHORS' file.
*/
package jalview.analysis;
import jalview.math.MatrixI;
import jalview.schemes.ResidueProperties;
import jalview.schemes.ScoreMatrix;
import java.io.PrintStream;
/**
* Performs Principal Component Analysis on given sequences
*/
public class PCA implements Runnable
{
boolean jvCalcMode = true;
MatrixI symm;
double[] eigenvalue;
MatrixI eigenvector;
StringBuilder details = new StringBuilder(1024);
private String[] seqs;
private ScoreMatrix scoreMatrix;
/**
* Creates a new PCA object. By default, uses blosum62 matrix to generate
* sequence similarity matrices
*
* @param s
* Set of amino acid sequences to perform PCA on
*/
public PCA(String[] s)
{
this(s, false);
}
/**
* Creates a new PCA object. By default, uses blosum62 matrix to generate
* sequence similarity matrices
*
* @param s
* Set of sequences to perform PCA on
* @param nucleotides
* if true, uses standard DNA/RNA matrix for sequence similarity
* calculation.
*/
public PCA(String[] s, boolean nucleotides)
{
this(s, nucleotides, null);
}
public PCA(String[] s, boolean nucleotides, String s_m)
{
this.seqs = s;
scoreMatrix = null;
String sm = s_m;
if (sm != null)
{
scoreMatrix = ResidueProperties.getScoreMatrix(sm);
}
if (scoreMatrix == null)
{
// either we were given a non-existent score matrix or a scoremodel that
// isn't based on a pairwise symbol score matrix
scoreMatrix = ResidueProperties
.getScoreMatrix(sm = (nucleotides ? "DNA" : "BLOSUM62"));
}
details.append("PCA calculation using " + sm
+ " sequence similarity matrix\n========\n\n");
}
/**
* Returns Eigenvalue
*
* @param i
* Index of diagonal within matrix
*
* @return Returns value of diagonal from matrix
*/
public double getEigenvalue(int i)
{
return eigenvector.getD()[i];
}
/**
* DOCUMENT ME!
*
* @param l
* DOCUMENT ME!
* @param n
* DOCUMENT ME!
* @param mm
* DOCUMENT ME!
* @param factor
* DOCUMENT ME!
*
* @return DOCUMENT ME!
*/
public float[][] getComponents(int l, int n, int mm, float factor)
{
float[][] out = new float[getHeight()][3];
for (int i = 0; i < getHeight(); i++)
{
out[i][0] = (float) component(i, l) * factor;
out[i][1] = (float) component(i, n) * factor;
out[i][2] = (float) component(i, mm) * factor;
}
return out;
}
/**
* DOCUMENT ME!
*
* @param n
* DOCUMENT ME!
*
* @return DOCUMENT ME!
*/
public double[] component(int n)
{
// n = index of eigenvector
double[] out = new double[getHeight()];
for (int i = 0; i < out.length; i++)
{
out[i] = component(i, n);
}
return out;
}
/**
* DOCUMENT ME!
*
* @param row
* DOCUMENT ME!
* @param n
* DOCUMENT ME!
*
* @return DOCUMENT ME!
*/
double component(int row, int n)
{
double out = 0.0;
for (int i = 0; i < symm.width(); i++)
{
out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
}
return out / eigenvector.getD()[n];
}
public String getDetails()
{
return details.toString();
}
/**
* DOCUMENT ME!
*/
@Override
public void run()
{
PrintStream ps = new PrintStream(System.out)
{
@Override
public void print(String x)
{
details.append(x);
}
@Override
public void println()
{
details.append("\n");
}
};
// long now = System.currentTimeMillis();
try
{
details.append("PCA Calculation Mode is "
+ (jvCalcMode ? "Jalview variant" : "Original SeqSpace")
+ "\n");
eigenvector = scoreMatrix.computePairwiseScores(seqs);
details.append(" --- OrigT * Orig ---- \n");
eigenvector.print(ps, "%8.2f");
symm = eigenvector.copy();
eigenvector.tred();
details.append(" ---Tridiag transform matrix ---\n");
details.append(" --- D vector ---\n");
eigenvector.printD(ps, "%15.4e");
ps.println();
details.append("--- E vector ---\n");
eigenvector.printE(ps, "%15.4e");
ps.println();
// Now produce the diagonalization matrix
eigenvector.tqli();
} catch (Exception q)
{
q.printStackTrace();
details.append("\n*** Unexpected exception when performing PCA ***\n"
+ q.getLocalizedMessage());
details.append("*** Matrices below may not be fully diagonalised. ***\n");
}
details.append(" --- New diagonalization matrix ---\n");
eigenvector.print(ps, "%8.2f");
details.append(" --- Eigenvalues ---\n");
eigenvector.printD(ps, "%15.4e");
ps.println();
/*
* for (int seq=0;seq