X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Futil%2FComparison.java;h=9fea705bade9be1b2442d357808c7f8292a986c9;hb=0de37d1a047209510bab82225109ae2a47931f79;hp=d4fc233b901a57835f85c8a3705bb797072fa254;hpb=be762d8d9c71a7aa3121e845c45911c7192b7827;p=jalview.git diff --git a/src/jalview/util/Comparison.java b/src/jalview/util/Comparison.java index d4fc233..9fea705 100644 --- a/src/jalview/util/Comparison.java +++ b/src/jalview/util/Comparison.java @@ -20,11 +20,13 @@ */ package jalview.util; -import jalview.datamodel.SequenceI; - import java.util.ArrayList; import java.util.List; +import jalview.bin.Cache; +import jalview.bin.Console; +import jalview.datamodel.SequenceI; + /** * Assorted methods for analysing or comparing sequences. */ @@ -32,7 +34,15 @@ public class Comparison { private static final int EIGHTY_FIVE = 85; - private static final int TO_UPPER_CASE = 'a' - 'A'; + private static final int NUCLEOTIDE_COUNT_PERCENT; + + private static final int NUCLEOTIDE_COUNT_LONG_SEQUENCE_AMBIGUITY_PERCENT; + + private static final int NUCLEOTIDE_COUNT_SHORT_SEQUENCE; + + private static final int NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE; + + private static final boolean NUCLEOTIDE_AMBIGUITY_DETECTION; public static final char GAP_SPACE = ' '; @@ -44,6 +54,21 @@ public class Comparison new char[] { GAP_SPACE, GAP_DOT, GAP_DASH }); + static + { + // these options read only at start of session + NUCLEOTIDE_COUNT_PERCENT = Cache.getDefault("NUCLEOTIDE_COUNT_PERCENT", + 55); + NUCLEOTIDE_COUNT_LONG_SEQUENCE_AMBIGUITY_PERCENT = Cache.getDefault( + "NUCLEOTIDE_COUNT_LONG_SEQUENCE_AMBIGUITY_PERCENT", 95); + NUCLEOTIDE_COUNT_SHORT_SEQUENCE = Cache + .getDefault("NUCLEOTIDE_COUNT_SHORT", 100); + NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE = Cache + .getDefault("NUCLEOTIDE_COUNT_VERY_SHORT", 4); + NUCLEOTIDE_AMBIGUITY_DETECTION = Cache + .getDefault("NUCLEOTIDE_AMBIGUITY_DETECTION", true); + } + /** * DOCUMENT ME! * @@ -91,7 +116,6 @@ public class Comparison jlen--; } - int count = 0; int match = 0; float pid = -1; @@ -104,8 +128,6 @@ public class Comparison { match++; } - - count++; } pid = (float) match / (float) ilen * 100; @@ -119,8 +141,6 @@ public class Comparison { match++; } - - count++; } pid = (float) match / (float) jlen * 100; @@ -256,19 +276,150 @@ public class Comparison */ public static final boolean isGap(char c) { - return (c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE) ? true : false; + return c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE; } /** * Overloaded method signature to test whether a single sequence is nucleotide - * (that is, more than 85% CGTA) + * (that is, more than 85% CGTAUNX) * * @param seq * @return */ public static final boolean isNucleotide(SequenceI seq) { - return isNucleotide(new SequenceI[] { seq }); + if (seq == null || seq.getLength() == 0) + { + return false; + } + long ntCount = 0; // nucleotide symbol count (does not include ntaCount) + long aaCount = 0; // non-nucleotide, non-gap symbol count (includes nCount + // and ntaCount) + long nCount = 0; // "Unknown" (N) symbol count + long xCount = 0; // Also used as "Unknown" (X) symbol count + long ntaCount = 0; // nucleotide ambiguity symbol count + + int len = seq.getLength(); + for (int i = 0; i < len; i++) + { + char c = seq.getCharAt(i); + if (isNucleotide(c)) + { + ntCount++; + } + else if (!isGap(c)) + { + aaCount++; + if (isN(c)) + { + nCount++; + } + else + { + if (isX(c)) + { + xCount++; + } + if (isNucleotideAmbiguity(c)) + { + ntaCount++; + } + } + } + } + long allCount = ntCount + aaCount; + + if (NUCLEOTIDE_AMBIGUITY_DETECTION) + { + Console.debug("Performing new nucleotide detection routine"); + if (allCount > NUCLEOTIDE_COUNT_SHORT_SEQUENCE) + { + // a long sequence. + // check for at least 55% nucleotide, and nucleotide and ambiguity codes + // (including N) must make up 95% + return ntCount * 100 > NUCLEOTIDE_COUNT_PERCENT * allCount + && 100 * (ntCount + nCount + + ntaCount) > NUCLEOTIDE_COUNT_LONG_SEQUENCE_AMBIGUITY_PERCENT + * allCount; + } + else if (allCount > NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE) + { + // a short sequence. + // check if a short sequence is at least 55% nucleotide and the rest of + // the symbols are all X or all N + if (ntCount * 100 > NUCLEOTIDE_COUNT_PERCENT * allCount + && (nCount == aaCount || xCount == aaCount)) + { + return true; + } + + // a short sequence. + // check for at least x% nucleotide and all the rest nucleotide + // ambiguity codes (including N), where x slides from 75% for sequences + // of length 4 (i.e. only one non-nucleotide) to 55% for sequences of + // length 100 + return myShortSequenceNucleotideProportionCount(ntCount, allCount) + && nCount + ntaCount == aaCount; + } + else + { + // a very short sequence. (<4) + // all bases must be nucleotide + return ntCount > 0 && ntCount == allCount; + } + } + else + { + Console.debug("Performing old nucleotide detection routine"); + /* + * Check for nucleotide count > 85% of total count (in a form that evades + * int / float conversion or divide by zero). + */ + if ((ntCount + nCount) * 100 > EIGHTY_FIVE * allCount) + { + return ntCount > 0; // all N is considered protein. Could use a + // threshold here too + } + } + return false; + } + + protected static boolean myShortSequenceNucleotideProportionCount( + long ntCount, long allCount) + { + /** + * this method is valid only for NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE <= + * allCount <= NUCLEOTIDE_COUNT_SHORT_SEQUENCE + */ + // the following is a simplified integer version of: + // + // a := allCount # the number of bases in the sequence + // n : = ntCount # the number of definite nucleotide bases + // vs := NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE + // s := NUCLEOTIDE_COUNT_SHORT_SEQUENCE + // lp := NUCLEOTIDE_COUNT_LOWER_PERCENT + // vsp := 1 - (1/a) # this is the proportion of required definite + // nucleotides + // # in a VERY_SHORT Sequence (4 bases). + // # This should be equivalent to all but one base in the sequence. + // p := (a - vs)/(s - vs) # proportion of the way between + // # VERY_SHORT and SHORT thresholds. + // tp := vsp + p * (lp/100 - vsp) # the proportion of definite nucleotides + // # required for this length of sequence. + // minNt := tp * a # the minimum number of definite nucleotide bases + // # required for this length of sequence. + // + // We are then essentially returning: + // # ntCount >= 55% of allCount and the rest are all nucleotide ambiguity: + // ntCount >= tp * allCount && nCount + ntaCount == aaCount + // but without going into float/double land + long LHS = 100 * allCount + * (NUCLEOTIDE_COUNT_SHORT_SEQUENCE + - NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE) + * (ntCount - allCount + 1); + long RHS = allCount * (allCount - NUCLEOTIDE_COUNT_VERY_SHORT_SEQUENCE) + * (allCount * NUCLEOTIDE_COUNT_PERCENT - 100 * allCount + 100); + return LHS >= RHS; } /** @@ -285,45 +436,24 @@ public class Comparison { return false; } - - int ntCount = 0; - int aaCount = 0; + // true if we have seen a nucleotide sequence + boolean na = false; for (SequenceI seq : seqs) { if (seq == null) { continue; } + na = true; // TODO could possibly make an informed guess just from the first sequence // to save a lengthy calculation - int len = seq.getLength(); - for (int i = 0; i < len; i++) + if (seq.isProtein()) { - char c = seq.getCharAt(i); - if (isNucleotide(c)) - { - ntCount++; - } - else if (!isGap(c)) - { - aaCount++; - } + // if even one looks like protein, the alignment is protein + return false; } } - - /* - * Check for nucleotide count > 85% of total count (in a form that evades - * int / float conversion or divide by zero). - */ - if (ntCount * 100 > EIGHTY_FIVE * (ntCount + aaCount)) - { - return true; - } - else - { - return false; - } - + return na; } /** @@ -334,12 +464,16 @@ public class Comparison */ public static boolean isNucleotide(char c) { - if ('a' <= c && c <= 'z') - { - c -= TO_UPPER_CASE; - } + return isNucleotide(c, false); + } - switch (c) + /** + * includeAmbiguity = true also includes Nucleotide Ambiguity codes + */ + public static boolean isNucleotide(char c, boolean includeAmbiguity) + { + char C = Character.toUpperCase(c); + switch (C) { case 'A': case 'C': @@ -348,9 +482,50 @@ public class Comparison case 'U': return true; } + if (includeAmbiguity) + { + boolean ambiguity = isNucleotideAmbiguity(C); + if (ambiguity) + return true; + } + return false; + } + + /** + * Tests *only* nucleotide ambiguity codes (and not definite nucleotide codes) + */ + public static boolean isNucleotideAmbiguity(char c) + { + switch (Character.toUpperCase(c)) + { + case 'I': + case 'X': + case 'R': + case 'Y': + case 'W': + case 'S': + case 'M': + case 'K': + case 'B': + case 'H': + case 'D': + case 'V': + return true; + case 'N': // not counting N as nucleotide + } return false; } + public static boolean isN(char c) + { + return 'n' == Character.toLowerCase(c); + } + + public static boolean isX(char c) + { + return 'x' == Character.toLowerCase(c); + } + /** * Answers true if every character in the string is one of aAcCgGtTuU, or * (optionally) a gap character (dot, dash, space), else false @@ -361,6 +536,12 @@ public class Comparison */ public static boolean isNucleotideSequence(String s, boolean allowGaps) { + return isNucleotideSequence(s, allowGaps, false); + } + + public static boolean isNucleotideSequence(String s, boolean allowGaps, + boolean includeAmbiguous) + { if (s == null) { return false; @@ -368,7 +549,7 @@ public class Comparison for (int i = 0; i < s.length(); i++) { char c = s.charAt(i); - if (!isNucleotide(c)) + if (!isNucleotide(c, includeAmbiguous)) { if (!allowGaps || !isGap(c)) { @@ -419,13 +600,7 @@ public class Comparison public static boolean isSameResidue(char c1, char c2, boolean caseSensitive) { - if (caseSensitive) - { - return (c1 == c2); - } - else - { - return Character.toUpperCase(c1) == Character.toUpperCase(c2); - } + return caseSensitive ? c1 == c2 + : Character.toUpperCase(c1) == Character.toUpperCase(c2); } }