X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FConservation.java;h=7e72dcd638cab9e327cbb4c5a01d502c452cc3fb;hb=0a5ce6145bb76fc7eb8a5cc2670e20453fbedd29;hp=2a065f6cb37fd8adff3d8ad6cccec467af3463a8;hpb=a45774ee31d9f35d4eff46d54d7deab719afb092;p=jalview.git diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index 2a065f6..7e72dcd 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -1,26 +1,40 @@ /* - * Jalview - A Sequence Alignment Editor and Viewer (Version 2.7) - * Copyright (C) 2011 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle + * 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. - * + * 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 . + * 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 java.awt.Color; -import java.util.*; +import jalview.datamodel.AlignmentAnnotation; +import jalview.datamodel.Annotation; +import jalview.datamodel.ResidueCount; +import jalview.datamodel.Sequence; +import jalview.datamodel.SequenceI; +import jalview.datamodel.ResidueCount.SymbolCounts; +import jalview.schemes.ResidueProperties; +import jalview.util.Comparison; -import jalview.datamodel.*; +import java.awt.Color; +import java.util.List; +import java.util.Map; +import java.util.Map.Entry; +import java.util.TreeMap; +import java.util.Vector; /** * Calculates conservation values for a given set of sequences @@ -30,19 +44,27 @@ import jalview.datamodel.*; */ public class Conservation { + private static final int TOUPPERCASE = 'a' - 'A'; + SequenceI[] sequences; int start; int end; - Vector seqNums; // vector of int vectors where first is sequence checksum + Vector seqNums; // vector of int vectors where first is sequence + // checksum int maxLength = 0; // used by quality calcs boolean seqNumsChanged = false; // updated after any change via calcSeqNum; - Hashtable[] total; + /* + * a map per column with {property, conservation} where conservation value is + * 1 (property is conserved), 0 (absence of property is conserved) or -1 + * (property is not conserved i.e. column has residues with and without it) + */ + Map[] total; boolean canonicaliseAa = true; // if true then conservation calculation will @@ -51,30 +73,29 @@ public class Conservation // symbol /** Stores calculated quality values */ - public Vector quality; + private Vector quality; /** Stores maximum and minimum values of quality values */ - public Double[] qualityRange = new Double[2]; - - String consString = ""; + private double[] qualityRange = new double[2]; - Sequence consSequence; + private Sequence consSequence; - Hashtable propHash; + /* + * percentage of residues in a column to qualify for counting conservation + */ + private int threshold; - int threshold; + private String name = ""; - String name = ""; + private int[][] cons2; - int[][] cons2; + private String[] consSymbs; /** * Creates a new Conservation object. * * @param name * Name of conservation - * @param propHash - * hash of properties for each symbol * @param threshold * to count the residues in residueHash(). commonly used value is 3 * @param sequences @@ -84,12 +105,10 @@ public class Conservation * @param end * end residue position */ - public Conservation(String name, Hashtable propHash, int threshold, - Vector sequences, int start, int end) + public Conservation(String name, int threshold, + List sequences, int start, int end) { - this.name = name; - this.propHash = propHash; this.threshold = threshold; this.start = start; this.end = end; @@ -100,20 +119,22 @@ public class Conservation int s, sSize = sequences.size(); SequenceI[] sarray = new SequenceI[sSize]; this.sequences = sarray; - try { - for (s = 0; s < sSize; s++) + try { - sarray[s] = (SequenceI) sequences.elementAt(s); - if (sarray[s].getLength() > maxLength) + for (s = 0; s < sSize; s++) { - maxLength = sarray[s].getLength(); + sarray[s] = sequences.get(s); + if (sarray[s].getLength() > maxLength) + { + maxLength = sarray[s].getLength(); + } } - } } catch (ArrayIndexOutOfBoundsException ex) { - // bail - another thread has modified the sequence array, so the current calculation is probably invalid. - this.sequences=new SequenceI[0]; - maxLength=0; + // bail - another thread has modified the sequence array, so the current + // calculation is probably invalid. + this.sequences = new SequenceI[0]; + maxLength = 0; } } @@ -139,7 +160,7 @@ public class Conservation seqNums.addElement(new int[sq.length() + 1]); } - if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0]) + if (sq.hashCode() != seqNums.elementAt(i)[0]) { int j; int len; @@ -181,178 +202,219 @@ public class Conservation */ public void calculate() { - Hashtable resultHash, ht; - int thresh, j, jSize = sequences.length; - int[] values; // Replaces residueHash - String type, res = null; - char c; - Enumeration enumeration2; + int height = sequences.length; - total = new Hashtable[maxLength]; + total = new Map[maxLength]; - for (int i = start; i <= end; i++) + for (int column = start; column <= end; column++) { - values = new int[255]; + ResidueCount values = countResidues(column); + + // TODO is threshold a percentage or count value? + int thresh = (threshold * height) / 100; - for (j = 0; j < jSize; j++) + /* + * check observed residues in column and record whether each + * physico-chemical property is conserved (+1), absence conserved (0), + * or not conserved (-1) + * Using TreeMap means properties are displayed in alphabetical order + */ + Map resultHash = new TreeMap(); + SymbolCounts symbolCounts = values.getSymbolCounts(); + char[] symbols = symbolCounts.symbols; + int[] counts = symbolCounts.values; + for (int j = 0; j < symbols.length; j++) { - if (sequences[j].getLength() > i) + char c = symbols[j]; + if (counts[j] > thresh) { - c = sequences[j].getCharAt(i); - - if (canonicaliseAa) - { // lookup the base aa code symbol - c = (char) jalview.schemes.ResidueProperties.aaIndex[sequences[j] - .getCharAt(i)]; - if (c > 20) - { - c = '-'; - } - else - { - // recover canonical aa symbol - c = jalview.schemes.ResidueProperties.aa[c].charAt(0); - } - } - else - { - // original behaviour - operate on ascii symbols directly - // No need to check if its a '-' - if (c == '.' || c == ' ') - { - c = '-'; - } - - if (!canonicaliseAa && 'a' <= c && c <= 'z') - { - c -= (32); // 32 = 'a' - 'A' - } - } - values[c]++; + recordConservation(resultHash, String.valueOf(c)); + } + } + if (values.getGapCount() > thresh) + { + recordConservation(resultHash, "-"); + } + + if (total.length > 0) + { + total[column - start] = resultHash; + } + } + } + + /** + * Updates the conservation results for an observed residue + * + * @param resultMap + * a map of {property, conservation} where conservation value is +1 + * (all residues have the property), 0 (no residue has the property) + * or -1 (some do, some don't) + * @param res + */ + protected static void recordConservation(Map resultMap, + String res) + { + res = res.toUpperCase(); + for (Entry> property : ResidueProperties.propHash + .entrySet()) + { + String propertyName = property.getKey(); + Integer residuePropertyValue = property.getValue().get(res); + + if (!resultMap.containsKey(propertyName)) + { + /* + * first time we've seen this residue - note whether it has this property + */ + if (residuePropertyValue != null) + { + resultMap.put(propertyName, residuePropertyValue); } else { - values['-']++; + /* + * unrecognised residue - use default value for property + */ + resultMap.put(propertyName, property.getValue().get("-")); } } - - // What is the count threshold to count the residues in residueHash() - thresh = (threshold * (jSize)) / 100; - - // loop over all the found residues - resultHash = new Hashtable(); - for (char v = '-'; v < 'Z'; v++) + else { - - if (values[v] > thresh) + Integer currentResult = resultMap.get(propertyName); + if (currentResult.intValue() != -1 + && !currentResult.equals(residuePropertyValue)) { - res = String.valueOf(v); + /* + * property is unconserved - residues seen both with and without it + */ + resultMap.put(propertyName, Integer.valueOf(-1)); + } + } + } + } - // Now loop over the properties - enumeration2 = propHash.keys(); + /** + * Counts residues (upper-cased) and gaps in the given column + * + * @param column + * @return + */ + protected ResidueCount countResidues(int column) + { + ResidueCount values = new ResidueCount(false); - while (enumeration2.hasMoreElements()) - { - type = (String) enumeration2.nextElement(); - ht = (Hashtable) propHash.get(type); - - // Have we ticked this before? - if (!resultHash.containsKey(type)) - { - if (ht.containsKey(res)) - { - resultHash.put(type, ht.get(res)); - } - else - { - resultHash.put(type, ht.get("-")); - } - } - else if (((Integer) resultHash.get(type)).equals((Integer) ht - .get(res)) == false) - { - resultHash.put(type, new Integer(-1)); - } - } + for (int row = 0; row < sequences.length; row++) + { + if (sequences[row].getLength() > column) + { + char c = sequences[row].getCharAt(column); + if (canonicaliseAa) + { + int index = ResidueProperties.aaIndex[c]; + c = index > 20 ? '-' : ResidueProperties.aa[index].charAt(0); + } + else + { + c = toUpperCase(c); + } + if (Comparison.isGap(c)) + { + values.addGap(); + } + else + { + values.add(c); } } - - if (total.length>0) { - total[i - start] = resultHash; + else + { + values.addGap(); } } + return values; } - /***************************************************************************** - * count conservation for the j'th column of the alignment + /** + * Counts conservation and gaps for a column of the alignment * - * @return { gap count, conserved residue count} + * @return { 1 if fully conserved, else 0, gap count } */ - public int[] countConsNGaps(int j) + public int[] countConservationAndGaps(int column) { - int count = 0; - int cons = 0; - int nres = 0; - int[] r = new int[2]; - char f = '$'; - int i, iSize = sequences.length; - char c; + int gapCount = 0; + boolean fullyConserved = true; + int iSize = sequences.length; - for (i = 0; i < iSize; i++) + if (iSize == 0) + { + return new int[] { 0, 0 }; + } + + char lastRes = '0'; + for (int i = 0; i < iSize; i++) { - if (j >= sequences[i].getLength()) + if (column >= sequences[i].getLength()) { - count++; + gapCount++; continue; } - c = sequences[i].getCharAt(j); // gaps do not have upper/lower case + char c = sequences[i].getCharAt(column); // gaps do not have upper/lower case - if (jalview.util.Comparison.isGap((c))) + if (Comparison.isGap((c))) { - count++; + gapCount++; } else { - nres++; - - if (nres == 1) + c = toUpperCase(c); + if (lastRes == '0') { - f = c; - cons++; + lastRes = c; } - else if (f == c) + if (c != lastRes) { - cons++; + fullyConserved = false; } } } - r[0] = (nres == cons) ? 1 : 0; - r[1] = count; - + int[] r = new int[] { fullyConserved ? 1 : 0, gapCount }; return r; } /** + * Returns the upper-cased character if between 'a' and 'z', else the + * unchanged value + * + * @param c + * @return + */ + char toUpperCase(char c) + { + if ('a' <= c && c <= 'z') + { + c -= TOUPPERCASE; + } + return c; + } + + /** * Calculates the conservation sequence * - * @param consflag - * if true, poitiveve conservation; false calculates negative - * conservation - * @param percentageGaps - * commonly used value is 25 + * @param positiveOnly + * if true, calculate positive conservation; else calculate both + * positive and negative conservation + * @param maxPercentageGaps + * the percentage of gaps in a column, at or above which no + * conservation is asserted */ - public void verdict(boolean consflag, float percentageGaps) + public void verdict(boolean positiveOnly, float maxPercentageGaps) { - StringBuffer consString = new StringBuffer(); - String type; - Integer result; - int[] gapcons; - int totGaps, count; - float pgaps; - Hashtable resultHash; - Enumeration enumeration; + // TODO call this at the end of calculate(), should not be a public method + + StringBuilder consString = new StringBuilder(end); // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE @@ -361,42 +423,53 @@ public class Conservation { consString.append('-'); } - + consSymbs = new String[end - start + 1]; for (int i = start; i <= end; i++) { - gapcons = countConsNGaps(i); - totGaps = gapcons[1]; - pgaps = ((float) totGaps * 100) / (float) sequences.length; + int[] gapcons = countConservationAndGaps(i); + boolean fullyConserved = gapcons[0] == 1; + int totGaps = gapcons[1]; + float pgaps = (totGaps * 100f) / sequences.length; - if (percentageGaps > pgaps) + if (maxPercentageGaps > pgaps) { - resultHash = total[i - start]; - - // Now find the verdict - count = 0; - enumeration = resultHash.keys(); - - while (enumeration.hasMoreElements()) + Map resultHash = total[i - start]; + int count = 0; + StringBuilder positives = new StringBuilder(64); + StringBuilder negatives = new StringBuilder(32); + for (String type : resultHash.keySet()) { - type = (String) enumeration.nextElement(); - result = (Integer) resultHash.get(type); - - // Do we want to count +ve conservation or +ve and -ve cons.? - if (consflag) + int result = resultHash.get(type).intValue(); + if (result == -1) + { + /* + * not conserved (present or absent) + */ + continue; + } + count++; + if (result == 1) { - if (result.intValue() == 1) - { - count++; - } + /* + * positively conserved property (all residues have it) + */ + positives.append(positives.length() == 0 ? "" : " "); + positives.append(type); } - else + if (result == 0 && !positiveOnly) { - if (result.intValue() != -1) - { - count++; - } + /* + * absense of property is conserved (all residues lack it) + */ + negatives.append(negatives.length() == 0 ? "" : " "); + negatives.append("!").append(type); } } + if (negatives.length() > 0) + { + positives.append(" ").append(negatives); + } + consSymbs[i - start] = positives.toString(); if (count < 10) { @@ -404,7 +477,7 @@ public class Conservation } else { - consString.append((gapcons[0] == 1) ? "*" : "+"); + consString.append(fullyConserved ? "*" : "+"); } } else @@ -437,7 +510,7 @@ public class Conservation */ private void percentIdentity2() { - seqNums = new Vector(); + seqNums = new Vector(); // calcSeqNum(s); int i = 0, iSize = sequences.length; // Do we need to calculate this again? @@ -464,7 +537,7 @@ public class Conservation while (j < sequences.length) { - sqnum = (int[]) seqNums.elementAt(j); + sqnum = seqNums.elementAt(j); for (i = 1; i < sqnum.length; i++) { @@ -494,17 +567,17 @@ public class Conservation /** * Calculates the quality of the set of sequences * - * @param start + * @param startRes * Start residue - * @param end + * @param endRes * End residue */ - public void findQuality(int start, int end) + public void findQuality(int startRes, int endRes) { - quality = new Vector(); + quality = new Vector(); double max = -10000; - int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62(); + int[][] BLOSUM62 = ResidueProperties.getBLOSUM62(); // Loop over columns // JBPNote Profiling info // long ts = System.currentTimeMillis(); @@ -519,10 +592,10 @@ public class Conservation for (l = 0; l < size; l++) { - lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1; + lengths[l] = seqNums.elementAt(l).length - 1; } - for (j = start; j <= end; j++) + for (j = startRes; j <= endRes; j++) { bigtot = 0; @@ -546,8 +619,10 @@ public class Conservation { tot = 0; xx = new double[24]; - seqNum = (j < lengths[k]) ? ((int[]) seqNums.elementAt(k))[j + 1] - : 23; // Sequence, or gap at the end + 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++) @@ -580,9 +655,9 @@ public class Conservation double newmax = -10000; - for (j = start; j <= end; j++) + for (j = startRes; j <= endRes; j++) { - tmp = ((Double) quality.elementAt(j)).doubleValue(); + tmp = quality.elementAt(j).doubleValue(); tmp = ((max - tmp) * (size - cons2[j][23])) / size; // System.out.println(tmp+ " " + j); @@ -595,12 +670,12 @@ public class Conservation } // System.out.println("Quality " + s); - qualityRange[0] = new Double(0); - qualityRange[1] = new Double(newmax); + qualityRange[0] = 0D; + qualityRange[1] = newmax; } /** - * complete the given consensus and quuality annotation rows. Note: currently + * 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. * @@ -638,7 +713,7 @@ public class Conservation char c; - if (conservation.annotations != null + if (conservation != null && conservation.annotations != null && conservation.annotations.length < alWidth) { conservation.annotations = new Annotation[alWidth]; @@ -646,17 +721,17 @@ public class Conservation if (quality2 != null) { - quality2.graphMax = qualityRange[1].floatValue(); + quality2.graphMax = (float) qualityRange[1]; if (quality2.annotations != null && quality2.annotations.length < alWidth) { quality2.annotations = new Annotation[alWidth]; } - qmin = qualityRange[0].floatValue(); - qmax = qualityRange[1].floatValue(); + qmin = (float) qualityRange[0]; + qmax = (float) qualityRange[1]; } - for (int i = 0; i < alWidth; i++) + for (int i = istart; i < alWidth; i++) { float value = 0; @@ -664,7 +739,7 @@ public class Conservation if (Character.isDigit(c)) { - value = (int) (c - '0'); + value = c - '0'; } else if (c == '*') { @@ -675,18 +750,23 @@ public class Conservation value = 10; } - float vprop = value - min; - vprop /= max; - conservation.annotations[i] = new Annotation(String.valueOf(c), - String.valueOf(value), ' ', value, new Color(minR - + (maxR * vprop), minG + (maxG * vprop), minB - + (maxB * vprop))); + if (conservation != null) + { + float vprop = value - min; + vprop /= max; + int consp = i - start; + 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))); + } // Quality calc if (quality2 != null) { - value = ((Double) quality.elementAt(i)).floatValue(); - vprop = value - qmin; + value = quality.elementAt(i).floatValue(); + float vprop = value - qmin; vprop /= qmax; quality2.annotations[i] = new Annotation(" ", String.valueOf(value), ' ', value, new Color(minR @@ -695,4 +775,62 @@ public class Conservation } } } + + /** + * construct and call the calculation methods on a new Conservation object + * + * @param name + * - name of conservation + * @param threshold + * - minimum number of conserved residues needed to indicate + * conservation (typically 3) + * @param seqs + * @param start + * first column in calculation window + * @param end + * last column in calculation window + * @param positiveOnly + * calculate positive (true) or positive and negative (false) + * conservation + * @param maxPercentGaps + * percentage of gaps tolerated in column + * @param calcQuality + * flag indicating if alignment quality should be calculated + * @return Conservation object ready for use in visualization + */ + public static Conservation calculateConservation(String name, + int threshold, List seqs, int start, int end, + boolean positiveOnly, int maxPercentGaps, boolean calcQuality) + { + Conservation cons = new Conservation(name, threshold, seqs, start, end); + cons.calculate(); + cons.verdict(positiveOnly, maxPercentGaps); + + if (calcQuality) + { + cons.findQuality(); + } + + return cons; + } + + /** + * Returns the computed tooltip (annotation description) for a given column. + * The tip is empty if the conservation score is zero, otherwise holds the + * conserved properties (and, optionally, properties whose absence is + * conserved). + * + * @param column + * @return + */ + String getTooltip(int column) + { + char[] sequence = getConsSequence().getSequence(); + char val = column < sequence.length ? sequence[column] : '-'; + boolean hasConservation = val != '-' && val != '0'; + int consp = column - start; + String tip = (hasConservation && consp > -1 && consp < consSymbs.length) ? consSymbs[consp] + : ""; + return tip; + } }