X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAAFrequency.java;h=796625a46ad1f84849dcd73f7b416cf736df83f7;hb=b254ce17e47c43d19804efa129367291a79b1315;hp=3636b5e0f79c816ce6691d05d82134286223d3d7;hpb=0af11fd54aa58f2b17e5f940391cb880e61bc57f;p=jalview.git diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java index 3636b5e..796625a 100755 --- a/src/jalview/analysis/AAFrequency.java +++ b/src/jalview/analysis/AAFrequency.java @@ -30,6 +30,7 @@ import jalview.datamodel.Profiles; import jalview.datamodel.ProfilesI; import jalview.datamodel.ResidueCount; import jalview.datamodel.ResidueCount.SymbolCounts; +import jalview.datamodel.SecondaryStructureCount; import jalview.datamodel.SequenceI; import jalview.ext.android.SparseIntArray; import jalview.util.Comparison; @@ -54,6 +55,8 @@ import java.util.List; public class AAFrequency { public static final String PROFILE = "P"; + private static final String SS_ANNOTATION_LABEL = "Secondary Structure"; + private static final char COIL = 'C'; /* * Quick look-up of String value of char 'A' to 'Z' @@ -147,14 +150,13 @@ public class AAFrequency { if (sequences[row] == null) { - System.err - .println("WARNING: Consensus skipping null sequence - possible race condition."); + jalview.bin.Console.errPrintln( + "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)) { @@ -189,7 +191,131 @@ public class AAFrequency } return new Profiles(result); // long elapsed = System.currentTimeMillis() - now; - // System.out.println(elapsed); + // jalview.bin.Console.outPrintln(elapsed); + } + + + public static final ProfilesI calculateSS(List list, int start, + int end) + { + return calculateSS(list, start, end, false); + } + + public static final ProfilesI calculateSS(List sequences, + int start, int end, boolean profile) + { + SequenceI[] seqs = new SequenceI[sequences.size()]; + int width = 0; + synchronized (sequences) + { + for (int i = 0; i < sequences.size(); i++) + { + seqs[i] = sequences.get(i); + int length = seqs[i].getLength(); + if (length > width) + { + width = length; + } + } + + if (end >= width) + { + end = width; + } + + ProfilesI reply = calculateSS(seqs, width, start, end, profile); + return reply; + } + } + + public static final ProfilesI calculateSS(final SequenceI[] sequences, + int width, int start, int end, boolean saveFullProfile) + { + // long now = System.currentTimeMillis(); + int seqCount = sequences.length; + + ProfileI[] result = new ProfileI[width]; + + for (int column = start; column < end; column++) + { + /* + * Apply a heuristic to detect nucleotide data (which can + * be counted in more compact arrays); here we test for + * more than 90% nucleotide; recheck every 10 columns in case + * of misleading data e.g. highly conserved Alanine in peptide! + * Mistakenly guessing nucleotide has a small performance cost, + * as it will result in counting in sparse arrays. + * Mistakenly guessing peptide has a small space cost, + * as it will use a larger than necessary array to hold counts. + */ + + int ssCount = 0; + + SecondaryStructureCount ssCounts = new SecondaryStructureCount(); + + for (int row = 0; row < seqCount; row++) + { + if (sequences[row] == null) + { + jalview.bin.Console.errPrintln( + "WARNING: Consensus skipping null sequence - possible race condition."); + continue; + } + + char c = sequences[row].getCharAt(column); + + if (sequences[row].getLength() > column && !Comparison.isGap(c)) + { + + AlignmentAnnotation[] aa = sequences[row].getAnnotation(SS_ANNOTATION_LABEL); + if(aa == null) { + continue; + } + int seqPosition = sequences[row].findPosition(column); + char ss; + if (aa[0].getAnnotationForPosition(seqPosition) != null) { + ss = aa[0].getAnnotationForPosition(seqPosition).secondaryStructure; + + //There is no representation for coil and it can be either ' ' or null. + if (ss == ' ') { + ss = COIL; + } + } + else { + ss = COIL; + } + + //secondaryStructures[row][column] = ss; + + ssCounts.add(ss); + ssCount++; + + } + else + { + /* + * count a gap if the sequence doesn't reach this column + */ + ssCounts.addGap(); + } + } + + int maxSSCount = ssCounts.getModalCount(); + String maxSS = ssCounts.getSSForCount(maxSSCount); + int gapCount = ssCounts.getGapCount(); + ProfileI profile = new Profile(maxSS, ssCount, gapCount, + maxSSCount); + + if (saveFullProfile) + { + profile.setSSCounts(ssCounts); + } + + result[column] = profile; + } + return new Profiles(result); + // long elapsed = System.currentTimeMillis() - now; + // jalview.bin.Console.outPrintln(elapsed); } /** @@ -286,8 +412,60 @@ public class AAFrequency ' ', value); } // long elapsed = System.currentTimeMillis() - now; - // System.out.println(-elapsed); + // jalview.bin.Console.outPrintln(-elapsed); } + + + public static void completeSSConsensus(AlignmentAnnotation ssConsensus, + ProfilesI profiles, int startCol, int endCol, boolean ignoreGaps, + boolean showSequenceLogo, long nseq) + { + // long now = System.currentTimeMillis(); + if (ssConsensus == null || ssConsensus.annotations == null + || ssConsensus.annotations.length < endCol) + { + /* + * called with a bad alignment annotation row + * wait for it to be initialised properly + */ + return; + } + + 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 + */ + ssConsensus.annotations[i] = null; + return; + } + + final int dp = getPercentageDp(nseq); + + float value = profile.getSSPercentageIdentity(ignoreGaps); + + String description = getSSTooltip(profile, value, showSequenceLogo, + ignoreGaps, dp); + + String modalSS = profile.getModalSS(); + if ("".equals(modalSS)) + { + modalSS = "-"; + } + else if (modalSS.length() > 1) + { + modalSS = "+"; + } + ssConsensus.annotations[i] = new Annotation(modalSS, description, + ' ', value); + } + // long elapsed = System.currentTimeMillis() - now; + // jalview.bin.Console.outPrintln(-elapsed); + } /** * Derive the gap count annotation row. @@ -316,7 +494,7 @@ public class AAFrequency // always set ranges again gaprow.graphMax = nseq; gaprow.graphMin = 0; - double scale = 0.8/nseq; + double scale = 0.8 / nseq; for (int i = startCol; i < endCol; i++) { ProfileI profile = profiles.get(i); @@ -332,11 +510,11 @@ public class AAFrequency final int gapped = profile.getNonGapped(); - String description = ""; + String description = "" + gapped; - gaprow.annotations[i] = new Annotation(description, description, - '\0', gapped, jalview.util.ColorUtils.bleachColour( - Color.DARK_GRAY, (float) scale * gapped)); + gaprow.annotations[i] = new Annotation("", description, '\0', gapped, + jalview.util.ColorUtils.bleachColour(Color.DARK_GRAY, + (float) scale * gapped)); } } @@ -345,7 +523,8 @@ public class AAFrequency * * Percentages are as a fraction of all sequence, or only ungapped sequences * if ignoreGaps is true. @@ -366,8 +545,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 @@ -392,13 +571,48 @@ public class AAFrequency } return description; } + + static String getSSTooltip(ProfileI profile, float pid, + boolean showSequenceLogo, boolean ignoreGaps, int dp) + { + SecondaryStructureCount counts = profile.getSSCounts(); + + String description = null; + if (counts != null && showSequenceLogo) + { + int normaliseBy = ignoreGaps ? profile.getNonGapped() + : profile.getHeight(); + description = counts.getTooltip(normaliseBy, dp); + } + else + { + StringBuilder sb = new StringBuilder(64); + String maxSS = profile.getModalSS(); + if (maxSS.length() > 1) + { + sb.append("[").append(maxSS).append("]"); + } + else + { + sb.append(maxSS); + } + if (maxSS.length() > 0) + { + sb.append(" "); + Format.appendPercentage(sb, pid, dp); + sb.append("%"); + } + description = sb.toString(); + } + return description; + } /** * Returns the sorted profile for the given consensus data. The returned array * contains * * 
-   *    [profileType, numberOfValues, nonGapCount, charValue1, percentage1, charValue2, percentage2, ...]
+   *    [profileType, numberOfValues, totalPercent, charValue1, percentage1, charValue2, percentage2, ...]
    * in descending order of percentage value
    *
* @@ -411,40 +625,76 @@ public class AAFrequency */ public static int[] extractProfile(ProfileI profile, boolean ignoreGaps) { - int[] rtnval = new int[64]; - ResidueCount counts = profile.getCounts(); - if (counts == null) + char[] symbols; + int[] values; + + if (profile.getCounts() != null) + { + ResidueCount counts = profile.getCounts(); + SymbolCounts symbolCounts = counts.getSymbolCounts(); + symbols = symbolCounts.symbols; + values = symbolCounts.values; + + } + else if(profile.getSSCounts() != null) { + SecondaryStructureCount counts = profile.getSSCounts(); + // to do + SecondaryStructureCount.SymbolCounts symbolCounts = counts.getSymbolCounts(); + symbols = symbolCounts.symbols; + values = symbolCounts.values; + } + else { return null; } + - SymbolCounts symbolCounts = counts.getSymbolCounts(); - 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; } @@ -454,15 +704,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 @@ -492,9 +742,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; } @@ -518,7 +775,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(); @@ -531,7 +788,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(); @@ -542,8 +799,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); @@ -551,6 +808,7 @@ public class AAFrequency { codonCounts[codonEncoded + 2]++; ungappedCount++; + break; } } } @@ -576,7 +834,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 @@ -591,7 +850,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? @@ -623,10 +882,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 @@ -685,8 +944,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);