X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Futil%2FComparison.java;h=9fea705bade9be1b2442d357808c7f8292a986c9;hb=8479d86d490c841168c9eed1ac80960ce4191c78;hp=cd98ee74bf560bf426c948950cc1e094809ff3a4;hpb=57738a1f3c19b1c3a00bd3ac5108f8cd0af32f99;p=jalview.git diff --git a/src/jalview/util/Comparison.java b/src/jalview/util/Comparison.java index cd98ee7..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,7 +276,7 @@ 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; } /** @@ -268,19 +288,22 @@ public class Comparison */ public static final boolean isNucleotide(SequenceI seq) { - if (seq == null) + if (seq == null || seq.getLength() == 0) { return false; } - long ntCount = 0; - long aaCount = 0; - long nCount = 0; + 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) || isX(c)) + if (isNucleotide(c)) { ntCount++; } @@ -291,21 +314,112 @@ public class Comparison { nCount++; } + else + { + if (isX(c)) + { + xCount++; + } + if (isNucleotideAmbiguity(c)) + { + ntaCount++; + } + } } } - /* - * 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 * (ntCount + aaCount)) + long allCount = ntCount + aaCount; + + if (NUCLEOTIDE_AMBIGUITY_DETECTION) { - return ntCount > 0; // all N is considered protein. Could use a threshold - // here too + 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 { - return false; + 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; } /** @@ -350,11 +464,16 @@ public class Comparison */ public static boolean isNucleotide(char c) { - if ('a' <= c && c <= 'z') - { - c -= TO_UPPER_CASE; - } - switch (c) + return isNucleotide(c, false); + } + + /** + * 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': @@ -363,29 +482,48 @@ public class Comparison case 'U': return true; } + if (includeAmbiguity) + { + boolean ambiguity = isNucleotideAmbiguity(C); + if (ambiguity) + return true; + } return false; } - public static boolean isN(char c) + /** + * Tests *only* nucleotide ambiguity codes (and not definite nucleotide codes) + */ + public static boolean isNucleotideAmbiguity(char c) { - switch (c) + switch (Character.toUpperCase(c)) { - case 'N': - case 'n': + 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) { - switch (c) - { - case 'X': - case 'x': - return true; - } - return false; + return 'x' == Character.toLowerCase(c); } /** @@ -398,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; @@ -405,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)) { @@ -456,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); } }