*/
package jalview.analysis;
+import jalview.analysis.scoremodels.ScoreMatrix;
+import jalview.analysis.scoremodels.ScoreModels;
import jalview.datamodel.AlignmentAnnotation;
import jalview.datamodel.Annotation;
import jalview.datamodel.ResidueCount;
private static final int TOUPPERCASE = 'a' - 'A';
+ private static final int GAP_INDEX = -1;
+
SequenceI[] sequences;
int start;
int end;
- Vector<int[]> seqNums; // vector of int vectors where first is sequence
- // checksum
+ /*
+ * a list whose i'th element is an array whose first entry is the checksum
+ * of the i'th sequence, followed by residues encoded to score matrix index
+ */
+ Vector<int[]> seqNums;
int maxLength = 0; // used by quality calcs
*/
Map<String, Integer>[] total;
- boolean canonicaliseAa = true; // if true then conservation calculation will
-
- // map all symbols to canonical aa numbering
- // rather than consider conservation of that
- // symbol
+ /*
+ * if true then conservation calculation will map all symbols to canonical aa
+ * numbering rather than consider conservation of that symbol
+ */
+ boolean canonicaliseAa = true;
- /** Stores calculated quality values */
private Vector<Double> quality;
- /** Stores maximum and minimum values of quality values */
- private double[] qualityRange = new double[2];
+ private double qualityMinimum;
+
+ private double qualityMaximum;
private Sequence consSequence;
private String name = "";
+ /*
+ * an array, for each column, of counts of symbols (by score matrix index)
+ */
private int[][] cons2;
+ /*
+ * gap counts for each column
+ */
+ private int[] cons2GapCounts;
+
private String[] consSymbs;
/**
}
/**
- * Translate sequence i into a numerical representation and store it in the
- * i'th position of the seqNums array.
+ * Translate sequence i into score matrix indices and store it in the i'th
+ * position of the seqNums array.
*
* @param i
+ * @param sm
*/
- private void calcSeqNum(int i)
+ private void calcSeqNum(int i, ScoreMatrix sm)
{
- String sq = null; // for dumb jbuilder not-inited exception warning
- int[] sqnum = null;
-
int sSize = sequences.length;
if ((i > -1) && (i < sSize))
{
- sq = sequences[i].getSequenceAsString();
+ String sq = sequences[i].getSequenceAsString();
if (seqNums.size() <= i)
{
seqNums.addElement(new int[sq.length() + 1]);
}
+ /*
+ * the first entry in the array is the sequence's hashcode,
+ * following entries are matrix indices of sequence characters
+ */
if (sq.hashCode() != seqNums.elementAt(i)[0])
{
int j;
maxLength = len;
}
- sqnum = new int[len + 1]; // better to always make a new array -
+ int[] sqnum = new int[len + 1]; // better to always make a new array -
// sequence can change its length
sqnum[0] = sq.hashCode();
for (j = 1; j <= len; j++)
{
- sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
- .charAt(j - 1)];
+ // sqnum[j] = ResidueProperties.aaIndex[sq.charAt(j - 1)];
+ char residue = sq.charAt(j - 1);
+ if (Comparison.isGap(residue))
+ {
+ sqnum[j] = GAP_INDEX;
+ }
+ else
+ {
+ sqnum[j] = sm.getMatrixIndex(residue);
+ if (sqnum[j] == -1)
+ {
+ sqnum[j] = GAP_INDEX;
+ }
+ }
}
seqNums.setElementAt(sqnum, i);
// From Alignment.java in jalview118
public void findQuality()
{
- findQuality(0, maxLength - 1);
+ findQuality(0, maxLength - 1, ScoreModels.getInstance().getBlosum62());
}
/**
* DOCUMENT ME!
+ *
+ * @param sm
*/
- private void percentIdentity2()
+ private void percentIdentity(ScoreMatrix sm)
{
seqNums = new Vector<int[]>();
- // calcSeqNum(s);
int i = 0, iSize = sequences.length;
// Do we need to calculate this again?
for (i = 0; i < iSize; i++)
{
- calcSeqNum(i);
+ calcSeqNum(i, sm);
}
if ((cons2 == null) || seqNumsChanged)
{
+ // FIXME remove magic number 24 without changing calc
+ // sm.getSize() returns 25 so doesn't quite do it...
cons2 = new int[maxLength][24];
+ cons2GapCounts = new int[maxLength];
- // Initialize the array
- for (int j = 0; j < 24; j++)
- {
- for (i = 0; i < maxLength; i++)
- {
- cons2[i][j] = 0;
- }
- }
-
- int[] sqnum;
int j = 0;
while (j < sequences.length)
{
- sqnum = seqNums.elementAt(j);
+ int[] sqnum = seqNums.elementAt(j);
for (i = 1; i < sqnum.length; i++)
{
- cons2[i - 1][sqnum[i]]++;
+ int index = sqnum[i];
+ if (index == GAP_INDEX)
+ {
+ cons2GapCounts[i - 1]++;
+ }
+ else
+ {
+ cons2[i - 1][index]++;
+ }
}
+ // TODO should this start from sqnum.length?
for (i = sqnum.length - 1; i < maxLength; i++)
{
- cons2[i][23]++; // gap count
+ cons2GapCounts[i]++;
}
-
j++;
}
-
- // unnecessary ?
-
- /*
- * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
- * maxj = -1;
- *
- * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
- * maxi = i; maxj = j; } } }
- */
}
}
/**
- * Calculates the quality of the set of sequences
+ * Calculates the quality of the set of sequences over the given inclusive
+ * column range, using the specified substitution score matrix
*
- * @param startRes
- * Start residue
- * @param endRes
- * End residue
+ * @param startCol
+ * @param endCol
+ * @param scoreMatrix
*/
- public void findQuality(int startRes, int endRes)
+ protected void findQuality(int startCol, int endCol, ScoreMatrix scoreMatrix)
{
quality = new Vector<Double>();
- double max = -10000;
- int[][] BLOSUM62 = ResidueProperties.getBLOSUM62();
+ double max = -Double.MAX_VALUE;
+ float[][] scores = scoreMatrix.getMatrix();
- // Loop over columns // JBPNote Profiling info
- // long ts = System.currentTimeMillis();
- // long te = System.currentTimeMillis();
- percentIdentity2();
+ percentIdentity(scoreMatrix);
int size = seqNums.size();
int[] lengths = new int[size];
- double tot, bigtot, sr, tmp;
- double[] x, xx;
- int l, j, i, ii, i2, k, seqNum;
- for (l = 0; l < size; l++)
+ for (int l = 0; l < size; l++)
{
lengths[l] = seqNums.elementAt(l).length - 1;
}
- for (j = startRes; j <= endRes; j++)
+ final int symbolCount = scoreMatrix.getSize();
+
+ for (int j = startCol; j <= endCol; j++)
{
- bigtot = 0;
+ double bigtot = 0;
// First Xr = depends on column only
- x = new double[24];
+ double[] x = new double[symbolCount];
- for (ii = 0; ii < 24; ii++)
+ for (int ii = 0; ii < symbolCount; ii++)
{
x[ii] = 0;
- for (i2 = 0; i2 < 24; i2++)
+ /*
+ * todo JAL-728 currently assuming last symbol in matrix is * for gap
+ * (which we ignore as counted separately); true for BLOSUM62 but may
+ * not be once alternative matrices are supported
+ */
+ for (int i2 = 0; i2 < symbolCount - 1; i2++)
{
- x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
+ x[ii] += (((double) cons2[j][i2] * scores[ii][i2]) + 4D);
}
+ x[ii] += 4D + cons2GapCounts[j] * scoreMatrix.getMinimumScore();
x[ii] /= size;
}
// Now calculate D for each position and sum
- for (k = 0; k < size; k++)
+ for (int k = 0; k < size; k++)
{
- tot = 0;
- xx = new double[24];
- seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : 23; // Sequence,
- // or gap
- // at the
- // end
-
- // This is a loop over r
- for (i = 0; i < 23; i++)
- {
- sr = 0;
+ double tot = 0;
+ double[] xx = new double[symbolCount];
+ // sequence character index, or implied gap if sequence too short
+ int seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1]
+ : GAP_INDEX;
- sr = (double) BLOSUM62[i][seqNum] + 4;
+ for (int i = 0; i < symbolCount - 1; i++)
+ {
+ double sr = 4D;
+ if (seqNum == GAP_INDEX)
+ {
+ sr += scoreMatrix.getMinimumScore();
+ }
+ else
+ {
+ sr += scores[i][seqNum];
+ }
- // Calculate X with another loop over residues
- // System.out.println("Xi " + i + " " + x[i] + " " + sr);
xx[i] = x[i] - sr;
tot += (xx[i] * xx[i]);
bigtot += Math.sqrt(tot);
}
- // This is the quality for one column
- if (max < bigtot)
- {
- max = bigtot;
- }
+ max = Math.max(max, bigtot);
- // bigtot = bigtot * (size-cons2[j][23])/size;
quality.addElement(new Double(bigtot));
-
- // Need to normalize by gaps
}
- double newmax = -10000;
+ double newmax = -Double.MAX_VALUE;
- for (j = startRes; j <= endRes; j++)
+ for (int j = startCol; j <= endCol; j++)
{
- tmp = quality.elementAt(j).doubleValue();
- tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+ double tmp = quality.elementAt(j).doubleValue();
+ // tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+ tmp = ((max - tmp) * (size - cons2GapCounts[j])) / size;
// System.out.println(tmp+ " " + j);
quality.setElementAt(new Double(tmp), j);
}
}
- // System.out.println("Quality " + s);
- qualityRange[0] = 0D;
- qualityRange[1] = newmax;
+ qualityMinimum = 0D;
+ qualityMaximum = newmax;
}
/**
if (quality2 != null)
{
- quality2.graphMax = (float) qualityRange[1];
+ quality2.graphMax = (float) qualityMaximum;
if (quality2.annotations != null
&& quality2.annotations.length < alWidth)
{
quality2.annotations = new Annotation[alWidth];
}
- qmin = (float) qualityRange[0];
- qmax = (float) qualityRange[1];
+ qmin = (float) qualityMinimum;
+ qmax = (float) qualityMaximum;
}
for (int i = istart; i < alWidth; i++)
git://source.jalview.org / jalview.git / blobdiff