Merge branch 'develop' into features/hmmer
[jalview.git] / src / jalview / analysis / AAFrequency.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
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
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.analysis;
22
23 import jalview.datamodel.AlignedCodonFrame;
24 import jalview.datamodel.AlignmentAnnotation;
25 import jalview.datamodel.AlignmentI;
26 import jalview.datamodel.Annotation;
27 import jalview.datamodel.HiddenMarkovModel;
28 import jalview.datamodel.Profile;
29 import jalview.datamodel.ProfileI;
30 import jalview.datamodel.Profiles;
31 import jalview.datamodel.ProfilesI;
32 import jalview.datamodel.ResidueCount;
33 import jalview.datamodel.ResidueCount.SymbolCounts;
34 import jalview.datamodel.SequenceI;
35 import jalview.ext.android.SparseIntArray;
36 import jalview.schemes.ResidueProperties;
37 import jalview.util.Comparison;
38 import jalview.util.Format;
39 import jalview.util.MappingUtils;
40 import jalview.util.QuickSort;
41
42 import java.awt.Color;
43 import java.util.Arrays;
44 import java.util.Hashtable;
45 import java.util.List;
46
47 /**
48  * Takes in a vector or array of sequences and column start and column end and
49  * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied.
50  * This class is used extensively in calculating alignment colourschemes that
51  * depend on the amount of conservation in each alignment column.
52  * 
53  * @author $author$
54  * @version $Revision$
55  */
56 public class AAFrequency
57 {
58   public static final String PROFILE = "P";
59
60   private static final String AMINO = "amino";
61
62   private static final String DNA = "DNA";
63
64   private static final String RNA = "RNA";
65
66   /*
67    * Quick look-up of String value of char 'A' to 'Z'
68    */
69   private static final String[] CHARS = new String['Z' - 'A' + 1];
70
71   static
72   {
73     for (char c = 'A'; c <= 'Z'; c++)
74     {
75       CHARS[c - 'A'] = String.valueOf(c);
76     }
77   }
78
79   public static final ProfilesI calculate(List<SequenceI> list, int start,
80           int end)
81   {
82     return calculate(list, start, end, false);
83   }
84
85   public static final ProfilesI calculate(List<SequenceI> sequences,
86           int start, int end, boolean profile)
87   {
88     SequenceI[] seqs = new SequenceI[sequences.size()];
89     int width = 0;
90     synchronized (sequences)
91     {
92       for (int i = 0; i < sequences.size(); i++)
93       {
94         seqs[i] = sequences.get(i);
95         int length = seqs[i].getLength();
96         if (length > width)
97         {
98           width = length;
99         }
100       }
101
102       if (end >= width)
103       {
104         end = width;
105       }
106
107       ProfilesI reply = calculate(seqs, width, start, end, profile);
108       return reply;
109     }
110   }
111
112
113
114   /**
115    * Calculate the consensus symbol(s) for each column in the given range.
116    * 
117    * @param sequences
118    * @param width
119    *          the full width of the alignment
120    * @param start
121    *          start column (inclusive, base zero)
122    * @param end
123    *          end column (exclusive)
124    * @param saveFullProfile
125    *          if true, store all symbol counts
126    */
127   public static final ProfilesI calculate(final SequenceI[] sequences,
128           int width, int start, int end, boolean saveFullProfile)
129   {
130     // long now = System.currentTimeMillis();
131     int seqCount = sequences.length;
132     boolean nucleotide = false;
133     int nucleotideCount = 0;
134     int peptideCount = 0;
135
136     ProfileI[] result = new ProfileI[width];
137
138     for (int column = start; column < end; column++)
139     {
140       /*
141        * Apply a heuristic to detect nucleotide data (which can
142        * be counted in more compact arrays); here we test for
143        * more than 90% nucleotide; recheck every 10 columns in case
144        * of misleading data e.g. highly conserved Alanine in peptide!
145        * Mistakenly guessing nucleotide has a small performance cost,
146        * as it will result in counting in sparse arrays.
147        * Mistakenly guessing peptide has a small space cost, 
148        * as it will use a larger than necessary array to hold counts. 
149        */
150       if (nucleotideCount > 100 && column % 10 == 0)
151       {
152         nucleotide = (9 * peptideCount < nucleotideCount);
153       }
154       ResidueCount residueCounts = new ResidueCount(nucleotide);
155
156       for (int row = 0; row < seqCount; row++)
157       {
158         if (sequences[row] == null)
159         {
160           System.err.println(
161                   "WARNING: Consensus skipping null sequence - possible race condition.");
162           continue;
163         }
164         char[] seq = sequences[row].getSequence();
165         if (seq.length > column)
166         {
167           char c = seq[column];
168           residueCounts.add(c);
169           if (Comparison.isNucleotide(c))
170           {
171             nucleotideCount++;
172           }
173           else if (!Comparison.isGap(c))
174           {
175             peptideCount++;
176           }
177         }
178         else
179         {
180           /*
181            * count a gap if the sequence doesn't reach this column
182            */
183           residueCounts.addGap();
184         }
185       }
186
187       int maxCount = residueCounts.getModalCount();
188       String maxResidue = residueCounts.getResiduesForCount(maxCount);
189       int gapCount = residueCounts.getGapCount();
190       ProfileI profile = new Profile(seqCount, gapCount, maxCount,
191               maxResidue);
192
193       if (saveFullProfile)
194       {
195         profile.setCounts(residueCounts);
196       }
197
198       result[column] = profile;
199     }
200     return new Profiles(result);
201     // long elapsed = System.currentTimeMillis() - now;
202     // System.out.println(elapsed);
203   }
204
205   /**
206    * Returns the full set of profiles for a hidden Markov model. The underlying
207    * data is the raw probabilities of a residue being emitted at each node,
208    * however the profiles returned by this function contain the percentage
209    * chance of a residue emission.
210    * 
211    * @param hmm
212    * @param width
213    *          The width of the Profile array (Profiles) to be returned.
214    * @param start
215    *          The alignment column on which the first profile is based.
216    * @param end
217    *          The alignment column on which the last profile is based.
218    * @param saveFullProfile
219    *          Flag for saving the counts for each profile
220    * @param removeBelowBackground
221    *          Flag for removing any characters with a match emission probability
222    *          less than its background frequency
223    * @return
224    */
225   public static ProfilesI calculateHMMProfiles(final HiddenMarkovModel hmm,
226           int width, int start, int end, boolean saveFullProfile,
227           boolean removeBelowBackground)
228   {
229     ProfileI[] result = new ProfileI[width];
230     int symbolCount = hmm.getNumberOfSymbols();
231     for (int column = start; column < end; column++)
232     {
233       ResidueCount counts = new ResidueCount();
234       for (char symbol : hmm.getSymbols())
235       {
236         int value = getAnalogueCount(hmm, column, symbol,
237                 removeBelowBackground);
238         counts.put(symbol, value);
239       }
240       int maxCount = counts.getModalCount();
241       String maxResidue = counts.getResiduesForCount(maxCount);
242       int gapCount = counts.getGapCount();
243       ProfileI profile = new Profile(symbolCount, gapCount, maxCount,
244               maxResidue);
245
246       if (saveFullProfile)
247       {
248         profile.setCounts(counts);
249       }
250
251       result[column] = profile;
252     }
253     return new Profiles(result);
254   }
255
256   /**
257    * Make an estimate of the profile size we are going to compute i.e. how many
258    * different characters may be present in it. Overestimating has a cost of
259    * using more memory than necessary. Underestimating has a cost of needing to
260    * extend the SparseIntArray holding the profile counts.
261    * 
262    * @param profileSizes
263    *          counts of sizes of profiles so far encountered
264    * @return
265    */
266   static int estimateProfileSize(SparseIntArray profileSizes)
267   {
268     if (profileSizes.size() == 0)
269     {
270       return 4;
271     }
272
273     /*
274      * could do a statistical heuristic here e.g. 75%ile
275      * for now just return the largest value
276      */
277     return profileSizes.keyAt(profileSizes.size() - 1);
278   }
279
280   /**
281    * Derive the consensus annotations to be added to the alignment for display.
282    * This does not recompute the raw data, but may be called on a change in
283    * display options, such as 'ignore gaps', which may in turn result in a
284    * change in the derived values.
285    * 
286    * @param consensus
287    *          the annotation row to add annotations to
288    * @param profiles
289    *          the source consensus data
290    * @param startCol
291    *          start column (inclusive)
292    * @param endCol
293    *          end column (exclusive)
294    * @param ignoreGaps
295    *          if true, normalise residue percentages ignoring gaps
296    * @param showSequenceLogo
297    *          if true include all consensus symbols, else just show modal
298    *          residue
299    * @param nseq
300    *          number of sequences
301    */
302   public static void completeConsensus(AlignmentAnnotation consensus,
303           ProfilesI profiles, int startCol, int endCol, boolean ignoreGaps,
304           boolean showSequenceLogo, long nseq)
305   {
306     // long now = System.currentTimeMillis();
307     if (consensus == null || consensus.annotations == null
308             || consensus.annotations.length < endCol)
309     {
310       /*
311        * called with a bad alignment annotation row 
312        * wait for it to be initialised properly
313        */
314       return;
315     }
316
317     for (int i = startCol; i < endCol; i++)
318     {
319       ProfileI profile = profiles.get(i);
320       if (profile == null)
321       {
322         /*
323          * happens if sequences calculated over were 
324          * shorter than alignment width
325          */
326         consensus.annotations[i] = null;
327         return;
328       }
329
330       final int dp = getPercentageDp(nseq);
331
332       float value = profile.getPercentageIdentity(ignoreGaps);
333
334       String description = getTooltip(profile, value, showSequenceLogo,
335               ignoreGaps, dp);
336
337       String modalResidue = profile.getModalResidue();
338       if ("".equals(modalResidue))
339       {
340         modalResidue = "-";
341       }
342       else if (modalResidue.length() > 1)
343       {
344         modalResidue = "+";
345       }
346       consensus.annotations[i] = new Annotation(modalResidue, description,
347               ' ', value);
348     }
349     // long elapsed = System.currentTimeMillis() - now;
350     // System.out.println(-elapsed);
351   }
352
353   /**
354    * Derive the information annotations to be added to the alignment for
355    * display. This does not recompute the raw data, but may be called on a
356    * change in display options, such as 'ignore below background frequency',
357    * which may in turn result in a change in the derived values.
358    * 
359    * @param information
360    *          the annotation row to add annotations to
361    * @param profiles
362    *          the source information data
363    * @param startCol
364    *          start column (inclusive)
365    * @param endCol
366    *          end column (exclusive)
367    * @param ignoreGaps
368    *          if true, normalise residue percentages 
369    * @param showSequenceLogo
370    *          if true include all information symbols, else just show modal
371    *          residue
372    * @param nseq
373    *          number of sequences
374    */
375   public static void completeInformation(AlignmentAnnotation information,
376           ProfilesI profiles, int startCol, int endCol,
377           boolean ignoreBelowBackground,
378           boolean showSequenceLogo, long nseq)
379   {
380     // long now = System.currentTimeMillis();
381     if (information == null || information.annotations == null
382             || information.annotations.length < endCol)
383     {
384       /*
385        * called with a bad alignment annotation row 
386        * wait for it to be initialised properly
387        */
388       return;
389     }
390
391     Float max = 0f;
392
393     for (int i = startCol; i < endCol; i++)
394     {
395       ProfileI profile = profiles.get(i);
396       if (profile == null)
397       {
398         /*
399          * happens if sequences calculated over were 
400          * shorter than alignment width
401          */
402         information.annotations[i] = null;
403         return;
404       }
405
406       HiddenMarkovModel hmm;
407       
408       SequenceI hmmSeq = information.sequenceRef;
409       
410       hmm = hmmSeq.getHMM();
411       
412       Float value = getInformationContent(i, hmm);
413
414       if (value > max)
415       {
416         max = value;
417       }
418
419       String description = value + " bits";
420       information.annotations[i] = new Annotation(
421               Character.toString(hmm.getConsensusAtAlignColumn(i)),
422               description, ' ', value);
423     }
424     information.graphMax = max;
425     // long elapsed = System.currentTimeMillis() - now;
426     // System.out.println(-elapsed);
427   }
428
429   /**
430    * Derive the gap count annotation row.
431    * 
432    * @param gaprow
433    *          the annotation row to add annotations to
434    * @param profiles
435    *          the source consensus data
436    * @param startCol
437    *          start column (inclusive)
438    * @param endCol
439    *          end column (exclusive)
440    */
441   public static void completeGapAnnot(AlignmentAnnotation gaprow,
442           ProfilesI profiles, int startCol, int endCol, long nseq)
443   {
444     if (gaprow == null || gaprow.annotations == null
445             || gaprow.annotations.length < endCol)
446     {
447       /*
448        * called with a bad alignment annotation row 
449        * wait for it to be initialised properly
450        */
451       return;
452     }
453     // always set ranges again
454     gaprow.graphMax = nseq;
455     gaprow.graphMin = 0;
456     double scale = 0.8 / nseq;
457     for (int i = startCol; i < endCol; i++)
458     {
459       ProfileI profile = profiles.get(i);
460       if (profile == null)
461       {
462         /*
463          * happens if sequences calculated over were 
464          * shorter than alignment width
465          */
466         gaprow.annotations[i] = null;
467         return;
468       }
469
470       final int gapped = profile.getNonGapped();
471
472       String description = "" + gapped;
473
474       gaprow.annotations[i] = new Annotation("", description, '\0', gapped,
475               jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY,
476                       (float) scale * gapped));
477     }
478   }
479
480   /**
481    * Returns a tooltip showing either
482    * <ul>
483    * <li>the full profile (percentages of all residues present), if
484    * showSequenceLogo is true, or</li>
485    * <li>just the modal (most common) residue(s), if showSequenceLogo is
486    * false</li>
487    * </ul>
488    * Percentages are as a fraction of all sequence, or only ungapped sequences
489    * if ignoreGaps is true.
490    * 
491    * @param profile
492    * @param pid
493    * @param showSequenceLogo
494    * @param ignoreGaps
495    * @param dp
496    *          the number of decimal places to format percentages to
497    * @return
498    */
499   static String getTooltip(ProfileI profile, float pid,
500           boolean showSequenceLogo, boolean ignoreGaps, int dp)
501   {
502     ResidueCount counts = profile.getCounts();
503
504     String description = null;
505     if (counts != null && showSequenceLogo)
506     {
507       int normaliseBy = ignoreGaps ? profile.getNonGapped()
508               : profile.getHeight();
509       description = counts.getTooltip(normaliseBy, dp);
510     }
511     else
512     {
513       StringBuilder sb = new StringBuilder(64);
514       String maxRes = profile.getModalResidue();
515       if (maxRes.length() > 1)
516       {
517         sb.append("[").append(maxRes).append("]");
518       }
519       else
520       {
521         sb.append(maxRes);
522       }
523       if (maxRes.length() > 0)
524       {
525         sb.append(" ");
526         Format.appendPercentage(sb, pid, dp);
527         sb.append("%");
528       }
529       description = sb.toString();
530     }
531     return description;
532   }
533
534   /**
535    * Returns the sorted profile for the given consensus data. The returned array
536    * contains
537    * 
538    * <pre>
539    *    [profileType, numberOfValues, nonGapCount, charValue1, percentage1, charValue2, percentage2, ...]
540    * in descending order of percentage value
541    * </pre>
542    * 
543    * @param profile
544    *          the data object from which to extract and sort values
545    * @param ignoreGaps
546    *          if true, only non-gapped values are included in percentage
547    *          calculations
548    * @return
549    */
550   public static int[] extractProfile(ProfileI profile, boolean ignoreGaps)
551   {
552     int[] rtnval = new int[64];
553     ResidueCount counts = profile.getCounts();
554     if (counts == null)
555     {
556       return null;
557     }
558
559     SymbolCounts symbolCounts = counts.getSymbolCounts();
560     char[] symbols = symbolCounts.symbols;
561     int[] values = symbolCounts.values;
562     QuickSort.sort(values, symbols);
563     int nextArrayPos = 2;
564     int totalPercentage = 0;
565     final int divisor = ignoreGaps ? profile.getNonGapped()
566             : profile.getHeight();
567
568     /*
569      * traverse the arrays in reverse order (highest counts first)
570      */
571     for (int i = symbols.length - 1; i >= 0; i--)
572     {
573       int theChar = symbols[i];
574       int charCount = values[i];
575
576       rtnval[nextArrayPos++] = theChar;
577       final int percentage = (charCount * 100) / divisor;
578       rtnval[nextArrayPos++] = percentage;
579       totalPercentage += percentage;
580     }
581     rtnval[0] = symbols.length;
582     rtnval[1] = totalPercentage;
583     int[] result = new int[rtnval.length + 1];
584     result[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
585     System.arraycopy(rtnval, 0, result, 1, rtnval.length);
586
587     return result;
588   }
589
590
591   /**
592    * Extract a sorted extract of cDNA codon profile data. The returned array
593    * contains
594    * 
595    * <pre>
596    *    [profileType, numberOfValues, totalCount, charValue1, percentage1, charValue2, percentage2, ...]
597    * in descending order of percentage value, where the character values encode codon triplets
598    * </pre>
599    * 
600    * @param hashtable
601    * @return
602    */
603   public static int[] extractCdnaProfile(Hashtable hashtable,
604           boolean ignoreGaps)
605   {
606     // this holds #seqs, #ungapped, and then codon count, indexed by encoded
607     // codon triplet
608     int[] codonCounts = (int[]) hashtable.get(PROFILE);
609     int[] sortedCounts = new int[codonCounts.length - 2];
610     System.arraycopy(codonCounts, 2, sortedCounts, 0,
611             codonCounts.length - 2);
612
613     int[] result = new int[3 + 2 * sortedCounts.length];
614     // first value is just the type of profile data
615     result[0] = AlignmentAnnotation.CDNA_PROFILE;
616
617     char[] codons = new char[sortedCounts.length];
618     for (int i = 0; i < codons.length; i++)
619     {
620       codons[i] = (char) i;
621     }
622     QuickSort.sort(sortedCounts, codons);
623     int totalPercentage = 0;
624     int distinctValuesCount = 0;
625     int j = 3;
626     int divisor = ignoreGaps ? codonCounts[1] : codonCounts[0];
627     for (int i = codons.length - 1; i >= 0; i--)
628     {
629       final int codonCount = sortedCounts[i];
630       if (codonCount == 0)
631       {
632         break; // nothing else of interest here
633       }
634       distinctValuesCount++;
635       result[j++] = codons[i];
636       final int percentage = codonCount * 100 / divisor;
637       result[j++] = percentage;
638       totalPercentage += percentage;
639     }
640     result[2] = totalPercentage;
641
642     /*
643      * Just return the non-zero values
644      */
645     // todo next value is redundant if we limit the array to non-zero counts
646     result[1] = distinctValuesCount;
647     return Arrays.copyOfRange(result, 0, j);
648   }
649
650   /**
651    * Compute a consensus for the cDNA coding for a protein alignment.
652    * 
653    * @param alignment
654    *          the protein alignment (which should hold mappings to cDNA
655    *          sequences)
656    * @param hconsensus
657    *          the consensus data stores to be populated (one per column)
658    */
659   public static void calculateCdna(AlignmentI alignment,
660           Hashtable[] hconsensus)
661   {
662     final char gapCharacter = alignment.getGapCharacter();
663     List<AlignedCodonFrame> mappings = alignment.getCodonFrames();
664     if (mappings == null || mappings.isEmpty())
665     {
666       return;
667     }
668
669     int cols = alignment.getWidth();
670     for (int col = 0; col < cols; col++)
671     {
672       // todo would prefer a Java bean for consensus data
673       Hashtable<String, int[]> columnHash = new Hashtable<>();
674       // #seqs, #ungapped seqs, counts indexed by (codon encoded + 1)
675       int[] codonCounts = new int[66];
676       codonCounts[0] = alignment.getSequences().size();
677       int ungappedCount = 0;
678       for (SequenceI seq : alignment.getSequences())
679       {
680         if (seq.getCharAt(col) == gapCharacter)
681         {
682           continue;
683         }
684         List<char[]> codons = MappingUtils.findCodonsFor(seq, col,
685                 mappings);
686         for (char[] codon : codons)
687         {
688           int codonEncoded = CodingUtils.encodeCodon(codon);
689           if (codonEncoded >= 0)
690           {
691             codonCounts[codonEncoded + 2]++;
692             ungappedCount++;
693           }
694         }
695       }
696       codonCounts[1] = ungappedCount;
697       // todo: sort values here, save counts and codons?
698       columnHash.put(PROFILE, codonCounts);
699       hconsensus[col] = columnHash;
700     }
701   }
702
703   /**
704    * Derive displayable cDNA consensus annotation from computed consensus data.
705    * 
706    * @param consensusAnnotation
707    *          the annotation row to be populated for display
708    * @param consensusData
709    *          the computed consensus data
710    * @param showProfileLogo
711    *          if true show all symbols present at each position, else only the
712    *          modal value
713    * @param nseqs
714    *          the number of sequences in the alignment
715    */
716   public static void completeCdnaConsensus(
717           AlignmentAnnotation consensusAnnotation,
718           Hashtable[] consensusData, boolean showProfileLogo, int nseqs)
719   {
720     if (consensusAnnotation == null
721             || consensusAnnotation.annotations == null
722             || consensusAnnotation.annotations.length < consensusData.length)
723     {
724       // called with a bad alignment annotation row - wait for it to be
725       // initialised properly
726       return;
727     }
728
729     // ensure codon triplet scales with font size
730     consensusAnnotation.scaleColLabel = true;
731     for (int col = 0; col < consensusData.length; col++)
732     {
733       Hashtable hci = consensusData[col];
734       if (hci == null)
735       {
736         // gapped protein column?
737         continue;
738       }
739       // array holds #seqs, #ungapped, then codon counts indexed by codon
740       final int[] codonCounts = (int[]) hci.get(PROFILE);
741       int totalCount = 0;
742
743       /*
744        * First pass - get total count and find the highest
745        */
746       final char[] codons = new char[codonCounts.length - 2];
747       for (int j = 2; j < codonCounts.length; j++)
748       {
749         final int codonCount = codonCounts[j];
750         codons[j - 2] = (char) (j - 2);
751         totalCount += codonCount;
752       }
753
754       /*
755        * Sort array of encoded codons by count ascending - so the modal value
756        * goes to the end; start by copying the count (dropping the first value)
757        */
758       int[] sortedCodonCounts = new int[codonCounts.length - 2];
759       System.arraycopy(codonCounts, 2, sortedCodonCounts, 0,
760               codonCounts.length - 2);
761       QuickSort.sort(sortedCodonCounts, codons);
762
763       int modalCodonEncoded = codons[codons.length - 1];
764       int modalCodonCount = sortedCodonCounts[codons.length - 1];
765       String modalCodon = String
766               .valueOf(CodingUtils.decodeCodon(modalCodonEncoded));
767       if (sortedCodonCounts.length > 1 && sortedCodonCounts[codons.length
768               - 2] == sortedCodonCounts[codons.length - 1])
769       {
770         /*
771          * two or more codons share the modal count
772          */
773         modalCodon = "+";
774       }
775       float pid = sortedCodonCounts[sortedCodonCounts.length - 1] * 100
776               / (float) totalCount;
777
778       /*
779        * todo ? Replace consensus hashtable with sorted arrays of codons and
780        * counts (non-zero only). Include total count in count array [0].
781        */
782
783       /*
784        * Scan sorted array backwards for most frequent values first. Show
785        * repeated values compactly.
786        */
787       StringBuilder mouseOver = new StringBuilder(32);
788       StringBuilder samePercent = new StringBuilder();
789       String percent = null;
790       String lastPercent = null;
791       int percentDecPl = getPercentageDp(nseqs);
792
793       for (int j = codons.length - 1; j >= 0; j--)
794       {
795         int codonCount = sortedCodonCounts[j];
796         if (codonCount == 0)
797         {
798           /*
799            * remaining codons are 0% - ignore, but finish off the last one if
800            * necessary
801            */
802           if (samePercent.length() > 0)
803           {
804             mouseOver.append(samePercent).append(": ").append(percent)
805                     .append("% ");
806           }
807           break;
808         }
809         int codonEncoded = codons[j];
810         final int pct = codonCount * 100 / totalCount;
811         String codon = String
812                 .valueOf(CodingUtils.decodeCodon(codonEncoded));
813         StringBuilder sb = new StringBuilder();
814         Format.appendPercentage(sb, pct, percentDecPl);
815         percent = sb.toString();
816         if (showProfileLogo || codonCount == modalCodonCount)
817         {
818           if (percent.equals(lastPercent) && j > 0)
819           {
820             samePercent.append(samePercent.length() == 0 ? "" : ", ");
821             samePercent.append(codon);
822           }
823           else
824           {
825             if (samePercent.length() > 0)
826             {
827               mouseOver.append(samePercent).append(": ").append(lastPercent)
828                       .append("% ");
829             }
830             samePercent.setLength(0);
831             samePercent.append(codon);
832           }
833           lastPercent = percent;
834         }
835       }
836
837       consensusAnnotation.annotations[col] = new Annotation(modalCodon,
838               mouseOver.toString(), ' ', pid);
839     }
840   }
841
842   /**
843    * Returns the number of decimal places to show for profile percentages. For
844    * less than 100 sequences, returns zero (the integer percentage value will be
845    * displayed). For 100-999 sequences, returns 1, for 1000-9999 returns 2, etc.
846    * 
847    * @param nseq
848    * @return
849    */
850   protected static int getPercentageDp(long nseq)
851   {
852     int scale = 0;
853     while (nseq >= 100)
854     {
855       scale++;
856       nseq /= 10;
857     }
858     return scale;
859   }
860
861   /**
862    * Returns the information content at a specified column.
863    * 
864    * @param column
865    *          Index of the column, starting from 0.
866    * @return
867    */
868   public static float getInformationContent(int column,
869           HiddenMarkovModel hmm)
870   {
871     float informationContent = 0f;
872
873     for (char symbol : hmm.getSymbols())
874     {
875       float freq = 0f;
876       freq = ResidueProperties.backgroundFrequencies
877               .get(hmm.getAlphabetType()).get(symbol);
878       Double hmmProb = hmm.getMatchEmissionProbability(column, symbol);
879       float prob = hmmProb.floatValue();
880       informationContent += prob * (Math.log(prob / freq) / Math.log(2));
881
882     }
883
884     return informationContent;
885   }
886
887   /**
888    * Produces a HMM profile for a column in an alignment
889    * 
890    * @param aa
891    *          Alignment annotation for which the profile is being calculated.
892    * @param column
893    *          Column in the alignment the profile is being made for.
894    * @param removeBelowBackground
895    *          Boolean indicating whether to ignore residues with probabilities
896    *          less than their background frequencies.
897    * @return
898    */
899   public static int[] extractHMMProfile(HiddenMarkovModel hmm, int column,
900           boolean removeBelowBackground)
901   {
902
903     if (hmm != null)
904     {
905       int size = hmm.getNumberOfSymbols();
906       char symbols[] = new char[size];
907       int values[] = new int[size];
908       List<Character> charList = hmm.getSymbols();
909       Integer totalCount = 0;
910
911       for (int i = 0; i < size; i++)
912       {
913         char symbol = charList.get(i);
914         symbols[i] = symbol;
915         int value = getAnalogueCount(hmm, column, symbol,
916                 removeBelowBackground);
917         values[i] = value;
918         totalCount += value;
919       }
920
921       QuickSort.sort(values, symbols);
922
923       int[] profile = new int[3 + size * 2];
924
925       profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
926       profile[1] = size;
927       profile[2] = 100;
928
929       if (totalCount != 0)
930       {
931         int arrayPos = 3;
932         for (int k = size - 1; k >= 0; k--)
933         {
934           Float percentage;
935           Integer value = values[k];
936           if (removeBelowBackground)
937           {
938             percentage = (value.floatValue() / totalCount.floatValue())
939                     * 100;
940           }
941           else
942           {
943             percentage = value.floatValue() / 100f;
944           }
945           int intPercent = Math.round(percentage);
946           profile[arrayPos] = symbols[k];
947           profile[arrayPos + 1] = intPercent;
948           arrayPos += 2;
949         }
950       }
951       return profile;
952     }
953     return null;
954   }
955
956   /**
957    * Converts the emission probability of a residue at a column in the alignment
958    * to a 'count' to allow for processing by the annotation renderer.
959    * 
960    * @param hmm
961    * @param column
962    * @param removeBelowBackground
963    *          When true, this method returns 0 for any symbols with a match
964    *          emission probability less than the background frequency.
965    * @param symbol
966    * @return
967    */
968   static int getAnalogueCount(HiddenMarkovModel hmm, int column,
969           char symbol, boolean removeBelowBackground)
970   {
971     Double value;
972
973     value = hmm.getMatchEmissionProbability(column, symbol);
974     double freq;
975
976     freq = ResidueProperties.backgroundFrequencies
977             .get(hmm.getAlphabetType()).get(symbol);
978     if (value < freq && removeBelowBackground)
979     {
980       return 0;
981     }
982
983     value = value * 10000;
984     return Math.round(value.floatValue());
985   }
986 }