From 7652d1d5889b737d52dbd15f3919d8b0dbdbc48b Mon Sep 17 00:00:00 2001 From: jprocter Date: Thu, 10 Apr 2008 13:01:31 +0000 Subject: [PATCH] quick fix for conservation calculation when residue ambiguity codes are used. --- src/jalview/analysis/Conservation.java | 215 +++++++++++++++++++------------- 1 file changed, 129 insertions(+), 86 deletions(-) diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index 4eb1218..c216471 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -24,44 +24,68 @@ import jalview.datamodel.*; /** * Calculates conservation values for a given set of sequences - * + * * @author $author$ * @version $Revision$ */ public class Conservation { SequenceI[] sequences; + int start; + int end; + Vector seqNums; // vector of int vectors where first is sequence checksum - int maxLength = 0; // used by quality calcs + + int maxLength = 0; // used by quality calcs + boolean seqNumsChanged = false; // updated after any change via calcSeqNum; + Hashtable[] total; + boolean canonicaliseAa = true; // if true then conservation calculation will + // map all symbols to canonical aa numbering + // rather than consider conservation of that + // symbol + /** Stores calculated quality values */ public Vector quality; - /** Stores maximum and minimum values of quality values */ + /** Stores maximum and minimum values of quality values */ public Double[] qualityRange = new Double[2]; + String consString = ""; + Sequence consSequence; + Hashtable propHash; + int threshold; + String name = ""; + int[][] cons2; /** * Creates a new Conservation object. - * - * @param name Name of conservation - * @param propHash DOCUMENT ME! - * @param threshold to count the residues in residueHash(). commonly used value is 3 - * @param sequences sequences to be used in calculation - * @param start start residue position - * @param end end residue position + * + * @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 + * sequences to be used in calculation + * @param start + * start residue position + * @param end + * end residue position */ public Conservation(String name, Hashtable propHash, int threshold, - Vector sequences, int start, int end) + Vector sequences, int start, int end) { this.name = name; @@ -70,7 +94,8 @@ public class Conservation this.start = start; this.end = end; - maxLength = end - start + 1; // default width includes bounds of calculation + maxLength = end - start + 1; // default width includes bounds of + // calculation int s, sSize = sequences.size(); SequenceI[] sarray = new SequenceI[sSize]; @@ -87,9 +112,10 @@ public class Conservation } /** - * DOCUMENT ME! - * - * @param i DOCUMENT ME! + * Translate sequence i into a numerical + * representation and store it in the i'th position of the seqNums array. + * + * @param i */ private void calcSeqNum(int i) { @@ -98,7 +124,7 @@ public class Conservation int sSize = sequences.length; - if ( (i > -1) && (i < sSize)) + if ((i > -1) && (i < sSize)) { sq = sequences[i].getSequenceAsString(); @@ -107,7 +133,7 @@ public class Conservation seqNums.addElement(new int[sq.length() + 1]); } - if (sq.hashCode() != ( (int[]) seqNums.elementAt(i))[0]) + if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0]) { int j; int len; @@ -119,12 +145,14 @@ public class Conservation maxLength = len; } - sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length + 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] = jalview.schemes.ResidueProperties.aaIndex[sq + .charAt(j - 1)]; } seqNums.setElementAt(sqnum, i); @@ -137,8 +165,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"); } } @@ -166,17 +194,34 @@ public class Conservation { c = sequences[j].getCharAt(i); - // No need to check if its a '-' - if (c == '.' || c == ' ') - { - c = '-'; + 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); + } } - - if ('a' <= c && c <= 'z') + else { - c -= (32); // 32 = 'a' - 'A' - } + // 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]++; } else @@ -185,10 +230,10 @@ public class Conservation } } - //What is the count threshold to count the residues in residueHash() + // What is the count threshold to count the residues in residueHash() thresh = (threshold * (jSize)) / 100; - //loop over all the found residues + // loop over all the found residues resultHash = new Hashtable(); for (char v = '-'; v < 'Z'; v++) { @@ -197,7 +242,7 @@ public class Conservation { res = String.valueOf(v); - //Now loop over the properties + // Now loop over the properties enumeration2 = propHash.keys(); while (enumeration2.hasMoreElements()) @@ -205,7 +250,7 @@ public class Conservation type = (String) enumeration2.nextElement(); ht = (Hashtable) propHash.get(type); - //Have we ticked this before? + // Have we ticked this before? if (!resultHash.containsKey(type)) { if (ht.containsKey(res)) @@ -217,10 +262,10 @@ public class Conservation resultHash.put(type, ht.get("-")); } } - else if ( ( (Integer) resultHash.get(type)).equals( - (Integer) ht.get(res)) == false) + else if (((Integer) resultHash.get(type)).equals((Integer) ht + .get(res)) == false) { - resultHash.put(type, new Integer( -1)); + resultHash.put(type, new Integer(-1)); } } } @@ -230,9 +275,9 @@ public class Conservation } } - /*** - * countConsNGaps - * returns gap count in int[0], and conserved residue count in int[1] + /***************************************************************************** + * count conservation for the j'th column of the alignment + * @return { gap count, conserved residue count} */ public int[] countConsNGaps(int j) { @@ -254,7 +299,7 @@ public class Conservation c = sequences[i].getCharAt(j); // gaps do not have upper/lower case - if (jalview.util.Comparison.isGap( (c))) + if (jalview.util.Comparison.isGap((c))) { count++; } @@ -282,9 +327,12 @@ public class Conservation /** * Calculates the conservation sequence - * - * @param consflag if true, poitiveve conservation; false calculates negative conservation - * @param percentageGaps commonly used value is 25 + * + * @param consflag + * if true, poitiveve conservation; false calculates negative + * conservation + * @param percentageGaps + * commonly used value is 25 */ public void verdict(boolean consflag, float percentageGaps) { @@ -297,9 +345,9 @@ public class Conservation Hashtable resultHash; Enumeration enumeration; - //NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY - //EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE - //DOES NOT EXIST IN JALVIEW 2.1.2 + // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY + // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE + // DOES NOT EXIST IN JALVIEW 2.1.2 for (int i = 0; i < start; i++) { consString.append('-'); @@ -309,13 +357,13 @@ public class Conservation { gapcons = countConsNGaps(i); totGaps = gapcons[1]; - pgaps = ( (float) totGaps * 100) / (float) sequences.length; + pgaps = ((float) totGaps * 100) / (float) sequences.length; if (percentageGaps > pgaps) { resultHash = total[i - start]; - //Now find the verdict + // Now find the verdict count = 0; enumeration = resultHash.keys(); @@ -324,7 +372,7 @@ public class Conservation type = (String) enumeration.nextElement(); result = (Integer) resultHash.get(type); - //Do we want to count +ve conservation or +ve and -ve cons.? + // Do we want to count +ve conservation or +ve and -ve cons.? if (consflag) { if (result.intValue() == 1) @@ -347,7 +395,7 @@ public class Conservation } else { - consString.append( (gapcons[0] == 1) ? "*" : "+"); + consString.append((gapcons[0] == 1) ? "*" : "+"); } } else @@ -360,8 +408,8 @@ public class Conservation } /** - * - * + * + * * @return Conservation sequence */ public Sequence getConsSequence() @@ -383,13 +431,13 @@ public class Conservation seqNums = new Vector(); // calcSeqNum(s); int i = 0, iSize = sequences.length; - //Do we need to calculate this again? + // Do we need to calculate this again? for (i = 0; i < iSize; i++) { calcSeqNum(i); } - if ( (cons2 == null) || seqNumsChanged) + if ((cons2 == null) || seqNumsChanged) { cons2 = new int[maxLength][24]; @@ -424,28 +472,24 @@ public class Conservation // 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; - } - - } - } */ + /* + * 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 - * - * @param start Start residue - * @param end End residue + * + * @param start + * Start residue + * @param end + * End residue */ public void findQuality(int start, int end) { @@ -454,9 +498,9 @@ public class Conservation double max = -10000; int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62(); - //Loop over columns // JBPNote Profiling info - //long ts = System.currentTimeMillis(); - //long te = System.currentTimeMillis(); + // Loop over columns // JBPNote Profiling info + // long ts = System.currentTimeMillis(); + // long te = System.currentTimeMillis(); percentIdentity2(); int size = seqNums.size(); @@ -467,7 +511,7 @@ public class Conservation for (l = 0; l < size; l++) { - lengths[l] = ( (int[]) seqNums.elementAt(l)).length - 1; + lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1; } for (j = start; j <= end; j++) @@ -483,8 +527,7 @@ public class Conservation for (i2 = 0; i2 < 24; i2++) { - x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) + - 4); + x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4); } x[ii] /= size; @@ -495,8 +538,8 @@ 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]) ? ((int[]) seqNums.elementAt(k))[j + 1] + : 23; // Sequence, or gap at the end // This is a loop over r for (i = 0; i < 23; i++) @@ -505,8 +548,8 @@ public class Conservation sr = (double) BLOSUM62[i][seqNum] + 4; - //Calculate X with another loop over residues - // System.out.println("Xi " + i + " " + x[i] + " " + sr); + // 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]); @@ -521,7 +564,7 @@ public class Conservation max = bigtot; } - // bigtot = bigtot * (size-cons2[j][23])/size; + // bigtot = bigtot * (size-cons2[j][23])/size; quality.addElement(new Double(bigtot)); // Need to normalize by gaps @@ -531,10 +574,10 @@ public class Conservation for (j = start; j <= end; j++) { - tmp = ( (Double) quality.elementAt(j)).doubleValue(); - tmp = ( (max - tmp) * (size - cons2[j][23])) / size; + tmp = ((Double) quality.elementAt(j)).doubleValue(); + tmp = ((max - tmp) * (size - cons2[j][23])) / size; - // System.out.println(tmp+ " " + j); + // System.out.println(tmp+ " " + j); quality.setElementAt(new Double(tmp), j); if (tmp > newmax) @@ -543,7 +586,7 @@ public class Conservation } } - // System.out.println("Quality " + s); + // System.out.println("Quality " + s); qualityRange[0] = new Double(0); qualityRange[1] = new Double(newmax); } -- 1.7.10.2