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