1 /* Jalview - a java multiple alignment editor
\r
2 * Copyright (C) 1998 Michele Clamp
\r
4 * This program is free software; you can redistribute it and/or
\r
5 * modify it under the terms of the GNU General Public License
\r
6 * as published by the Free Software Foundation; either version 2
\r
7 * of the License, or (at your option) any later version.
\r
9 * This program is distributed in the hope that it will be useful,
\r
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
\r
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
\r
12 * GNU General Public License for more details.
\r
14 * You should have received a copy of the GNU General Public License
\r
15 * along with this program; if not, write to the Free Software
\r
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
\r
18 package jalview.analysis;
\r
22 import jalview.gui.*;
\r
23 import jalview.datamodel.*;
\r
26 public class Conservation {
\r
30 Vector seqNums; // vector of int vectors where first is sequence checksum
\r
31 int maxLength=0; // used by quality calcs
\r
32 boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
\r
33 private void calcSeqNums() {
\r
34 for (int i=0; i<sequences.size(); i++) {
\r
38 private void calcSeqNum(int i) {
\r
39 String sq=null; // for dumb jbuilder not-inited exception warning
\r
41 if (i>-1 && i<sequences.size()) {
\r
42 sq = ((SequenceI)sequences.elementAt(i)).getSequence();
\r
43 if (seqNums.size()<=i) {
\r
44 seqNums.addElement(new int[sq.length()+1]);
\r
46 if (sq.hashCode()!=((int[])seqNums.elementAt(i))[0])
\r
49 seqNumsChanged = true;
\r
50 sq = ((SequenceI) sequences.elementAt(i)).getSequence();
\r
52 if (maxLength < len)
\r
54 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
\r
55 sqnum[0] = sq.hashCode();
\r
56 for (j = 1; j <= len; j++)
\r
58 sqnum[j] = ( (Integer) jalview.schemes.ResidueProperties.aaHash.get(new
\r
59 String(sq.substring(j - 1, j)))).intValue(); // yuk
\r
61 seqNums.setElementAt(sqnum, i);
\r
64 // JBPNote INFO level debug
\r
65 System.err.println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
\r
68 Vector total = new Vector();
\r
69 public Vector quality;
\r
70 public Double[] qualityRange = new Double[2];
\r
71 String consString = "";
\r
73 Sequence consSequence;
\r
79 public Conservation(String name,Hashtable propHash, int threshold, Vector sequences, int start, int end) {
\r
81 this.propHash = propHash;
\r
82 this.threshold = threshold;
\r
83 this.sequences = sequences;
\r
86 seqNums = new Vector(sequences.size());
\r
91 public void calculate() {
\r
93 for (int i = start;i <= end; i++) {
\r
94 Hashtable resultHash = null;
\r
95 Hashtable residueHash = null;
\r
97 resultHash = new Hashtable();
\r
98 residueHash = new Hashtable();
\r
100 for (int j=0; j < sequences.size(); j++) {
\r
101 // JBPNote - have to make sure elements of the sequences vector
\r
102 // are tested like this everywhere...
\r
103 if (sequences.elementAt(j) instanceof Sequence) {
\r
104 Sequence s = (Sequence)sequences.elementAt(j);
\r
106 if (s.getLength() > i) {
\r
107 String res = s.getSequence().substring(i,i+1);
\r
109 if (residueHash.containsKey(res)) {
\r
110 int count = ((Integer)residueHash.get(res)).intValue() ;
\r
112 residueHash.put(res,new Integer(count));
\r
114 residueHash.put(res,new Integer(1));
\r
117 if (residueHash.containsKey("-")) {
\r
118 int count = ((Integer)residueHash.get("-")).intValue() ;
\r
120 residueHash.put("-",new Integer(count));
\r
122 residueHash.put("-",new Integer(1));
\r
128 //What is the count threshold to count the residues in residueHash()
\r
129 int thresh = threshold*(sequences.size())/100;
\r
131 //loop over all the found residues
\r
132 Enumeration e = residueHash.keys();
\r
134 while (e.hasMoreElements()) {
\r
136 String res = (String)e.nextElement();
\r
137 if (((Integer)residueHash.get(res)).intValue() > thresh) {
\r
139 //Now loop over the properties
\r
140 Enumeration e2 = propHash.keys();
\r
142 while (e2.hasMoreElements()) {
\r
143 String type = (String)e2.nextElement();
\r
144 Hashtable ht = (Hashtable)propHash.get(type);
\r
146 //Have we ticked this before?
\r
147 if (! resultHash.containsKey(type)) {
\r
148 if (ht.containsKey(res)) {
\r
149 resultHash.put(type,ht.get(res));
\r
151 resultHash.put(type,ht.get("-"));
\r
153 } else if ( ((Integer)resultHash.get(type)).equals((Integer)ht.get(res)) == false) {
\r
154 resultHash.put(type,new Integer(-1));
\r
160 total.addElement(resultHash);
\r
164 public int countGaps(int j)
\r
168 for (int i = 0; i < sequences.size();i++)
\r
170 if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length())
\r
171 { count++; continue;}
\r
173 char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);
\r
174 if (jalview.util.Comparison.isGap((c)))
\r
182 * returns gap count in int[0], and conserved residue count in int[1]
\r
184 public int[] countConsNGaps(int j)
\r
189 int[] r = new int[2];
\r
191 for (int i = 0; i < sequences.size();i++)
\r
193 if( j >= ((Sequence)sequences.elementAt(i)).getSequence().length())
\r
197 char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);
\r
198 if (jalview.util.Comparison.isGap((c)))
\r
211 r[0] = (nres==cons) ? 1 : 0;
\r
218 public void verdict(boolean consflag, float percentageGaps) {
\r
219 String consString = "";
\r
221 for (int i=start; i <= end; i++) {
\r
222 int[] gapcons = countConsNGaps(i);
\r
223 boolean cons = (gapcons[0]==1) ? true : false;
\r
224 int totGaps = gapcons[1];
\r
225 float pgaps = (float)totGaps*100/(float)sequences.size();
\r
226 // System.out.println("percentage gaps = "+pgaps+"\n");
\r
227 if (percentageGaps > pgaps)
\r
229 Hashtable resultHash = (Hashtable)total.elementAt(i-start);
\r
231 //Now find the verdict
\r
233 Enumeration e3 = resultHash.keys();
\r
235 while (e3.hasMoreElements())
\r
237 String type = (String)e3.nextElement();
\r
238 Integer result = (Integer)resultHash.get(type);
\r
240 //Do we want to count +ve conservation or +ve and -ve cons.?
\r
244 if (result.intValue() == 1)
\r
249 if (result.intValue() != -1)
\r
255 consString = consString + String.valueOf(count);// Conserved props!=Identity
\r
257 consString = consString + ((gapcons[0]==1) ? "*" : "+");
\r
262 consString = consString + "-";
\r
266 consSequence = new Sequence(name,consString,start,end);
\r
269 public Sequence getConsSequence() {
\r
270 return consSequence;
\r
273 // From Alignment.java in jalview118
\r
275 public void findQuality() {
\r
276 findQuality(0,maxLength-1);
\r
281 private void percentIdentity2() {
\r
282 calcSeqNums(); // updates maxLength, too.
\r
283 if (cons2==null || seqNumsChanged) {
\r
284 cons2 = new int[maxLength][24];
\r
285 // Initialize the array
\r
286 for (int j=0;j<24;j++) {
\r
287 for (int i=0; i < maxLength;i++) {
\r
294 while(j < sequences.size()) {
\r
295 sqnum=(int[])seqNums.elementAt(j);
\r
296 for (int i = 1; i < sqnum.length; i++) {
\r
297 cons2[i-1][sqnum[i]]++;
\r
299 for (int i=sqnum.length-1; i<maxLength; i++) {
\r
300 cons2[i][23]++; // gap count
\r
306 /* for (int i=start; i <= end; i++) {
\r
311 for (int j=0;j<24;j++) {
\r
312 if (cons2[i][j] > max) {
\r
325 public void findQuality(int start, int end) {
\r
326 quality = new Vector();
\r
327 double max = -10000;
\r
329 int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
\r
330 //Loop over columns // JBPNote Profiling info
\r
331 // long ts = System.currentTimeMillis();
\r
332 //long te = System.currentTimeMillis();
\r
333 percentIdentity2();
\r
335 int size = seqNums.size();
\r
336 int[] lengths = new int[size];
\r
338 for (int l = 0; l < size; l++)
\r
339 lengths[l] = ((int[]) seqNums.elementAt(l)).length-1;
\r
341 for (int j=start; j <= end; j++) {
\r
344 // First Xr = depends on column only
\r
345 double x[] = new double[24];
\r
347 for (int ii=0; ii < 24; ii++) {
\r
350 for (int i2=0; i2 < 24; i2++) {
\r
351 x[ii] += (double)cons2[j][i2] * BLOSUM62[ii][i2]+4;
\r
353 } catch (Exception e) {
\r
354 System.err.println("Exception during quality calculation.");
\r
355 e.printStackTrace();
\r
357 //System.out.println("X " + ii + " " + x[ii]);
\r
359 //System.out.println("X " + ii + " " + x[ii]);
\r
361 // Now calculate D for each position and sum
\r
362 for (int k=0; k < size; k++) {
\r
364 double[] xx = new double[24];
\r
367 ? ((int[]) seqNums.elementAt(k))[j+1]
\r
368 : 23; // Sequence, or gap at the end
\r
370 // This is a loop over r
\r
371 for (int i=0; i < 23; i++) {
\r
374 sr = (double)BLOSUM62[i][seqNum]+4;
\r
375 } catch (Exception e) {
\r
376 System.out.println("Exception in sr: " + e);
\r
377 e.printStackTrace();
\r
379 //Calculate X with another loop over residues
\r
381 // System.out.println("Xi " + i + " " + x[i] + " " + sr);
\r
384 tot += xx[i]*xx[i];
\r
386 bigtot += Math.sqrt(tot);
\r
388 // This is the quality for one column
\r
389 if (max < bigtot) {max = bigtot;}
\r
390 // bigtot = bigtot * (size-cons2[j][23])/size;
\r
392 quality.addElement(new Double(bigtot));
\r
396 // Need to normalize by gaps
\r
398 double newmax=-10000;
\r
399 for (int j=start; j <= end; j++) {
\r
400 double tmp = ((Double)quality.elementAt(j)).doubleValue();
\r
401 tmp = (max - tmp)*(size-cons2[j][23])/size;
\r
402 // System.out.println(tmp+ " " + j);
\r
403 quality.setElementAt(new Double(tmp),j);
\r
407 // System.out.println("Quality " + s);
\r
408 qualityRange[0] = new Double(0);
\r
409 qualityRange[1] = new Double(newmax);
\r