src/jalview/datamodel/ResidueCount.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.datamodel;
  22
  23 import jalview.util.Comparison;
  24 import jalview.util.Format;
  25 import jalview.util.QuickSort;
  26 import jalview.util.SparseCount;
  27
  28 /**
  29  * A class to count occurrences of residues in a profile, optimised for speed
  30  * and memory footprint.
  31  *
  32  * @author gmcarstairs
  33  *
  34  */
  35 public class ResidueCount
  36 {
  37   /**
  38    * A data bean to hold the results of counting symbols
  39    */
  40   public class SymbolCounts
  41   {
  42     /**
  43      * the symbols seen (as char values), in no particular order
  44      */
  45     public final char[] symbols;
  46
  47     /**
  48      * the counts for each symbol, in the same order as the symbols
  49      */
  50     public final int[] values;
  51
  52     SymbolCounts(char[] s, int[] v)
  53     {
  54       symbols = s;
  55       values = v;
  56     }
  57   }
  58
  59   private static final int TOUPPERCASE = 'A' - 'a';
  60
  61   /*
  62    * nucleotide symbols to count (including N unknown)
  63    */
  64   private static final String NUCS = "ACGNTU";
  65
  66   /*
  67    * amino acid symbols to count (including X unknown)
  68    * NB we also include U so as to support counting of RNA bases
  69    * in the "don't know" case of nucleotide / peptide
  70    */
  71   private static final String AAS = "ACDEFGHIKLMNPQRSTUVWXY";
  72
  73   static final int GAP_COUNT = 0;
  74
  75   /*
  76    * fast lookup tables holding the index into our count
  77    * arrays of each symbol; index 0 is reserved for gap counting
  78    */
  79   private static int[] NUC_INDEX = new int[26];
  80
  81   private static int[] AA_INDEX = new int[26];
  82   static
  83   {
  84     for (int i = 0; i < NUCS.length(); i++)
  85     {
  86       NUC_INDEX[NUCS.charAt(i) - 'A'] = i + 1;
  87     }
  88     for (int i = 0; i < AAS.length(); i++)
  89     {
  90       AA_INDEX[AAS.charAt(i) - 'A'] = i + 1;
  91     }
  92   }
  93
  94   /*
  95    * counts array, just big enough for the nucleotide or peptide
  96    * character set (plus gap counts in position 0)
  97    */
  98   private short[] counts;
  99
 100   /*
 101    * alternative array of int counts for use if any count
 102    * exceeds the maximum value of short (32767)
 103    */
 104   private int[] intCounts;
 105
 106   /*
 107    * flag set if we switch from short to int counts
 108    */
 109   private boolean useIntCounts;
 110
 111   /*
 112    * general-purpose counter, only for use for characters
 113    * that are not in the expected alphabet
 114    */
 115   private SparseCount otherData;
 116
 117   /*
 118    * keeps track of the maximum count value recorded
 119    * (if this class ever allows decrements, would need to
 120    * calculate this on request instead)
 121    */
 122   int maxCount;
 123
 124   /*
 125    * if we think we are counting nucleotide, can get by with smaller
 126    * array to hold counts
 127    */
 128   private boolean isNucleotide;
 129
 130   /**
 131    * Default constructor allocates arrays able to count either nucleotide or
 132    * peptide bases. Use this constructor if not sure which the data is.
 133    */
 134   public ResidueCount()
 135   {
 136     this(false);
 137   }
 138
 139   /**
 140    * Constructor that allocates an array just big enough for the anticipated
 141    * characters, plus one position to count gaps
 142    */
 143   public ResidueCount(boolean nucleotide)
 144   {
 145     isNucleotide = nucleotide;
 146     int charsToCount = nucleotide ? NUCS.length() : AAS.length();
 147     counts = new short[charsToCount + 1];
 148   }
 149
 150   /**
 151    * Increments the count for the given character. The supplied character may be
 152    * upper or lower case but counts are for the upper case only. Gap characters
 153    * (space, ., -) are all counted together.
 154    *
 155    * @param c
 156    * @return the new value of the count for the character
 157    */
 158   public int add(final char c)
 159   {
 160     char u = toUpperCase(c);
 161     int newValue = 0;
 162     int offset = getOffset(u);
 163
 164     /*
 165      * offset 0 is reserved for gap counting, so 0 here means either
 166      * an unexpected character, or a gap character passed in error
 167      */
 168     if (offset == 0)
 169     {
 170       if (Comparison.isGap(u))
 171       {
 172         newValue = addGap();
 173       }
 174       else
 175       {
 176         newValue = addOtherCharacter(u);
 177       }
 178     }
 179     else
 180     {
 181       newValue = increment(offset);
 182     }
 183     return newValue;
 184   }
 185
 186   /**
 187    * Increment the count at the specified offset. If this would result in short
 188    * overflow, promote to counting int values instead.
 189    *
 190    * @param offset
 191    * @return the new value of the count at this offset
 192    */
 193   int increment(int offset)
 194   {
 195     int newValue = 0;
 196     if (useIntCounts)
 197     {
 198       newValue = intCounts[offset];
 199       intCounts[offset] = ++newValue;
 200     }
 201     else
 202     {
 203       if (counts[offset] == Short.MAX_VALUE)
 204       {
 205         handleOverflow();
 206         newValue = intCounts[offset];
 207         intCounts[offset] = ++newValue;
 208       }
 209       else
 210       {
 211         newValue = counts[offset];
 212         counts[offset] = (short) ++newValue;
 213       }
 214     }
 215     maxCount = Math.max(maxCount, newValue);
 216     return newValue;
 217   }
 218
 219   /**
 220    * Switch from counting in short to counting in int
 221    */
 222   synchronized void handleOverflow()
 223   {
 224     intCounts = new int[counts.length];
 225     for (int i = 0; i < counts.length; i++)
 226     {
 227       intCounts[i] = counts[i];
 228     }
 229     counts = null;
 230     useIntCounts = true;
 231   }
 232
 233   /**
 234    * Returns this character's offset in the count array
 235    *
 236    * @param c
 237    * @return
 238    */
 239   int getOffset(char c)
 240   {
 241     int offset = 0;
 242     if ('A' <= c && c <= 'Z')
 243     {
 244       offset = isNucleotide ? NUC_INDEX[c - 'A'] : AA_INDEX[c - 'A'];
 245     }
 246     return offset;
 247   }
 248
 249   /**
 250    * @param c
 251    * @return
 252    */
 253   protected char toUpperCase(final char c)
 254   {
 255     char u = c;
 256     if ('a' <= c && c <= 'z')
 257     {
 258       u = (char) (c + TOUPPERCASE);
 259     }
 260     return u;
 261   }
 262
 263   /**
 264    * Increment count for some unanticipated character. The first time this
 265    * called, a SparseCount is instantiated to hold these 'extra' counts.
 266    *
 267    * @param c
 268    * @return the new value of the count for the character
 269    */
 270   int addOtherCharacter(char c)
 271   {
 272     if (otherData == null)
 273     {
 274       otherData = new SparseCount();
 275     }
 276     int newValue = otherData.add(c, 1);
 277     maxCount = Math.max(maxCount, newValue);
 278     return newValue;
 279   }
 280
 281   /**
 282    * Set count for some unanticipated character. The first time this called, a
 283    * SparseCount is instantiated to hold these 'extra' counts.
 284    *
 285    * @param c
 286    * @param value
 287    */
 288   void setOtherCharacter(char c, int value)
 289   {
 290     if (otherData == null)
 291     {
 292       otherData = new SparseCount();
 293     }
 294     otherData.put(c, value);
 295   }
 296
 297   /**
 298    * Increment count of gap characters
 299    *
 300    * @return the new count of gaps
 301    */
 302   public int addGap()
 303   {
 304     int newValue = increment(GAP_COUNT);
 305     return newValue;
 306   }
 307
 308   /**
 309    * Answers true if we are counting ints (only after overflow of short counts)
 310    *
 311    * @return
 312    */
 313   boolean isCountingInts()
 314   {
 315     return useIntCounts;
 316   }
 317
 318   /**
 319    * Sets the count for the given character. The supplied character may be upper
 320    * or lower case but counts are for the upper case only.
 321    *
 322    * @param c
 323    * @param count
 324    */
 325   public void put(char c, int count)
 326   {
 327     char u = toUpperCase(c);
 328     int offset = getOffset(u);
 329
 330     /*
 331      * offset 0 is reserved for gap counting, so 0 here means either
 332      * an unexpected character, or a gap character passed in error
 333      */
 334     if (offset == 0)
 335     {
 336       if (Comparison.isGap(u))
 337       {
 338         set(0, count);
 339       }
 340       else
 341       {
 342         setOtherCharacter(u, count);
 343         maxCount = Math.max(maxCount, count);
 344       }
 345     }
 346     else
 347     {
 348       set(offset, count);
 349       maxCount = Math.max(maxCount, count);
 350     }
 351   }
 352
 353   /**
 354    * Sets the count at the specified offset. If this would result in short
 355    * overflow, promote to counting int values instead.
 356    *
 357    * @param offset
 358    * @param value
 359    */
 360   void set(int offset, int value)
 361   {
 362     if (useIntCounts)
 363     {
 364       intCounts[offset] = value;
 365     }
 366     else
 367     {
 368       if (value > Short.MAX_VALUE || value < Short.MIN_VALUE)
 369       {
 370         handleOverflow();
 371         intCounts[offset] = value;
 372       }
 373       else
 374       {
 375         counts[offset] = (short) value;
 376       }
 377     }
 378   }
 379
 380   /**
 381    * Returns the count for the given character, or zero if no count held
 382    *
 383    * @param c
 384    * @return
 385    */
 386   public int getCount(char c)
 387   {
 388     char u = toUpperCase(c);
 389     int offset = getOffset(u);
 390     if (offset == 0)
 391     {
 392       if (!Comparison.isGap(u))
 393       {
 394         // should have called getGapCount()
 395         return otherData == null ? 0 : otherData.get(u);
 396       }
 397     }
 398     return useIntCounts ? intCounts[offset] : counts[offset];
 399   }
 400
 401   public int getGapCount()
 402   {
 403     return useIntCounts ? intCounts[0] : counts[0];
 404   }
 405
 406   /**
 407    * Answers true if this object wraps a counter for unexpected characters
 408    *
 409    * @return
 410    */
 411   boolean isUsingOtherData()
 412   {
 413     return otherData != null;
 414   }
 415
 416   /**
 417    * Returns the character (or concatenated characters) for the symbol(s) with
 418    * the given count in the profile. Can be used to get the modal residue by
 419    * supplying the modal count value. Returns an empty string if no symbol has
 420    * the given count. The symbols are in alphabetic order of standard peptide or
 421    * nucleotide characters, followed by 'other' symbols if any.
 422    *
 423    * @return
 424    */
 425   public String getResiduesForCount(int count)
 426   {
 427     if (count == 0)
 428     {
 429       return "";
 430     }
 431
 432     /*
 433      * find counts for the given value and append the
 434      * corresponding symbol
 435      */
 436     StringBuilder modal = new StringBuilder();
 437     if (useIntCounts)
 438     {
 439       for (int i = 1; i < intCounts.length; i++)
 440       {
 441         if (intCounts[i] == count)
 442         {
 443           modal.append(
 444                   isNucleotide ? NUCS.charAt(i - 1) : AAS.charAt(i - 1));
 445         }
 446       }
 447     }
 448     else
 449     {
 450       for (int i = 1; i < counts.length; i++)
 451       {
 452         if (counts[i] == count)
 453         {
 454           modal.append(
 455                   isNucleotide ? NUCS.charAt(i - 1) : AAS.charAt(i - 1));
 456         }
 457       }
 458     }
 459     if (otherData != null)
 460     {
 461       for (int i = 0; i < otherData.size(); i++)
 462       {
 463         if (otherData.valueAt(i) == count)
 464         {
 465           modal.append((char) otherData.keyAt(i));
 466         }
 467       }
 468     }
 469     return modal.toString();
 470   }
 471
 472   /**
 473    * Returns the highest count for any symbol(s) in the profile (excluding gap)
 474    *
 475    * @return
 476    */
 477   public int getModalCount()
 478   {
 479     return maxCount;
 480   }
 481
 482   /**
 483    * Returns the number of distinct symbols with a non-zero count (excluding the
 484    * gap symbol)
 485    *
 486    * @return
 487    */
 488   public int size()
 489   {
 490     int size = 0;
 491     if (useIntCounts)
 492     {
 493       for (int i = 1; i < intCounts.length; i++)
 494       {
 495         if (intCounts[i] > 0)
 496         {
 497           size++;
 498         }
 499       }
 500     }
 501     else
 502     {
 503       for (int i = 1; i < counts.length; i++)
 504       {
 505         if (counts[i] > 0)
 506         {
 507           size++;
 508         }
 509       }
 510     }
 511
 512     /*
 513      * include 'other' characters recorded (even if count is zero
 514      * though that would be a strange use case)
 515      */
 516     if (otherData != null)
 517     {
 518       size += otherData.size();
 519     }
 520
 521     return size;
 522   }
 523
 524   /**
 525    * Returns a data bean holding those symbols that have a non-zero count
 526    * (excluding the gap symbol), with their counts.
 527    *
 528    * @return
 529    */
 530   public SymbolCounts getSymbolCounts()
 531   {
 532     int size = size();
 533     char[] symbols = new char[size];
 534     int[] values = new int[size];
 535     int j = 0;
 536
 537     if (useIntCounts)
 538     {
 539       for (int i = 1; i < intCounts.length; i++)
 540       {
 541         if (intCounts[i] > 0)
 542         {
 543           char symbol = isNucleotide ? NUCS.charAt(i - 1)
 544                   : AAS.charAt(i - 1);
 545           symbols[j] = symbol;
 546           values[j] = intCounts[i];
 547           j++;
 548         }
 549       }
 550     }
 551     else
 552     {
 553       for (int i = 1; i < counts.length; i++)
 554       {
 555         if (counts[i] > 0)
 556         {
 557           char symbol = isNucleotide ? NUCS.charAt(i - 1)
 558                   : AAS.charAt(i - 1);
 559           symbols[j] = symbol;
 560           values[j] = counts[i];
 561           j++;
 562         }
 563       }
 564     }
 565     if (otherData != null)
 566     {
 567       for (int i = 0; i < otherData.size(); i++)
 568       {
 569         symbols[j] = (char) otherData.keyAt(i);
 570         values[j] = otherData.valueAt(i);
 571         j++;
 572       }
 573     }
 574
 575     return new SymbolCounts(symbols, values);
 576   }
 577
 578   /**
 579    * Returns a tooltip string showing residues in descending order of their
 580    * percentage frequency in the profile
 581    *
 582    * @param normaliseBy
 583    *          the divisor for residue counts (may or may not include gapped
 584    *          sequence count)
 585    * @param percentageDecPl
 586    *          the number of decimal places to show in percentages
 587    * @return
 588    */
 589   public String getTooltip(int normaliseBy, int percentageDecPl)
 590   {
 591     SymbolCounts symbolCounts = getSymbolCounts();
 592     char[] ca = symbolCounts.symbols;
 593     int[] vl = symbolCounts.values;
 594
 595     /*
 596      * sort characters into ascending order of their counts
 597      */
 598     QuickSort.sort(vl, ca);
 599
 600     /*
 601      * traverse in reverse order (highest count first) to build tooltip
 602      */
 603     boolean first = true;
 604     StringBuilder sb = new StringBuilder(64);
 605     for (int c = ca.length - 1; c >= 0; c--)
 606     {
 607       final char residue = ca[c];
 608       // TODO combine residues which share a percentage
 609       // (see AAFrequency.completeCdnaConsensus)
 610       float tval = (vl[c] * 100f) / normaliseBy;
 611       sb.append(first ? "" : "; ").append(residue).append(" ");
 612       Format.appendPercentage(sb, tval, percentageDecPl);
 613       sb.append("%");
 614       first = false;
 615     }
 616     return sb.toString();
 617   }
 618
 619   /**
 620    * Returns a string representation of the symbol counts, for debug purposes.
 621    */
 622   @Override
 623   public String toString()
 624   {
 625     StringBuilder sb = new StringBuilder();
 626     sb.append("[ ");
 627     SymbolCounts sc = getSymbolCounts();
 628     for (int i = 0; i < sc.symbols.length; i++)
 629     {
 630       sb.append(sc.symbols[i]).append(":").append(sc.values[i]).append(" ");
 631     }
 632     sb.append("]");
 633     return sb.toString();
 634   }
 635 }