X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAAFrequency.java;h=85a9bb03904dedfef151ba5f743448c46b446419;hb=bd6ce8f5f9fc8e5bc8a6188d15987ce0ffd2c1ee;hp=30d53737cda0309f04a4ed591563656f30177356;hpb=0a5ce6145bb76fc7eb8a5cc2670e20453fbedd29;p=jalview.git diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java index 30d5373..85a9bb0 100755 --- a/src/jalview/analysis/AAFrequency.java +++ b/src/jalview/analysis/AAFrequency.java @@ -24,16 +24,22 @@ import jalview.datamodel.AlignedCodonFrame; import jalview.datamodel.AlignmentAnnotation; import jalview.datamodel.AlignmentI; import jalview.datamodel.Annotation; +import jalview.datamodel.HiddenMarkovModel; import jalview.datamodel.Profile; +import jalview.datamodel.ProfileI; +import jalview.datamodel.Profiles; +import jalview.datamodel.ProfilesI; import jalview.datamodel.ResidueCount; -import jalview.datamodel.SequenceI; import jalview.datamodel.ResidueCount.SymbolCounts; +import jalview.datamodel.SequenceI; import jalview.ext.android.SparseIntArray; +import jalview.schemes.ResidueProperties; import jalview.util.Comparison; import jalview.util.Format; import jalview.util.MappingUtils; import jalview.util.QuickSort; +import java.awt.Color; import java.util.Arrays; import java.util.Hashtable; import java.util.List; @@ -51,6 +57,12 @@ 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"; + /* * Quick look-up of String value of char 'A' to 'Z' */ @@ -64,13 +76,13 @@ public class AAFrequency } } - public static final Profile[] calculate(List list, - int start, int end) + public static final ProfilesI calculate(List list, int start, + int end) { return calculate(list, start, end, false); } - public static final Profile[] calculate(List sequences, + public static final ProfilesI calculate(List sequences, int start, int end, boolean profile) { SequenceI[] seqs = new SequenceI[sequences.size()]; @@ -80,39 +92,40 @@ public class AAFrequency for (int i = 0; i < sequences.size(); i++) { seqs[i] = sequences.get(i); - if (seqs[i].getLength() > width) + int length = seqs[i].getLength(); + if (length > width) { - width = seqs[i].getLength(); + width = length; } } - Profile[] reply = new Profile[width]; - if (end >= width) { end = width; } - calculate(seqs, start, end, reply, profile); + ProfilesI reply = calculate(seqs, width, start, end, profile); return reply; } } + + /** * Calculate the consensus symbol(s) for each column in the given range. * * @param sequences + * @param width + * the full width of the alignment * @param start * start column (inclusive, base zero) * @param end * end column (exclusive) - * @param result - * array in which to store profile per column * @param saveFullProfile * if true, store all symbol counts */ - public static final void calculate(final SequenceI[] sequences, - int start, int end, Profile[] result, boolean saveFullProfile) + public static final ProfilesI calculate(final SequenceI[] sequences, + int width, int start, int end, boolean saveFullProfile) { // long now = System.currentTimeMillis(); int seqCount = sequences.length; @@ -120,6 +133,8 @@ public class AAFrequency int nucleotideCount = 0; int peptideCount = 0; + ProfileI[] result = new ProfileI[width]; + for (int column = start; column < end; column++) { /* @@ -142,14 +157,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)) { @@ -172,7 +186,7 @@ public class AAFrequency int maxCount = residueCounts.getModalCount(); String maxResidue = residueCounts.getResiduesForCount(maxCount); int gapCount = residueCounts.getGapCount(); - Profile profile = new Profile(seqCount, gapCount, maxCount, + ProfileI profile = new Profile(seqCount, gapCount, maxCount, maxResidue); if (saveFullProfile) @@ -182,11 +196,63 @@ public class AAFrequency result[column] = profile; } + return new Profiles(result); // long elapsed = System.currentTimeMillis() - now; // System.out.println(elapsed); } /** + * 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 saveFullProfile + * Flag for saving the counts for each profile + * @param removeBelowBackground + * Flag for removing any characters with a match emission probability + * less than its background frequency + * @return + */ + public static ProfilesI calculateHMMProfiles(final HiddenMarkovModel hmm, + int width, int start, int end, boolean saveFullProfile, + boolean removeBelowBackground, boolean infoLetterHeight) + { + ProfileI[] result = new ProfileI[width]; + int symbolCount = hmm.getNumberOfSymbols(); + for (int column = start; column < end; column++) + { + ResidueCount counts = new ResidueCount(); + for (char symbol : hmm.getSymbols()) + { + int value = getAnalogueCount(hmm, column, symbol, + removeBelowBackground, infoLetterHeight); + counts.put(symbol, value); + } + int maxCount = counts.getModalCount(); + String maxResidue = counts.getResiduesForCount(maxCount); + int gapCount = counts.getGapCount(); + ProfileI profile = new Profile(symbolCount, gapCount, maxCount, + maxResidue); + + if (saveFullProfile) + { + profile.setCounts(counts); + } + + result[column] = profile; + } + return new Profiles(result); + } + + /** * Make an estimate of the profile size we are going to compute i.e. how many * different characters may be present in it. Overestimating has a cost of * using more memory than necessary. Underestimating has a cost of needing to @@ -220,10 +286,10 @@ public class AAFrequency * the annotation row to add annotations to * @param profiles * the source consensus data - * @param iStart - * start column - * @param width - * end column + * @param startCol + * start column (inclusive) + * @param endCol + * end column (exclusive) * @param ignoreGaps * if true, normalise residue percentages ignoring gaps * @param showSequenceLogo @@ -233,12 +299,12 @@ public class AAFrequency * number of sequences */ public static void completeConsensus(AlignmentAnnotation consensus, - Profile[] profiles, int iStart, int width, boolean ignoreGaps, + ProfilesI profiles, int startCol, int endCol, boolean ignoreGaps, boolean showSequenceLogo, long nseq) { // long now = System.currentTimeMillis(); if (consensus == null || consensus.annotations == null - || consensus.annotations.length < width) + || consensus.annotations.length < endCol) { /* * called with a bad alignment annotation row @@ -247,21 +313,21 @@ public class AAFrequency return; } - final int dp = getPercentageDp(nseq); - - for (int i = iStart; i < width; i++) + for (int i = startCol; i < endCol; i++) { - Profile profile; - if (i >= profiles.length || ((profile = profiles[i]) == null)) + ProfileI profile = profiles.get(i); + if (profile == null) { /* * happens if sequences calculated over were * shorter than alignment width */ consensus.annotations[i] = null; - continue; + return; } + final int dp = getPercentageDp(nseq); + float value = profile.getPercentageIdentity(ignoreGaps); String description = getTooltip(profile, value, showSequenceLogo, @@ -276,19 +342,154 @@ public class AAFrequency { modalResidue = "+"; } - consensus.annotations[i] = new Annotation(modalResidue, - description, ' ', value); + consensus.annotations[i] = new Annotation(modalResidue, description, + ' ', value); } // long elapsed = System.currentTimeMillis() - now; // System.out.println(-elapsed); } /** + * Derive the information annotations to be added to the alignment for + * display. This does not recompute the raw data, but may be called on a + * change in display options, such as 'ignore below background frequency', + * which may in turn result in a change in the derived values. + * + * @param information + * the annotation row to add annotations to + * @param profiles + * the source information data + * @param startCol + * start column (inclusive) + * @param endCol + * end column (exclusive) + * @param ignoreGaps + * 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 float completeInformation(AlignmentAnnotation information, + ProfilesI profiles, int startCol, int endCol, long nseq, + Float currentMax) + { + // long now = System.currentTimeMillis(); + if (information == null || information.annotations == null + || information.annotations.length < endCol) + { + /* + * called with a bad alignment annotation row + * wait for it to be initialised properly + */ + return 0; + } + + Float max = 0f; + + for (int i = startCol; i < endCol; i++) + { + ProfileI profile = profiles.get(i); + if (profile == null) + { + /* + * happens if sequences calculated over were + * shorter than alignment width + */ + information.annotations[i] = null; + return 0; + } + + HiddenMarkovModel hmm; + + SequenceI hmmSeq = information.sequenceRef; + + hmm = hmmSeq.getHMM(); + + Float value = getInformationContent(i, hmm); + + if (value > max) + { + max = value; + } + + String description = value + " bits"; + information.annotations[i] = new Annotation( + Character.toString(Character + .toUpperCase(hmm.getConsensusAtAlignColumn(i))), + description, ' ', value); + } + if (max > currentMax) + { + information.graphMax = max; + return max; + } + else + { + information.graphMax = currentMax; + return currentMax; + } + } + + /** + * Derive the gap count annotation row. + * + * @param gaprow + * the annotation row to add annotations to + * @param profiles + * the source consensus data + * @param startCol + * start column (inclusive) + * @param endCol + * end column (exclusive) + */ + public static void completeGapAnnot(AlignmentAnnotation gaprow, + ProfilesI profiles, int startCol, int endCol, long nseq) + { + if (gaprow == null || gaprow.annotations == null + || gaprow.annotations.length < endCol) + { + /* + * called with a bad alignment annotation row + * wait for it to be initialised properly + */ + return; + } + // always set ranges again + gaprow.graphMax = nseq; + gaprow.graphMin = 0; + double scale = 0.8 / nseq; + for (int i = startCol; i < endCol; i++) + { + ProfileI profile = profiles.get(i); + if (profile == null) + { + /* + * happens if sequences calculated over were + * shorter than alignment width + */ + gaprow.annotations[i] = null; + return; + } + + final int gapped = profile.getNonGapped(); + + String description = "" + gapped; + + gaprow.annotations[i] = new Annotation("", description, '\0', gapped, + jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY, + (float) scale * gapped)); + } + } + + /** * Returns a tooltip showing either *
    *
  • the full profile (percentages of all residues present), if * showSequenceLogo is true, or
    • - *
  • just the modal (most common) residue(s), if showSequenceLogo is false
    • + *
  • just the modal (most common) residue(s), if showSequenceLogo is + * false
    • *
