formatting
[jalview.git] / src / jalview / analysis / StructureFrequency.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.7)
3  * Copyright (C) 2011 J Procter, AM Waterhouse, J Engelhardt, LM Lui, 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
19 package jalview.analysis;
20
21 import java.util.*;
22
23 import jalview.datamodel.*;
24
25 /**
26  * Takes in a vector or array of sequences and column start and column end and
27  * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied.
28  * This class is used extensively in calculating alignment colourschemes that
29  * depend on the amount of conservation in each alignment column.
30  * 
31  * @author $author$
32  * @version $Revision$
33  */
34 public class StructureFrequency
35 {
36   // No need to store 1000s of strings which are not
37   // visible to the user.
38   public static final String MAXCOUNT = "C";
39
40   public static final String MAXRESIDUE = "R";
41
42   public static final String PID_GAPS = "G";
43
44   public static final String PID_NOGAPS = "N";
45
46   public static final String PROFILE = "P";
47
48   public static final String PAIRPROFILE = "B";
49
50   /**
51    * Returns the 3' position of a base pair
52    * 
53    * @param pairs
54    *          Secondary structure annotation
55    * @param indice
56    *          5' position of a base pair
57    * @return 3' position of a base pair
58    */
59   public static int findPair(SequenceFeature[] pairs, int indice)
60   {
61     for (int i = 0; i < pairs.length; i++)
62     {
63       if (pairs[i].getBegin() == indice)
64       {
65         return pairs[i].getEnd();
66       }
67     }
68     return -1;
69   }
70
71   /**
72    * Method to calculate a 'base pair consensus row', very similar to nucleotide
73    * consensus but takes into account a given structure
74    * 
75    * @param sequences
76    * @param start
77    * @param end
78    * @param result
79    * @param profile
80    * @param rnaStruc
81    */
82   public static final void calculate(SequenceI[] sequences, int start,
83           int end, Hashtable[] result, boolean profile,
84           AlignmentAnnotation rnaStruc)
85   {
86     Hashtable residueHash;
87     String maxResidue;
88     char[] seq, struc = rnaStruc.getRNAStruc().toCharArray();
89     SequenceFeature[] rna = rnaStruc._rnasecstr;
90     char c, s, cEnd;
91     int count, nonGap = 0, i, bpEnd = -1, j, jSize = sequences.length;
92     int[] values;
93     int[][] pairs;
94     float percentage;
95
96     for (i = start; i < end; i++) // foreach column
97     {
98       residueHash = new Hashtable();
99       maxResidue = "-";
100       values = new int[255];
101       pairs = new int[255][255];
102       bpEnd = -1;
103       if (i < struc.length)
104       {
105         s = struc[i];
106       }
107       else
108       {
109         s = '-';
110       }
111       if (s == '.' || s == ' ')
112       {
113         s = '-';
114       }
115
116       if (s != '(')
117       {
118         if (s == '-')
119         {
120           values['-']++;
121         }
122       }
123       else
124       {
125         bpEnd = findPair(rna, i);
126         if (bpEnd>-1)
127         {
128         for (j = 0; j < jSize; j++) // foreach row
129         {
130           if (sequences[j] == null)
131           {
132             System.err
133                     .println("WARNING: Consensus skipping null sequence - possible race condition.");
134             continue;
135           }
136           c = sequences[j].getCharAt(i);
137           {
138
139             // standard representation for gaps in sequence and structure
140             if (c == '.' || c == ' ')
141             {
142               c = '-';
143             }
144
145             if (c == '-')
146             {
147               values['-']++;
148               continue;
149             }
150             cEnd = sequences[j].getCharAt(bpEnd);
151             if (checkBpType(c, cEnd))
152             {
153               values['(']++; // H means it's a helix (structured)
154             }
155             pairs[c][cEnd]++;
156
157             maxResidue = "(";
158           }
159         }
160         }
161         // nonGap++;
162       }
163       // UPDATE this for new values
164       if (profile)
165       {
166         residueHash.put(PROFILE, new int[][]
167         { values, new int[]
168         { jSize, (jSize - values['-']) } });
169
170         residueHash.put(PAIRPROFILE, pairs);
171       }
172
173       count = values['('];
174
175       residueHash.put(MAXCOUNT, new Integer(count));
176       residueHash.put(MAXRESIDUE, maxResidue);
177
178       percentage = ((float) count * 100) / (float) jSize;
179       residueHash.put(PID_GAPS, new Float(percentage));
180
181       // percentage = ((float) count * 100) / (float) nongap;
182       // residueHash.put(PID_NOGAPS, new Float(percentage));
183       if (result[i] == null)
184       {
185         result[i] = residueHash;
186       }
187       if (bpEnd > 0)
188       {
189         values[')'] = values['('];
190         values['('] = 0;
191
192         residueHash = new Hashtable();
193         maxResidue = ")";
194
195         if (profile)
196         {
197           residueHash.put(PROFILE, new int[][]
198           { values, new int[]
199           { jSize, (jSize - values['-']) } });
200
201           residueHash.put(PAIRPROFILE, pairs);
202         }
203
204         residueHash.put(MAXCOUNT, new Integer(count));
205         residueHash.put(MAXRESIDUE, maxResidue);
206
207         percentage = ((float) count * 100) / (float) jSize;
208         residueHash.put(PID_GAPS, new Float(percentage));
209
210         result[bpEnd] = residueHash;
211       }
212     }
213   }
214
215   /**
216    * Method to check if a base-pair is a canonical or a wobble bp
217    * 
218    * @param up
219    *          5' base
220    * @param down
221    *          3' base
222    * @return True if it is a canonical/wobble bp
223    */
224   public static boolean checkBpType(char up, char down)
225   {
226     if (up > 'Z')
227     {
228       up -= 32;
229     }
230     if (down > 'Z')
231     {
232       down -= 32;
233     }
234
235     switch (up)
236     {
237     case 'A':
238       switch (down)
239       {
240       case 'T':
241         return true;
242       case 'U':
243         return true;
244       }
245       break;
246     case 'C':
247       switch (down)
248       {
249       case 'G':
250         return true;
251       }
252       break;
253     case 'T':
254       switch (down)
255       {
256       case 'A':
257         return true;
258       case 'G':
259         return true;
260       }
261       break;
262     case 'G':
263       switch (down)
264       {
265       case 'C':
266         return true;
267       case 'T':
268         return true;
269       case 'U':
270         return true;
271       }
272       break;
273     case 'U':
274       switch (down)
275       {
276       case 'A':
277         return true;
278       case 'G':
279         return true;
280       }
281       break;
282     }
283     return false;
284   }
285
286   /**
287    * Compute all or part of the annotation row from the given consensus
288    * hashtable
289    * 
290    * @param consensus
291    *          - pre-allocated annotation row
292    * @param hconsensus
293    * @param iStart
294    * @param width
295    * @param ignoreGapsInConsensusCalculation
296    * @param includeAllConsSymbols
297    */
298   public static void completeConsensus(AlignmentAnnotation consensus,
299           Hashtable[] hconsensus, int iStart, int width,
300           boolean ignoreGapsInConsensusCalculation,
301           boolean includeAllConsSymbols)
302   {
303     float tval, value;
304     if (consensus == null || consensus.annotations == null
305             || consensus.annotations.length < width)
306     {
307       // called with a bad alignment annotation row - wait for it to be
308       // initialised properly
309       return;
310     }
311     for (int i = iStart; i < width; i++)
312     {
313       if (i >= hconsensus.length)
314       {
315         // happens if sequences calculated over were shorter than alignment
316         // width
317         consensus.annotations[i] = null;
318         continue;
319       }
320       value = 0;
321       if (ignoreGapsInConsensusCalculation)
322       {
323         value = ((Float) hconsensus[i].get(StructureFrequency.PID_NOGAPS))
324                 .floatValue();
325       }
326       else
327       {
328         value = ((Float) hconsensus[i].get(StructureFrequency.PID_GAPS))
329                 .floatValue();
330       }
331
332       String maxRes = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
333               .toString();
334       String mouseOver = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
335               + " ";
336       if (maxRes.length() > 1)
337       {
338         mouseOver = "[" + maxRes + "] ";
339         maxRes = "+";
340       }
341       int[][] profile = (int[][]) hconsensus[i]
342               .get(StructureFrequency.PROFILE);
343       int[][] pairs = (int[][]) hconsensus[i]
344               .get(StructureFrequency.PAIRPROFILE);
345
346       if (pairs != null && includeAllConsSymbols) // Just responsible for the
347       // tooltip
348       // TODO Update tooltips for Structure row
349       {
350         mouseOver = "";
351
352         /* TODO It's not sure what is the purpose of the alphabet and wheter it is useful for structure?
353          * 
354          * if (alphabet != null) { for (int c = 0; c < alphabet.length; c++) {
355          * tval = ((float) profile[0][alphabet[c]]) 100f / (float)
356          * profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; mouseOver +=
357          * ((c == 0) ? "" : "; ") + alphabet[c] + " " + ((int) tval) + "%"; } }
358          * else {
359          */
360         Object[] ca = new Object[625];
361         float[] vl = new float[625];
362         int x = 0;
363         for (int c = 65; c < 90; c++)
364         {
365           for (int d = 65; d < 90; d++)
366           {
367             ca[x] = new int[]
368             { c, d };
369             vl[x] = (float) pairs[c][d];
370             x++;
371           }
372         }
373         jalview.util.QuickSort.sort(vl, ca);
374         int p = 0;
375
376         for (int c = 624; c > 0; c--)
377         {
378           if (vl[c] > 0)
379           {
380             tval = ((float) vl[c] * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
381                     : 0]);
382             mouseOver += ((p == 0) ? "" : "; ") + (char) ((int[]) ca[c])[0]
383                     + (char) ((int[]) ca[c])[1] + " " + ((int) tval) + "%";
384             p++;
385
386           }
387         }
388
389         // }
390       }
391       else
392       {
393         mouseOver += ((int) value + "%");
394       }
395       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
396               value);
397     }
398   }
399
400   /**
401    * get the sorted base-pair profile for the given position of the consensus
402    * 
403    * @param hconsensus
404    * @return profile of the given column
405    */
406   public static int[] extractProfile(Hashtable hconsensus,
407           boolean ignoreGapsInConsensusCalculation)
408   {
409     int[] rtnval = new int[74]; // 2*(5*5)+2
410     int[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
411     int[][] pairs = (int[][]) hconsensus
412             .get(StructureFrequency.PAIRPROFILE);
413
414     if (profile == null)
415       return null;
416
417     // TODO fix the object length, also do it in completeConsensus
418     Object[] ca = new Object[625];
419     float[] vl = new float[625];
420     int x = 0;
421     for (int c = 65; c < 90; c++)
422     {
423       for (int d = 65; d < 90; d++)
424       {
425         ca[x] = new int[]
426         { c, d };
427         vl[x] = (float) pairs[c][d];
428         x++;
429       }
430     }
431     jalview.util.QuickSort.sort(vl, ca);
432
433     rtnval[0] = 2;
434     rtnval[1] = 0;
435     for (int c = 624; c > 0; c--)
436     {
437       if (vl[c] > 0)
438       {
439         rtnval[rtnval[0]++] = ((int[]) ca[c])[0];
440         rtnval[rtnval[0]++] = ((int[]) ca[c])[1];
441         rtnval[rtnval[0]] = (int) ((float) vl[c] * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
442                 : 0]);
443         rtnval[1]+=rtnval[rtnval[0]++];
444       }
445     }
446
447     return rtnval;
448   }
449
450   public static void main(String args[])
451   {
452     // Short test to see if checkBpType works
453     ArrayList<String> test = new ArrayList<String>();
454     test.add("A");
455     test.add("c");
456     test.add("g");
457     test.add("T");
458     test.add("U");
459     for (String i : test)
460     {
461       for (String j : test)
462       {
463         System.out.println(i + "-" + j + ": "
464                 + StructureFrequency.checkBpType(i.charAt(0), j.charAt(0)));
465       }
466     }
467   }
468 }