src/jalview/util/Comparison.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2)
   3  * Copyright (C) 2014 The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.util;
  22
  23 import jalview.datamodel.SequenceI;
  24
  25 /**
  26  * Assorted methods for analysing or comparing sequences.
  27  */
  28 public class Comparison
  29 {
  30   private static final int EIGHTY_FIVE = 85;
  31
  32   private static final int TO_UPPER_CASE = 'a' - 'A';
  33
  34   public static final String GapChars = " .-";
  35
  36   /**
  37    * DOCUMENT ME!
  38    *
  39    * @param ii
  40    *          DOCUMENT ME!
  41    * @param jj
  42    *          DOCUMENT ME!
  43    *
  44    * @return DOCUMENT ME!
  45    */
  46   public static final float compare(SequenceI ii, SequenceI jj)
  47   {
  48     return Comparison.compare(ii, jj, 0, ii.getLength() - 1);
  49   }
  50
  51   /**
  52    * this was supposed to be an ungapped pid calculation
  53    *
  54    * @param ii
  55    *          SequenceI
  56    * @param jj
  57    *          SequenceI
  58    * @param start
  59    *          int
  60    * @param end
  61    *          int
  62    * @return float
  63    */
  64   public static float compare(SequenceI ii, SequenceI jj, int start, int end)
  65   {
  66     String si = ii.getSequenceAsString();
  67     String sj = jj.getSequenceAsString();
  68
  69     int ilen = si.length() - 1;
  70     int jlen = sj.length() - 1;
  71
  72     while (Comparison.isGap(si.charAt(start + ilen)))
  73     {
  74       ilen--;
  75     }
  76
  77     while (Comparison.isGap(sj.charAt(start + jlen)))
  78     {
  79       jlen--;
  80     }
  81
  82     int count = 0;
  83     int match = 0;
  84     float pid = -1;
  85
  86     if (ilen > jlen)
  87     {
  88       for (int j = 0; j < jlen; j++)
  89       {
  90         if (si.substring(start + j, start + j + 1).equals(
  91                 sj.substring(start + j, start + j + 1)))
  92         {
  93           match++;
  94         }
  95
  96         count++;
  97       }
  98
  99       pid = (float) match / (float) ilen * 100;
 100     }
 101     else
 102     {
 103       for (int j = 0; j < jlen; j++)
 104       {
 105         if (si.substring(start + j, start + j + 1).equals(
 106                 sj.substring(start + j, start + j + 1)))
 107         {
 108           match++;
 109         }
 110
 111         count++;
 112       }
 113
 114       pid = (float) match / (float) jlen * 100;
 115     }
 116
 117     return pid;
 118   }
 119
 120   /**
 121    * this is a gapped PID calculation
 122    *
 123    * @param s1
 124    *          SequenceI
 125    * @param s2
 126    *          SequenceI
 127    * @return float
 128    */
 129   public final static float PID(String seq1, String seq2)
 130   {
 131     return PID(seq1, seq2, 0, seq1.length());
 132   }
 133
 134   static final int caseShift = 'a' - 'A';
 135
 136   // Another pid with region specification
 137   public final static float PID(String seq1, String seq2, int start, int end)
 138   {
 139     return PID(seq1, seq2, start, end, true, false);
 140   }
 141
 142   /**
 143    * Calculate percent identity for a pair of sequences over a particular range,
 144    * with different options for ignoring gaps.
 145    *
 146    * @param seq1
 147    * @param seq2
 148    * @param start
 149    *          - position in seqs
 150    * @param end
 151    *          - position in seqs
 152    * @param wcGaps
 153    *          - if true - gaps match any character, if false, do not match
 154    *          anything
 155    * @param ungappedOnly
 156    *          - if true - only count PID over ungapped columns
 157    * @return
 158    */
 159   public final static float PID(String seq1, String seq2, int start,
 160           int end, boolean wcGaps, boolean ungappedOnly)
 161   {
 162     int s1len = seq1.length();
 163     int s2len = seq2.length();
 164
 165     int len = Math.min(s1len, s2len);
 166
 167     if (end < len)
 168     {
 169       len = end;
 170     }
 171
 172     if (len < start)
 173     {
 174       start = len - 1; // we just use a single residue for the difference
 175     }
 176
 177     int elen = len - start, bad = 0;
 178     char chr1;
 179     char chr2;
 180     boolean agap;
 181     for (int i = start; i < len; i++)
 182     {
 183       chr1 = seq1.charAt(i);
 184
 185       chr2 = seq2.charAt(i);
 186       agap = isGap(chr1) || isGap(chr2);
 187       if ('a' <= chr1 && chr1 <= 'z')
 188       {
 189         // TO UPPERCASE !!!
 190         // Faster than toUpperCase
 191         chr1 -= caseShift;
 192       }
 193       if ('a' <= chr2 && chr2 <= 'z')
 194       {
 195         // TO UPPERCASE !!!
 196         // Faster than toUpperCase
 197         chr2 -= caseShift;
 198       }
 199
 200       if (chr1 != chr2)
 201       {
 202         if (agap)
 203         {
 204           if (ungappedOnly)
 205           {
 206             elen--;
 207           }
 208           else if (!wcGaps)
 209           {
 210             bad++;
 211           }
 212         }
 213         else
 214         {
 215           bad++;
 216         }
 217       }
 218
 219     }
 220     if (elen < 1)
 221     {
 222       return 0f;
 223     }
 224     return ((float) 100 * (elen - bad)) / elen;
 225   }
 226
 227   /**
 228    * Answers true if the supplied character is a recognised gap character, else
 229    * false. Currently hard-coded to recognise '-', '-' or ' ' (hyphen / dot /
 230    * space).
 231    *
 232    * @param c
 233    *
 234    * @return
 235    */
 236   public static final boolean isGap(char c)
 237   {
 238     return (c == '-' || c == '.' || c == ' ') ? true : false;
 239   }
 240
 241   /**
 242    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
 243    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
 244    * wrong answer (as AGCT are also animo acid codes).
 245    *
 246    * @param seqs
 247    * @return
 248    */
 249   public static final boolean isNucleotide(SequenceI[] seqs)
 250   {
 251     if (seqs == null)
 252     {
 253       return false;
 254     }
 255     int ntCount = 0;
 256     int aaCount = 0;
 257     for (SequenceI seq : seqs)
 258     {
 259       for (char c : seq.getSequence())
 260       {
 261         if ('a' <= c && c <= 'z')
 262         {
 263           c -= TO_UPPER_CASE;
 264         }
 265
 266         if (c == 'A' || c == 'G' || c == 'C' || c == 'T' || c == 'U')
 267         {
 268           ntCount++;
 269         }
 270         else if (!Comparison.isGap(c))
 271         {
 272           aaCount++;
 273         }
 274       }
 275     }
 276
 277     /*
 278      * Check for nucleotide count > 85% of total count (in a form that evades
 279      * int / float conversion or divide by zero).
 280      */
 281     if (ntCount * 100 > EIGHTY_FIVE * (ntCount + aaCount))
 282     {
 283       return true;
 284     }
 285     else
 286     {
 287       return false;
 288     }
 289
 290   }
 291 }