javadoc
[jalview.git] / src / jalview / datamodel / SeqCigar.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer
3  * Copyright (C) 2007 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
18  */
19 package jalview.datamodel;
20
21 import jalview.analysis.*;
22 import jalview.util.*;
23
24 public class SeqCigar
25     extends CigarSimple
26 {
27   /**
28    * start(inclusive) and end(exclusive) of subsequence on refseq
29    */
30   private int start, end;
31   private SequenceI refseq = null;
32   /**
33    * Reference dataset sequence for the cigar string
34    * @return SequenceI
35    */
36   public SequenceI getRefSeq()
37   {
38     return refseq;
39   }
40
41   /**
42    *
43    * @return int start index of cigar ops on refSeq
44    */
45   public int getStart()
46   {
47     return start;
48   }
49
50   /**
51    *
52    * @return int end index (exclusive) of cigar ops on refSeq
53    */
54   public int getEnd()
55   {
56     return end;
57   }
58
59   /**
60    * Returns sequence as a string with cigar operations applied to it
61    * @return String
62    */
63   public String getSequenceString(char GapChar)
64   {
65     return (length == 0) ? "" :
66         (String) getSequenceAndDeletions(refseq.getSequenceAsString(start, end),
67                                          GapChar)[0];
68   }
69
70   /**
71    * recreates a gapped and edited version of RefSeq or null for an empty cigar string
72    * @return SequenceI
73    */
74   public SequenceI getSeq(char GapChar)
75   {
76     Sequence seq;
77     if (refseq == null || length == 0)
78     {
79       return null;
80     }
81     Object[] edit_result = getSequenceAndDeletions(refseq.getSequenceAsString(
82         start, end),
83         GapChar);
84     if (edit_result == null)
85     {
86       throw new Error(
87           "Implementation Error - unexpected null from getSequenceAndDeletions");
88     }
89     int bounds[] = (int[]) edit_result[1];
90     seq = new Sequence(refseq.getName(), (String) edit_result[0],
91                        refseq.getStart() + start + bounds[0],
92                        refseq.getStart() + start +
93                        ( (bounds[2] == 0) ? -1 : bounds[2]));
94     // seq.checkValidRange(); probably not needed
95     seq.setDatasetSequence(refseq);
96     seq.setDescription(refseq.getDescription());
97     return seq;
98   }
99
100   /*
101      We don't allow this - refseq is given at construction time only
102    public void setSeq(SequenceI seq) {
103     this.seq = seq;
104      }
105    */
106   /**
107    * internal constructor - sets seq to a gapless sequence derived from seq
108    * and prepends any 'D' operations needed to get to the first residue of seq.
109    * @param seq SequenceI
110    * @param initialDeletion true to mark initial dataset sequence residues as deleted in subsequence
111    * @param _s index of first position in seq
112    * @param _e index after last position in (possibly gapped) seq
113    * @return true if gaps are present in seq
114    */
115   private boolean _setSeq(SequenceI seq, boolean initialDeletion, int _s,
116                           int _e)
117   {
118     boolean hasgaps = false;
119     if (seq == null)
120     {
121       throw new Error("Implementation Error - _setSeq(null,...)");
122     }
123     if (_s < 0)
124     {
125       throw new Error("Implementation Error: _s=" + _s);
126     }
127     String seq_string = seq.getSequenceAsString();
128     if (_e == 0 || _e < _s || _e > seq_string.length())
129     {
130       _e = seq_string.length();
131     }
132     // resolve start and end positions relative to ungapped reference sequence
133     start = seq.findPosition(_s) - seq.getStart();
134     end = seq.findPosition(_e) - seq.getStart();
135     int l_ungapped = end - start;
136     // Find correct sequence to reference and correct start and end - if necessary
137     SequenceI ds = seq.getDatasetSequence();
138     if (ds == null)
139     {
140       // make a new dataset sequence
141       String ungapped = AlignSeq.extractGaps(jalview.util.Comparison.GapChars,
142                                              new String(seq_string));
143       l_ungapped = ungapped.length();
144       // check that we haven't just duplicated an ungapped sequence.
145       if (l_ungapped == seq.getLength())
146       {
147         ds = seq;
148       }
149       else
150       {
151         ds = new Sequence(seq.getName(), ungapped,
152                           seq.getStart(),
153                           seq.getStart() + ungapped.length() - 1);
154         // JBPNote: this would be consistent but may not be useful
155         //        seq.setDatasetSequence(ds);
156       }
157     }
158     // add in offset between seq and the dataset sequence
159     if (ds.getStart() < seq.getStart())
160     {
161       int offset = seq.getStart() - ds.getStart();
162       if (initialDeletion)
163       {
164         // absolute cigar string
165         addDeleted(_s + offset);
166         start = 0;
167         end += offset;
168       }
169       else
170       {
171         // normal behaviour - just mark start and end subsequence
172         start += offset;
173         end += offset;
174
175       }
176
177     }
178
179     // any gaps to process ?
180     if (l_ungapped != (_e - _s))
181     {
182       hasgaps = true;
183     }
184
185     this.refseq = ds;
186
187     // Check  offsets
188     if (end > ds.getLength())
189     {
190       throw new Error("SeqCigar: Possible implementation error: sequence is longer than dataset sequence");
191 //      end = ds.getLength();
192     }
193
194     return hasgaps;
195   }
196
197   /**
198    * directly initialise a cigar object with a sequence of range, operation pairs and a sequence to apply it to.
199    * operation and range should be relative to the seq.getStart()'th residue of the dataset seq resolved from seq.
200    * @param seq SequenceI
201    * @param operation char[]
202    * @param range int[]
203    */
204   public SeqCigar(SequenceI seq, char operation[], int range[])
205   {
206     super();
207     if (seq == null)
208     {
209       throw new Error("Implementation Bug. Null seq !");
210     }
211     if (operation.length != range.length)
212     {
213       throw new Error("Implementation Bug. Cigar Operation list!= range list");
214     }
215
216     if (operation != null)
217     {
218       this.operation = new char[operation.length + _inc_length];
219       this.range = new int[operation.length + _inc_length];
220
221       if (_setSeq(seq, false, 0, 0))
222       {
223         throw new Error("NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
224       }
225       for (int i = this.length, j = 0; j < operation.length; i++, j++)
226       {
227         char op = operation[j];
228         if (op != M && op != I && op != D)
229         {
230           throw new Error(
231               "Implementation Bug. Cigar Operation '" + j + "' '" + op +
232               "' not one of '" + M + "', '" + I + "', or '" + D + "'.");
233         }
234         this.operation[i] = op;
235         this.range[i] = range[j];
236       }
237       this.length += operation.length;
238     }
239     else
240     {
241       this.operation = null;
242       this.range = null;
243       this.length = 0;
244       if (_setSeq(seq, false, 0, 0))
245       {
246         throw new Error("NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
247       }
248     }
249   }
250
251   /**
252    * add range matched residues to cigar string
253    * @param range int
254    */
255   public void addMatch(int range)
256   {
257     this.addOperation(M, range);
258   }
259
260   /**
261    * Adds
262    * insertion and match operations based on seq to the cigar up to
263    * the endpos column of seq.
264    *
265    * @param cigar CigarBase
266    * @param seq SequenceI
267    * @param startpos int
268    * @param endpos int
269    * @param initialDeletions if true then initial deletions will be added from start of seq to startpos
270    */
271   protected static void addSequenceOps(CigarBase cigar, SequenceI seq,
272                                        int startpos, int endpos,
273                                        boolean initialDeletions)
274   {
275     char op = '\0';
276     int range = 0;
277     int p = 0, res = seq.getLength();
278
279     if (!initialDeletions)
280     {
281       p = startpos;
282     }
283
284     while (p <= endpos)
285     {
286       boolean isGap = (p < res) ? jalview.util.Comparison.isGap(seq.getCharAt(p)) : true;
287       if ( (startpos <= p) && (p <= endpos))
288       {
289         if (isGap)
290         {
291           if (range > 0 && op != I)
292           {
293             cigar.addOperation(op, range);
294             range = 0;
295           }
296           op = I;
297           range++;
298         }
299         else
300         {
301           if (range > 0 && op != M)
302           {
303             cigar.addOperation(op, range);
304             range = 0;
305           }
306           op = M;
307           range++;
308         }
309       }
310       else
311       {
312         if (!isGap)
313         {
314           if (range > 0 && op != D)
315           {
316             cigar.addOperation(op, range);
317             range = 0;
318           }
319           op = D;
320           range++;
321         }
322         else
323         {
324           // do nothing - insertions are not made in flanking regions
325         }
326       }
327       p++;
328     }
329     if (range > 0)
330     {
331       cigar.addOperation(op, range);
332     }
333   }
334
335   /**
336    * create a cigar string for given sequence
337    * @param seq SequenceI
338    */
339   public SeqCigar(SequenceI seq)
340   {
341     super();
342     if (seq == null)
343     {
344       throw new Error("Implementation error for new Cigar(SequenceI)");
345     }
346     _setSeq(seq, false, 0, 0);
347     // there is still work to do
348     addSequenceOps(this, seq, 0, seq.getLength() - 1, false);
349   }
350
351   /**
352    * Create Cigar from a range of gaps and residues on a sequence object
353    * @param seq SequenceI
354    * @param start int - first column in range
355    * @param end int - last column in range
356    */
357   public SeqCigar(SequenceI seq, int start, int end)
358   {
359     super();
360     if (seq == null)
361     {
362       throw new Error("Implementation error for new Cigar(SequenceI)");
363     }
364     _setSeq(seq, false, start, end + 1);
365     // there is still work to do
366     addSequenceOps(this, seq, start, end, false);
367   }
368
369   /**
370    * Create a cigar object from a cigar string like '[<I|D|M><range>]+'
371    * Will fail if the given seq already contains gaps (JBPNote: future implementation will fix)
372    * @param seq SequenceI object resolvable to a dataset sequence
373    * @param cigarString String
374    * @return Cigar
375    */
376   public static SeqCigar parseCigar(SequenceI seq, String cigarString)
377       throws Exception
378   {
379     Object[] opsandrange = parseCigarString(cigarString);
380     return new SeqCigar(seq, (char[]) opsandrange[0], (int[]) opsandrange[1]);
381   }
382
383   /**
384    * createAlignment
385    *
386    * @param alseqs SeqCigar[]
387    * @param gapCharacter char
388    * @return SequenceI[]
389    */
390   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
391       char gapCharacter, ColumnSelection colsel, int[] segments)
392   {
393     SequenceI[] seqs = new SequenceI[alseqs.length];
394     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
395     String[] alseqs_string = new String[alseqs.length];
396     Object[] gs_regions = new Object[alseqs.length];
397     for (int i = 0; i < alseqs.length; i++)
398     {
399       alseqs_string[i] = alseqs[i].getRefSeq().
400           getSequenceAsString(alseqs[i].start, alseqs[i].end);
401       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
402           gapCharacter); // gapped sequence, {start, start col, end. endcol}, hidden regions {{start, end, col}})
403       if (gs_regions[i] == null)
404       {
405         throw new Error("Implementation error: " + i +
406                         "'th sequence Cigar has no operations.");
407       }
408       g_seqs[i] = new StringBuffer( (String) ( (Object[]) gs_regions[i])[0]); // the visible gapped sequence
409     }
410     // Now account for insertions. (well - deletions)
411     // this is complicated because we must keep track of shifted positions in each sequence
412     ShiftList shifts = new ShiftList();
413     for (int i = 0; i < alseqs.length; i++)
414     {
415       Object[] gs_region = ( (Object[]) ( (Object[]) gs_regions[i])[2]);
416       if (gs_region != null)
417
418       {
419         for (int hr = 0; hr < gs_region.length; hr++)
420         {
421           int[] region = (int[]) gs_region[hr];
422           char[] insert = new char[region[1] - region[0] + 1];
423           for (int s = 0; s < insert.length; s++)
424           {
425             insert[s] = gapCharacter;
426           }
427           int inspos = shifts.shift(region[2]); // resolve insertion position in current alignment frame of reference
428           for (int s = 0; s < alseqs.length; s++)
429           {
430             if (s != i)
431             {
432               if (g_seqs[s].length() <= inspos)
433               {
434                 // prefix insertion with more gaps.
435                 for (int l = inspos - g_seqs[s].length(); l > 0; l--)
436                 {
437                   g_seqs[s].append(gapCharacter); // to debug - use a diffferent gap character here
438                 }
439               }
440               g_seqs[s].insert(inspos, insert);
441             }
442             else
443             {
444               g_seqs[s].insert(inspos,
445                                alseqs_string[i].substring(region[0],
446                   region[1] + 1));
447             }
448           }
449           shifts.addShift(region[2], insert.length); // update shift in alignment frame of reference
450           if (segments == null)
451           {
452             // add a hidden column for this deletion
453             colsel.hideColumns(inspos, inspos + insert.length - 1);
454           }
455         }
456       }
457     }
458     for (int i = 0; i < alseqs.length; i++)
459     {
460       int[] bounds = ( (int[]) ( (Object[]) gs_regions[i])[1]);
461       SequenceI ref = alseqs[i].getRefSeq();
462       seqs[i] = new Sequence(ref.getName(), g_seqs[i].toString(),
463                              ref.getStart() + alseqs[i].start + bounds[0],
464                              ref.getStart() + alseqs[i].start +
465                              (bounds[2] == 0 ? -1 : bounds[2]));
466       seqs[i].setDatasetSequence(ref);
467       seqs[i].setDescription(ref.getDescription());
468     }
469     if (segments != null)
470     {
471       for (int i = 0; i < segments.length; i += 3)
472       {
473         //int start=shifts.shift(segments[i]-1)+1;
474         //int end=shifts.shift(segments[i]+segments[i+1]-1)-1;
475         colsel.hideColumns(segments[i + 1],
476                            segments[i + 1] + segments[i + 2] - 1);
477       }
478     }
479     return seqs;
480   }
481
482   /**
483    * non rigorous testing
484    */
485   /**
486    *
487    * @param seq Sequence
488    * @param ex_cs_gapped String
489    * @return String
490    */
491   public static String testCigar_string(Sequence seq, String ex_cs_gapped)
492   {
493     SeqCigar c_sgapped = new SeqCigar(seq);
494     String cs_gapped = c_sgapped.getCigarstring();
495     if (!cs_gapped.equals(ex_cs_gapped))
496     {
497       System.err.println("Failed getCigarstring: incorect string '" + cs_gapped +
498                          "' != " + ex_cs_gapped);
499     }
500     return cs_gapped;
501   }
502
503   public static boolean testSeqRecovery(SeqCigar gen_sgapped,
504                                         SequenceI s_gapped)
505   {
506     // this is non-rigorous - start and end  recovery is not tested.
507     SequenceI gen_sgapped_s = gen_sgapped.getSeq('-');
508     if (!gen_sgapped_s.getSequence().equals(s_gapped.getSequence()))
509     {
510       System.err.println("Couldn't reconstruct sequence.\n" +
511                          gen_sgapped_s.getSequenceAsString() + "\n" +
512                          s_gapped.getSequenceAsString());
513       return false;
514     }
515     return true;
516   }
517
518   public static void main(String argv[])
519       throws Exception
520   {
521     String o_seq;
522     Sequence s = new Sequence("MySeq",
523                               o_seq =
524                               "asdfktryasdtqwrtsaslldddptyipqqwaslchvhttt",
525                               39, 80);
526     String orig_gapped;
527     Sequence s_gapped = new Sequence("MySeq",
528                                      orig_gapped =
529         "----asdf------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhttt",
530                                      39, 80);
531     String ex_cs_gapped = "4I4M6I6M3I11M4I12M4I9M";
532     s_gapped.setDatasetSequence(s);
533     String sub_gapped_s;
534     Sequence s_subsequence_gapped = new Sequence("MySeq",
535                                                  sub_gapped_s =
536         "------ktryas---dtqwrtsasll----dddptyipqqwa----slchvh",
537                                                  43, 77);
538
539     s_subsequence_gapped.setDatasetSequence(s);
540     SeqCigar c_null = new SeqCigar(s);
541     String cs_null = c_null.getCigarstring();
542     if (!cs_null.equals("42M"))
543     {
544       System.err.println(
545           "Failed to recover ungapped sequence cigar operations:" +
546           ( (cs_null == "") ? "empty string" : cs_null));
547     }
548     testCigar_string(s_gapped, ex_cs_gapped);
549     SeqCigar gen_sgapped = SeqCigar.parseCigar(s, ex_cs_gapped);
550     if (!gen_sgapped.getCigarstring().equals(ex_cs_gapped))
551     {
552       System.err.println("Failed parseCigar(" + ex_cs_gapped +
553                          ")->getCigarString()->'" + gen_sgapped.getCigarstring() +
554                          "'");
555     }
556     testSeqRecovery(gen_sgapped, s_gapped);
557     // Test dataset resolution
558     SeqCigar sub_gapped = new SeqCigar(s_subsequence_gapped);
559     if (!testSeqRecovery(sub_gapped, s_subsequence_gapped))
560     {
561       System.err.println("Failed recovery for subsequence of dataset sequence");
562     }
563     // width functions
564     if (sub_gapped.getWidth() != sub_gapped_s.length())
565     {
566       System.err.println("Failed getWidth()");
567     }
568
569     sub_gapped.getFullWidth();
570     if (sub_gapped.hasDeletedRegions())
571     {
572       System.err.println("hasDeletedRegions is incorrect.");
573     }
574     // Test start-end region SeqCigar
575     SeqCigar sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
576     if (sub_se_gp.getWidth() != 41)
577     {
578       System.err.println(
579           "SeqCigar(seq, start, end) not properly clipped alignsequence.");
580     }
581     System.out.println("Original sequence align:\n" + sub_gapped_s +
582                        "\nReconstructed window from 8 to 48\n"
583                        + "XXXXXXXX" + sub_se_gp.getSequenceString('-') + "..."
584                        + "\nCigar String:" + sub_se_gp.getCigarstring() + "\n"
585         );
586     SequenceI ssgp = sub_se_gp.getSeq('-');
587     System.out.println("\t " + ssgp.getSequenceAsString());
588     for (int r = 0; r < 10; r++)
589     {
590       sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
591       int sl = sub_se_gp.getWidth();
592       int st = sl - 1 - r;
593       for (int rs = 0; rs < 10; rs++)
594       {
595         int e = st + rs;
596         sub_se_gp.deleteRange(st, e);
597         String ssgapedseq = sub_se_gp.getSeq('-').getSequenceAsString();
598         System.out.println(st + "," + e + "\t:" + ssgapedseq);
599         st -= 3;
600       }
601     }
602     {
603       SeqCigar[] set = new SeqCigar[]
604           {
605           new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
606           new SeqCigar(s_gapped)};
607       Alignment al = new Alignment(set);
608       for (int i = 0; i < al.getHeight(); i++)
609       {
610         System.out.println("" + al.getSequenceAt(i).getName() + "\t" +
611                            al.getSequenceAt(i).getStart() + "\t" +
612                            al.getSequenceAt(i).getEnd() + "\t" +
613                            al.getSequenceAt(i).getSequenceAsString());
614       }
615     }
616     {
617       System.out.println("Gapped.");
618       SeqCigar[] set = new SeqCigar[]
619           {
620           new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
621           new SeqCigar(s_gapped)};
622       set[0].deleteRange(20, 25);
623       Alignment al = new Alignment(set);
624       for (int i = 0; i < al.getHeight(); i++)
625       {
626         System.out.println("" + al.getSequenceAt(i).getName() + "\t" +
627                            al.getSequenceAt(i).getStart() + "\t" +
628                            al.getSequenceAt(i).getEnd() + "\t" +
629                            al.getSequenceAt(i).getSequenceAsString());
630       }
631     }
632 //    if (!ssgapedseq.equals("ryas---dtqqwa----slchvh"))
633 //      System.err.println("Subseqgaped\n------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhryas---dtqwrtsasll--qwa----slchvh\n"+ssgapedseq+"\n"+sub_se_gp.getCigarstring());
634   }
635
636 }