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
23 import jalview.gui.*;
\r
28 public class Conservation {
\r
32 Vector seqNums; // vector of int vectors where first is sequence checksum
\r
33 int maxLength = 0; // used by quality calcs
\r
34 boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
\r
35 Vector total = new Vector();
\r
36 public Vector quality;
\r
37 public Double[] qualityRange = new Double[2];
\r
38 String consString = "";
\r
39 Sequence consSequence;
\r
45 public Conservation(String name, Hashtable propHash, int threshold,
\r
46 Vector sequences, int start, int end) {
\r
48 this.propHash = propHash;
\r
49 this.threshold = threshold;
\r
50 this.sequences = sequences;
\r
53 seqNums = new Vector(sequences.size());
\r
57 private void calcSeqNums() {
\r
58 for (int i = 0; i < sequences.size(); i++) {
\r
63 private void calcSeqNum(int i) {
\r
64 String sq = null; // for dumb jbuilder not-inited exception warning
\r
67 if ((i > -1) && (i < sequences.size())) {
\r
68 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
\r
70 if (seqNums.size() <= i) {
\r
71 seqNums.addElement(new int[sq.length() + 1]);
\r
74 if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0]) {
\r
77 seqNumsChanged = true;
\r
78 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
\r
81 if (maxLength < len) {
\r
85 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
\r
86 sqnum[0] = sq.hashCode();
\r
88 for (j = 1; j <= len; j++) {
\r
89 sqnum[j] = ((Integer) jalview.schemes.ResidueProperties.aaHash.get(new String(
\r
90 sq.substring(j - 1, j)))).intValue(); // yuk
\r
93 seqNums.setElementAt(sqnum, i);
\r
96 // JBPNote INFO level debug
\r
98 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
\r
102 public void calculate() {
\r
103 for (int i = start; i <= end; i++) {
\r
104 Hashtable resultHash = null;
\r
105 Hashtable residueHash = null;
\r
107 resultHash = new Hashtable();
\r
108 residueHash = new Hashtable();
\r
110 for (int j = 0; j < sequences.size(); j++) {
\r
111 // JBPNote - have to make sure elements of the sequences vector
\r
112 // are tested like this everywhere...
\r
113 if (sequences.elementAt(j) instanceof Sequence) {
\r
114 Sequence s = (Sequence) sequences.elementAt(j);
\r
116 if (s.getLength() > i) {
\r
117 String res = s.getSequence().substring(i, i + 1);
\r
119 if (residueHash.containsKey(res)) {
\r
120 int count = ((Integer) residueHash.get(res)).intValue();
\r
122 residueHash.put(res, new Integer(count));
\r
124 residueHash.put(res, new Integer(1));
\r
127 if (residueHash.containsKey("-")) {
\r
128 int count = ((Integer) residueHash.get("-")).intValue();
\r
130 residueHash.put("-", new Integer(count));
\r
132 residueHash.put("-", new Integer(1));
\r
138 //What is the count threshold to count the residues in residueHash()
\r
139 int thresh = (threshold * (sequences.size())) / 100;
\r
141 //loop over all the found residues
\r
142 Enumeration e = residueHash.keys();
\r
144 while (e.hasMoreElements()) {
\r
145 String res = (String) e.nextElement();
\r
147 if (((Integer) residueHash.get(res)).intValue() > thresh) {
\r
148 //Now loop over the properties
\r
149 Enumeration e2 = propHash.keys();
\r
151 while (e2.hasMoreElements()) {
\r
152 String type = (String) e2.nextElement();
\r
153 Hashtable ht = (Hashtable) propHash.get(type);
\r
155 //Have we ticked this before?
\r
156 if (!resultHash.containsKey(type)) {
\r
157 if (ht.containsKey(res)) {
\r
158 resultHash.put(type, ht.get(res));
\r
160 resultHash.put(type, ht.get("-"));
\r
162 } else if (((Integer) resultHash.get(type)).equals(
\r
163 (Integer) ht.get(res)) == false) {
\r
164 resultHash.put(type, new Integer(-1));
\r
170 total.addElement(resultHash);
\r
174 public int countGaps(int j) {
\r
177 for (int i = 0; i < sequences.size(); i++) {
\r
178 if ((j + 1) > ((Sequence) sequences.elementAt(i)).getSequence()
\r
185 char c = ((Sequence) sequences.elementAt(i)).getSequence().charAt(j);
\r
187 if (jalview.util.Comparison.isGap((c))) {
\r
197 * returns gap count in int[0], and conserved residue count in int[1]
\r
199 public int[] countConsNGaps(int j) {
\r
203 int[] r = new int[2];
\r
206 for (int i = 0; i < sequences.size(); i++) {
\r
207 if (j >= ((Sequence) sequences.elementAt(i)).getSequence().length()) {
\r
213 char c = ((Sequence) sequences.elementAt(i)).getSequence().charAt(j);
\r
215 if (jalview.util.Comparison.isGap((c))) {
\r
223 } else if (f == c) {
\r
229 r[0] = (nres == cons) ? 1 : 0;
\r
235 public void verdict(boolean consflag, float percentageGaps) {
\r
236 String consString = "";
\r
238 for (int i = start; i <= end; i++) {
\r
239 int[] gapcons = countConsNGaps(i);
\r
240 boolean cons = (gapcons[0] == 1) ? true : false;
\r
241 int totGaps = gapcons[1];
\r
242 float pgaps = ((float) totGaps * 100) / (float) sequences.size();
\r
244 // System.out.println("percentage gaps = "+pgaps+"\n");
\r
245 if (percentageGaps > pgaps) {
\r
246 Hashtable resultHash = (Hashtable) total.elementAt(i - start);
\r
248 //Now find the verdict
\r
250 Enumeration e3 = resultHash.keys();
\r
252 while (e3.hasMoreElements()) {
\r
253 String type = (String) e3.nextElement();
\r
254 Integer result = (Integer) resultHash.get(type);
\r
256 //Do we want to count +ve conservation or +ve and -ve cons.?
\r
258 if (result.intValue() == 1) {
\r
262 if (result.intValue() != -1) {
\r
269 consString = consString + String.valueOf(count); // Conserved props!=Identity
\r
271 consString = consString + ((gapcons[0] == 1) ? "*" : "+");
\r
274 consString = consString + "-";
\r
278 consSequence = new Sequence(name, consString, start, end);
\r
281 public Sequence getConsSequence() {
\r
282 return consSequence;
\r
285 // From Alignment.java in jalview118
\r
286 public void findQuality() {
\r
287 findQuality(0, maxLength - 1);
\r
290 private void percentIdentity2() {
\r
291 calcSeqNums(); // updates maxLength, too.
\r
293 if ((cons2 == null) || seqNumsChanged) {
\r
294 cons2 = new int[maxLength][24];
\r
296 // Initialize the array
\r
297 for (int j = 0; j < 24; j++) {
\r
298 for (int i = 0; i < maxLength; i++) {
\r
306 while (j < sequences.size()) {
\r
307 sqnum = (int[]) seqNums.elementAt(j);
\r
309 for (int i = 1; i < sqnum.length; i++) {
\r
310 cons2[i - 1][sqnum[i]]++;
\r
313 for (int i = sqnum.length - 1; i < maxLength; i++) {
\r
314 cons2[i][23]++; // gap count
\r
322 /* for (int i=start; i <= end; i++) {
\r
327 for (int j=0;j<24;j++) {
\r
328 if (cons2[i][j] > max) {
\r
339 public void findQuality(int start, int end) {
\r
340 quality = new Vector();
\r
342 double max = -10000;
\r
344 int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
\r
346 //Loop over columns // JBPNote Profiling info
\r
347 // long ts = System.currentTimeMillis();
\r
348 //long te = System.currentTimeMillis();
\r
349 percentIdentity2();
\r
351 int size = seqNums.size();
\r
352 int[] lengths = new int[size];
\r
354 for (int l = 0; l < size; l++)
\r
355 lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
\r
357 for (int j = start; j <= end; j++) {
\r
360 // First Xr = depends on column only
\r
361 double[] x = new double[24];
\r
363 for (int ii = 0; ii < 24; ii++) {
\r
367 for (int i2 = 0; i2 < 24; i2++) {
\r
368 x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) +
\r
371 } catch (Exception e) {
\r
372 System.err.println("Exception during quality calculation.");
\r
373 e.printStackTrace();
\r
376 //System.out.println("X " + ii + " " + x[ii]);
\r
379 //System.out.println("X " + ii + " " + x[ii]);
\r
382 // Now calculate D for each position and sum
\r
383 for (int k = 0; k < size; k++) {
\r
385 double[] xx = new double[24];
\r
386 int seqNum = (j < lengths[k])
\r
387 ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
\r
389 // This is a loop over r
\r
390 for (int i = 0; i < 23; i++) {
\r
394 sr = (double) BLOSUM62[i][seqNum] + 4;
\r
395 } catch (Exception e) {
\r
396 System.out.println("Exception in sr: " + e);
\r
397 e.printStackTrace();
\r
400 //Calculate X with another loop over residues
\r
401 // System.out.println("Xi " + i + " " + x[i] + " " + sr);
\r
404 tot += (xx[i] * xx[i]);
\r
407 bigtot += Math.sqrt(tot);
\r
410 // This is the quality for one column
\r
411 if (max < bigtot) {
\r
415 // bigtot = bigtot * (size-cons2[j][23])/size;
\r
416 quality.addElement(new Double(bigtot));
\r
420 // Need to normalize by gaps
\r
423 double newmax = -10000;
\r
425 for (int j = start; j <= end; j++) {
\r
426 double tmp = ((Double) quality.elementAt(j)).doubleValue();
\r
427 tmp = ((max - tmp) * (size - cons2[j][23])) / size;
\r
429 // System.out.println(tmp+ " " + j);
\r
430 quality.setElementAt(new Double(tmp), j);
\r
432 if (tmp > newmax) {
\r
437 // System.out.println("Quality " + s);
\r
438 qualityRange[0] = new Double(0);
\r
439 qualityRange[1] = new Double(newmax);
\r