X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FConservation.java;h=ca5bbb3b399e94ea97ff55792355fe184c52ca3d;hb=256f99fa9de8ee5a2f206163a6a55be69b35e00c;hp=deeb2e9151acd4385d018af9b01f53862beb4d4d;hpb=e5127e5f27f02c8f328539eaac68eabfbee1d135;p=jalview.git diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java index deeb2e9..ca5bbb3 100755 --- a/src/jalview/analysis/Conservation.java +++ b/src/jalview/analysis/Conservation.java @@ -1,413 +1,543 @@ -/* 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. - */ -package jalview.analysis; - - -import java.util.*; -import jalview.gui.*; -import jalview.datamodel.*; - - -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; - private void calcSeqNums() { - for (int i=0; i-1 && i i) { - String res = s.getSequence().substring(i,i+1); - - 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)); - } - } else { - if (residueHash.containsKey("-")) { - int 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() - int thresh = threshold*(sequences.size())/100; - - //loop over all the found residues - Enumeration e = residueHash.keys(); - - while (e.hasMoreElements()) { - - String res = (String)e.nextElement(); - if (((Integer)residueHash.get(res)).intValue() > thresh) { - - //Now loop over the properties - Enumeration e2 = propHash.keys(); - - while (e2.hasMoreElements()) { - String type = (String)e2.nextElement(); - Hashtable 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 int countGaps(int j) - { - int count = 0; - - for (int i = 0; i < sequences.size();i++) - { - if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length()) - { count++; continue;} - - char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j); - if (jalview.util.Comparison.isGap((c))) - count++; - - } - return count; - } - /*** - * 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='$'; - for (int i = 0; i < sequences.size();i++) - { - if( j >= ((Sequence)sequences.elementAt(i)).getSequence().length()) - { count++; - continue;} - - char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j); - if (jalview.util.Comparison.isGap((c))) - count++; - else { - nres++; - if (nres==1) { - f = c; - cons++; - } else - if (f == c) { - cons++; - } - } - } - r[0] = (nres==cons) ? 1 : 0; - r[1] = count; - - - return r; - } - - public void verdict(boolean consflag, float percentageGaps) { - String consString = ""; - - for (int i=start; i <= end; i++) { - int[] gapcons = countConsNGaps(i); - boolean cons = (gapcons[0]==1) ? true : false; - int totGaps = gapcons[1]; - float pgaps = (float)totGaps*100/(float)sequences.size(); - // System.out.println("percentage gaps = "+pgaps+"\n"); - if (percentageGaps > pgaps) - { - Hashtable resultHash = (Hashtable)total.elementAt(i-start); - - //Now find the verdict - int count = 0; - Enumeration e3 = resultHash.keys(); - - while (e3.hasMoreElements()) - { - 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++; - } - } - - if (count < 10) - consString = consString + String.valueOf(count);// Conserved props!=Identity - else - consString = consString + ((gapcons[0]==1) ? "*" : "+"); - - } - else - { - consString = consString + "-"; - } - } - - consSequence = new Sequence(name,consString,start,end); - } - - public Sequence getConsSequence() { - return consSequence; - } - - // From Alignment.java in jalview118 - - public void findQuality() { - findQuality(0,maxLength-1); - } - - int[][] cons2; - - private void percentIdentity2() { - calcSeqNums(); // updates maxLength, too. - 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; - } - } - - int sqnum[]; - int j = 0; - while(j < sequences.size()) { - sqnum=(int[])seqNums.elementAt(j); - for (int i = 1; i < sqnum.length; i++) { - cons2[i-1][sqnum[i]]++; - } - for (int i=sqnum.length-1; i max) { - max = cons2[i][j]; - maxi = i; - maxj = j; - } - - } - } */ - } - -} - - -public void findQuality(int start, int end) { - quality = new Vector(); - double max = -10000; - String s = ""; - int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62(); - //Loop over columns // JBPNote Profiling info - // long ts = System.currentTimeMillis(); - //long te = System.currentTimeMillis(); - percentIdentity2(); - - int size = seqNums.size(); - int[] lengths = new int[size]; - - for (int l = 0; l < size; l++) - lengths[l] = ((int[]) seqNums.elementAt(l)).length-1; - - for (int j=start; j <= end; j++) { - double bigtot = 0; - - // First Xr = depends on column only - double x[] = new double[24]; - - for (int 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++) { - double tot = 0; - double[] xx = new double[24]; - int seqNum = - (jnewmax) - newmax = tmp; - } - // System.out.println("Quality " + s); - qualityRange[0] = new Double(0); - qualityRange[1] = new Double(newmax); - } - - -} +/* +* 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.*; + + +/** + * 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 + boolean seqNumsChanged = false; // updated after any change via calcSeqNum; + Hashtable [] total; + + /** 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.start = start; + this.end = end; + + maxLength=end-start+1; // default width includes bounds of calculation + + int s, sSize = sequences.size(); + SequenceI[] sarray = new SequenceI[sSize]; + this.sequences = sarray; + + for (s = 0; s < sSize; s++) + { + sarray[s] = (SequenceI) sequences.elementAt(s); + if(sarray[s].getLength()>maxLength) + maxLength = sarray[s].getLength(); + } + } + + + /** + * DOCUMENT ME! + * + * @param i DOCUMENT ME! + */ + private void calcSeqNum(int i) + { + String sq = null; // for dumb jbuilder not-inited exception warning + int[] sqnum = null; + + int sSize = sequences.length; + + if ((i > -1) && (i < sSize)) + { + sq = sequences[i].getSequence(); + + if (seqNums.size() <= i) + { + seqNums.addElement(new int[sq.length() + 1]); + } + + if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0]) + { + int j; + int len; + seqNumsChanged = true; + 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] = jalview.schemes.ResidueProperties.aaIndex[sq.charAt(j-1)]; + } + + + seqNums.setElementAt(sqnum, i); + } + else + System.out.println("SEQUENCE HAS BEEN DELETED!!!"); + } + else + { + // JBPNote INFO level debug + System.err.println( + "ERROR: calcSeqNum called with out of range sequence index for Alignment\n"); + } + } + + /** + * Calculates the conservation values for given set of sequences + */ + public void calculate() + { + Hashtable resultHash, ht; + int thresh, j, jSize = sequences.length; + int[] values; // Replaces residueHash + String type, res=null; + char c; + Enumeration enumeration2; + + total = new Hashtable[maxLength]; + + for (int i = start; i <= end; i++) + { + values = new int[132]; + + for (j = 0; j < jSize; j++) + { + if (sequences[j].getLength() > i) + { + c = sequences[j].getCharAt(i); + + // No need to check if its a '-' + if (c == '.' || c == ' ') + c = '-'; + + if ('a' <= c && c <= 'z') + { + c -= (32);// 32 = 'a' - 'A' + } + + values[c]++; + } + else + { + values['-']++; + } + } + + //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 (int v = '-'; v < 'Z'; v++) + { + + if (values[v] > thresh) + { + res = String.valueOf( (char) v); + + //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[i-start] = resultHash; + } + } + + + /*** + * 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.length; + char c; + + for (i = 0; i < iSize; i++) + { + if (j >= sequences[i].getLength()) + { + count++; + continue; + } + + c = sequences[i].getCharAt(j); // gaps do not have upper/lower case + + if (jalview.util.Comparison.isGap((c))) + { + count++; + } + else + { + nres++; + + if (nres == 1) + { + f = c; + cons++; + } + else if (f == c) + { + cons++; + } + } + } + + r[0] = (nres == cons) ? 1 : 0; + r[1] = count; + + return r; + } + + /** + * 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=0; i pgaps) + { + resultHash = total[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 + { + consString.append('-'); + } + } + + 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() + { + 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); + } + + + if ((cons2 == null) || seqNumsChanged) + { + 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 = (int[]) seqNums.elementAt(j); + + for (i = 1; i < sqnum.length; i++) + { + cons2[i - 1][sqnum[i]]++; + } + + for (i = sqnum.length - 1; i < maxLength; i++) + { + cons2[i][23]++; // 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; + } + + } + } */ + } + } + + /** + * 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(); + + double max = -10000; + int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62(); + + //Loop over columns // JBPNote Profiling info + //long ts = System.currentTimeMillis(); + //long te = System.currentTimeMillis(); + percentIdentity2(); + + 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++) + lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1; + + for (j = start; j <= end; j++) + { + bigtot = 0; + + // First Xr = depends on column only + x = new double[24]; + + for (ii = 0; ii < 24; ii++) + { + x[ii] = 0; + + for (i2 = 0; i2 < 24; i2++) + { + x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) + + 4); + } + + x[ii] /= size; + } + + // Now calculate D for each position and sum + for (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; + + sr = (double) BLOSUM62[i][seqNum] + 4; + + //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 + } + + double newmax = -10000; + + for (j = start; j <= end; j++) + { + tmp = ((Double) quality.elementAt(j)).doubleValue(); + tmp = ((max - tmp) * (size - cons2[j][23])) / size; + + // System.out.println(tmp+ " " + j); + quality.setElementAt(new Double(tmp), j); + + if (tmp > newmax) + { + newmax = tmp; + } + } + + // System.out.println("Quality " + s); + qualityRange[0] = new Double(0); + qualityRange[1] = new Double(newmax); + } +}