From c650abb5e4572557e1d175ecc72a0f7c8892202e Mon Sep 17 00:00:00 2001 From: gmungoc Date: Tue, 21 Feb 2017 10:40:02 +0000 Subject: [PATCH] JAL-2416 quality calculation now only depends on loaded score matrix --- src/jalview/analysis/Conservation.java | 102 +++++++++++++++----------------- 1 file changed, 49 insertions(+), 53 deletions(-) diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index 79fa3d6..d1a3c37 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -20,6 +20,7 @@ */ package jalview.analysis; +import jalview.analysis.scoremodels.ScoreMatrix; import jalview.analysis.scoremodels.ScoreModels; import jalview.datamodel.AlignmentAnnotation; import jalview.datamodel.Annotation; @@ -163,27 +164,30 @@ public class Conservation } /** - * 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 gapIndex = sm.getMatrixIndex(' '); 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; @@ -196,14 +200,18 @@ public class Conservation 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)]; + sqnum[j] = sm.getMatrixIndex(sq.charAt(j - 1)); + if (sqnum[j] == -1) + { + sqnum[j] = gapIndex; + } } seqNums.setElementAt(sqnum, i); @@ -533,37 +541,32 @@ public class Conservation /** * DOCUMENT ME! + * + * @param sm */ - private void percentIdentity2() + private void percentIdentity(ScoreMatrix sm) { seqNums = new Vector(); - // 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); } + int gapIndex = sm.getMatrixIndex(' '); + 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]; - // 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++) { @@ -572,21 +575,10 @@ public class Conservation for (i = sqnum.length - 1; i < maxLength; i++) { - cons2[i][23]++; // gap count + cons2[i][gapIndex]++; // gap count } - 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; } } } - */ } } @@ -602,14 +594,15 @@ public class Conservation { quality = new Vector(); - double max = -10000; - float[][] BLOSUM62 = ScoreModels.getInstance().getBlosum62() - .getMatrix(); + double max = -Double.MAX_VALUE; + ScoreMatrix blosum62 = ScoreModels.getInstance().getBlosum62(); + float[][] blosumScores = blosum62.getMatrix(); + int gapIndex = blosum62.getMatrixIndex(' '); // Loop over columns // JBPNote Profiling info // long ts = System.currentTimeMillis(); // long te = System.currentTimeMillis(); - percentIdentity2(); + percentIdentity(blosum62); int size = seqNums.size(); int[] lengths = new int[size]; @@ -622,20 +615,24 @@ public class Conservation lengths[l] = seqNums.elementAt(l).length - 1; } + // todo ? remove '*' (unused?) from score matrix and + // use getSize() here instead of getSize() - 1 ?? + final int symbolCount = blosum62.getSize() - 1; // 24; + for (j = startRes; j <= endRes; j++) { bigtot = 0; // First Xr = depends on column only - x = new double[24]; + x = new double[symbolCount]; - for (ii = 0; ii < 24; ii++) + for (ii = 0; ii < symbolCount; ii++) { x[ii] = 0; - for (i2 = 0; i2 < 24; i2++) + for (i2 = 0; i2 < symbolCount; i2++) { - x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4); + x[ii] += (((double) cons2[j][i2] * blosumScores[ii][i2]) + 4); } x[ii] /= size; @@ -645,18 +642,16 @@ public class Conservation for (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 + xx = new double[symbolCount]; + seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : gapIndex; + // Sequence, or gap at the end // This is a loop over r - for (i = 0; i < 23; i++) + for (i = 0; i < symbolCount - 1; i++) { sr = 0; - sr = (double) BLOSUM62[i][seqNum] + 4; + sr = (double) blosumScores[i][seqNum] + 4; // Calculate X with another loop over residues // System.out.println("Xi " + i + " " + x[i] + " " + sr); @@ -680,12 +675,13 @@ public class Conservation // Need to normalize by gaps } - double newmax = -10000; + double newmax = -Double.MAX_VALUE; for (j = startRes; j <= endRes; j++) { tmp = quality.elementAt(j).doubleValue(); - tmp = ((max - tmp) * (size - cons2[j][23])) / size; + // tmp = ((max - tmp) * (size - cons2[j][23])) / size; + tmp = ((max - tmp) * (size - cons2[j][gapIndex])) / size; // System.out.println(tmp+ " " + j); quality.setElementAt(new Double(tmp), j); -- 1.7.10.2