src/jalview/datamodel/CigarBase.java

   1 /*\r
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)\r
   3  * Copyright (C) 2010 J Procter, AM Waterhouse, 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       throw new Error(\r
 274               "Invalid range string (must be non-zero positive number)");\r
 275     }\r
 276     int lngth = 0;\r
 277     if (operation == null)\r
 278     {\r
 279       this.operation = new char[_inc_length];\r
 280       this.range = new int[_inc_length];\r
 281     }\r
 282     if (length + 1 == operation.length)\r
 283     {\r
 284       char[] ops = this.operation;\r
 285       this.operation = new char[length + 1 + _inc_length];\r
 286       System.arraycopy(ops, 0, this.operation, 0, length);\r
 287       ops = null;\r
 288       int[] rng = this.range;\r
 289       this.range = new int[length + 1 + _inc_length];\r
 290       System.arraycopy(rng, 0, this.range, 0, length);\r
 291       rng = null;\r
 292     }\r
 293     if ((length > 0) && (operation[length - 1] == op))\r
 294     {\r
 295       length--; // modify existing operation.\r
 296     }\r
 297     else\r
 298     {\r
 299       this.range[length] = 0; // reset range\r
 300     }\r
 301     this.operation[length] = op;\r
 302     this.range[length++] += range;\r
 303   }\r
 304 \r
 305   /**\r
 306    * semi-efficient insert an operation on the current cigar string set at\r
 307    * column pos (from 1) NOTE: Insertion operations simply extend width of cigar\r
 308    * result - affecting registration of alignment Deletion ops will shorten\r
 309    * length of result - and affect registration of alignment Match ops will also\r
 310    * affect length of result - affecting registration of alignment (ie\r
 311    * "10M".insert(4,I,3)->"4M3I3M") - (replace?) (ie\r
 312    * "10M".insert(4,D,3)->"4M3D3M") - (shortens alignment) (ie\r
 313    * "5I5M".insert(4,I,3)->"8I5M") - real insertion (ie\r
 314    * "5I5M".insert(4,D,3)->"4I2D3M") - shortens aligment - I's are removed, Ms\r
 315    * changed to Ds (ie "10M".insert(4,M,3)->"13M") - lengthens - Is changed to\r
 316    * M, Ds changed to M. (ie "5I5M".insert(4,M,3)->"4I8M") - effectively shifts\r
 317    * sequence left by 1 residue and extends it by 3 (\r
 318    * "10D5M".insert(-1,M,3)->"3M7D5M") ( "10D5M".insert(0,M,3)->"7D8M") (\r
 319    * "10D5M".insert(1,M,3)->"10D8M") ( "1M10D5M".insert(0,M,3)->"1M10D8M") (\r
 320    * "1M10D5M".insert(1,M,3)->"\r
 321    * \r
 322    * if pos is beyond width - I operations are added before the operation\r
 323    * \r
 324    * @param pos\r
 325    *          int -1, 0-length of visible region, or greater to append new ops\r
 326    *          (with insertions in between)\r
 327    * @param op\r
 328    *          char\r
 329    * @param range\r
 330    *          int public void addOperationAt(int pos, char op, int range) { int\r
 331    *          cursor = -1; // mark the position for the current operation being\r
 332    *          edited. int o = 0; boolean last_d = false; // previous op was a\r
 333    *          deletion. if (pos < -1) throw new\r
 334    *          Error("pos<-1 is not supported."); while (o<length) { if\r
 335    *          (operation[o] != D) { if ( (cursor + this.range[o]) < pos) {\r
 336    *          cursor += this.range[o]; o++; last_d=false; } else { break; } }\r
 337    *          else { last_d=true; o++; } } if (o==length) { // must insert more\r
 338    *          operations before pos if (pos-cursor>0) addInsertion(pos-cursor);\r
 339    *          // then just add the new operation. Regardless of what it is.\r
 340    *          addOperation(op, range); } else { int diff = pos - cursor;\r
 341    * \r
 342    *          int e_length = length-o; // new edit operation array length. //\r
 343    *          diff<0 - can only happen before first insertion or match. -\r
 344    *          affects op and all following // dif==0 - only when at first\r
 345    *          position of existing op - // diff>0 - must preserve some existing\r
 346    *          operations int[] e_range = new int[e_length];\r
 347    *          System.arraycopy(this.range, o, e_range, 0, e_length); char[] e_op\r
 348    *          = new char[e_length]; System.arraycopy(this.operation, o, e_op, 0,\r
 349    *          e_length); length = o; // can now use add_operation to extend\r
 350    *          list. int e_o=0; // current operation being edited. switch (op) {\r
 351    *          case M: switch (e_op[e_o]) { case M: if (last_d && diff <= 0) { //\r
 352    *          reduce D's, if possible if (range<=this.range[o-1]) { this.range[o\r
 353    *          - 1] -= range; } else { this.range[o-1]=0; } if\r
 354    *          (this.range[o-1]==0) o--; // lose this op. } e_range[e_o] +=\r
 355    *          range; // just add more matched residues break; case I: // change\r
 356    *          from insertion to match if (last_d && diff<=0) { // reduce D's, if\r
 357    *          possible if (range<=this.range[o-1]) { this.range[o - 1] -= range;\r
 358    *          } else { this.range[o-1]=0; } if (this.range[o-1]==0) o--; // lose\r
 359    *          this op. } e_range[e_o] break; default: throw new Inp }\r
 360    * \r
 361    *          break; case I: break; case D: } break; default: throw new\r
 362    *          Error("Implementation Error: Unknown operation in addOperation!");\r
 363    *          } // finally, add remaining ops. while (e_o<e_length) {\r
 364    *          addOperation(e_op[e_o], e_range[e_o]); e_o++; } } }\r
 365    */\r
 366   /**\r
 367    * Mark residues from start to end (inclusive) as deleted from the alignment,\r
 368    * and removes any insertions.\r
 369    * \r
 370    * @param start\r
 371    *          int\r
 372    * @param end\r
 373    *          int\r
 374    * @return deleted int - number of symbols marked as deleted\r
 375    */\r
 376   public int deleteRange(int start, int end)\r
 377   {\r
 378     int deleted = 0;\r
 379     if (length == 0)\r
 380     {\r
 381       // nothing to do here\r
 382       return deleted;\r
 383     }\r
 384     if (start < 0 || start > end)\r
 385     {\r
 386       throw new Error(\r
 387               "Implementation Error: deleteRange out of bounds: start must be non-negative and less than end.");\r
 388     }\r
 389     // find beginning\r
 390     int cursor = 0; // mark the position for the current operation being edited.\r
 391     int rlength = 1 + end - start; // number of positions to delete\r
 392     int oldlen = length;\r
 393     int o = 0;\r
 394     boolean editing = false;\r
 395     char[] oldops = operation;\r
 396     int[] oldrange = range;\r
 397     length = 0;\r
 398     operation = null;\r
 399     range = null;\r
 400     compact_operations();\r
 401     while (o < oldlen && cursor <= end && rlength > 0)\r
 402     {\r
 403       if (oldops[o] == D)\r
 404       {\r
 405         // absorbed into new deleted region.\r
 406         addDeleted(oldrange[o++]);\r
 407         continue;\r
 408       }\r
 409 \r
 410       int remain = oldrange[o]; // number of op characters left to edit\r
 411       if (!editing)\r
 412       {\r
 413         if ((cursor + remain) <= start)\r
 414         {\r
 415           addOperation(oldops[o], oldrange[o]);\r
 416           cursor += oldrange[o++];\r
 417           continue; // next operation\r
 418         }\r
 419         editing = true;\r
 420         // add operations before hidden region\r
 421         if (start - cursor > 0)\r
 422         {\r
 423           addOperation(oldops[o], start - cursor);\r
 424           remain -= start - cursor;\r
 425         }\r
 426       }\r
 427       // start inserting new ops\r
 428       if (o < oldlen && editing && rlength > 0 && remain > 0)\r
 429       {\r
 430         switch (oldops[o])\r
 431         {\r
 432         case M:\r
 433           if (rlength > remain)\r
 434           {\r
 435             addDeleted(remain);\r
 436             deleted += remain;\r
 437           }\r
 438           else\r
 439           {\r
 440             deleted += rlength;\r
 441             addDeleted(rlength);\r
 442             if (remain - rlength > 0)\r
 443             {\r
 444               this.addOperation(M, remain - rlength); // add remaining back.\r
 445             }\r
 446             rlength = 0;\r
 447             remain = 0;\r
 448           }\r
 449           break;\r
 450         case I:\r
 451           if (remain - rlength > 0)\r
 452           {\r
 453             // only remove some gaps\r
 454             addInsertion(remain - rlength);\r
 455             rlength = 0;\r
 456           }\r
 457           break;\r
 458         case D:\r
 459           throw new Error("Implementation error."); // do nothing;\r
 460         default:\r
 461           throw new Error("Implementation Error! Unknown operation '"\r
 462                   + oldops[o] + "'");\r
 463         }\r
 464         rlength -= remain;\r
 465         remain = oldrange[++o]; // number of op characters left to edit\r
 466       }\r
 467     }\r
 468     // add remaining\r
 469     while (o < oldlen)\r
 470     {\r
 471       addOperation(oldops[o], oldrange[o++]);\r
 472     }\r
 473     // if (cursor<(start+1)) {\r
 474     // ran out of ops - nothing to do here ?\r
 475     // addInsertion(start-cursor);\r
 476     // }\r
 477     return deleted;\r
 478   }\r
 479 \r
 480   /**\r
 481    * Deleted regions mean that there will be discontinuous sequence numbering in\r
 482    * the sequence returned by getSeq(char).\r
 483    * \r
 484    * @return true if there deletions\r
 485    */\r
 486   public boolean hasDeletedRegions()\r
 487   {\r
 488     for (int i = 0; i < length; i++)\r
 489     {\r
 490       if (operation[i] == D)\r
 491       {\r
 492         return true;\r
 493       }\r
 494     }\r
 495     return false;\r
 496   }\r
 497 \r
 498   /**\r
 499    * enumerate the ranges on seq that are marked as deleted in this cigar\r
 500    * \r
 501    * @return int[] { vis_start, sym_start, length }\r
 502    */\r
 503   public int[] getDeletedRegions()\r
 504   {\r
 505     if (length == 0)\r
 506     {\r
 507       return null;\r
 508     }\r
 509     Vector dr = new Vector();\r
 510     int cursor = 0, vcursor = 0;\r
 511     for (int i = 0; i < length; i++)\r
 512     {\r
 513       switch (operation[i])\r
 514       {\r
 515       case M:\r
 516         cursor += range[i];\r
 517       case I:\r
 518         vcursor += range[i];\r
 519         break;\r
 520       case D:\r
 521         dr.addElement(new int[]\r
 522         { vcursor, cursor, range[i] });\r
 523         cursor += range[i];\r
 524       }\r
 525     }\r
 526     if (dr.size() == 0)\r
 527     {\r
 528       return null;\r
 529     }\r
 530     int[] delregions = new int[dr.size() * 3];\r
 531     for (int i = 0, l = dr.size(); i < l; i++)\r
 532     {\r
 533       int[] reg = (int[]) dr.elementAt(i);\r
 534       delregions[i * 3] = reg[0];\r
 535       delregions[i * 3 + 1] = reg[1];\r
 536       delregions[i * 3 + 2] = reg[2];\r
 537     }\r
 538     return delregions;\r
 539   }\r
 540 \r
 541   /**\r
 542    * sum of ranges in cigar string\r
 543    * \r
 544    * @return int number of residues hidden, matched, or gaps inserted into\r
 545    *         sequence\r
 546    */\r
 547   public int getFullWidth()\r
 548   {\r
 549     int w = 0;\r
 550     if (range != null)\r
 551     {\r
 552       for (int i = 0; i < length; i++)\r
 553       {\r
 554         w += range[i];\r
 555       }\r
 556     }\r
 557     return w;\r
 558   }\r
 559 \r
 560   /**\r
 561    * Visible length of aligned sequence\r
 562    * \r
 563    * @return int length of including gaps and less hidden regions\r
 564    */\r
 565   public int getWidth()\r
 566   {\r
 567     int w = 0;\r
 568     if (range != null)\r
 569     {\r
 570       for (int i = 0; i < length; i++)\r
 571       {\r
 572         if (operation[i] == M || operation[i] == I)\r
 573         {\r
 574           w += range[i];\r
 575         }\r
 576       }\r
 577     }\r
 578     return w;\r
 579   }\r
 580 \r
 581   /**\r
 582    * mark a range of inserted residues\r
 583    * \r
 584    * @param range\r
 585    *          int\r
 586    */\r
 587   public void addInsertion(int range)\r
 588   {\r
 589     this.addOperation(I, range);\r
 590   }\r
 591 \r
 592   /**\r
 593    * mark the next range residues as hidden (not aligned) or deleted\r
 594    * \r
 595    * @param range\r
 596    *          int\r
 597    */\r
 598   public void addDeleted(int range)\r
 599   {\r
 600     this.addOperation(D, range);\r
 601   }\r
 602 \r
 603   /**\r
 604    * Modifies operation list to delete columns from start to end (inclusive)\r
 605    * editing will remove insertion operations, and convert matches to deletions\r
 606    * \r
 607    * @param start\r
 608    *          alignment column\r
 609    * @param end\r
 610    *          alignment column\r
 611    * @return boolean true if residues were marked as deleted. public boolean\r
 612    *         deleteRange(int start, int end) { boolean deleted = false; int op =\r
 613    *         0, prevop = -1, firstm = -1, lastm = -1, postop = -1; int width =\r
 614    *         0; // zero'th column if (length > 0) { // find operation bracketing\r
 615    *         start of the range do { if (operation[op] != D) { width +=\r
 616    *         range[prevop = op]; } op++; } while (op < length && width < start);\r
 617    *         } if (width < start) { // run off end - add more operations up to\r
 618    *         deletion. addInsertion(start - width); } else { // edit existing\r
 619    *         operations. op = prevop; width -= range[prevop]; int[] oldrange =\r
 620    *         range; char[] oldops = operation; range = new int[oldrange.length];\r
 621    *         operation = new char[oldops.length]; if (op < length) { do { if\r
 622    *         (operation[op] != D) { width += range[postop = op]; } op++; } while\r
 623    *         (op < length && width <= end); } } if (deleted == true) {\r
 624    *         addDeleted(end - start + 1); } return deleted; }\r
 625    */\r
 626   /**\r
 627    * Return an ENSEMBL style cigar string where D may indicates excluded parts\r
 628    * of seq\r
 629    * \r
 630    * @return String of form ([0-9]+[IMD])+\r
 631    */\r
 632   public String getCigarstring()\r
 633   {\r
 634     StringBuffer cigarString = new StringBuffer();\r
 635     for (int i = 0; i < length; i++)\r
 636     {\r
 637       cigarString.append("" + range[i]);\r
 638       cigarString.append(operation[i]);\r
 639     }\r
 640     return cigarString.toString();\r
 641   }\r
 642 }\r