2 * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3 * Copyright (C) $$Year-Rel$$ The Jalview Authors
5 * This file is part of Jalview.
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.
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.
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.
21 package jalview.analysis;
23 import jalview.datamodel.AlignedCodonFrame;
24 import jalview.datamodel.AlignmentAnnotation;
25 import jalview.datamodel.AlignmentI;
26 import jalview.datamodel.Annotation;
27 import jalview.datamodel.SequenceI;
28 import jalview.util.Format;
29 import jalview.util.MappingUtils;
30 import jalview.util.QuickSort;
32 import java.util.Arrays;
33 import java.util.Hashtable;
34 import java.util.List;
38 * Takes in a vector or array of sequences and column start and column end and
39 * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied.
40 * This class is used extensively in calculating alignment colourschemes that
41 * depend on the amount of conservation in each alignment column.
46 public class AAFrequency
48 private static final int TO_UPPER_CASE = 'A' - 'a'; // -32
50 public static final String MAXCOUNT = "C";
52 public static final String MAXRESIDUE = "R";
54 public static final String PID_GAPS = "G";
56 public static final String PID_NOGAPS = "N";
58 public static final String PROFILE = "P";
60 public static final String ENCODED_CHARS = "E";
63 * Quick look-up of String value of char 'A' to 'Z'
65 private static final String[] CHARS = new String['Z' - 'A' + 1];
69 for (char c = 'A'; c <= 'Z'; c++)
71 CHARS[c - 'A'] = String.valueOf(c);
75 public static final Hashtable[] calculate(List<SequenceI> list,
78 return calculate(list, start, end, false);
81 public static final Hashtable[] calculate(List<SequenceI> sequences,
82 int start, int end, boolean profile)
84 SequenceI[] seqs = new SequenceI[sequences.size()];
86 synchronized (sequences)
88 for (int i = 0; i < sequences.size(); i++)
90 seqs[i] = sequences.get(i);
91 if (seqs[i].getLength() > width)
93 width = seqs[i].getLength();
97 Hashtable[] reply = new Hashtable[width];
104 calculate(seqs, start, end, reply, profile);
109 public static final void calculate(SequenceI[] sequences, int start,
110 int end, Hashtable[] result, boolean profile)
112 Hashtable residueHash;
113 int maxCount, nongap, i, j, v;
114 int jSize = sequences.length;
119 int[] values = new int[255];
123 for (i = start; i < end; i++)
125 residueHash = new Hashtable();
129 values = new int[255];
131 for (j = 0; j < jSize; j++)
133 if (sequences[j] == null)
136 .println("WARNING: Consensus skipping null sequence - possible race condition.");
139 seq = sequences[j].getSequence();
144 if (c == '.' || c == ' ')
154 else if ('a' <= c && c <= 'z')
170 maxResidue = String.valueOf(c);
175 for (v = 'A'; v <= 'Z'; v++)
177 // TODO why ignore values[v] == 1?
178 if (values[v] < 1 /* 2 */|| values[v] < maxCount)
183 if (values[v] > maxCount)
185 maxResidue = CHARS[v - 'A'];
187 else if (values[v] == maxCount)
189 maxResidue += CHARS[v - 'A'];
191 maxCount = values[v];
194 if (maxResidue.length() == 0)
200 // TODO use a 1-dimensional array with jSize, nongap in [0] and [1]
201 residueHash.put(PROFILE, new int[][] { values,
202 new int[] { jSize, nongap } });
204 residueHash.put(MAXCOUNT, new Integer(maxCount));
205 residueHash.put(MAXRESIDUE, maxResidue);
207 percentage = ((float) maxCount * 100) / jSize;
208 residueHash.put(PID_GAPS, new Float(percentage));
212 // calculate for non-gapped too
213 percentage = ((float) maxCount * 100) / nongap;
215 residueHash.put(PID_NOGAPS, new Float(percentage));
217 result[i] = residueHash;
222 * Compute all or part of the annotation row from the given consensus
226 * - pre-allocated annotation row
230 * @param ignoreGapsInConsensusCalculation
231 * @param includeAllConsSymbols
234 public static void completeConsensus(AlignmentAnnotation consensus,
235 Hashtable[] hconsensus, int iStart, int width,
236 boolean ignoreGapsInConsensusCalculation,
237 boolean includeAllConsSymbols, long nseq)
239 completeConsensus(consensus, hconsensus, iStart, width,
240 ignoreGapsInConsensusCalculation, includeAllConsSymbols, null,
245 * Derive the consensus annotations to be added to the alignment for display.
246 * This does not recompute the raw data, but may be called on a change in
247 * display options, such as 'show logo', which may in turn result in a change
248 * in the derived values.
251 * the annotation row to add annotations to
253 * the source consensus data
258 * @param ignoreGapsInConsensusCalculation
259 * if true, use the consensus calculated ignoring gaps
260 * @param includeAllConsSymbols
261 * if true include all consensus symbols, else just show modal
265 * number of sequences
267 public static void completeConsensus(AlignmentAnnotation consensus,
268 Hashtable[] hconsensus, int iStart, int width,
269 boolean ignoreGapsInConsensusCalculation,
270 boolean includeAllConsSymbols, char[] alphabet, long nseq)
272 if (consensus == null || consensus.annotations == null
273 || consensus.annotations.length < width)
275 // called with a bad alignment annotation row - wait for it to be
276 // initialised properly
280 final Format fmt = getPercentageFormat(nseq);
282 for (int i = iStart; i < width; i++)
285 if (i >= hconsensus.length || ((hci = hconsensus[i]) == null))
287 // happens if sequences calculated over were shorter than alignment
289 consensus.annotations[i] = null;
292 Float fv = (Float) hci
293 .get(ignoreGapsInConsensusCalculation ? PID_NOGAPS : PID_GAPS);
296 consensus.annotations[i] = null;
297 // data has changed below us .. give up and
300 float value = fv.floatValue();
301 String maxRes = hci.get(AAFrequency.MAXRESIDUE).toString();
302 StringBuilder mouseOver = new StringBuilder(64);
303 if (maxRes.length() > 1)
305 mouseOver.append("[").append(maxRes).append("] ");
310 mouseOver.append(hci.get(AAFrequency.MAXRESIDUE) + " ");
312 int[][] profile = (int[][]) hci.get(AAFrequency.PROFILE);
313 if (profile != null && includeAllConsSymbols)
315 int sequenceCount = profile[1][0];
316 int nonGappedCount = profile[1][1];
317 int normalisedBy = ignoreGapsInConsensusCalculation ? nonGappedCount
319 mouseOver.setLength(0);
320 if (alphabet != null)
322 for (int c = 0; c < alphabet.length; c++)
324 float tval = profile[0][alphabet[c]] * 100f / normalisedBy;
326 .append(((c == 0) ? "" : "; "))
329 .append(((fmt != null) ? fmt.form(tval) : ((int) tval)))
335 // TODO do this sort once only in calculate()?
336 // char[][] ca = new char[profile[0].length][];
337 char[] ca = new char[profile[0].length];
338 float[] vl = new float[profile[0].length];
339 for (int c = 0; c < ca.length; c++)
342 // ca[c] = new char[]
344 vl[c] = profile[0][c];
346 QuickSort.sort(vl, ca);
347 for (int p = 0, c = ca.length - 1; profile[0][ca[c]] > 0; c--)
349 final char residue = ca[c];
352 float tval = profile[0][residue] * 100f / normalisedBy;
354 .append((((p == 0) ? "" : "; ")))
357 .append(((fmt != null) ? fmt.form(tval)
358 : ((int) tval))).append("%");
367 (((fmt != null) ? fmt.form(value) : ((int) value))))
370 consensus.annotations[i] = new Annotation(maxRes,
371 mouseOver.toString(), ' ', value);
376 * Returns a Format designed to show all significant figures for profile
377 * percentages. For less than 100 sequences, returns null (the integer
378 * percentage value will be displayed). For 100-999 sequences, returns "%3.1f"
383 protected static Format getPercentageFormat(long nseq)
391 return scale <= 1 ? null : new Format("%3." + (scale - 1) + "f");
395 * Returns the sorted profile for the given consensus data. The returned array
399 * [profileType, numberOfValues, nonGapCount, charValue1, percentage1, charValue2, percentage2, ...]
400 * in descending order of percentage value
404 * the data table from which to extract and sort values
406 * if true, only non-gapped values are included in percentage
410 public static int[] extractProfile(Hashtable hconsensus,
413 int[] rtnval = new int[64];
414 int[][] profile = (int[][]) hconsensus.get(AAFrequency.PROFILE);
419 char[] ca = new char[profile[0].length];
420 float[] vl = new float[profile[0].length];
421 for (int c = 0; c < ca.length; c++)
424 vl[c] = profile[0][c];
426 QuickSort.sort(vl, ca);
427 int nextArrayPos = 2;
428 int totalPercentage = 0;
429 int distinctValuesCount = 0;
430 final int divisor = profile[1][ignoreGaps ? 1 : 0];
431 for (int c = ca.length - 1; profile[0][ca[c]] > 0; c--)
435 rtnval[nextArrayPos++] = ca[c];
436 final int percentage = (int) (profile[0][ca[c]] * 100f / divisor);
437 rtnval[nextArrayPos++] = percentage;
438 totalPercentage += percentage;
439 distinctValuesCount++;
442 rtnval[0] = distinctValuesCount;
443 rtnval[1] = totalPercentage;
444 int[] result = new int[rtnval.length + 1];
445 result[0] = AlignmentAnnotation.SEQUENCE_PROFILE;
446 System.arraycopy(rtnval, 0, result, 1, rtnval.length);
452 * Extract a sorted extract of cDNA codon profile data. The returned array
456 * [profileType, numberOfValues, totalCount, charValue1, percentage1, charValue2, percentage2, ...]
457 * in descending order of percentage value, where the character values encode codon triplets
463 public static int[] extractCdnaProfile(Hashtable hashtable,
466 // this holds #seqs, #ungapped, and then codon count, indexed by encoded
468 int[] codonCounts = (int[]) hashtable.get(PROFILE);
469 int[] sortedCounts = new int[codonCounts.length - 2];
470 System.arraycopy(codonCounts, 2, sortedCounts, 0,
471 codonCounts.length - 2);
473 int[] result = new int[3 + 2 * sortedCounts.length];
474 // first value is just the type of profile data
475 result[0] = AlignmentAnnotation.CDNA_PROFILE;
477 char[] codons = new char[sortedCounts.length];
478 for (int i = 0; i < codons.length; i++)
480 codons[i] = (char) i;
482 QuickSort.sort(sortedCounts, codons);
483 int totalPercentage = 0;
484 int distinctValuesCount = 0;
486 int divisor = ignoreGaps ? codonCounts[1] : codonCounts[0];
487 for (int i = codons.length - 1; i >= 0; i--)
489 final int codonCount = sortedCounts[i];
492 break; // nothing else of interest here
494 distinctValuesCount++;
495 result[j++] = codons[i];
496 final int percentage = codonCount * 100 / divisor;
497 result[j++] = percentage;
498 totalPercentage += percentage;
500 result[2] = totalPercentage;
503 * Just return the non-zero values
505 // todo next value is redundant if we limit the array to non-zero counts
506 result[1] = distinctValuesCount;
507 return Arrays.copyOfRange(result, 0, j);
511 * Compute a consensus for the cDNA coding for a protein alignment.
514 * the protein alignment (which should hold mappings to cDNA
517 * the consensus data stores to be populated (one per column)
519 public static void calculateCdna(AlignmentI alignment,
520 Hashtable[] hconsensus)
522 final char gapCharacter = alignment.getGapCharacter();
523 Set<AlignedCodonFrame> mappings = alignment.getCodonFrames();
524 if (mappings == null || mappings.isEmpty())
529 int cols = alignment.getWidth();
530 for (int col = 0; col < cols; col++)
532 // todo would prefer a Java bean for consensus data
533 Hashtable<String, int[]> columnHash = new Hashtable<String, int[]>();
534 // #seqs, #ungapped seqs, counts indexed by (codon encoded + 1)
535 int[] codonCounts = new int[66];
536 codonCounts[0] = alignment.getSequences().size();
537 int ungappedCount = 0;
538 for (SequenceI seq : alignment.getSequences())
540 if (seq.getCharAt(col) == gapCharacter)
544 char[] codon = MappingUtils.findCodonFor(seq, col, mappings);
545 int codonEncoded = CodingUtils.encodeCodon(codon);
546 if (codonEncoded >= 0)
548 codonCounts[codonEncoded + 2]++;
552 codonCounts[1] = ungappedCount;
553 // todo: sort values here, save counts and codons?
554 columnHash.put(PROFILE, codonCounts);
555 hconsensus[col] = columnHash;
560 * Derive displayable cDNA consensus annotation from computed consensus data.
562 * @param consensusAnnotation
563 * the annotation row to be populated for display
564 * @param consensusData
565 * the computed consensus data
566 * @param showProfileLogo
567 * if true show all symbols present at each position, else only the
570 * the number of sequences in the alignment
572 public static void completeCdnaConsensus(
573 AlignmentAnnotation consensusAnnotation,
574 Hashtable[] consensusData, boolean showProfileLogo, int nseqs)
576 if (consensusAnnotation == null
577 || consensusAnnotation.annotations == null
578 || consensusAnnotation.annotations.length < consensusData.length)
580 // called with a bad alignment annotation row - wait for it to be
581 // initialised properly
585 // ensure codon triplet scales with font size
586 consensusAnnotation.scaleColLabel = true;
587 for (int col = 0; col < consensusData.length; col++)
589 Hashtable hci = consensusData[col];
592 // gapped protein column?
595 // array holds #seqs, #ungapped, then codon counts indexed by codon
596 final int[] codonCounts = (int[]) hci.get(PROFILE);
600 * First pass - get total count and find the highest
602 final char[] codons = new char[codonCounts.length - 2];
603 for (int j = 2; j < codonCounts.length; j++)
605 final int codonCount = codonCounts[j];
606 codons[j - 2] = (char) (j - 2);
607 totalCount += codonCount;
611 * Sort array of encoded codons by count ascending - so the modal value
612 * goes to the end; start by copying the count (dropping the first value)
614 int[] sortedCodonCounts = new int[codonCounts.length - 2];
615 System.arraycopy(codonCounts, 2, sortedCodonCounts, 0,
616 codonCounts.length - 2);
617 QuickSort.sort(sortedCodonCounts, codons);
619 int modalCodonEncoded = codons[codons.length - 1];
620 int modalCodonCount = sortedCodonCounts[codons.length - 1];
621 String modalCodon = String.valueOf(CodingUtils
622 .decodeCodon(modalCodonEncoded));
623 if (sortedCodonCounts.length > 1
624 && sortedCodonCounts[codons.length - 2] == modalCodonEncoded)
628 float pid = sortedCodonCounts[sortedCodonCounts.length - 1] * 100
629 / (float) totalCount;
632 * todo ? Replace consensus hashtable with sorted arrays of codons and
633 * counts (non-zero only). Include total count in count array [0].
637 * Scan sorted array backwards for most frequent values first. Show
638 * repeated values compactly.
640 StringBuilder mouseOver = new StringBuilder(32);
641 StringBuilder samePercent = new StringBuilder();
642 String percent = null;
643 String lastPercent = null;
644 Format fmt = getPercentageFormat(nseqs);
646 for (int j = codons.length - 1; j >= 0; j--)
648 int codonCount = sortedCodonCounts[j];
652 * remaining codons are 0% - ignore, but finish off the last one if
655 if (samePercent.length() > 0)
657 mouseOver.append(samePercent).append(": ").append(percent)
662 int codonEncoded = codons[j];
663 final int pct = codonCount * 100 / totalCount;
664 String codon = String
665 .valueOf(CodingUtils.decodeCodon(codonEncoded));
666 percent = fmt == null ? Integer.toString(pct) : fmt.form(pct);
667 if (showProfileLogo || codonCount == modalCodonCount)
669 if (percent.equals(lastPercent) && j > 0)
671 samePercent.append(samePercent.length() == 0 ? "" : ", ");
672 samePercent.append(codon);
676 if (samePercent.length() > 0)
678 mouseOver.append(samePercent).append(": ")
679 .append(lastPercent).append("% ");
681 samePercent.setLength(0);
682 samePercent.append(codon);
684 lastPercent = percent;
688 consensusAnnotation.annotations[col] = new Annotation(modalCodon,
689 mouseOver.toString(), ' ', pid);