2 * Jalview - A Sequence Alignment Editor and Viewer
\r
3 * Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
\r
5 * This program is free software; you can redistribute it and/or
\r
6 * modify it under the terms of the GNU General Public License
\r
7 * as published by the Free Software Foundation; either version 2
\r
8 * of the License, or (at your option) any later version.
\r
10 * This program is distributed in the hope that it will be useful,
\r
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
\r
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
\r
13 * GNU General Public License for more details.
\r
15 * You should have received a copy of the GNU General Public License
\r
16 * along with this program; if not, write to the Free Software
\r
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
\r
19 package jalview.analysis;
\r
21 import jalview.datamodel.*;
\r
27 * Calculates conservation values for a given set of sequences
\r
30 * @version $Revision$
\r
32 public class Conservation
\r
37 Vector seqNums; // vector of int vectors where first is sequence checksum
\r
38 int maxLength = 0; // used by quality calcs
\r
39 boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
\r
40 Vector total = new Vector();
\r
42 /** Stores calculated quality values */
\r
43 public Vector quality;
\r
45 /** Stores maximum and minimum values of quality values */
\r
46 public Double[] qualityRange = new Double[2];
\r
47 String consString = "";
\r
48 Sequence consSequence;
\r
55 * Creates a new Conservation object.
\r
57 * @param name Name of conservation
\r
58 * @param propHash DOCUMENT ME!
\r
59 * @param threshold to count the residues in residueHash(). commonly used value is 3
\r
60 * @param sequences sequences to be used in calculation
\r
61 * @param start start residue position
\r
62 * @param end end residue position
\r
64 public Conservation(String name, Hashtable propHash, int threshold,
\r
65 Vector sequences, int start, int end)
\r
68 this.propHash = propHash;
\r
69 this.threshold = threshold;
\r
70 this.sequences = sequences;
\r
73 seqNums = new Vector(sequences.size());
\r
80 private void calcSeqNums()
\r
82 int i=0, iSize=sequences.size();
\r
83 for (i=0; i < iSize; i++)
\r
92 * @param i DOCUMENT ME!
\r
94 private void calcSeqNum(int i)
\r
96 String sq = null; // for dumb jbuilder not-inited exception warning
\r
99 if ((i > -1) && (i < sequences.size()))
\r
101 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
\r
103 if (seqNums.size() <= i)
\r
105 seqNums.addElement(new int[sq.length() + 1]);
\r
108 if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
\r
112 seqNumsChanged = true;
\r
113 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
\r
116 if (maxLength < len)
\r
121 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
\r
122 sqnum[0] = sq.hashCode();
\r
124 for (j = 1; j <= len; j++)
\r
126 sqnum[j] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(String.valueOf(
\r
127 sq.charAt(j - 1)))).intValue(); // yuk - JBPNote - case taken care of in aaHash
\r
130 seqNums.setElementAt(sqnum, i);
\r
135 // JBPNote INFO level debug
\r
136 System.err.println(
\r
137 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
\r
142 * Calculates the conservation values for given set of sequences
\r
144 public void calculate()
\r
146 Hashtable resultHash, residueHash, ht;
\r
147 int count, thresh, j, jSize = sequences.size();
\r
148 String sequence, res, type;
\r
149 Enumeration enumeration, enumeration2;
\r
151 for (int i = start; i <= end; i++)
\r
153 resultHash = new Hashtable();
\r
154 residueHash = new Hashtable();
\r
156 for (j = 0; j < jSize; j++)
\r
158 // JBPNote - have to make sure elements of the sequences vector
\r
159 // are tested like this everywhere...
\r
160 if (sequences.elementAt(j) instanceof Sequence)
\r
162 sequence = ((Sequence) sequences.elementAt(j)).getSequence();
\r
164 if (sequence.length() > i)
\r
166 res = String.valueOf(Character.toUpperCase(sequence.charAt(i)));
\r
168 if (residueHash.containsKey(res))
\r
170 count = ((Integer) residueHash.get(res)).intValue();
\r
172 residueHash.put(res, new Integer(count));
\r
176 residueHash.put(res, new Integer(1));
\r
181 if (residueHash.containsKey("-"))
\r
183 count = ((Integer) residueHash.get("-")).intValue();
\r
185 residueHash.put("-", new Integer(count));
\r
189 residueHash.put("-", new Integer(1));
\r
195 //What is the count threshold to count the residues in residueHash()
\r
196 thresh = (threshold * (sequences.size())) / 100;
\r
198 //loop over all the found residues
\r
199 enumeration = residueHash.keys();
\r
201 while (enumeration.hasMoreElements())
\r
203 res = (String) enumeration.nextElement();
\r
205 if (((Integer) residueHash.get(res)).intValue() > thresh)
\r
207 //Now loop over the properties
\r
208 enumeration2 = propHash.keys();
\r
210 while (enumeration2.hasMoreElements())
\r
212 type = (String) enumeration2.nextElement();
\r
213 ht = (Hashtable) propHash.get(type);
\r
215 //Have we ticked this before?
\r
216 if (!resultHash.containsKey(type))
\r
218 if (ht.containsKey(res))
\r
220 resultHash.put(type, ht.get(res));
\r
224 resultHash.put(type, ht.get("-"));
\r
227 else if (((Integer) resultHash.get(type)).equals(
\r
228 (Integer) ht.get(res)) == false)
\r
230 resultHash.put(type, new Integer(-1));
\r
236 total.addElement(resultHash);
\r
243 * returns gap count in int[0], and conserved residue count in int[1]
\r
245 public int[] countConsNGaps(int j)
\r
250 int[] r = new int[2];
\r
252 int i, iSize = sequences.size();
\r
255 for (i = 0; i < iSize; i++)
\r
257 if (j >= ((Sequence) sequences.elementAt(i)).getLength())
\r
263 c = ((Sequence) sequences.elementAt(i)).getCharAt(j); // gaps do not have upper/lower case
\r
265 if (jalview.util.Comparison.isGap((c)))
\r
285 r[0] = (nres == cons) ? 1 : 0;
\r
292 * Calculates the conservation sequence
\r
294 * @param consflag if true, poitiveve conservation; false calculates negative conservation
\r
295 * @param percentageGaps commonly used value is 25
\r
297 public void verdict(boolean consflag, float percentageGaps)
\r
299 StringBuffer consString = new StringBuffer();
\r
303 int totGaps, count;
\r
305 Hashtable resultHash ;
\r
306 Enumeration enumeration;
\r
309 for (int i = start; i <= end; i++)
\r
311 gapcons = countConsNGaps(i);
\r
312 totGaps = gapcons[1];
\r
313 pgaps = ((float) totGaps * 100) / (float) sequences.size();
\r
315 if (percentageGaps > pgaps)
\r
317 resultHash = (Hashtable) total.elementAt(i - start);
\r
319 //Now find the verdict
\r
321 enumeration = resultHash.keys();
\r
323 while (enumeration.hasMoreElements())
\r
325 type = (String) enumeration.nextElement();
\r
326 result = (Integer) resultHash.get(type);
\r
328 //Do we want to count +ve conservation or +ve and -ve cons.?
\r
331 if (result.intValue() == 1)
\r
338 if (result.intValue() != -1)
\r
347 consString.append(count); // Conserved props!=Identity
\r
351 consString.append((gapcons[0] == 1) ? "*" : "+");
\r
356 consString.append("-");
\r
360 consSequence = new Sequence(name, consString.toString(), start, end);
\r
366 * @return Conservation sequence
\r
368 public Sequence getConsSequence()
\r
370 return consSequence;
\r
373 // From Alignment.java in jalview118
\r
374 public void findQuality()
\r
376 findQuality(0, maxLength - 1);
\r
382 private void percentIdentity2()
\r
384 calcSeqNums(); // updates maxLength, too.
\r
386 if ((cons2 == null) || seqNumsChanged)
\r
388 cons2 = new int[maxLength][24];
\r
390 // Initialize the array
\r
391 for (int j = 0; j < 24; j++)
\r
393 for (int i = 0; i < maxLength; i++)
\r
402 while (j < sequences.size())
\r
404 sqnum = (int[]) seqNums.elementAt(j);
\r
406 for (int i = 1; i < sqnum.length; i++)
\r
408 cons2[i - 1][sqnum[i]]++;
\r
411 for (int i = sqnum.length - 1; i < maxLength; i++)
\r
413 cons2[i][23]++; // gap count
\r
421 /* for (int i=start; i <= end; i++) {
\r
426 for (int j=0;j<24;j++) {
\r
427 if (cons2[i][j] > max) {
\r
439 * Calculates the quality of the set of sequences
\r
441 * @param start Start residue
\r
442 * @param end End residue
\r
444 public void findQuality(int start, int end)
\r
446 quality = new Vector();
\r
448 double max = -10000;
\r
449 int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
\r
451 //Loop over columns // JBPNote Profiling info
\r
452 // long ts = System.currentTimeMillis();
\r
453 //long te = System.currentTimeMillis();
\r
454 percentIdentity2();
\r
456 int size = seqNums.size();
\r
457 int[] lengths = new int[size];
\r
458 double tot, bigtot, sr, tmp;
\r
460 int l, j, i, ii, seqNum;
\r
462 for (l = 0; l < size; l++)
\r
463 lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
\r
466 for (j = start; j <= end; j++)
\r
470 // First Xr = depends on column only
\r
471 x = new double[24];
\r
473 for (ii = 0; ii < 24; ii++)
\r
479 for (int i2 = 0; i2 < 24; i2++)
\r
481 x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) +
\r
485 catch (Exception e)
\r
487 System.err.println("Exception during quality calculation.");
\r
488 e.printStackTrace();
\r
491 //System.out.println("X " + ii + " " + x[ii]);
\r
494 //System.out.println("X " + ii + " " + x[ii]);
\r
497 // Now calculate D for each position and sum
\r
498 for (int k = 0; k < size; k++)
\r
501 xx = new double[24];
\r
502 seqNum = (j < lengths[k])
\r
503 ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
\r
505 // This is a loop over r
\r
506 for (i = 0; i < 23; i++)
\r
512 sr = (double) BLOSUM62[i][seqNum] + 4;
\r
514 catch (Exception e)
\r
516 System.out.println("Exception in sr: " + e);
\r
517 e.printStackTrace();
\r
520 //Calculate X with another loop over residues
\r
521 // System.out.println("Xi " + i + " " + x[i] + " " + sr);
\r
524 tot += (xx[i] * xx[i]);
\r
527 bigtot += Math.sqrt(tot);
\r
530 // This is the quality for one column
\r
536 // bigtot = bigtot * (size-cons2[j][23])/size;
\r
537 quality.addElement(new Double(bigtot));
\r
540 // Need to normalize by gaps
\r
543 double newmax = -10000;
\r
545 for (j = start; j <= end; j++)
\r
547 tmp = ((Double) quality.elementAt(j)).doubleValue();
\r
548 tmp = ((max - tmp) * (size - cons2[j][23])) / size;
\r
550 // System.out.println(tmp+ " " + j);
\r
551 quality.setElementAt(new Double(tmp), j);
\r
559 // System.out.println("Quality " + s);
\r
560 qualityRange[0] = new Double(0);
\r
561 qualityRange[1] = new Double(newmax);
\r