X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FConservation.java;h=f65420c311cd8ea5305447299ebb48eb41ac9e16;hb=174230b4233d9ce80f94527768d2cd2f76da11ab;hp=4a93eef2dedb82a35b69a6ca40d596a9464a279d;hpb=597b36936d89f172923e3321c8c5ea322cc5e5e5;p=jalview.git diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index 4a93eef..f65420c 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -1,192 +1,576 @@ -/* Jalview - a java multiple alignment editor - * Copyright (C) 1998 Michele Clamp - * - * This program 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 2 - * of the License, or (at your option) any later version. - * - * This program 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 this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - */ +/* +* Jalview - A Sequence Alignment Editor and Viewer +* Copyright (C) 2006 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle +* +* This program 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 2 +* of the License, or (at your option) any later version. +* +* This program 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 this program; if not, write to the Free Software +* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA +*/ package jalview.analysis; +import jalview.datamodel.*; import java.util.*; -import jalview.gui.*; -import jalview.datamodel.*; -public class Conservation { - Vector sequences; - int start; - int end; +/** + * Calculates conservation values for a given set of sequences + * + * @author $author$ + * @version $Revision$ + */ +public class Conservation +{ + Vector sequences; + int start; + int end; + 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; + Vector total = new Vector(); + + /** Stores calculated quality values */ + public Vector quality; + + /** 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 + */ + public Conservation(String name, Hashtable propHash, int threshold, + Vector sequences, int start, int end) + { + this.name = name; + this.propHash = propHash; + this.threshold = threshold; + this.sequences = sequences; + this.start = start; + this.end = end; + seqNums = new Vector(sequences.size()); + calcSeqNums(); + } - Vector total = new Vector(); + /** + * DOCUMENT ME! + */ + private void calcSeqNums() + { + int i=0, iSize=sequences.size(); + for (i=0; i < iSize; i++) + { + calcSeqNum(i); + } + } - String consString = ""; + /** + * DOCUMENT ME! + * + * @param i DOCUMENT ME! + */ + private void calcSeqNum(int i) + { + String sq = null; // for dumb jbuilder not-inited exception warning + int[] sqnum = null; + + if ((i > -1) && (i < sequences.size())) + { + sq = ((SequenceI) sequences.elementAt(i)).getSequence(); - DrawableSequence consSequence; - Hashtable propHash; - int threshold; - Hashtable[] freqs; + if (seqNums.size() <= i) + { + seqNums.addElement(new int[sq.length() + 1]); + } - String name = ""; + if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0]) + { + int j; + int len; + seqNumsChanged = true; + sq = ((SequenceI) sequences.elementAt(i)).getSequence(); + len = sq.length(); + + if (maxLength < len) + { + maxLength = len; + } + + 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] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(String.valueOf( + sq.charAt(j - 1)))).intValue(); // yuk - JBPNote - case taken care of in aaHash + } + + seqNums.setElementAt(sqnum, i); + } + } + else + { + // JBPNote INFO level debug + System.err.println( + "ERROR: calcSeqNum called with out of range sequence index for Alignment\n"); + } + } - public Conservation(String name,Hashtable[] freqs,Hashtable propHash, int threshold, Vector sequences, int start, int end) { - this.name = name; - this.freqs = freqs; - this.propHash = propHash; - this.threshold = threshold; - this.sequences = sequences; - this.start = start; - this.end = end; - } + /** + * Calculates the conservation values for given set of sequences + */ + public void calculate() + { + Hashtable resultHash, residueHash, ht; + int count, thresh, j, jSize = sequences.size(); + String type, res=null; + SequenceI sequence; + char c; + Enumeration enumeration, enumeration2; + + for (int i = start; i <= end; i++) + { + resultHash = new Hashtable(); + residueHash = new Hashtable(); + + for (j = 0; j < jSize; j++) + { + // JBPNote - have to make sure elements of the sequences vector + // are tested like this everywhere... + sequence = (Sequence) sequences.elementAt(j); + + if (sequence.getLength() > i) + { + c = sequence.getCharAt(i); + + // No need to check if its a '-' + if(c == '.' || c==' ') + c = '-'; + + if ('a' <= c && c <= 'z') + { + // TO UPPERCASE !!! + //Faster than toUpperCase + c -= ('a' - 'A') ; + } + + res = String.valueOf( c ); + + + if (residueHash.containsKey(res)) + { + count = ((Integer) residueHash.get(res)).intValue(); + count++; + residueHash.put(res, new Integer(count)); + } + else + { + residueHash.put(res, new Integer(1)); + } + } + else + { + if (residueHash.containsKey("-")) + { + count = ((Integer) residueHash.get("-")).intValue(); + count++; + residueHash.put("-", new Integer(count)); + } + else + { + residueHash.put("-", new Integer(1)); + } + } + } + + //What is the count threshold to count the residues in residueHash() + thresh = (threshold * (sequences.size())) / 100; + + //loop over all the found residues + enumeration = residueHash.keys(); + + while (enumeration.hasMoreElements()) + { + res = (String) enumeration.nextElement(); + + if (((Integer) residueHash.get(res)).intValue() > thresh) + { + //Now loop over the properties + enumeration2 = propHash.keys(); + + 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)); + } + } + } + } + total.addElement(resultHash); + } + } - public void calculate() { - for (int i = start;i <= end; i++) { - Hashtable resultHash = null; - Hashtable residueHash = null; + /*** + * countConsNGaps + * returns gap count in int[0], and conserved residue count in int[1] + */ + public int[] countConsNGaps(int j) + { + int count = 0; + int cons = 0; + int nres = 0; + int[] r = new int[2]; + char f = '$'; + int i, iSize = sequences.size(); + char c; + + for (i = 0; i < iSize; i++) + { + if (j >= ((Sequence) sequences.elementAt(i)).getLength()) + { + count++; + continue; + } - resultHash = new Hashtable(); - residueHash = new Hashtable(); + c = ((Sequence) sequences.elementAt(i)).getCharAt(j); // gaps do not have upper/lower case - for (int j=0; j < sequences.size(); j++) { + if (jalview.util.Comparison.isGap((c))) + { + count++; + } + else + { + nres++; + + if (nres == 1) + { + f = c; + cons++; + } + else if (f == c) + { + cons++; + } + } + } - if (sequences.elementAt(j) instanceof Sequence) { - Sequence s = (Sequence)sequences.elementAt(j); + r[0] = (nres == cons) ? 1 : 0; + r[1] = count; - if (s.getLength() > i) { - String res = s.getSequence().substring(i,i+1); + return r; + } - if (residueHash.containsKey(res)) { - int count = ((Integer)residueHash.get(res)).intValue() ; - count++; - residueHash.put(res,new Integer(count)); - } else { - residueHash.put(res,new Integer(1)); + /** + * Calculates the conservation sequence + * + * @param consflag if true, poitiveve conservation; false calculates negative conservation + * @param percentageGaps commonly used value is 25 + */ + public void verdict(boolean consflag, float percentageGaps) + { + StringBuffer consString = new StringBuffer(); + String type; + Integer result; + int[] gapcons; + int totGaps, count; + float pgaps; + Hashtable resultHash ; + Enumeration enumeration; + + + for (int i = start; i <= end; i++) + { + gapcons = countConsNGaps(i); + totGaps = gapcons[1]; + pgaps = ((float) totGaps * 100) / (float) sequences.size(); + + if (percentageGaps > pgaps) + { + resultHash = (Hashtable) total.elementAt(i - start); + + //Now find the verdict + count = 0; + enumeration = resultHash.keys(); + + while (enumeration.hasMoreElements()) + { + type = (String) enumeration.nextElement(); + result = (Integer) resultHash.get(type); + + //Do we want to count +ve conservation or +ve and -ve cons.? + if (consflag) + { + if (result.intValue() == 1) + { + count++; + } + } + else + { + if (result.intValue() != -1) + { + count++; + } + } + } + + if (count < 10) + { + consString.append(count); // Conserved props!=Identity + } + else + { + consString.append((gapcons[0] == 1) ? "*" : "+"); + } } - } else { - if (residueHash.containsKey("-")) { - int count = ((Integer)residueHash.get("-")).intValue() ; - count++; - residueHash.put("-",new Integer(count)); - } else { - residueHash.put("-",new Integer(1)); + else + { + consString.append("-"); } - } } - } - //What is the count threshold to count the residues in residueHash() - int thresh = threshold*(sequences.size())/100; + consSequence = new Sequence(name, consString.toString(), start, end); + } + + /** + * + * + * @return Conservation sequence + */ + public Sequence getConsSequence() + { + return consSequence; + } + + // From Alignment.java in jalview118 + public void findQuality() + { + findQuality(0, maxLength - 1); + } + + /** + * DOCUMENT ME! + */ + private void percentIdentity2() + { + calcSeqNums(); // updates maxLength, too. - //loop over all the found residues - Enumeration e = residueHash.keys(); + if ((cons2 == null) || seqNumsChanged) + { + cons2 = new int[maxLength][24]; + + // Initialize the array + for (int j = 0; j < 24; j++) + { + for (int i = 0; i < maxLength; i++) + { + cons2[i][j] = 0; + } + } - while (e.hasMoreElements()) { + int[] sqnum; + int j = 0; - String res = (String)e.nextElement(); - if (((Integer)residueHash.get(res)).intValue() > thresh) { + while (j < sequences.size()) + { + sqnum = (int[]) seqNums.elementAt(j); - //Now loop over the properties - Enumeration e2 = propHash.keys(); + for (int i = 1; i < sqnum.length; i++) + { + cons2[i - 1][sqnum[i]]++; + } - while (e2.hasMoreElements()) { - String type = (String)e2.nextElement(); - Hashtable ht = (Hashtable)propHash.get(type); + for (int i = sqnum.length - 1; i < maxLength; i++) + { + cons2[i][23]++; // gap count + } - //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)); + j++; } - } - } - } - total.addElement(resultHash); - } - } - public int countGaps(int j) { - int count = 0; + // unnecessary ? + + /* for (int i=start; i <= end; i++) { + int max = -1000; + int maxi = -1; + int maxj = -1; - for (int i = 0; i < sequences.size();i++) { - if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length()) - { count++; continue;} + for (int j=0;j<24;j++) { + if (cons2[i][j] > max) { + max = cons2[i][j]; + maxi = i; + maxj = j; + } - String tmp = ((Sequence)sequences.elementAt(i)).getSequence().substring(j,j+1); - if (tmp.equals(" ") || tmp.equals(".") || tmp.equals("-")) { - count++; - } + } + } */ + } } - return count; - } - public void verdict(boolean consflag, float percentageGaps) { - String consString = ""; + /** + * Calculates the quality of the set of sequences + * + * @param start Start residue + * @param end End residue + */ + public void findQuality(int start, int end) + { + quality = new Vector(); - for (int i=start; i <= end; i++) { - int totGaps = countGaps(i); - float pgaps = (float)totGaps*100/(float)sequences.size(); + double max = -10000; + int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62(); - if (percentageGaps > pgaps) - { - Hashtable resultHash = (Hashtable)total.elementAt(i); + //Loop over columns // JBPNote Profiling info + // long ts = System.currentTimeMillis(); + //long te = System.currentTimeMillis(); + percentIdentity2(); - //Now find the verdict - int count = 0; - Enumeration e3 = resultHash.keys(); + int size = seqNums.size(); + int[] lengths = new int[size]; + double tot, bigtot, sr, tmp; + double [] x, xx; + int l, j, i, ii, seqNum; - while (e3.hasMoreElements()) + for (l = 0; l < size; l++) + lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1; + + + for (j = start; j <= end; j++) { - String type = (String)e3.nextElement(); - Integer result = (Integer)resultHash.get(type); - - //Do we want to count +ve conservation or +ve and -ve cons.? - - if (consflag) - { - if (result.intValue() == 1) - count++; - } - else - { - if (result.intValue() != -1) - count++; - } + bigtot = 0; + + // First Xr = depends on column only + x = new double[24]; + + for (ii = 0; ii < 24; ii++) + { + x[ii] = 0; + + try + { + for (int i2 = 0; i2 < 24; i2++) + { + x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + + 4); + } + } + catch (Exception e) + { + System.err.println("Exception during quality calculation."); + e.printStackTrace(); + } + + //System.out.println("X " + ii + " " + x[ii]); + x[ii] /= (size); + + //System.out.println("X " + ii + " " + x[ii]); + } + + // Now calculate D for each position and sum + for (int k = 0; k < size; k++) + { + tot = 0; + xx = new double[24]; + 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++) + { + sr = 0; + + try + { + sr = (double) BLOSUM62[i][seqNum] + 4; + } + catch (Exception e) + { + System.out.println("Exception in sr: " + e); + e.printStackTrace(); + } + + //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; + } + + // bigtot = bigtot * (size-cons2[j][23])/size; + quality.addElement(new Double(bigtot)); + + + // Need to normalize by gaps } - if (count < 10) - consString = consString + String.valueOf(count); - else - consString = consString + "*"; + double newmax = -10000; - } - else - { - consString = consString + "-"; - } - } + for (j = start; j <= end; j++) + { + tmp = ((Double) quality.elementAt(j)).doubleValue(); + tmp = ((max - tmp) * (size - cons2[j][23])) / size; - consSequence = new DrawableSequence(name,consString,start,end); - } + // System.out.println(tmp+ " " + j); + quality.setElementAt(new Double(tmp), j); - public jalview.gui.DrawableSequence getConsSequence() { - return consSequence; - } + if (tmp > newmax) + { + newmax = tmp; + } + } + // System.out.println("Quality " + s); + qualityRange[0] = new Double(0); + qualityRange[1] = new Double(newmax); + } }