src/jalview/datamodel/CigarBase.java

   1 /*\r
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.7)\r
   3  * Copyright (C) 2011 J Procter, AM Waterhouse, J Engelhardt, LM Lui, G Barton, M Clamp, S Searle\r
   4  * \r
   5  * This file is part of Jalview.\r
   6  * \r
   7  * Jalview is free software: you can redistribute it and/or\r
   8  * modify it under the terms of the GNU General Public License \r
   9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.\r
  10  * \r
  11  * Jalview is distributed in the hope that it will be useful, but \r
  12  * WITHOUT ANY WARRANTY; without even the implied warranty \r
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR \r
  14  * PURPOSE.  See the GNU General Public License for more details.\r
  15  * \r
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.\r
  17  */\r
  18 package jalview.datamodel;\r
  19 \r
  20 import java.util.*;\r
  21 \r
  22 public abstract class CigarBase\r
  23 {\r
  24   /**\r
  25    * Base class for compact idiosyncratic representation of gaps and aligned\r
  26    * residues Regards to Tom Oldfield for his DynamicArray class. 17th July 2006\r
  27    * Not thread safe.\r
  28    */\r
  29   public CigarBase()\r
  30   {\r
  31     // nothing to be done (probably)\r
  32   }\r
  33 \r
  34   protected int length = 0;\r
  35 \r
  36   protected int _inc_length = 10; // extension range for addition of new\r
  37 \r
  38   // operations\r
  39 \r
  40   protected char[] operation = null;\r
  41 \r
  42   protected int[] range = null;\r
  43 \r
  44   /**\r
  45    * Range of Hidden residues in seq (translated as deleted in seq)\r
  46    */\r
  47   public static final char D = 'D';\r
  48 \r
  49   /**\r
  50    * Range of insertions to seq\r
  51    */\r
  52   public static final char I = 'I';\r
  53 \r
  54   /**\r
  55    * Range of aligned residues\r
  56    */\r
  57   public static final char M = 'M';\r
  58 \r
  59   static protected final char _case_shift = 'a' - 'A';\r
  60 \r
  61   /**\r
  62    * Ugly function to get edited sequence string, start and end symbol positions\r
  63    * and the deletion regions as an array of int pairs May return null for an\r
  64    * empty cigar string. May return null for deletion ranges if there are none.\r
  65    * \r
  66    * @param reference\r
  67    *          - the symbol sequence to apply the cigar operations to (or null if\r
  68    *          no sequence)\r
  69    * @param GapChar\r
  70    *          - the symbol to use for Insert operations\r
  71    * @return Object[] { String, int[] {start, startcol, end, endcol}, int[][3]\r
  72    *         {start, end, col} or null} the gapped sequence, first and last\r
  73    *         residue index, and the deletion ranges on the reference sequence\r
  74    */\r
  75   public Object[] getSequenceAndDeletions(String reference, char GapChar)\r
  76   {\r
  77     int rlength = 0;\r
  78     int[][] deletions = new int[length][];\r
  79     int[][] trunc_deletions = null;\r
  80     StringBuffer sq = new StringBuffer();\r
  81     int cursor = 0, alcursor = 0, start = 0, startpos = 0, end = 0, endpos = 0, delcount = -1;\r
  82     boolean consecutive_del = false;\r
  83     if (length == 0)\r
  84     {\r
  85       return null;\r
  86     }\r
  87     if (reference != null)\r
  88     {\r
  89       rlength = reference.length();\r
  90     }\r
  91     boolean modstart = true;\r
  92     for (int i = 0; i < length; i++)\r
  93     {\r
  94       switch (operation[i])\r
  95       {\r
  96       case D:\r
  97         if (!consecutive_del)\r
  98         {\r
  99           deletions[++delcount] = new int[]\r
 100           { cursor, 0, alcursor };\r
 101         }\r
 102         cursor += range[i];\r
 103         deletions[delcount][1] = cursor - 1;\r
 104         consecutive_del = true;\r
 105         break;\r
 106       case I:\r
 107         consecutive_del = false;\r
 108         for (int r = 0; r < range[i]; r++)\r
 109         {\r
 110           sq.append(GapChar);\r
 111           alcursor++;\r
 112         }\r
 113         break;\r
 114       case M:\r
 115         consecutive_del = false;\r
 116         if (modstart)\r
 117         {\r
 118           start = cursor;\r
 119           startpos = alcursor;\r
 120           modstart = false;\r
 121         }\r
 122         if (reference != null)\r
 123         {\r
 124           int sbend = cursor + range[i];\r
 125           if (sbend > rlength)\r
 126           {\r
 127             sq.append(reference.substring(cursor, rlength));\r
 128             while (sbend-- >= rlength)\r
 129             {\r
 130               sq.append(GapChar);\r
 131             }\r
 132           }\r
 133           else\r
 134           {\r
 135             sq.append(reference.substring(cursor, sbend));\r
 136           }\r
 137         }\r
 138         alcursor += range[i];\r
 139         cursor += range[i];\r
 140         end = cursor - 1;\r
 141         endpos = alcursor;\r
 142         break;\r
 143       default:\r
 144         throw new Error("Unknown SeqCigar operation '" + operation[i] + "'");\r
 145       }\r
 146     }\r
 147     if (++delcount > 0)\r
 148     {\r
 149       trunc_deletions = new int[delcount][];\r
 150       System.arraycopy(deletions, 0, trunc_deletions, 0, delcount);\r
 151     }\r
 152     deletions = null;\r
 153     return new Object[]\r
 154     { ((reference != null) ? sq.toString() : null), new int[]\r
 155     { start, startpos, end, endpos }, trunc_deletions };\r
 156   }\r
 157 \r
 158   protected void compact_operations()\r
 159   {\r
 160     int i = 1;\r
 161     if (operation == null)\r
 162     {\r
 163       return;\r
 164     }\r
 165     char last = operation[0];\r
 166     while (i < length)\r
 167     {\r
 168       if (last == operation[i])\r
 169       {\r
 170         range[i - 1] += range[i];\r
 171         int r = length - i;\r
 172         if (r > 0)\r
 173         {\r
 174           System.arraycopy(range, i + 1, range, i, r);\r
 175           System.arraycopy(operation, i + 1, operation, i, r);\r
 176         }\r
 177         length--;\r
 178       }\r
 179       else\r
 180       {\r
 181         last = operation[i++];\r
 182       }\r
 183     }\r
 184   }\r
 185 \r
 186   /**\r
 187    * turn a cigar string into a series of operation range pairs\r
 188    * \r
 189    * @param cigarString\r
 190    *          String\r
 191    * @return object[] {char[] operation, int[] range}\r
 192    * @throws java.lang.Exception\r
 193    *           for improperly formated cigar strings or ones with unknown\r
 194    *           operations\r
 195    */\r
 196   public static Object[] parseCigarString(String cigarString)\r
 197           throws Exception\r
 198   {\r
 199     int ops = 0;\r
 200     for (int i = 0, l = cigarString.length(); i < l; i++)\r
 201     {\r
 202       char c = cigarString.charAt(i);\r
 203       if (c == M || c == (M - _case_shift) || c == I\r
 204               || c == (I - _case_shift) || c == D || c == (D - _case_shift))\r
 205       {\r
 206         ops++;\r
 207       }\r
 208     }\r
 209     char[] operation = new char[ops];\r
 210     int[] range = new int[ops];\r
 211     int op = 0;\r
 212     int i = 0, l = cigarString.length();\r
 213     while (i < l)\r
 214     {\r
 215       char c;\r
 216       int j = i;\r
 217       do\r
 218       {\r
 219         c = cigarString.charAt(j++);\r
 220       } while (c >= '0' && c <= '9' && j < l);\r
 221       if (j >= l && c >= '0' && c <= '9')\r
 222       {\r
 223         throw new Exception("Unterminated cigar string.");\r
 224       }\r
 225       try\r
 226       {\r
 227         String rangeint = cigarString.substring(i, j - 1);\r
 228         range[op] = Integer.parseInt(rangeint);\r
 229         i = j;\r
 230       } catch (Exception e)\r
 231       {\r
 232         throw new Error("Implementation bug in parseCigarString");\r
 233       }\r
 234       if (c >= 'a' && c <= 'z')\r
 235       {\r
 236         c -= _case_shift;\r
 237       }\r
 238       if ((c == M || c == I || c == D))\r
 239       {\r
 240         operation[op++] = c;\r
 241       }\r
 242       else\r
 243       {\r
 244         throw new Exception("Unexpected operation '" + c\r
 245                 + "' in cigar string (position " + i + " in '"\r
 246                 + cigarString + "'");\r
 247       }\r
 248     }\r
 249     return new Object[]\r
 250     { operation, range };\r
 251   }\r
 252 \r
 253   /**\r
 254    * add an operation to cigar string\r
 255    * \r
 256    * @param op\r
 257    *          char\r
 258    * @param range\r
 259    *          int\r
 260    */\r
 261   public void addOperation(char op, int range)\r
 262   {\r
 263     if (op >= 'a' && op <= 'z')\r
 264     {\r
 265       op -= _case_shift;\r
 266     }\r
 267     if (op != M && op != D && op != I)\r
 268     {\r
 269       throw new Error("Implementation error. Invalid operation string.");\r
 270     }\r
 271     if (range==0)\r
 272     {\r
 273       return; // No Operation to add.\r
 274     }\r
 275     if (range < 0)\r
 276     {\r
 277       throw new Error(\r
 278               "Invalid range string (must be zero or positive number)");\r
 279     }\r
 280     int lngth = 0;\r
 281     if (operation == null)\r
 282     {\r
 283       this.operation = new char[_inc_length];\r
 284       this.range = new int[_inc_length];\r
 285     }\r
 286     if (length + 1 == operation.length)\r
 287     {\r
 288       char[] ops = this.operation;\r
 289       this.operation = new char[length + 1 + _inc_length];\r
 290       System.arraycopy(ops, 0, this.operation, 0, length);\r
 291       ops = null;\r
 292       int[] rng = this.range;\r
 293       this.range = new int[length + 1 + _inc_length];\r
 294       System.arraycopy(rng, 0, this.range, 0, length);\r
 295       rng = null;\r
 296     }\r
 297     if ((length > 0) && (operation[length - 1] == op))\r
 298     {\r
 299       length--; // modify existing operation.\r
 300     }\r
 301     else\r
 302     {\r
 303       this.range[length] = 0; // reset range\r
 304     }\r
 305     this.operation[length] = op;\r
 306     this.range[length++] += range;\r
 307   }\r
 308 \r
 309   /**\r
 310    * semi-efficient insert an operation on the current cigar string set at\r
 311    * column pos (from 1) NOTE: Insertion operations simply extend width of cigar\r
 312    * result - affecting registration of alignment Deletion ops will shorten\r
 313    * length of result - and affect registration of alignment Match ops will also\r
 314    * affect length of result - affecting registration of alignment (ie\r
 315    * "10M".insert(4,I,3)->"4M3I3M") - (replace?) (ie\r
 316    * "10M".insert(4,D,3)->"4M3D3M") - (shortens alignment) (ie\r
 317    * "5I5M".insert(4,I,3)->"8I5M") - real insertion (ie\r
 318    * "5I5M".insert(4,D,3)->"4I2D3M") - shortens aligment - I's are removed, Ms\r
 319    * changed to Ds (ie "10M".insert(4,M,3)->"13M") - lengthens - Is changed to\r
 320    * M, Ds changed to M. (ie "5I5M".insert(4,M,3)->"4I8M") - effectively shifts\r
 321    * sequence left by 1 residue and extends it by 3 (\r
 322    * "10D5M".insert(-1,M,3)->"3M7D5M") ( "10D5M".insert(0,M,3)->"7D8M") (\r
 323    * "10D5M".insert(1,M,3)->"10D8M") ( "1M10D5M".insert(0,M,3)->"1M10D8M") (\r
 324    * "1M10D5M".insert(1,M,3)->"\r
 325    * \r
 326    * if pos is beyond width - I operations are added before the operation\r
 327    * \r
 328    * @param pos\r
 329    *          int -1, 0-length of visible region, or greater to append new ops\r
 330    *          (with insertions in between)\r
 331    * @param op\r
 332    *          char\r
 333    * @param range\r
 334    *          int public void addOperationAt(int pos, char op, int range) { int\r
 335    *          cursor = -1; // mark the position for the current operation being\r
 336    *          edited. int o = 0; boolean last_d = false; // previous op was a\r
 337    *          deletion. if (pos < -1) throw new\r
 338    *          Error("pos<-1 is not supported."); while (o<length) { if\r
 339    *          (operation[o] != D) { if ( (cursor + this.range[o]) < pos) {\r
 340    *          cursor += this.range[o]; o++; last_d=false; } else { break; } }\r
 341    *          else { last_d=true; o++; } } if (o==length) { // must insert more\r
 342    *          operations before pos if (pos-cursor>0) addInsertion(pos-cursor);\r
 343    *          // then just add the new operation. Regardless of what it is.\r
 344    *          addOperation(op, range); } else { int diff = pos - cursor;\r
 345    * \r
 346    *          int e_length = length-o; // new edit operation array length. //\r
 347    *          diff<0 - can only happen before first insertion or match. -\r
 348    *          affects op and all following // dif==0 - only when at first\r
 349    *          position of existing op - // diff>0 - must preserve some existing\r
 350    *          operations int[] e_range = new int[e_length];\r
 351    *          System.arraycopy(this.range, o, e_range, 0, e_length); char[] e_op\r
 352    *          = new char[e_length]; System.arraycopy(this.operation, o, e_op, 0,\r
 353    *          e_length); length = o; // can now use add_operation to extend\r
 354    *          list. int e_o=0; // current operation being edited. switch (op) {\r
 355    *          case M: switch (e_op[e_o]) { case M: if (last_d && diff <= 0) { //\r
 356    *          reduce D's, if possible if (range<=this.range[o-1]) { this.range[o\r
 357    *          - 1] -= range; } else { this.range[o-1]=0; } if\r
 358    *          (this.range[o-1]==0) o--; // lose this op. } e_range[e_o] +=\r
 359    *          range; // just add more matched residues break; case I: // change\r
 360    *          from insertion to match if (last_d && diff<=0) { // reduce D's, if\r
 361    *          possible if (range<=this.range[o-1]) { this.range[o - 1] -= range;\r
 362    *          } else { this.range[o-1]=0; } if (this.range[o-1]==0) o--; // lose\r
 363    *          this op. } e_range[e_o] break; default: throw new Inp }\r
 364    * \r
 365    *          break; case I: break; case D: } break; default: throw new\r
 366    *          Error("Implementation Error: Unknown operation in addOperation!");\r
 367    *          } // finally, add remaining ops. while (e_o<e_length) {\r
 368    *          addOperation(e_op[e_o], e_range[e_o]); e_o++; } } }\r
 369    */\r
 370   /**\r
 371    * Mark residues from start to end (inclusive) as deleted from the alignment,\r
 372    * and removes any insertions.\r
 373    * \r
 374    * @param start\r
 375    *          int\r
 376    * @param end\r
 377    *          int\r
 378    * @return deleted int - number of symbols marked as deleted\r
 379    */\r
 380   public int deleteRange(int start, int end)\r
 381   {\r
 382     int deleted = 0;\r
 383     if (length == 0)\r
 384     {\r
 385       // nothing to do here\r
 386       return deleted;\r
 387     }\r
 388     if (start < 0 || start > end)\r
 389     {\r
 390       throw new Error(\r
 391               "Implementation Error: deleteRange out of bounds: start must be non-negative and less than end.");\r
 392     }\r
 393     // find beginning\r
 394     int cursor = 0; // mark the position for the current operation being edited.\r
 395     int rlength = 1 + end - start; // number of positions to delete\r
 396     int oldlen = length;\r
 397     int o = 0;\r
 398     boolean editing = false;\r
 399     char[] oldops = operation;\r
 400     int[] oldrange = range;\r
 401     length = 0;\r
 402     operation = null;\r
 403     range = null;\r
 404     compact_operations();\r
 405     while (o < oldlen && cursor <= end && rlength > 0)\r
 406     {\r
 407       if (oldops[o] == D)\r
 408       {\r
 409         // absorbed into new deleted region.\r
 410         addDeleted(oldrange[o++]);\r
 411         continue;\r
 412       }\r
 413 \r
 414       int remain = oldrange[o]; // number of op characters left to edit\r
 415       if (!editing)\r
 416       {\r
 417         if ((cursor + remain) <= start)\r
 418         {\r
 419           addOperation(oldops[o], oldrange[o]);\r
 420           cursor += oldrange[o++];\r
 421           continue; // next operation\r
 422         }\r
 423         editing = true;\r
 424         // add operations before hidden region\r
 425         if (start - cursor > 0)\r
 426         {\r
 427           addOperation(oldops[o], start - cursor);\r
 428           remain -= start - cursor;\r
 429         }\r
 430       }\r
 431       // start inserting new ops\r
 432       if (o < oldlen && editing && rlength > 0 && remain > 0)\r
 433       {\r
 434         switch (oldops[o])\r
 435         {\r
 436         case M:\r
 437           if (rlength > remain)\r
 438           {\r
 439             addDeleted(remain);\r
 440             deleted += remain;\r
 441           }\r
 442           else\r
 443           {\r
 444             deleted += rlength;\r
 445             addDeleted(rlength);\r
 446             if (remain - rlength > 0)\r
 447             {\r
 448               this.addOperation(M, remain - rlength); // add remaining back.\r
 449             }\r
 450             rlength = 0;\r
 451             remain = 0;\r
 452           }\r
 453           break;\r
 454         case I:\r
 455           if (remain - rlength > 0)\r
 456           {\r
 457             // only remove some gaps\r
 458             addInsertion(remain - rlength);\r
 459             rlength = 0;\r
 460           }\r
 461           break;\r
 462         case D:\r
 463           throw new Error("Implementation error."); // do nothing;\r
 464         default:\r
 465           throw new Error("Implementation Error! Unknown operation '"\r
 466                   + oldops[o] + "'");\r
 467         }\r
 468         rlength -= remain;\r
 469         remain = oldrange[++o]; // number of op characters left to edit\r
 470       }\r
 471     }\r
 472     // add remaining\r
 473     while (o < oldlen)\r
 474     {\r
 475       addOperation(oldops[o], oldrange[o++]);\r
 476     }\r
 477     // if (cursor<(start+1)) {\r
 478     // ran out of ops - nothing to do here ?\r
 479     // addInsertion(start-cursor);\r
 480     // }\r
 481     return deleted;\r
 482   }\r
 483 \r
 484   /**\r
 485    * Deleted regions mean that there will be discontinuous sequence numbering in\r
 486    * the sequence returned by getSeq(char).\r
 487    * \r
 488    * @return true if there deletions\r
 489    */\r
 490   public boolean hasDeletedRegions()\r
 491   {\r
 492     for (int i = 0; i < length; i++)\r
 493     {\r
 494       if (operation[i] == D)\r
 495       {\r
 496         return true;\r
 497       }\r
 498     }\r
 499     return false;\r
 500   }\r
 501 \r
 502   /**\r
 503    * enumerate the ranges on seq that are marked as deleted in this cigar\r
 504    * \r
 505    * @return int[] { vis_start, sym_start, length }\r
 506    */\r
 507   public int[] getDeletedRegions()\r
 508   {\r
 509     if (length == 0)\r
 510     {\r
 511       return null;\r
 512     }\r
 513     Vector dr = new Vector();\r
 514     int cursor = 0, vcursor = 0;\r
 515     for (int i = 0; i < length; i++)\r
 516     {\r
 517       switch (operation[i])\r
 518       {\r
 519       case M:\r
 520         cursor += range[i];\r
 521       case I:\r
 522         vcursor += range[i];\r
 523         break;\r
 524       case D:\r
 525         dr.addElement(new int[]\r
 526         { vcursor, cursor, range[i] });\r
 527         cursor += range[i];\r
 528       }\r
 529     }\r
 530     if (dr.size() == 0)\r
 531     {\r
 532       return null;\r
 533     }\r
 534     int[] delregions = new int[dr.size() * 3];\r
 535     for (int i = 0, l = dr.size(); i < l; i++)\r
 536     {\r
 537       int[] reg = (int[]) dr.elementAt(i);\r
 538       delregions[i * 3] = reg[0];\r
 539       delregions[i * 3 + 1] = reg[1];\r
 540       delregions[i * 3 + 2] = reg[2];\r
 541     }\r
 542     return delregions;\r
 543   }\r
 544 \r
 545   /**\r
 546    * sum of ranges in cigar string\r
 547    * \r
 548    * @return int number of residues hidden, matched, or gaps inserted into\r
 549    *         sequence\r
 550    */\r
 551   public int getFullWidth()\r
 552   {\r
 553     int w = 0;\r
 554     if (range != null)\r
 555     {\r
 556       for (int i = 0; i < length; i++)\r
 557       {\r
 558         w += range[i];\r
 559       }\r
 560     }\r
 561     return w;\r
 562   }\r
 563 \r
 564   /**\r
 565    * Visible length of aligned sequence\r
 566    * \r
 567    * @return int length of including gaps and less hidden regions\r
 568    */\r
 569   public int getWidth()\r
 570   {\r
 571     int w = 0;\r
 572     if (range != null)\r
 573     {\r
 574       for (int i = 0; i < length; i++)\r
 575       {\r
 576         if (operation[i] == M || operation[i] == I)\r
 577         {\r
 578           w += range[i];\r
 579         }\r
 580       }\r
 581     }\r
 582     return w;\r
 583   }\r
 584 \r
 585   /**\r
 586    * mark a range of inserted residues\r
 587    * \r
 588    * @param range\r
 589    *          int\r
 590    */\r
 591   public void addInsertion(int range)\r
 592   {\r
 593     this.addOperation(I, range);\r
 594   }\r
 595 \r
 596   /**\r
 597    * mark the next range residues as hidden (not aligned) or deleted\r
 598    * \r
 599    * @param range\r
 600    *          int\r
 601    */\r
 602   public void addDeleted(int range)\r
 603   {\r
 604     this.addOperation(D, range);\r
 605   }\r
 606 \r
 607   /**\r
 608    * Modifies operation list to delete columns from start to end (inclusive)\r
 609    * editing will remove insertion operations, and convert matches to deletions\r
 610    * \r
 611    * @param start\r
 612    *          alignment column\r
 613    * @param end\r
 614    *          alignment column\r
 615    * @return boolean true if residues were marked as deleted. public boolean\r
 616    *         deleteRange(int start, int end) { boolean deleted = false; int op =\r
 617    *         0, prevop = -1, firstm = -1, lastm = -1, postop = -1; int width =\r
 618    *         0; // zero'th column if (length > 0) { // find operation bracketing\r
 619    *         start of the range do { if (operation[op] != D) { width +=\r
 620    *         range[prevop = op]; } op++; } while (op < length && width < start);\r
 621    *         } if (width < start) { // run off end - add more operations up to\r
 622    *         deletion. addInsertion(start - width); } else { // edit existing\r
 623    *         operations. op = prevop; width -= range[prevop]; int[] oldrange =\r
 624    *         range; char[] oldops = operation; range = new int[oldrange.length];\r
 625    *         operation = new char[oldops.length]; if (op < length) { do { if\r
 626    *         (operation[op] != D) { width += range[postop = op]; } op++; } while\r
 627    *         (op < length && width <= end); } } if (deleted == true) {\r
 628    *         addDeleted(end - start + 1); } return deleted; }\r
 629    */\r
 630   /**\r
 631    * Return an ENSEMBL style cigar string where D may indicates excluded parts\r
 632    * of seq\r
 633    * \r
 634    * @return String of form ([0-9]+[IMD])+\r
 635    */\r
 636   public String getCigarstring()\r
 637   {\r
 638     StringBuffer cigarString = new StringBuffer();\r
 639     for (int i = 0; i < length; i++)\r
 640     {\r
 641       cigarString.append("" + range[i]);\r
 642       cigarString.append(operation[i]);\r
 643     }\r
 644     return cigarString.toString();\r
 645   }\r
 646 }\r