X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAAFrequency.java;h=5ea892bd5e465f99e6f8b8f14095866fae2fdc0a;hb=3a7e301b545949dc3e5906c0e4343b4f6db56b0c;hp=e682e2ca2091f5e82a41f70737953503c042cbba;hpb=41c34ccf7054239554932deca6f64042231977b5;p=jalview.git diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java index e682e2c..5ea892b 100755 --- a/src/jalview/analysis/AAFrequency.java +++ b/src/jalview/analysis/AAFrequency.java @@ -24,6 +24,7 @@ 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; @@ -32,11 +33,13 @@ import jalview.datamodel.ResidueCount; 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; @@ -54,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' */ @@ -67,8 +76,8 @@ public class AAFrequency } } - public static final ProfilesI calculate(List list, - int start, int end) + public static final ProfilesI calculate(List list, int start, + int end) { return calculate(list, start, end, false); } @@ -100,6 +109,8 @@ public class AAFrequency } } + + /** * Calculate the consensus symbol(s) for each column in the given range. * @@ -191,6 +202,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) + { + ProfileI[] result = new ProfileI[width]; + int symbolCount = hmm.getNumberOfSymbols(); + String alph = hmm.getAlphabetType(); + for (int column = start; column < end; column++) + { + ResidueCount counts = new ResidueCount(); + for (char symbol : hmm.getSymbols()) + { + int value = getAnalogueCount(hmm, column, removeBelowBackground, + alph, symbol); + 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 @@ -289,9 +332,85 @@ public class AAFrequency } /** + * 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 void completeInformation(AlignmentAnnotation information, + ProfilesI profiles, int startCol, int endCol, + boolean ignoreBelowBackground, + boolean showSequenceLogo, long nseq) + { + // 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; + } + + 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; + } + + 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(" ", description, + ' ', value); + } + information.graphMax = max; + // long elapsed = System.currentTimeMillis() - now; + // System.out.println(-elapsed); + } + + /** * Derive the gap count annotation row. * - * @param consensus + * @param gaprow * the annotation row to add annotations to * @param profiles * the source consensus data @@ -300,12 +419,11 @@ public class AAFrequency * @param endCol * end column (exclusive) */ - public static void completeGapAnnot(AlignmentAnnotation consensus, + public static void completeGapAnnot(AlignmentAnnotation gaprow, ProfilesI profiles, int startCol, int endCol, long nseq) { - // long now = System.currentTimeMillis(); - if (consensus == null || consensus.annotations == null - || consensus.annotations.length < endCol) + if (gaprow == null || gaprow.annotations == null + || gaprow.annotations.length < endCol) { /* * called with a bad alignment annotation row @@ -314,8 +432,9 @@ public class AAFrequency return; } // always set ranges again - consensus.graphMax = nseq; - consensus.graphMin = 0; + gaprow.graphMax = nseq; + gaprow.graphMin = 0; + double scale = 0.8/nseq; for (int i = startCol; i < endCol; i++) { ProfileI profile = profiles.get(i); @@ -325,18 +444,18 @@ public class AAFrequency * happens if sequences calculated over were * shorter than alignment width */ - consensus.annotations[i] = null; + gaprow.annotations[i] = null; return; } - final int gapped = profile.getGapped(); + final int gapped = profile.getNonGapped(); String description = "" + gapped; - consensus.annotations[i] = new Annotation(gapped); + gaprow.annotations[i] = new Annotation("", description, + '\0', gapped, jalview.util.ColorUtils.bleachColour( + Color.DARK_GRAY, (float) scale * gapped)); } - // long elapsed = System.currentTimeMillis() - now; - // System.out.println(-elapsed); } /** @@ -408,8 +527,7 @@ public class AAFrequency * calculations * @return */ - public static int[] extractProfile(ProfileI profile, - boolean ignoreGaps) + public static int[] extractProfile(ProfileI profile, boolean ignoreGaps) { int[] rtnval = new int[64]; ResidueCount counts = profile.getCounts(); @@ -449,6 +567,7 @@ public class AAFrequency return result; } + /** * Extract a sorted extract of cDNA codon profile data. The returned array * contains @@ -531,7 +650,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(); @@ -718,4 +837,146 @@ 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; + if ("amino".equals(hmm.getAlphabetType())) + { + freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol); + } + if ("DNA".equals(hmm.getAlphabetType())) + { + freq = ResidueProperties.dnaBackgroundFrequencies.get(symbol); + } + if ("RNA".equals(hmm.getAlphabetType())) + { + freq = ResidueProperties.rnaBackgroundFrequencies.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) + { + + if (hmm != null) + { + 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; + } + + 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; + } + + private static int getAnalogueCount(HiddenMarkovModel hmm, int column, + boolean removeBelowBackground, String alph, char symbol) + { + Double value; + + value = hmm.getMatchEmissionProbability(column, symbol); + double freq; + + if (AMINO.equals(alph) && removeBelowBackground) + { + freq = ResidueProperties.aminoBackgroundFrequencies.get(symbol); + if (value < freq) + { + value = 0d; + } + } + else if (DNA.equals(alph) && removeBelowBackground) + { + freq = ResidueProperties.dnaBackgroundFrequencies.get(symbol); + if (value < freq) + { + value = 0d; + } + } + else if (RNA.equals(alph) && removeBelowBackground) + { + freq = ResidueProperties.rnaBackgroundFrequencies.get(symbol); + if (value < freq) + { + value = 0d; + } + } + value = value * 10000; + return Math.round(value.floatValue()); + } }