X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FConservation.java;h=6cc9dd344addced46f9e8f5a06f7184ce47a86a7;hb=9d2408483e451285fd555c3cd6e0273977acbaa7;hp=f4346c380a70b03549eb41275a9cade5ae9570e3;hpb=7733b7433ee8771d85ccfe2ed32119feac8056e0;p=jalview.git diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index f4346c3..6cc9dd3 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -20,6 +20,8 @@ */ package jalview.analysis; +import java.util.Locale; + import jalview.analysis.scoremodels.ScoreMatrix; import jalview.analysis.scoremodels.ScoreModels; import jalview.datamodel.AlignmentAnnotation; @@ -30,6 +32,7 @@ import jalview.datamodel.Sequence; import jalview.datamodel.SequenceI; import jalview.schemes.ResidueProperties; import jalview.util.Comparison; +import jalview.util.Format; import java.awt.Color; import java.util.List; @@ -52,14 +55,21 @@ public class Conservation private static final int TOUPPERCASE = 'a' - 'A'; + private static final int GAP_INDEX = -1; + + private static final Format FORMAT_3DP = new Format("%2.5f"); + SequenceI[] sequences; int start; int end; - Vector 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 seqNums; int maxLength = 0; // used by quality calcs @@ -72,17 +82,17 @@ public class Conservation */ Map[] 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 quality; - /** Stores maximum and minimum values of quality values */ - private double[] qualityRange = new double[2]; + private double qualityMinimum; + + private double qualityMaximum; private Sequence consSequence; @@ -93,8 +103,16 @@ public class Conservation 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; /** @@ -130,8 +148,8 @@ public class Conservation * @param end * end column position */ - public Conservation(String name, int threshold, - List sequences, int start, int end) + public Conservation(String name, int threshold, List sequences, + int start, int end) { this.name = name; this.threshold = threshold; @@ -172,7 +190,6 @@ public class Conservation */ private void calcSeqNum(int i, ScoreMatrix sm) { - int gapIndex = sm.getGapIndex(); int sSize = sequences.length; if ((i > -1) && (i < sSize)) @@ -207,10 +224,18 @@ public class Conservation for (j = 1; j <= len; j++) { // sqnum[j] = ResidueProperties.aaIndex[sq.charAt(j - 1)]; - sqnum[j] = sm.getMatrixIndex(sq.charAt(j - 1)); - if (sqnum[j] == -1) + char residue = sq.charAt(j - 1); + if (Comparison.isGap(residue)) + { + sqnum[j] = GAP_INDEX; + } + else { - sqnum[j] = gapIndex; + sqnum[j] = sm.getMatrixIndex(residue); + if (sqnum[j] == -1) + { + sqnum[j] = GAP_INDEX; + } } } @@ -224,8 +249,8 @@ public class Conservation else { // JBPNote INFO level debug - System.err - .println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n"); + System.err.println( + "ERROR: calcSeqNum called with out of range sequence index for Alignment\n"); } } @@ -253,7 +278,7 @@ public class Conservation * or not conserved (-1) * Using TreeMap means properties are displayed in alphabetical order */ - SortedMap resultHash = new TreeMap(); + SortedMap resultHash = new TreeMap<>(); SymbolCounts symbolCounts = values.getSymbolCounts(); char[] symbols = symbolCounts.symbols; int[] counts = symbolCounts.values; @@ -289,7 +314,7 @@ public class Conservation protected static void recordConservation(Map resultMap, String res) { - res = res.toUpperCase(); + res = res.toUpperCase(Locale.ROOT); for (Entry> property : ResidueProperties.propHash .entrySet()) { @@ -394,7 +419,8 @@ public class Conservation continue; } - char c = sequences[i].getCharAt(column); // gaps do not have upper/lower case + char c = sequences[i].getCharAt(column); // gaps do not have upper/lower + // case if (Comparison.isGap((c))) { @@ -536,7 +562,7 @@ public class Conservation // From Alignment.java in jalview118 public void findQuality() { - findQuality(0, maxLength - 1); + findQuality(0, maxLength - 1, ScoreModels.getInstance().getBlosum62()); } /** @@ -546,7 +572,7 @@ public class Conservation */ private void percentIdentity(ScoreMatrix sm) { - seqNums = new Vector(); + seqNums = new Vector<>(); int i = 0, iSize = sequences.length; // Do we need to calculate this again? for (i = 0; i < iSize; i++) @@ -554,13 +580,12 @@ public class Conservation calcSeqNum(i, sm); } - int gapIndex = sm.getGapIndex(); - 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]; int j = 0; @@ -570,12 +595,21 @@ public class Conservation 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][gapIndex]++; // gap count + cons2GapCounts[i]++; } j++; } @@ -583,76 +617,79 @@ public class Conservation } /** - * 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(); + quality = new Vector<>(); double max = -Double.MAX_VALUE; - ScoreMatrix blosum62 = ScoreModels.getInstance().getBlosum62(); - float[][] blosumScores = blosum62.getMatrix(); + float[][] scores = scoreMatrix.getMatrix(); - // Loop over columns // JBPNote Profiling info - // long ts = System.currentTimeMillis(); - // long te = System.currentTimeMillis(); - percentIdentity(blosum62); + 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; } - final int symbolCount = blosum62.getSize(); - int gapIndex = blosum62.getGapIndex(); + final int symbolCount = scoreMatrix.getSize(); - for (j = startRes; j <= endRes; j++) + for (int j = startCol; j <= endCol; j++) { - bigtot = 0; + double bigtot = 0; // First Xr = depends on column only - x = new double[symbolCount]; + double[] x = new double[symbolCount]; - for (ii = 0; ii < symbolCount; ii++) + for (int ii = 0; ii < symbolCount; ii++) { x[ii] = 0; - for (i2 = 0; i2 < symbolCount; 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] * blosumScores[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[symbolCount]; - seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : gapIndex; - // Sequence, or gap at the end + 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; - // This is a loop over r - for (i = 0; i < symbolCount - 1; i++) + for (int i = 0; i < symbolCount - 1; i++) { - sr = 0; - - sr = (double) blosumScores[i][seqNum] + 4; + 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]); @@ -661,28 +698,21 @@ public class Conservation 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 + quality.addElement(Double.valueOf(bigtot)); } double newmax = -Double.MAX_VALUE; - for (j = startRes; j <= endRes; j++) + for (int j = startCol; j <= endCol; j++) { - tmp = quality.elementAt(j).doubleValue(); + double tmp = quality.elementAt(j).doubleValue(); // tmp = ((max - tmp) * (size - cons2[j][23])) / size; - tmp = ((max - tmp) * (size - cons2[j][gapIndex])) / size; + tmp = ((max - tmp) * (size - cons2GapCounts[j])) / size; // System.out.println(tmp+ " " + j); - quality.setElementAt(new Double(tmp), j); + quality.setElementAt(Double.valueOf(tmp), j); if (tmp > newmax) { @@ -690,15 +720,14 @@ public class Conservation } } - // System.out.println("Quality " + s); - qualityRange[0] = 0D; - qualityRange[1] = newmax; + qualityMinimum = 0D; + qualityMaximum = newmax; } /** * Complete the given consensus and quuality annotation rows. Note: currently - * this method will enlarge the given annotation row if it is too small, - * otherwise will leave its length unchanged. + * this method will reallocate the given annotation row if it is different to + * the calculated width, otherwise will leave its length unchanged. * * @param conservation * conservation annotation row @@ -712,51 +741,46 @@ public class Conservation public void completeAnnotations(AlignmentAnnotation conservation, AlignmentAnnotation quality2, int istart, int alWidth) { - char[] sequence = getConsSequence().getSequence(); - float minR; - float minG; - float minB; - float maxR; - float maxG; - float maxB; - minR = 0.3f; - minG = 0.0f; - minB = 0f; - maxR = 1.0f - minR; - maxG = 0.9f - minG; - maxB = 0f - minB; // scalable range for colouring both Conservation and - // Quality + SequenceI cons = getConsSequence(); + + /* + * colour scale for Conservation and Quality; + */ + float minR = 0.3f; + float minG = 0.0f; + float minB = 0f; + float maxR = 1.0f - minR; + float maxG = 0.9f - minG; + float maxB = 0f - minB; float min = 0f; float max = 11f; float qmin = 0f; float qmax = 0f; - char c; - if (conservation != null && conservation.annotations != null - && conservation.annotations.length < alWidth) + && conservation.annotations.length != alWidth) { conservation.annotations = new Annotation[alWidth]; } if (quality2 != null) { - quality2.graphMax = (float) qualityRange[1]; + quality2.graphMax = (float) qualityMaximum; if (quality2.annotations != null - && quality2.annotations.length < alWidth) + && 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++) { float value = 0; - c = sequence[i]; + char c = cons.getCharAt(i); if (Character.isDigit(c)) { @@ -776,11 +800,11 @@ public class Conservation float vprop = value - min; vprop /= max; int consp = i - start; - String conssym = (value > 0 && consp > -1 && consp < consSymbs.length) ? consSymbs[consp] - : ""; + String conssym = (value > 0 && consp > -1 + && consp < consSymbs.length) ? consSymbs[consp] : ""; conservation.annotations[i] = new Annotation(String.valueOf(c), - conssym, ' ', value, new Color(minR + (maxR * vprop), minG - + (maxG * vprop), minB + (maxB * vprop))); + conssym, ' ', value, new Color(minR + (maxR * vprop), + minG + (maxG * vprop), minB + (maxB * vprop))); } // Quality calc @@ -789,10 +813,10 @@ public class Conservation value = quality.elementAt(i).floatValue(); float vprop = value - qmin; vprop /= qmax; - quality2.annotations[i] = new Annotation(" ", - String.valueOf(value), ' ', value, new Color(minR - + (maxR * vprop), minG + (maxG * vprop), minB - + (maxB * vprop))); + String description = FORMAT_3DP.form(value); + quality2.annotations[i] = new Annotation(" ", description, + ' ', value, new Color(minR + (maxR * vprop), + minG + (maxG * vprop), minB + (maxB * vprop))); } } } @@ -843,11 +867,12 @@ public class Conservation */ String getTooltip(int column) { - char[] sequence = getConsSequence().getSequence(); - char val = column < sequence.length ? sequence[column] : '-'; + SequenceI cons = getConsSequence(); + char val = column < cons.getLength() ? cons.getCharAt(column) : '-'; boolean hasConservation = val != '-' && val != '0'; int consp = column - start; - String tip = (hasConservation && consp > -1 && consp < consSymbs.length) ? consSymbs[consp] + String tip = (hasConservation && consp > -1 && consp < consSymbs.length) + ? consSymbs[consp] : ""; return tip; }