X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Futil%2FComparison.java;h=9fea705bade9be1b2442d357808c7f8292a986c9;hb=1f41372a2bda3efcc27abaed092495548135b6b4;hp=0beb45bcefca1ad9319236dac9d93bd18af68636;hpb=37de9310bec3501cbc6381e0c3dcb282fcaad812;p=jalview.git

diff --git a/src/jalview/util/Comparison.java b/src/jalview/util/Comparison.java
index 0beb45b..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,16 +34,40 @@ 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 char GAP_SPACE = ' ';
+  private static final boolean NUCLEOTIDE_AMBIGUITY_DETECTION;
 
-  private static final char GAP_DOT = '.';
+  public static final char GAP_SPACE = ' ';
 
-  private static final char GAP_DASH = '-';
+  public static final char GAP_DOT = '.';
 
-  public static final String GapChars = new String(new char[] { GAP_SPACE,
-      GAP_DOT, GAP_DASH });
+  public static final char GAP_DASH = '-';
+
+  public static final String GapChars = new String(
+          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!
@@ -71,7 +97,8 @@ public class Comparison
    *          int
    * @return float
    */
-  public static float compare(SequenceI ii, SequenceI jj, int start, int end)
+  public static float compare(SequenceI ii, SequenceI jj, int start,
+          int end)
   {
     String si = ii.getSequenceAsString();
     String sj = jj.getSequenceAsString();
@@ -89,7 +116,6 @@ public class Comparison
       jlen--;
     }
 
-    int count = 0;
     int match = 0;
     float pid = -1;
 
@@ -97,13 +123,11 @@ public class Comparison
     {
       for (int j = 0; j < jlen; j++)
       {
-        if (si.substring(start + j, start + j + 1).equals(
-                sj.substring(start + j, start + j + 1)))
+        if (si.substring(start + j, start + j + 1)
+                .equals(sj.substring(start + j, start + j + 1)))
         {
           match++;
         }
-
-        count++;
       }
 
       pid = (float) match / (float) ilen * 100;
@@ -112,13 +136,11 @@ public class Comparison
     {
       for (int j = 0; j < jlen; j++)
       {
-        if (si.substring(start + j, start + j + 1).equals(
-                sj.substring(start + j, start + j + 1)))
+        if (si.substring(start + j, start + j + 1)
+                .equals(sj.substring(start + j, start + j + 1)))
         {
           match++;
         }
-
-        count++;
       }
 
       pid = (float) match / (float) jlen * 100;
@@ -135,7 +157,9 @@ public class Comparison
    * @param s2
    *          SequenceI
    * @return float
+   * @deprecated use PIDModel.computePID()
    */
+  @Deprecated
   public final static float PID(String seq1, String seq2)
   {
     return PID(seq1, seq2, 0, seq1.length());
@@ -144,7 +168,12 @@ public class Comparison
   static final int caseShift = 'a' - 'A';
 
   // Another pid with region specification
-  public final static float PID(String seq1, String seq2, int start, int end)
+  /**
+   * @deprecated use PIDModel.computePID()
+   */
+  @Deprecated
+  public final static float PID(String seq1, String seq2, int start,
+          int end)
   {
     return PID(seq1, seq2, start, end, true, false);
   }
@@ -165,7 +194,9 @@ public class Comparison
    * @param ungappedOnly
    *          - if true - only count PID over ungapped columns
    * @return
+   * @deprecated use PIDModel.computePID()
    */
+  @Deprecated
   public final static float PID(String seq1, String seq2, int start,
           int end, boolean wcGaps, boolean ungappedOnly)
   {
@@ -245,49 +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 });
-  }
-
-  /**
-   * Answers true if more than 85% of the sequence residues (ignoring gaps) are
-   * A, G, C, T or U, else false. This is just a heuristic guess and may give a
-   * wrong answer (as AGCT are also amino acid codes).
-   * 
-   * @param seqs
-   * @return
-   */
-  public static final boolean isNucleotide(SequenceI[] seqs)
-  {
-    if (seqs == null)
+    if (seq == null || seq.getLength() == 0)
     {
       return false;
     }
-    char[][] letters = new char[seqs.length][];
-    for (int i = 0; i < seqs.length; i++)
+    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)
     {
-      if (seqs[i] != null)
+      Console.debug("Performing new nucleotide detection routine");
+      if (allCount > NUCLEOTIDE_COUNT_SHORT_SEQUENCE)
       {
-        char[] sequence = seqs[i].getSequence();
-        if (sequence != null)
+        // 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))
         {
-          letters[i] = sequence;
+          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;
+  }
 
-    return areNucleotide(letters);
+  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;
   }
 
   /**
@@ -295,47 +427,33 @@ public class Comparison
    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
    * wrong answer (as AGCT are also amino acid codes).
    * 
-   * @param letters
+   * @param seqs
    * @return
    */
-  static final boolean areNucleotide(char[][] letters)
+  public static final boolean isNucleotide(SequenceI[] seqs)
   {
-    int ntCount = 0;
-    int aaCount = 0;
-    for (char[] seq : letters)
+    if (seqs == null)
+    {
+      return false;
+    }
+    // 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
-      for (char c : seq)
+      if (seq.isProtein())
       {
-        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;
   }
 
   /**
@@ -346,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':
@@ -360,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
@@ -373,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;
@@ -380,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))
         {
@@ -411,8 +580,27 @@ public class Comparison
         flattened.add(s);
       }
     }
-    final SequenceI[] oneDArray = flattened.toArray(new SequenceI[flattened
-            .size()]);
+    final SequenceI[] oneDArray = flattened
+            .toArray(new SequenceI[flattened.size()]);
     return isNucleotide(oneDArray);
   }
+
+  /**
+   * Compares two residues either case sensitively or case insensitively
+   * depending on the caseSensitive flag
+   * 
+   * @param c1
+   *          first char
+   * @param c2
+   *          second char to compare with
+   * @param caseSensitive
+   *          if true comparison will be case sensitive otherwise its not
+   * @return
+   */
+  public static boolean isSameResidue(char c1, char c2,
+          boolean caseSensitive)
+  {
+    return caseSensitive ? c1 == c2
+            : Character.toUpperCase(c1) == Character.toUpperCase(c2);
+  }
 }