src/jalview/util/Comparison.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ 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   private static final char GAP_SPACE = ' ';
  35
  36   private static final char GAP_DOT = '.';
  37
  38   private static final char GAP_DASH = '-';
  39
  40   public static final String GapChars = new String(new char[]
  41   { GAP_SPACE, GAP_DOT, GAP_DASH });
  42
  43   /**
  44    * DOCUMENT ME!
  45    *
  46    * @param ii
  47    *          DOCUMENT ME!
  48    * @param jj
  49    *          DOCUMENT ME!
  50    *
  51    * @return DOCUMENT ME!
  52    */
  53   public static final float compare(SequenceI ii, SequenceI jj)
  54   {
  55     return Comparison.compare(ii, jj, 0, ii.getLength() - 1);
  56   }
  57
  58   /**
  59    * this was supposed to be an ungapped pid calculation
  60    *
  61    * @param ii
  62    *          SequenceI
  63    * @param jj
  64    *          SequenceI
  65    * @param start
  66    *          int
  67    * @param end
  68    *          int
  69    * @return float
  70    */
  71   public static float compare(SequenceI ii, SequenceI jj, int start, int end)
  72   {
  73     String si = ii.getSequenceAsString();
  74     String sj = jj.getSequenceAsString();
  75
  76     int ilen = si.length() - 1;
  77     int jlen = sj.length() - 1;
  78
  79     while (Comparison.isGap(si.charAt(start + ilen)))
  80     {
  81       ilen--;
  82     }
  83
  84     while (Comparison.isGap(sj.charAt(start + jlen)))
  85     {
  86       jlen--;
  87     }
  88
  89     int count = 0;
  90     int match = 0;
  91     float pid = -1;
  92
  93     if (ilen > jlen)
  94     {
  95       for (int j = 0; j < jlen; j++)
  96       {
  97         if (si.substring(start + j, start + j + 1).equals(
  98                 sj.substring(start + j, start + j + 1)))
  99         {
 100           match++;
 101         }
 102
 103         count++;
 104       }
 105
 106       pid = (float) match / (float) ilen * 100;
 107     }
 108     else
 109     {
 110       for (int j = 0; j < jlen; j++)
 111       {
 112         if (si.substring(start + j, start + j + 1).equals(
 113                 sj.substring(start + j, start + j + 1)))
 114         {
 115           match++;
 116         }
 117
 118         count++;
 119       }
 120
 121       pid = (float) match / (float) jlen * 100;
 122     }
 123
 124     return pid;
 125   }
 126
 127   /**
 128    * this is a gapped PID calculation
 129    *
 130    * @param s1
 131    *          SequenceI
 132    * @param s2
 133    *          SequenceI
 134    * @return float
 135    */
 136   public final static float PID(String seq1, String seq2)
 137   {
 138     return PID(seq1, seq2, 0, seq1.length());
 139   }
 140
 141   static final int caseShift = 'a' - 'A';
 142
 143   // Another pid with region specification
 144   public final static float PID(String seq1, String seq2, int start, int end)
 145   {
 146     return PID(seq1, seq2, start, end, true, false);
 147   }
 148
 149   /**
 150    * Calculate percent identity for a pair of sequences over a particular range,
 151    * with different options for ignoring gaps.
 152    *
 153    * @param seq1
 154    * @param seq2
 155    * @param start
 156    *          - position in seqs
 157    * @param end
 158    *          - position in seqs
 159    * @param wcGaps
 160    *          - if true - gaps match any character, if false, do not match
 161    *          anything
 162    * @param ungappedOnly
 163    *          - if true - only count PID over ungapped columns
 164    * @return
 165    */
 166   public final static float PID(String seq1, String seq2, int start,
 167           int end, boolean wcGaps, boolean ungappedOnly)
 168   {
 169     int s1len = seq1.length();
 170     int s2len = seq2.length();
 171
 172     int len = Math.min(s1len, s2len);
 173
 174     if (end < len)
 175     {
 176       len = end;
 177     }
 178
 179     if (len < start)
 180     {
 181       start = len - 1; // we just use a single residue for the difference
 182     }
 183
 184     int elen = len - start, bad = 0;
 185     char chr1;
 186     char chr2;
 187     boolean agap;
 188     for (int i = start; i < len; i++)
 189     {
 190       chr1 = seq1.charAt(i);
 191
 192       chr2 = seq2.charAt(i);
 193       agap = isGap(chr1) || isGap(chr2);
 194       if ('a' <= chr1 && chr1 <= 'z')
 195       {
 196         // TO UPPERCASE !!!
 197         // Faster than toUpperCase
 198         chr1 -= caseShift;
 199       }
 200       if ('a' <= chr2 && chr2 <= 'z')
 201       {
 202         // TO UPPERCASE !!!
 203         // Faster than toUpperCase
 204         chr2 -= caseShift;
 205       }
 206
 207       if (chr1 != chr2)
 208       {
 209         if (agap)
 210         {
 211           if (ungappedOnly)
 212           {
 213             elen--;
 214           }
 215           else if (!wcGaps)
 216           {
 217             bad++;
 218           }
 219         }
 220         else
 221         {
 222           bad++;
 223         }
 224       }
 225
 226     }
 227     if (elen < 1)
 228     {
 229       return 0f;
 230     }
 231     return ((float) 100 * (elen - bad)) / elen;
 232   }
 233
 234   /**
 235    * Answers true if the supplied character is a recognised gap character, else
 236    * false. Currently hard-coded to recognise '-', '-' or ' ' (hyphen / dot /
 237    * space).
 238    *
 239    * @param c
 240    *
 241    * @return
 242    */
 243   public static final boolean isGap(char c)
 244   {
 245     return (c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE) ? true : false;
 246   }
 247
 248   /**
 249    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
 250    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
 251    * wrong answer (as AGCT are also animo acid codes).
 252    *
 253    * @param seqs
 254    * @return
 255    */
 256   public static final boolean isNucleotide(SequenceI[] seqs)
 257   {
 258     if (seqs == null)
 259     {
 260       return false;
 261     }
 262     int ntCount = 0;
 263     int aaCount = 0;
 264     for (SequenceI seq : seqs)
 265     {
 266       for (char c : seq.getSequence())
 267       {
 268         if ('a' <= c && c <= 'z')
 269         {
 270           c -= TO_UPPER_CASE;
 271         }
 272
 273         if (c == 'A' || c == 'G' || c == 'C' || c == 'T' || c == 'U')
 274         {
 275           ntCount++;
 276         }
 277         else if (!Comparison.isGap(c))
 278         {
 279           aaCount++;
 280         }
 281       }
 282     }
 283
 284     /*
 285      * Check for nucleotide count > 85% of total count (in a form that evades
 286      * int / float conversion or divide by zero).
 287      */
 288     if (ntCount * 100 > EIGHTY_FIVE * (ntCount + aaCount))
 289     {
 290       return true;
 291     }
 292     else
 293     {
 294       return false;
 295     }
 296
 297   }
 298 }