2 * Jalview - A Sequence Alignment Editor and Viewer
3 * Copyright (C) 2006 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
19 package jalview.analysis;
21 import jalview.datamodel.*;
27 * Calculates conservation values for a given set of sequences
32 public class Conservation
37 Vector seqNums; // vector of int vectors where first is sequence checksum
38 int maxLength = 0; // used by quality calcs
39 boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
40 Vector total = new Vector();
42 /** Stores calculated quality values */
43 public Vector quality;
45 /** Stores maximum and minimum values of quality values */
46 public Double[] qualityRange = new Double[2];
47 String consString = "";
48 Sequence consSequence;
55 * Creates a new Conservation object.
57 * @param name Name of conservation
58 * @param propHash DOCUMENT ME!
59 * @param threshold to count the residues in residueHash(). commonly used value is 3
60 * @param sequences sequences to be used in calculation
61 * @param start start residue position
62 * @param end end residue position
64 public Conservation(String name, Hashtable propHash, int threshold,
65 Vector sequences, int start, int end)
68 this.propHash = propHash;
69 this.threshold = threshold;
70 this.sequences = sequences;
73 seqNums = new Vector(sequences.size());
80 private void calcSeqNums()
82 int i=0, iSize=sequences.size();
83 for (i=0; i < iSize; i++)
92 * @param i DOCUMENT ME!
94 private void calcSeqNum(int i)
96 String sq = null; // for dumb jbuilder not-inited exception warning
99 if ((i > -1) && (i < sequences.size()))
101 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
103 if (seqNums.size() <= i)
105 seqNums.addElement(new int[sq.length() + 1]);
108 if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
112 seqNumsChanged = true;
113 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
121 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
122 sqnum[0] = sq.hashCode();
124 for (j = 1; j <= len; j++)
126 sqnum[j] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(String.valueOf(
127 sq.charAt(j - 1)))).intValue(); // yuk - JBPNote - case taken care of in aaHash
130 seqNums.setElementAt(sqnum, i);
135 // JBPNote INFO level debug
137 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
142 * Calculates the conservation values for given set of sequences
144 public void calculate()
146 Hashtable resultHash, residueHash, ht;
147 int count, thresh, j, jSize = sequences.size();
148 String type, res=null;
151 Enumeration enumeration, enumeration2;
153 for (int i = start; i <= end; i++)
155 resultHash = new Hashtable();
156 residueHash = new Hashtable();
158 for (j = 0; j < jSize; j++)
160 // JBPNote - have to make sure elements of the sequences vector
161 // are tested like this everywhere...
162 sequence = (Sequence) sequences.elementAt(j);
164 if (sequence.getLength() > i)
166 c = sequence.getCharAt(i);
168 // No need to check if its a '-'
169 if(c == '.' || c==' ')
172 if ('a' <= c && c <= 'z')
175 //Faster than toUpperCase
179 res = String.valueOf( c );
182 if (residueHash.containsKey(res))
184 count = ((Integer) residueHash.get(res)).intValue();
186 residueHash.put(res, new Integer(count));
190 residueHash.put(res, new Integer(1));
195 if (residueHash.containsKey("-"))
197 count = ((Integer) residueHash.get("-")).intValue();
199 residueHash.put("-", new Integer(count));
203 residueHash.put("-", new Integer(1));
208 //What is the count threshold to count the residues in residueHash()
209 thresh = (threshold * (sequences.size())) / 100;
211 //loop over all the found residues
212 enumeration = residueHash.keys();
214 while (enumeration.hasMoreElements())
216 res = (String) enumeration.nextElement();
218 if (((Integer) residueHash.get(res)).intValue() > thresh)
220 //Now loop over the properties
221 enumeration2 = propHash.keys();
223 while (enumeration2.hasMoreElements())
225 type = (String) enumeration2.nextElement();
226 ht = (Hashtable) propHash.get(type);
228 //Have we ticked this before?
229 if (!resultHash.containsKey(type))
231 if (ht.containsKey(res))
233 resultHash.put(type, ht.get(res));
237 resultHash.put(type, ht.get("-"));
240 else if (((Integer) resultHash.get(type)).equals(
241 (Integer) ht.get(res)) == false)
243 resultHash.put(type, new Integer(-1));
249 total.addElement(resultHash);
256 * returns gap count in int[0], and conserved residue count in int[1]
258 public int[] countConsNGaps(int j)
263 int[] r = new int[2];
265 int i, iSize = sequences.size();
268 for (i = 0; i < iSize; i++)
270 if (j >= ((Sequence) sequences.elementAt(i)).getLength())
276 c = ((Sequence) sequences.elementAt(i)).getCharAt(j); // gaps do not have upper/lower case
278 if (jalview.util.Comparison.isGap((c)))
298 r[0] = (nres == cons) ? 1 : 0;
305 * Calculates the conservation sequence
307 * @param consflag if true, poitiveve conservation; false calculates negative conservation
308 * @param percentageGaps commonly used value is 25
310 public void verdict(boolean consflag, float percentageGaps)
312 StringBuffer consString = new StringBuffer();
318 Hashtable resultHash ;
319 Enumeration enumeration;
322 for (int i = start; i <= end; i++)
324 gapcons = countConsNGaps(i);
325 totGaps = gapcons[1];
326 pgaps = ((float) totGaps * 100) / (float) sequences.size();
328 if (percentageGaps > pgaps)
330 resultHash = (Hashtable) total.elementAt(i - start);
332 //Now find the verdict
334 enumeration = resultHash.keys();
336 while (enumeration.hasMoreElements())
338 type = (String) enumeration.nextElement();
339 result = (Integer) resultHash.get(type);
341 //Do we want to count +ve conservation or +ve and -ve cons.?
344 if (result.intValue() == 1)
351 if (result.intValue() != -1)
360 consString.append(count); // Conserved props!=Identity
364 consString.append((gapcons[0] == 1) ? "*" : "+");
369 consString.append("-");
373 consSequence = new Sequence(name, consString.toString(), start, end);
379 * @return Conservation sequence
381 public Sequence getConsSequence()
386 // From Alignment.java in jalview118
387 public void findQuality()
389 findQuality(0, maxLength - 1);
395 private void percentIdentity2()
397 calcSeqNums(); // updates maxLength, too.
399 if ((cons2 == null) || seqNumsChanged)
401 cons2 = new int[maxLength][24];
403 // Initialize the array
404 for (int j = 0; j < 24; j++)
406 for (int i = 0; i < maxLength; i++)
415 while (j < sequences.size())
417 sqnum = (int[]) seqNums.elementAt(j);
419 for (int i = 1; i < sqnum.length; i++)
421 cons2[i - 1][sqnum[i]]++;
424 for (int i = sqnum.length - 1; i < maxLength; i++)
426 cons2[i][23]++; // gap count
434 /* for (int i=start; i <= end; i++) {
439 for (int j=0;j<24;j++) {
440 if (cons2[i][j] > max) {
452 * Calculates the quality of the set of sequences
454 * @param start Start residue
455 * @param end End residue
457 public void findQuality(int start, int end)
459 quality = new Vector();
462 int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
464 //Loop over columns // JBPNote Profiling info
465 // long ts = System.currentTimeMillis();
466 //long te = System.currentTimeMillis();
469 int size = seqNums.size();
470 int[] lengths = new int[size];
471 double tot, bigtot, sr, tmp;
473 int l, j, i, ii, seqNum;
475 for (l = 0; l < size; l++)
476 lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
479 for (j = start; j <= end; j++)
483 // First Xr = depends on column only
486 for (ii = 0; ii < 24; ii++)
492 for (int i2 = 0; i2 < 24; i2++)
494 x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) +
500 System.err.println("Exception during quality calculation.");
504 //System.out.println("X " + ii + " " + x[ii]);
507 //System.out.println("X " + ii + " " + x[ii]);
510 // Now calculate D for each position and sum
511 for (int k = 0; k < size; k++)
515 seqNum = (j < lengths[k])
516 ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
518 // This is a loop over r
519 for (i = 0; i < 23; i++)
525 sr = (double) BLOSUM62[i][seqNum] + 4;
529 System.out.println("Exception in sr: " + e);
533 //Calculate X with another loop over residues
534 // System.out.println("Xi " + i + " " + x[i] + " " + sr);
537 tot += (xx[i] * xx[i]);
540 bigtot += Math.sqrt(tot);
543 // This is the quality for one column
549 // bigtot = bigtot * (size-cons2[j][23])/size;
550 quality.addElement(new Double(bigtot));
553 // Need to normalize by gaps
556 double newmax = -10000;
558 for (j = start; j <= end; j++)
560 tmp = ((Double) quality.elementAt(j)).doubleValue();
561 tmp = ((max - tmp) * (size - cons2[j][23])) / size;
563 // System.out.println(tmp+ " " + j);
564 quality.setElementAt(new Double(tmp), j);
572 // System.out.println("Quality " + s);
573 qualityRange[0] = new Double(0);
574 qualityRange[1] = new Double(newmax);