4bff6ef10941d5d30b1103bedadee91872e04c81
[jalview.git] / src / jalview / analysis / AAFrequency.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)
3  * Copyright (C) 2010 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
10  * 
11  * Jalview is distributed in the hope that it will be useful, but 
12  * WITHOUT ANY WARRANTY; without even the implied warranty 
13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14  * PURPOSE.  See the GNU General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 package jalview.analysis;
19
20 import java.util.*;
21
22 import jalview.datamodel.*;
23
24 /**
25  * Takes in a vector or array of sequences and column start and column end and
26  * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied.
27  * This class is used extensively in calculating alignment colourschemes that
28  * depend on the amount of conservation in each alignment column.
29  * 
30  * @author $author$
31  * @version $Revision$
32  */
33 public class AAFrequency
34 {
35   // No need to store 1000s of strings which are not
36   // visible to the user.
37   public static final String MAXCOUNT = "C";
38
39   public static final String MAXRESIDUE = "R";
40
41   public static final String PID_GAPS = "G";
42
43   public static final String PID_NOGAPS = "N";
44
45   public static final String PROFILE = "P";
46
47   public static final Hashtable[] calculate(Vector sequences, int start,
48           int end)
49   {
50     return calculate(sequences, start, end, false);
51   }
52
53   public static final Hashtable[] calculate(Vector sequences, int start,
54           int end, boolean profile)
55   {
56     SequenceI[] seqs = new SequenceI[sequences.size()];
57     int width = 0;
58     for (int i = 0; i < sequences.size(); i++)
59     {
60       seqs[i] = (SequenceI) sequences.elementAt(i);
61       if (seqs[i].getLength() > width)
62       {
63         width = seqs[i].getLength();
64       }
65     }
66
67     Hashtable[] reply = new Hashtable[width];
68
69     if (end >= width)
70     {
71       end = width;
72     }
73
74     calculate(seqs, start, end, reply, profile);
75
76     return reply;
77   }
78
79   public static final void calculate(SequenceI[] sequences, int start,
80           int end, Hashtable[] result)
81   {
82     calculate(sequences, start, end, result, false);
83   }
84
85   public static final void calculate(SequenceI[] sequences, int start,
86           int end, Hashtable[] result, boolean profile)
87   {
88     Hashtable residueHash;
89     int maxCount, nongap, i, j, v, jSize = sequences.length;
90     String maxResidue;
91     char c;
92     float percentage;
93
94     int[] values = new int[255];
95
96     char[] seq;
97
98     for (i = start; i < end; i++)
99     {
100       residueHash = new Hashtable();
101       maxCount = 0;
102       maxResidue = "";
103       nongap = 0;
104       values = new int[255];
105
106       for (j = 0; j < jSize; j++)
107       {
108         seq = sequences[j].getSequence();
109         if (seq.length > i)
110         {
111           c = seq[i];
112
113           if (c == '.' || c == ' ')
114           {
115             c = '-';
116           }
117
118           if (c == '-')
119           {
120             values['-']++;
121             continue;
122           }
123           else if ('a' <= c && c <= 'z')
124           {
125             c -= 32; // ('a' - 'A');
126           }
127
128           nongap++;
129           values[c]++;
130
131         }
132         else
133         {
134           values['-']++;
135         }
136       }
137
138       for (v = 'A'; v < 'Z'; v++)
139       {
140         if (values[v] < 2 || values[v] < maxCount)
141         {
142           continue;
143         }
144
145         if (values[v] > maxCount)
146         {
147           maxResidue = String.valueOf((char) v);
148         }
149         else if (values[v] == maxCount)
150         {
151           maxResidue += String.valueOf((char) v);
152         }
153         maxCount = values[v];
154       }
155
156       if (maxResidue.length() == 0)
157       {
158         maxResidue = "-";
159       }
160       if (profile)
161       {
162         residueHash.put(PROFILE, new int[][]
163         { values, new int[]
164         { jSize, nongap } });
165       }
166       residueHash.put(MAXCOUNT, new Integer(maxCount));
167       residueHash.put(MAXRESIDUE, maxResidue);
168
169       percentage = ((float) maxCount * 100) / (float) jSize;
170       residueHash.put(PID_GAPS, new Float(percentage));
171
172       percentage = ((float) maxCount * 100) / (float) nongap;
173       residueHash.put(PID_NOGAPS, new Float(percentage));
174       result[i] = residueHash;
175     }
176   }
177
178   /**
179    * Compute all or part of the annotation row from the given consensus
180    * hashtable
181    * 
182    * @param consensus
183    *          - pre-allocated annotation row
184    * @param hconsensus
185    * @param iStart
186    * @param width
187    * @param ignoreGapsInConsensusCalculation
188    * @param includeAllConsSymbols
189    */
190   public static void completeConsensus(AlignmentAnnotation consensus,
191           Hashtable[] hconsensus, int iStart, int width,
192           boolean ignoreGapsInConsensusCalculation,
193           boolean includeAllConsSymbols)
194   {
195     completeConsensus(consensus, hconsensus, iStart, width,
196             ignoreGapsInConsensusCalculation, includeAllConsSymbols, null); // new
197                                                                             // char[]
198     // { 'A', 'C', 'G', 'T', 'U' });
199   }
200
201   public static void completeConsensus(AlignmentAnnotation consensus,
202           Hashtable[] hconsensus, int iStart, int width,
203           boolean ignoreGapsInConsensusCalculation,
204           boolean includeAllConsSymbols, char[] alphabet)
205   {
206     float tval, value;
207     if (consensus == null || consensus.annotations == null
208             || consensus.annotations.length < width)
209     {
210       // called with a bad alignment annotation row - wait for it to be
211       // initialised properly
212       return;
213     }
214     for (int i = iStart; i < width; i++)
215     {
216       if (i >= hconsensus.length)
217       {
218         // happens if sequences calculated over were shorter than alignment
219         // width
220         consensus.annotations[i] = null;
221         continue;
222       }
223       value = 0;
224       if (ignoreGapsInConsensusCalculation)
225       {
226         value = ((Float) hconsensus[i].get(AAFrequency.PID_NOGAPS))
227                 .floatValue();
228       }
229       else
230       {
231         value = ((Float) hconsensus[i].get(AAFrequency.PID_GAPS))
232                 .floatValue();
233       }
234
235       String maxRes = hconsensus[i].get(AAFrequency.MAXRESIDUE).toString();
236       String mouseOver = hconsensus[i].get(AAFrequency.MAXRESIDUE) + " ";
237       if (maxRes.length() > 1)
238       {
239         mouseOver = "[" + maxRes + "] ";
240         maxRes = "+";
241       }
242       int[][] profile = (int[][]) hconsensus[i].get(AAFrequency.PROFILE);
243       if (profile != null && includeAllConsSymbols)
244       {
245         mouseOver = "";
246         if (alphabet != null)
247         {
248           for (int c = 0; c < alphabet.length; c++)
249           {
250             tval = ((float) profile[0][alphabet[c]])
251                     * 100f
252                     / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
253                             : 0];
254             mouseOver += ((c == 0) ? "" : "; ") + alphabet[c] + " "
255                     + ((int) tval) + "%";
256           }
257         }
258         else
259         {
260           Object[] ca = new Object[profile[0].length];
261           float[] vl = new float[profile[0].length];
262           for (int c = 0; c < ca.length; c++)
263           {
264             ca[c] = new char[]
265             { (char) c };
266             vl[c] = (float) profile[0][c];
267           }
268           ;
269           jalview.util.QuickSort.sort(vl, ca);
270           for (int p = 0, c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
271           {
272             if (((char[]) ca[c])[0] != '-')
273             {
274               tval = ((float) profile[0][((char[]) ca[c])[0]])
275                       * 100f
276                       / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
277                               : 0];
278               mouseOver += ((p == 0) ? "" : "; ") + ((char[]) ca[c])[0]
279                       + " " + ((int) tval) + "%";
280               p++;
281
282             }
283           }
284
285         }
286       }
287       else
288       {
289         mouseOver += ((int) value + "%");
290       }
291       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
292               value);
293     }
294   }
295
296   /**
297    * get the sorted profile for the given position of the consensus
298    * 
299    * @param hconsensus
300    * @return
301    */
302   public static int[] extractProfile(Hashtable hconsensus,
303           boolean ignoreGapsInConsensusCalculation)
304   {
305     int[] rtnval = new int[64];
306     int[][] profile = (int[][]) hconsensus.get(AAFrequency.PROFILE);
307     if (profile == null)
308       return null;
309     Object[] ca = new Object[profile[0].length];
310     float[] vl = new float[profile[0].length];
311     for (int c = 0; c < ca.length; c++)
312     {
313       ca[c] = new char[]
314       { (char) c };
315       vl[c] = (float) profile[0][c];
316     }
317     ;
318     jalview.util.QuickSort.sort(vl, ca);
319     rtnval[0] = 1;
320     for (int c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
321     {
322       if (((char[]) ca[c])[0] != '-')
323       {
324         rtnval[rtnval[0]++] = ((char[]) ca[c])[0];
325         rtnval[rtnval[0]++] = (int) (((float) profile[0][((char[]) ca[c])[0]]) * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
326                 : 0]);
327       }
328     }
329     return rtnval;
330   }
331 }