Jalview 2.6 source licence
[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 || consensus.annotations.length<width)
208     {
209       // called with a bad alignment annotation row - wait for it to be initialised properly
210       return;
211     }
212     for (int i = iStart; i < width; i++)
213     {
214       if (i>=hconsensus.length) {
215         // happens if sequences calculated over were shorter than alignment width
216         consensus.annotations[i]=null;
217         continue;
218       }
219       value = 0;
220       if (ignoreGapsInConsensusCalculation)
221       {
222         value = ((Float) hconsensus[i].get(AAFrequency.PID_NOGAPS))
223                 .floatValue();
224       }
225       else
226       {
227         value = ((Float) hconsensus[i].get(AAFrequency.PID_GAPS))
228                 .floatValue();
229       }
230
231       String maxRes = hconsensus[i].get(AAFrequency.MAXRESIDUE).toString();
232       String mouseOver = hconsensus[i].get(AAFrequency.MAXRESIDUE) + " ";
233       if (maxRes.length() > 1)
234       {
235         mouseOver = "[" + maxRes + "] ";
236         maxRes = "+";
237       }
238       int[][] profile = (int[][]) hconsensus[i].get(AAFrequency.PROFILE);
239       if (profile != null && includeAllConsSymbols)
240       {
241         mouseOver = "";
242         if (alphabet != null)
243         {
244           for (int c = 0; c < alphabet.length; c++)
245           {
246             tval = ((float) profile[0][alphabet[c]])
247                     * 100f
248                     / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
249                             : 0];
250             mouseOver += ((c == 0) ? "" : "; ") + alphabet[c] + " "
251                     + ((int) tval) + "%";
252           }
253         }
254         else
255         {
256           Object[] ca = new Object[profile[0].length];
257           float[] vl = new float[profile[0].length];
258           for (int c = 0; c < ca.length; c++)
259           {
260             ca[c] = new char[]
261             { (char) c };
262             vl[c] = (float) profile[0][c];
263           }
264           ;
265           jalview.util.QuickSort.sort(vl, ca);
266           for (int p = 0, c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
267           {
268             if (((char[]) ca[c])[0] != '-')
269             {
270               tval = ((float) profile[0][((char[]) ca[c])[0]])
271                       * 100f
272                       / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
273                               : 0];
274               mouseOver += ((p == 0) ? "" : "; ") + ((char[]) ca[c])[0]
275                       + " " + ((int) tval) + "%";
276               p++;
277
278             }
279           }
280
281         }
282       }
283       else
284       {
285         mouseOver += ((int) value + "%");
286       }
287       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
288               value);
289     }
290   }
291
292   /**
293    * get the sorted profile for the given position of the consensus
294    * 
295    * @param hconsensus
296    * @return
297    */
298   public static int[] extractProfile(Hashtable hconsensus,
299           boolean ignoreGapsInConsensusCalculation)
300   {
301     int[] rtnval = new int[64];
302     int[][] profile = (int[][]) hconsensus.get(AAFrequency.PROFILE);
303     if (profile == null)
304       return null;
305     Object[] ca = new Object[profile[0].length];
306     float[] vl = new float[profile[0].length];
307     for (int c = 0; c < ca.length; c++)
308     {
309       ca[c] = new char[]
310       { (char) c };
311       vl[c] = (float) profile[0][c];
312     }
313     ;
314     jalview.util.QuickSort.sort(vl, ca);
315     rtnval[0] = 1;
316     for (int c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--)
317     {
318       if (((char[]) ca[c])[0] != '-')
319       {
320         rtnval[rtnval[0]++] = ((char[]) ca[c])[0];
321         rtnval[rtnval[0]++] = (int) (((float) profile[0][((char[]) ca[c])[0]]) * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
322                 : 0]);
323       }
324     }
325     return rtnval;
326   }
327 }