X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAAFrequency.java;h=892b1b1ff8cb4efc79cdb7f3bac0ba3a6542bb7e;hb=d043ce47fc710d3eb2629ba926a8a7417bd67d8c;hp=046c589583251008e50dde78ed6ea21d552b725d;hpb=93b7c00c600a36f81c9a5f251366ed09ceae7597;p=jalview.git diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java index 046c589..892b1b1 100755 --- a/src/jalview/analysis/AAFrequency.java +++ b/src/jalview/analysis/AAFrequency.java @@ -50,31 +50,12 @@ import java.util.List; * This class is used extensively in calculating alignment colourschemes that * depend on the amount of conservation in each alignment column. * - * @author $author$ - * @version $Revision$ */ public class AAFrequency { - public static final String PROFILE = "P"; - - private static final String AMINO = "amino"; - - private static final String DNA = "DNA"; - - private static final String RNA = "RNA"; + private static final double LOG2 = Math.log(2); - /* - * Quick look-up of String value of char 'A' to 'Z' - */ - private static final String[] CHARS = new String['Z' - 'A' + 1]; - - static - { - for (char c = 'A'; c <= 'Z'; c++) - { - CHARS[c - 'A'] = String.valueOf(c); - } - } + public static final String PROFILE = "P"; public static final ProfilesI calculate(List list, int start, int end) @@ -109,8 +90,6 @@ public class AAFrequency } } - - /** * Calculate the consensus symbol(s) for each column in the given range. * @@ -157,14 +136,13 @@ public class AAFrequency { if (sequences[row] == null) { - System.err - .println("WARNING: Consensus skipping null sequence - possible race condition."); + System.err.println( + "WARNING: Consensus skipping null sequence - possible race condition."); continue; } - char[] seq = sequences[row].getSequence(); - if (seq.length > column) + if (sequences[row].getLength() > column) { - char c = seq[column]; + char c = sequences[row].getCharAt(column); residueCounts.add(c); if (Comparison.isNucleotide(c)) { @@ -202,20 +180,38 @@ public class AAFrequency // System.out.println(elapsed); } - public static ProfilesI calculateInformation(final HiddenMarkovModel hmm, - int width, int start, int end, boolean saveFullProfile, - boolean removeBelowBackground) + /** + * Returns the full set of profiles for a hidden Markov model. The underlying + * data is the raw probabilities of a residue being emitted at each node, + * however the profiles returned by this function contain the percentage + * chance of a residue emission. + * + * @param hmm + * @param width + * The width of the Profile array (Profiles) to be returned. + * @param start + * The alignment column on which the first profile is based. + * @param end + * The alignment column on which the last profile is based. + * @param removeBelowBackground + * if true, symbols with a match emission probability less than + * background frequency are ignored + * @return + */ + public static ProfilesI calculateHMMProfiles(final HiddenMarkovModel hmm, + int width, int start, int end, boolean removeBelowBackground, + boolean infoLetterHeight) { ProfileI[] result = new ProfileI[width]; - int symbolCount = hmm.getNumberOfSymbols(); - String alph = hmm.getAlphabetType(); + char[] symbols = hmm.getSymbols().toCharArray(); + int symbolCount = symbols.length; for (int column = start; column < end; column++) { ResidueCount counts = new ResidueCount(); - for (char symbol : hmm.getSymbols()) + for (char symbol : symbols) { - int value = getAnalogueCount(hmm, column, removeBelowBackground, - alph, symbol); + int value = getAnalogueCount(hmm, column, symbol, + removeBelowBackground, infoLetterHeight); counts.put(symbol, value); } int maxCount = counts.getModalCount(); @@ -223,11 +219,7 @@ public class AAFrequency int gapCount = counts.getGapCount(); ProfileI profile = new Profile(symbolCount, gapCount, maxCount, maxResidue); - - if (saveFullProfile) - { - profile.setCounts(counts); - } + profile.setCounts(counts); result[column] = profile; } @@ -346,71 +338,66 @@ public class AAFrequency * @param endCol * end column (exclusive) * @param ignoreGaps - * if true, normalise residue percentages + * if true, normalise residue percentages * @param showSequenceLogo * if true include all information symbols, else just show modal * residue - * @param nseq - * number of sequences */ - public static void completeInformation(AlignmentAnnotation information, - ProfilesI profiles, int startCol, int endCol, - boolean ignoreBelowBackground, - boolean showSequenceLogo, long nseq) + public static float completeInformation(AlignmentAnnotation information, + ProfilesI profiles, int startCol, int endCol) { // long now = System.currentTimeMillis(); - if (information == null || information.annotations == null - || information.annotations.length < endCol) + if (information == null || information.annotations == null) { /* * called with a bad alignment annotation row * wait for it to be initialised properly */ - return; + return 0; } - Float max = 0f; + float max = 0f; + SequenceI hmmSeq = information.sequenceRef; - for (int i = startCol; i < endCol; i++) + int seqLength = hmmSeq.getLength(); + if (information.annotations.length < seqLength) { - ProfileI profile = profiles.get(i); - if (profile == null) + return 0; + } + + HiddenMarkovModel hmm = hmmSeq.getHMM(); + + for (int column = startCol; column < endCol; column++) + { + if (column >= seqLength) { - /* - * happens if sequences calculated over were - * shorter than alignment width - */ - information.annotations[i] = null; - return; + // hmm consensus sequence is shorter than the alignment + break; } - - HiddenMarkovModel hmm; - - SequenceI hmmSeq = information.sequenceRef; - hmm = hmmSeq.getHMM(); - - Float value = getInformationContent(i, hmm); - - if (value > max) + float value = hmm.getInformationContent(column); + boolean isNaN = Float.isNaN(value); + if (!isNaN) { - max = value; + max = Math.max(max, value); } - String description = value + " bits"; - - information.annotations[i] = new Annotation(" ", description, - ' ', value); + String description = isNaN ? null + : String.format("%.4f bits", value); + information.annotations[column] = new Annotation( + Character.toString( + Character.toUpperCase(hmmSeq.getCharAt(column))), + description, ' ', value); } + information.graphMax = max; - // long elapsed = System.currentTimeMillis() - now; - // System.out.println(-elapsed); + return max; } /** - * Derive the gap count annotation row. + * Derive the occupancy count annotation * - * @param gaprow + * @param occupancy * the annotation row to add annotations to * @param profiles * the source consensus data @@ -419,11 +406,11 @@ public class AAFrequency * @param endCol * end column (exclusive) */ - public static void completeGapAnnot(AlignmentAnnotation gaprow, + public static void completeGapAnnot(AlignmentAnnotation occupancy, ProfilesI profiles, int startCol, int endCol, long nseq) { - if (gaprow == null || gaprow.annotations == null - || gaprow.annotations.length < endCol) + if (occupancy == null || occupancy.annotations == null + || occupancy.annotations.length < endCol) { /* * called with a bad alignment annotation row @@ -432,9 +419,9 @@ public class AAFrequency return; } // always set ranges again - gaprow.graphMax = nseq; - gaprow.graphMin = 0; - double scale = 0.8/nseq; + occupancy.graphMax = nseq; + occupancy.graphMin = 0; + double scale = 0.8 / nseq; for (int i = startCol; i < endCol; i++) { ProfileI profile = profiles.get(i); @@ -444,7 +431,7 @@ public class AAFrequency * happens if sequences calculated over were * shorter than alignment width */ - gaprow.annotations[i] = null; + occupancy.annotations[i] = null; return; } @@ -452,9 +439,10 @@ public class AAFrequency String description = "" + gapped; - gaprow.annotations[i] = new Annotation("", description, - '\0', gapped, jalview.util.ColorUtils.bleachColour( - Color.DARK_GRAY, (float) scale * gapped)); + occupancy.annotations[i] = new Annotation("", description, '\0', + gapped, + jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY, + (float) scale * gapped)); } } @@ -463,7 +451,8 @@ public class AAFrequency * * Percentages are as a fraction of all sequence, or only ungapped sequences * if ignoreGaps is true. @@ -484,8 +473,8 @@ public class AAFrequency String description = null; if (counts != null && showSequenceLogo) { - int normaliseBy = ignoreGaps ? profile.getNonGapped() : profile - .getHeight(); + int normaliseBy = ignoreGaps ? profile.getNonGapped() + : profile.getHeight(); description = counts.getTooltip(normaliseBy, dp); } else @@ -516,7 +505,7 @@ public class AAFrequency * contains * * 
-   *    [profileType, numberOfValues, nonGapCount, charValue1, percentage1, charValue2, percentage2, ...]
+   *    [profileType, numberOfValues, totalPercent, charValue1, percentage1, charValue2, percentage2, ...]
    * in descending order of percentage value
    *
* @@ -529,7 +518,6 @@ public class AAFrequency */ public static int[] extractProfile(ProfileI profile, boolean ignoreGaps) { - int[] rtnval = new int[64]; ResidueCount counts = profile.getCounts(); if (counts == null) { @@ -540,29 +528,51 @@ public class AAFrequency char[] symbols = symbolCounts.symbols; int[] values = symbolCounts.values; QuickSort.sort(values, symbols); - int nextArrayPos = 2; int totalPercentage = 0; - final int divisor = ignoreGaps ? profile.getNonGapped() : profile - .getHeight(); + final int divisor = ignoreGaps ? profile.getNonGapped() + : profile.getHeight(); /* * traverse the arrays in reverse order (highest counts first) */ + int[] result = new int[3 + 2 * symbols.length]; + int nextArrayPos = 3; + int nonZeroCount = 0; + for (int i = symbols.length - 1; i >= 0; i--) { int theChar = symbols[i]; int charCount = values[i]; - - rtnval[nextArrayPos++] = theChar; final int percentage = (charCount * 100) / divisor; - rtnval[nextArrayPos++] = percentage; + if (percentage == 0) + { + /* + * this count (and any remaining) round down to 0% - discard + */ + break; + } + nonZeroCount++; + result[nextArrayPos++] = theChar; + result[nextArrayPos++] = percentage; totalPercentage += percentage; } - rtnval[0] = symbols.length; - rtnval[1] = totalPercentage; - int[] result = new int[rtnval.length + 1]; + + /* + * truncate array if any zero values were discarded + */ + if (nonZeroCount < symbols.length) + { + int[] tmp = new int[3 + 2 * nonZeroCount]; + System.arraycopy(result, 0, tmp, 0, tmp.length); + result = tmp; + } + + /* + * fill in 'header' values + */ result[0] = AlignmentAnnotation.SEQUENCE_PROFILE; - System.arraycopy(rtnval, 0, result, 1, rtnval.length); + result[1] = nonZeroCount; + result[2] = totalPercentage; return result; } @@ -573,15 +583,15 @@ public class AAFrequency * contains * * 
-   *    [profileType, numberOfValues, totalCount, charValue1, percentage1, charValue2, percentage2, ...]
+   *    [profileType, numberOfValues, totalPercentage, charValue1, percentage1, charValue2, percentage2, ...]
    * in descending order of percentage value, where the character values encode codon triplets
    *
* * @param hashtable * @return */ - public static int[] extractCdnaProfile(Hashtable hashtable, - boolean ignoreGaps) + public static int[] extractCdnaProfile( + Hashtable hashtable, boolean ignoreGaps) { // this holds #seqs, #ungapped, and then codon count, indexed by encoded // codon triplet @@ -611,9 +621,16 @@ public class AAFrequency { break; // nothing else of interest here } + final int percentage = codonCount * 100 / divisor; + if (percentage == 0) + { + /* + * this (and any remaining) values rounded down to 0 - discard + */ + break; + } distinctValuesCount++; result[j++] = codons[i]; - final int percentage = codonCount * 100 / divisor; result[j++] = percentage; totalPercentage += percentage; } @@ -637,7 +654,7 @@ public class AAFrequency * the consensus data stores to be populated (one per column) */ public static void calculateCdna(AlignmentI alignment, - Hashtable[] hconsensus) + Hashtable[] hconsensus) { final char gapCharacter = alignment.getGapCharacter(); List mappings = alignment.getCodonFrames(); @@ -650,7 +667,7 @@ public class AAFrequency for (int col = 0; col < cols; col++) { // todo would prefer a Java bean for consensus data - Hashtable columnHash = new Hashtable<>(); + Hashtable columnHash = new Hashtable<>(); // #seqs, #ungapped seqs, counts indexed by (codon encoded + 1) int[] codonCounts = new int[66]; codonCounts[0] = alignment.getSequences().size(); @@ -661,8 +678,8 @@ public class AAFrequency { continue; } - List codons = MappingUtils - .findCodonsFor(seq, col, mappings); + List codons = MappingUtils.findCodonsFor(seq, col, + mappings); for (char[] codon : codons) { int codonEncoded = CodingUtils.encodeCodon(codon); @@ -670,6 +687,7 @@ public class AAFrequency { codonCounts[codonEncoded + 2]++; ungappedCount++; + break; } } } @@ -695,7 +713,8 @@ public class AAFrequency */ public static void completeCdnaConsensus( AlignmentAnnotation consensusAnnotation, - Hashtable[] consensusData, boolean showProfileLogo, int nseqs) + Hashtable[] consensusData, + boolean showProfileLogo, int nseqs) { if (consensusAnnotation == null || consensusAnnotation.annotations == null @@ -710,7 +729,7 @@ public class AAFrequency consensusAnnotation.scaleColLabel = true; for (int col = 0; col < consensusData.length; col++) { - Hashtable hci = consensusData[col]; + Hashtable hci = consensusData[col]; if (hci == null) { // gapped protein column? @@ -742,10 +761,10 @@ public class AAFrequency int modalCodonEncoded = codons[codons.length - 1]; int modalCodonCount = sortedCodonCounts[codons.length - 1]; - String modalCodon = String.valueOf(CodingUtils - .decodeCodon(modalCodonEncoded)); - if (sortedCodonCounts.length > 1 - && sortedCodonCounts[codons.length - 2] == sortedCodonCounts[codons.length - 1]) + String modalCodon = String + .valueOf(CodingUtils.decodeCodon(modalCodonEncoded)); + if (sortedCodonCounts.length > 1 && sortedCodonCounts[codons.length + - 2] == sortedCodonCounts[codons.length - 1]) { /* * two or more codons share the modal count @@ -804,8 +823,8 @@ public class AAFrequency { if (samePercent.length() > 0) { - mouseOver.append(samePercent).append(": ") - .append(lastPercent).append("% "); + mouseOver.append(samePercent).append(": ").append(lastPercent) + .append("% "); } samePercent.setLength(0); samePercent.append(codon); @@ -839,116 +858,116 @@ public class AAFrequency } /** - * Returns the information content at a specified column. + * Returns the sorted HMM profile for the given column of the alignment. The + * returned array contains * - * @param column - * Index of the column, starting from 0. - * @return - */ - public static float getInformationContent(int column, - HiddenMarkovModel hmm) - { - float informationContent = 0f; - - for (char symbol : hmm.getSymbols()) - { - float freq = 0f; - freq = ResidueProperties.backgroundFrequencies - .get(hmm.getAlphabetType()).get(symbol); - Double hmmProb = hmm.getMatchEmissionProbability(column, symbol); - float prob = hmmProb.floatValue(); - informationContent += prob * (Math.log(prob / freq) / Math.log(2)); - - } - - return informationContent; - } - - /** - * Produces a HMM profile for a column in an alignment + * 
+   *    [profileType=0, numberOfValues, 100, charValue1, percentage1, charValue2, percentage2, ...]
+   * in descending order of percentage value
+   *
* - * @param aa - * Alignment annotation for which the profile is being calculated. + * @param hmm * @param column - * Column in the alignment the profile is being made for. * @param removeBelowBackground - * Boolean indicating whether to ignore residues with probabilities - * less than their background frequencies. + * if true, ignores residues with probability less than their + * background frequency + * @param infoHeight + * if true, uses the log ratio 'information' measure to scale the + * value * @return */ public static int[] extractHMMProfile(HiddenMarkovModel hmm, int column, - boolean removeBelowBackground) + boolean removeBelowBackground, boolean infoHeight) { + if (hmm == null) + { + return null; + } + String alphabet = hmm.getSymbols(); + int size = alphabet.length(); + char symbols[] = new char[size]; + int values[] = new int[size]; + int totalCount = 0; - if (hmm != null) + for (int i = 0; i < size; i++) { - String alph = hmm.getAlphabetType(); - int size = hmm.getNumberOfSymbols(); - char symbols[] = new char[size]; - int values[] = new int[size]; - List charList = hmm.getSymbols(); - Integer totalCount = 0; - - for (int i = 0; i < size; i++) - { - char symbol = charList.get(i); - symbols[i] = symbol; - int value = getAnalogueCount(hmm, column, removeBelowBackground, - alph, symbol); - values[i] = value; - totalCount += value; - } + char symbol = alphabet.charAt(i); + symbols[i] = symbol; + int value = getAnalogueCount(hmm, column, symbol, + removeBelowBackground, infoHeight); + values[i] = value; + totalCount += value; + } - QuickSort.sort(values, symbols); + /* + * sort symbols by increasing emission probability + */ + QuickSort.sort(values, symbols); - int[] profile = new int[3 + size * 2]; + int[] profile = new int[3 + size * 2]; - profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE; - profile[1] = size; - profile[2] = 100; + profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE; + profile[1] = size; + profile[2] = 100; - if (totalCount != 0) + /* + * order symbol/count profile by decreasing emission probability + */ + if (totalCount != 0) + { + int arrayPos = 3; + for (int k = size - 1; k >= 0; k--) { - int arrayPos = 3; - for (int k = size - 1; k >= 0; k--) + Float percentage; + int value = values[k]; + if (removeBelowBackground) { - Float percentage; - Integer value = values[k]; - if (removeBelowBackground) - { - percentage = (value.floatValue() / totalCount.floatValue()) - * 100; - } - else - { - percentage = value.floatValue() / 100f; - } - int intPercent = Math.round(percentage); - profile[arrayPos] = symbols[k]; - profile[arrayPos + 1] = intPercent; - arrayPos += 2; + percentage = ((float) value) / totalCount * 100f; } + else + { + percentage = value / 100f; + } + int intPercent = Math.round(percentage); + profile[arrayPos] = symbols[k]; + profile[arrayPos + 1] = intPercent; + arrayPos += 2; } - return profile; } - return null; + return profile; } - private static int getAnalogueCount(HiddenMarkovModel hmm, int column, - boolean removeBelowBackground, String alph, char symbol) + /** + * Converts the emission probability of a residue at a column in the alignment + * to a 'count', suitable for rendering as an annotation value + * + * @param hmm + * @param column + * @param symbol + * @param removeBelowBackground + * if true, returns 0 for any symbol with a match emission + * probability less than the background frequency + * @infoHeight if true, uses the log ratio 'information content' to scale the + * value + * @return + */ + static int getAnalogueCount(HiddenMarkovModel hmm, int column, + char symbol, boolean removeBelowBackground, boolean infoHeight) { - Double value; - - value = hmm.getMatchEmissionProbability(column, symbol); - double freq; - - freq = ResidueProperties.backgroundFrequencies.get(alph).get(symbol); + double value = hmm.getMatchEmissionProbability(column, symbol); + double freq = ResidueProperties.backgroundFrequencies + .get(hmm.getAlphabetType()).get(symbol); if (value < freq && removeBelowBackground) { return 0; } - value = value * 10000; - return Math.round(value.floatValue()); + if (infoHeight) + { + value = value * (Math.log(value / freq) / LOG2); + } + + value = value * 10000d; + return Math.round((float) value); } }