src/jalview/datamodel/CigarBase.java

   1 /*\r
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.4)\r
   3  * Copyright (C) 2008 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
   4  * \r
   5  * This program is free software; you can redistribute it and/or\r
   6  * modify it under the terms of the GNU General Public License\r
   7  * as published by the Free Software Foundation; either version 2\r
   8  * of the License, or (at your option) any later version.\r
   9  * \r
  10  * This program is distributed in the hope that it will be useful,\r
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  13  * GNU General Public License for more details.\r
  14  * \r
  15  * You should have received a copy of the GNU General Public License\r
  16  * along with this program; if not, write to the Free Software\r
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA\r
  18  */\r
  19 package jalview.datamodel;\r
  20 \r
  21 import java.util.*;\r
  22 \r
  23 public abstract class CigarBase\r
  24 {\r
  25   /**\r
  26    * Base class for compact idiosyncratic representation of gaps and aligned\r
  27    * residues Regards to Tom Oldfield for his DynamicArray class. 17th July 2006\r
  28    * Not thread safe.\r
  29    */\r
  30   public CigarBase()\r
  31   {\r
  32     // nothing to be done (probably)\r
  33   }\r
  34 \r
  35   protected int length = 0;\r
  36 \r
  37   protected int _inc_length = 10; // extension range for addition of new\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\r
  68    *                null if 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")\r
 320    *  ( "1M10D5M".insert(0,M,3)->"1M10D8M") ( "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\r
 326    *                ops (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) {\r
 331    *                int cursor = -1; // mark the position for the current\r
 332    *                operation being edited. int o = 0; boolean last_d = false; //\r
 333    *                previous op was a 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\r
 338    *                more operations before pos if (pos-cursor>0)\r
 339    *                addInsertion(pos-cursor); // then just add the new\r
 340    *                operation. Regardless of what it is. addOperation(op,\r
 341    *                range); } else { int diff = pos - cursor;\r
 342    * \r
 343    * int e_length = length-o; // new edit operation array length. // diff<0 -\r
 344    * can only happen before first insertion or match. - affects op and all\r
 345    * following // dif==0 - only when at first position of existing op - //\r
 346    * diff>0 - must preserve some existing operations int[] e_range = new\r
 347    * int[e_length]; System.arraycopy(this.range, o, e_range, 0, e_length);\r
 348    * char[] e_op = new char[e_length]; System.arraycopy(this.operation, o, e_op,\r
 349    * 0, e_length); length = o; // can now use add_operation to extend list. int\r
 350    * e_o=0; // current operation being edited. switch (op) { case M: switch\r
 351    * (e_op[e_o]) { case M: if (last_d && diff <= 0) { // reduce D's, if possible\r
 352    * if (range<=this.range[o-1]) { this.range[o - 1] -= range; } else {\r
 353    * this.range[o-1]=0; } if (this.range[o-1]==0) o--; // lose this op. }\r
 354    * e_range[e_o] += range; // just add more matched residues break; case I: //\r
 355    * change from insertion to match if (last_d && diff<=0) { // reduce D's, if\r
 356    * possible if (range<=this.range[o-1]) { this.range[o - 1] -= range; } else {\r
 357    * this.range[o-1]=0; } if (this.range[o-1]==0) o--; // lose this op. }\r
 358    * e_range[e_o] break; default: throw new Inp }\r
 359    * \r
 360    * break; case I: break; case D: } break; default: throw new\r
 361    * Error("Implementation Error: Unknown operation in addOperation!"); } //\r
 362    * finally, add remaining ops. while (e_o<e_length) { addOperation(e_op[e_o],\r
 363    * e_range[e_o]); e_o++; } } }\r
 364    */\r
 365   /**\r
 366    * Mark residues from start to end (inclusive) as deleted from the alignment,\r
 367    * and removes any insertions.\r
 368    * \r
 369    * @param start\r
 370    *                int\r
 371    * @param end\r
 372    *                int\r
 373    * @return deleted int - number of symbols marked as deleted\r
 374    */\r
 375   public int deleteRange(int start, int end)\r
 376   {\r
 377     int deleted = 0;\r
 378     if (length == 0)\r
 379     {\r
 380       // nothing to do here\r
 381       return deleted;\r
 382     }\r
 383     if (start < 0 || start > end)\r
 384     {\r
 385       throw new Error(\r
 386               "Implementation Error: deleteRange out of bounds: start must be non-negative and less than end.");\r
 387     }\r
 388     // find beginning\r
 389     int cursor = 0; // mark the position for the current operation being edited.\r
 390     int rlength = 1 + end - start; // number of positions to delete\r
 391     int oldlen = length;\r
 392     int o = 0;\r
 393     boolean editing = false;\r
 394     char[] oldops = operation;\r
 395     int[] oldrange = range;\r
 396     length = 0;\r
 397     operation = null;\r
 398     range = null;\r
 399     compact_operations();\r
 400     while (o < oldlen && cursor <= end && rlength > 0)\r
 401     {\r
 402       if (oldops[o] == D)\r
 403       {\r
 404         // absorbed into new deleted region.\r
 405         addDeleted(oldrange[o++]);\r
 406         continue;\r
 407       }\r
 408 \r
 409       int remain = oldrange[o]; // number of op characters left to edit\r
 410       if (!editing)\r
 411       {\r
 412         if ((cursor + remain) <= start)\r
 413         {\r
 414           addOperation(oldops[o], oldrange[o]);\r
 415           cursor += oldrange[o++];\r
 416           continue; // next operation\r
 417         }\r
 418         editing = true;\r
 419         // add operations before hidden region\r
 420         if (start - cursor > 0)\r
 421         {\r
 422           addOperation(oldops[o], start - cursor);\r
 423           remain -= start - cursor;\r
 424         }\r
 425       }\r
 426       // start inserting new ops\r
 427       if (o < oldlen && editing && rlength > 0 && remain > 0)\r
 428       {\r
 429         switch (oldops[o])\r
 430         {\r
 431         case M:\r
 432           if (rlength > remain)\r
 433           {\r
 434             addDeleted(remain);\r
 435             deleted += remain;\r
 436           }\r
 437           else\r
 438           {\r
 439             deleted += rlength;\r
 440             addDeleted(rlength);\r
 441             if (remain - rlength > 0)\r
 442             {\r
 443               this.addOperation(M, remain - rlength); // add remaining back.\r
 444             }\r
 445             rlength = 0;\r
 446             remain = 0;\r
 447           }\r
 448           break;\r
 449         case I:\r
 450           if (remain - rlength > 0)\r
 451           {\r
 452             // only remove some gaps\r
 453             addInsertion(remain - rlength);\r
 454             rlength = 0;\r
 455           }\r
 456           break;\r
 457         case D:\r
 458           throw new Error("Implementation error."); // do nothing;\r
 459         default:\r
 460           throw new Error("Implementation Error! Unknown operation '"\r
 461                   + oldops[o] + "'");\r
 462         }\r
 463         rlength -= remain;\r
 464         remain = oldrange[++o]; // number of op characters left to edit\r
 465       }\r
 466     }\r
 467     // add remaining\r
 468     while (o < oldlen)\r
 469     {\r
 470       addOperation(oldops[o], oldrange[o++]);\r
 471     }\r
 472     // if (cursor<(start+1)) {\r
 473     // ran out of ops - nothing to do here ?\r
 474     // addInsertion(start-cursor);\r
 475     // }\r
 476     return deleted;\r
 477   }\r
 478 \r
 479   /**\r
 480    * Deleted regions mean that there will be discontinuous sequence numbering in\r
 481    * the sequence returned by getSeq(char).\r
 482    * \r
 483    * @return true if there deletions\r
 484    */\r
 485   public boolean hasDeletedRegions()\r
 486   {\r
 487     for (int i = 0; i < length; i++)\r
 488     {\r
 489       if (operation[i] == D)\r
 490       {\r
 491         return true;\r
 492       }\r
 493     }\r
 494     return false;\r
 495   }\r
 496 \r
 497   /**\r
 498    * enumerate the ranges on seq that are marked as deleted in this cigar\r
 499    * \r
 500    * @return int[] { vis_start, sym_start, length }\r
 501    */\r
 502   public int[] getDeletedRegions()\r
 503   {\r
 504     if (length == 0)\r
 505     {\r
 506       return null;\r
 507     }\r
 508     Vector dr = new Vector();\r
 509     int cursor = 0, vcursor = 0;\r
 510     for (int i = 0; i < length; i++)\r
 511     {\r
 512       switch (operation[i])\r
 513       {\r
 514       case M:\r
 515         cursor += range[i];\r
 516       case I:\r
 517         vcursor += range[i];\r
 518         break;\r
 519       case D:\r
 520         dr.addElement(new int[]\r
 521         { vcursor, cursor, range[i] });\r
 522         cursor += range[i];\r
 523       }\r
 524     }\r
 525     if (dr.size() == 0)\r
 526     {\r
 527       return null;\r
 528     }\r
 529     int[] delregions = new int[dr.size() * 3];\r
 530     for (int i = 0, l = dr.size(); i < l; i++)\r
 531     {\r
 532       int[] reg = (int[]) dr.elementAt(i);\r
 533       delregions[i * 3] = reg[0];\r
 534       delregions[i * 3 + 1] = reg[1];\r
 535       delregions[i * 3 + 2] = reg[2];\r
 536     }\r
 537     return delregions;\r
 538   }\r
 539 \r
 540   /**\r
 541    * sum of ranges in cigar string\r
 542    * \r
 543    * @return int number of residues hidden, matched, or gaps inserted into\r
 544    *         sequence\r
 545    */\r
 546   public int getFullWidth()\r
 547   {\r
 548     int w = 0;\r
 549     if (range != null)\r
 550     {\r
 551       for (int i = 0; i < length; i++)\r
 552       {\r
 553         w += range[i];\r
 554       }\r
 555     }\r
 556     return w;\r
 557   }\r
 558 \r
 559   /**\r
 560    * Visible length of aligned sequence\r
 561    * \r
 562    * @return int length of including gaps and less hidden regions\r
 563    */\r
 564   public int getWidth()\r
 565   {\r
 566     int w = 0;\r
 567     if (range != null)\r
 568     {\r
 569       for (int i = 0; i < length; i++)\r
 570       {\r
 571         if (operation[i] == M || operation[i] == I)\r
 572         {\r
 573           w += range[i];\r
 574         }\r
 575       }\r
 576     }\r
 577     return w;\r
 578   }\r
 579 \r
 580   /**\r
 581    * mark a range of inserted residues\r
 582    * \r
 583    * @param range\r
 584    *                int\r
 585    */\r
 586   public void addInsertion(int range)\r
 587   {\r
 588     this.addOperation(I, range);\r
 589   }\r
 590 \r
 591   /**\r
 592    * mark the next range residues as hidden (not aligned) or deleted\r
 593    * \r
 594    * @param range\r
 595    *                int\r
 596    */\r
 597   public void addDeleted(int range)\r
 598   {\r
 599     this.addOperation(D, range);\r
 600   }\r
 601 \r
 602   /**\r
 603    * Modifies operation list to delete columns from start to end (inclusive)\r
 604    * editing will remove insertion operations, and convert matches to deletions\r
 605    * \r
 606    * @param start\r
 607    *                alignment column\r
 608    * @param end\r
 609    *                alignment column\r
 610    * @return boolean true if residues were marked as deleted. public boolean\r
 611    *         deleteRange(int start, int end) { boolean deleted = false; int op =\r
 612    *         0, prevop = -1, firstm = -1, lastm = -1, postop = -1; int width =\r
 613    *         0; // zero'th column if (length > 0) { // find operation bracketing\r
 614    *         start of the range do { if (operation[op] != D) { width +=\r
 615    *         range[prevop = op]; } op++; } while (op < length && width < start); }\r
 616    *         if (width < start) { // run off end - add more operations up to\r
 617    *         deletion. addInsertion(start - width); } else { // edit existing\r
 618    *         operations. op = prevop; width -= range[prevop]; int[] oldrange =\r
 619    *         range; char[] oldops = operation; range = new int[oldrange.length];\r
 620    *         operation = new char[oldops.length]; if (op < length) { do { if\r
 621    *         (operation[op] != D) { width += range[postop = op]; } op++; } while\r
 622    *         (op < length && width <= end); } } if (deleted == true) {\r
 623    *         addDeleted(end - start + 1); } return deleted; }\r
 624    */\r
 625   /**\r
 626    * Return an ENSEMBL style cigar string where D may indicates excluded parts\r
 627    * of seq\r
 628    * \r
 629    * @return String of form ([0-9]+[IMD])+\r
 630    */\r
 631   public String getCigarstring()\r
 632   {\r
 633     StringBuffer cigarString = new StringBuffer();\r
 634     for (int i = 0; i < length; i++)\r
 635     {\r
 636       cigarString.append("" + range[i]);\r
 637       cigarString.append(operation[i]);\r
 638     }\r
 639     return cigarString.toString();\r
 640   }\r
 641 }\r