* Percentages are as a fraction of all sequence, or only ungapped sequences * if ignoreGaps is true. @@ -301,7 +502,7 @@ public class AAFrequency * the number of decimal places to format percentages to * @return */ - static String getTooltip(Profile profile, float pid, + static String getTooltip(ProfileI profile, float pid, boolean showSequenceLogo, boolean ignoreGaps, int dp) { ResidueCount counts = profile.getCounts(); @@ -309,8 +510,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 @@ -352,8 +553,7 @@ public class AAFrequency * calculations * @return */ - public static int[] extractProfile(Profile profile, - boolean ignoreGaps) + public static int[] extractProfile(ProfileI profile, boolean ignoreGaps) { int[] rtnval = new int[64]; ResidueCount counts = profile.getCounts(); @@ -368,8 +568,8 @@ public class AAFrequency 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) @@ -393,6 +593,7 @@ public class AAFrequency return result; } + /** * Extract a sorted extract of cDNA codon profile data. The returned array * contains @@ -475,7 +676,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(); @@ -486,8 +687,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); @@ -567,10 +768,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 @@ -629,8 +830,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); @@ -662,4 +863,135 @@ public class AAFrequency } return scale; } + + /** + * Returns the information content at a specified column. + * + * @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 + * + * @param aa + * Alignment annotation for which the profile is being calculated. + * @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. + * @return + */ + public static int[] extractHMMProfile(HiddenMarkovModel hmm, int column, + boolean removeBelowBackground, boolean infoHeight) + { + + if (hmm != null) + { + 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, symbol, + removeBelowBackground, infoHeight); + values[i] = value; + totalCount += value; + } + + QuickSort.sort(values, symbols); + + int[] profile = new int[3 + size * 2]; + + profile[0] = AlignmentAnnotation.SEQUENCE_PROFILE; + profile[1] = size; + profile[2] = 100; + + if (totalCount != 0) + { + int arrayPos = 3; + for (int k = size - 1; k >= 0; k--) + { + 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; + } + } + return profile; + } + return null; + } + + /** + * Converts the emission probability of a residue at a column in the alignment + * to a 'count' to allow for processing by the annotation renderer. + * + * @param hmm + * @param column + * @param removeBelowBackground + * When true, this method returns 0 for any symbols with a match + * emission probability less than the background frequency. + * @param symbol + * @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(hmm.getAlphabetType()).get(symbol); + if (value < freq && removeBelowBackground) + { + return 0; + } + + if (infoHeight) + { + value = value * (Math.log(value / freq) / Math.log(2)); + } + + value = value * 10000; + return Math.round(value.floatValue()); + } }