JAL-1632 add score model params to PCAModel and PCA constructors
[jalview.git] / src / jalview / datamodel / SeqCigar.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.datamodel;
22
23 import jalview.analysis.AlignSeq;
24 import jalview.analysis.SeqsetUtils;
25 import jalview.util.MessageManager;
26 import jalview.util.ShiftList;
27
28 import java.util.Enumeration;
29 import java.util.Hashtable;
30
31 public class SeqCigar extends CigarSimple
32 {
33   /**
34    * start(inclusive) and end(exclusive) of subsequence on refseq
35    */
36   private int start, end;
37
38   private SequenceI refseq = null;
39
40   private Hashtable seqProps;
41
42   /**
43    * Reference dataset sequence for the cigar string
44    * 
45    * @return SequenceI
46    */
47   public SequenceI getRefSeq()
48   {
49     return refseq;
50   }
51
52   /**
53    * 
54    * @return int start index of cigar ops on refSeq
55    */
56   public int getStart()
57   {
58     return start;
59   }
60
61   /**
62    * 
63    * @return int end index (exclusive) of cigar ops on refSeq
64    */
65   public int getEnd()
66   {
67     return end;
68   }
69
70   /**
71    * 
72    * @param column
73    * @return position in sequence for column (or -1 if no match state exists)
74    */
75   public int findPosition(int column)
76   {
77     int w = 0, ew, p = refseq.findPosition(start);
78     if (column < 0)
79     {
80       return -1;
81     }
82     if (range != null)
83     {
84       for (int i = 0; i < length; i++)
85       {
86         if (operation[i] == M || operation[i] == D)
87         {
88           p += range[i];
89         }
90         if (operation[i] == M || operation[i] == I)
91         {
92           ew = w + range[i];
93           if (column < ew)
94           {
95             if (operation[i] == I)
96             {
97               return -1;
98             }
99             return p - (ew - column);
100           }
101           w = ew;
102         }
103       }
104     }
105     return -1;
106   }
107
108   /**
109    * Returns sequence as a string with cigar operations applied to it
110    * 
111    * @return String
112    */
113   @Override
114   public String getSequenceString(char GapChar)
115   {
116     return (length == 0) ? "" : (String) getSequenceAndDeletions(
117             refseq.getSequenceAsString(start, end), GapChar)[0];
118   }
119
120   /**
121    * recreates a gapped and edited version of RefSeq or null for an empty cigar
122    * string
123    * 
124    * @return SequenceI
125    */
126   public SequenceI getSeq(char GapChar)
127   {
128     Sequence seq;
129     if (refseq == null || length == 0)
130     {
131       return null;
132     }
133     Object[] edit_result = getSequenceAndDeletions(
134             refseq.getSequenceAsString(start, end), GapChar);
135     if (edit_result == null)
136     {
137       throw new Error(
138               MessageManager
139                       .getString("error.implementation_error_unexpected_null_from_get_sequence_and_deletions"));
140     }
141     int bounds[] = (int[]) edit_result[1];
142     seq = new Sequence(refseq.getName(), (String) edit_result[0],
143             refseq.getStart() + start + bounds[0], refseq.getStart()
144                     + start + ((bounds[2] == 0) ? -1 : bounds[2]));
145     seq.setDescription(refseq.getDescription());
146     int sstart = seq.getStart(), send = seq.getEnd();
147     // seq.checkValidRange(); probably not needed
148     // recover local properties if present
149     if (seqProps != null)
150     {
151       // this recovers dataset sequence reference as well as local features,
152       // names, start/end settings.
153       SeqsetUtils.SeqCharacterUnhash(seq, seqProps);
154     }
155     // ensure dataset sequence is up to date from local reference
156     seq.setDatasetSequence(refseq);
157     seq.setStart(sstart);
158     seq.setEnd(send);
159     return seq;
160   }
161
162   /*
163    * We don't allow this - refseq is given at construction time only public void
164    * setSeq(SequenceI seq) { this.seq = seq; }
165    */
166   /**
167    * internal constructor - sets seq to a gapless sequence derived from seq and
168    * prepends any 'D' operations needed to get to the first residue of seq.
169    * 
170    * @param seq
171    *          SequenceI
172    * @param initialDeletion
173    *          true to mark initial dataset sequence residues as deleted in
174    *          subsequence
175    * @param _s
176    *          index of first position in seq
177    * @param _e
178    *          index after last position in (possibly gapped) seq
179    * @return true if gaps are present in seq
180    */
181   private boolean _setSeq(SequenceI seq, boolean initialDeletion, int _s,
182           int _e)
183   {
184     boolean hasgaps = false;
185     if (seq == null)
186     {
187       throw new Error(
188               MessageManager
189                       .getString("error.implementation_error_set_seq_null"));
190     }
191     if (_s < 0)
192     {
193       throw new Error(MessageManager.formatMessage(
194               "error.implementation_error_s", new String[] { Integer
195                       .valueOf(_s).toString() }));
196     }
197     String seq_string = seq.getSequenceAsString();
198     if (_e == 0 || _e < _s || _e > seq_string.length())
199     {
200       _e = seq_string.length();
201     }
202     // resolve start and end positions relative to ungapped reference sequence
203     start = seq.findPosition(_s) - seq.getStart();
204     end = seq.findPosition(_e) - seq.getStart();
205     int l_ungapped = end - start;
206     // Find correct sequence to reference and correct start and end - if
207     // necessary
208     SequenceI ds = seq.getDatasetSequence();
209     if (ds == null)
210     {
211       // make a new dataset sequence
212       String ungapped = AlignSeq.extractGaps(
213               jalview.util.Comparison.GapChars, new String(seq_string));
214       l_ungapped = ungapped.length();
215       // check that we haven't just duplicated an ungapped sequence.
216       if (l_ungapped == seq.getLength())
217       {
218         ds = seq;
219       }
220       else
221       {
222         ds = new Sequence(seq.getName(), ungapped, seq.getStart(),
223                 seq.getStart() + ungapped.length() - 1);
224         // JBPNote: this would be consistent but may not be useful
225         // seq.setDatasetSequence(ds);
226       }
227     }
228     // add in offset between seq and the dataset sequence
229     if (ds.getStart() < seq.getStart())
230     {
231       int offset = seq.getStart() - ds.getStart();
232       if (initialDeletion)
233       {
234         // absolute cigar string
235         addDeleted(_s + offset);
236         start = 0;
237         end += offset;
238       }
239       else
240       {
241         // normal behaviour - just mark start and end subsequence
242         start += offset;
243         end += offset;
244
245       }
246
247     }
248
249     // any gaps to process ?
250     if (l_ungapped != (_e - _s))
251     {
252       hasgaps = true;
253     }
254
255     refseq = ds;
256     // copy over local properties for the sequence instance of the refseq
257     seqProps = SeqsetUtils.SeqCharacterHash(seq);
258     // Check offsets
259     if (end > ds.getLength())
260     {
261       throw new Error(
262               MessageManager
263                       .getString("error.implementation_error_seqcigar_possible"));
264       // end = ds.getLength();
265     }
266
267     return hasgaps;
268   }
269
270   /**
271    * directly initialise a cigar object with a sequence of range, operation
272    * pairs and a sequence to apply it to. operation and range should be relative
273    * to the seq.getStart()'th residue of the dataset seq resolved from seq.
274    * 
275    * @param seq
276    *          SequenceI
277    * @param operation
278    *          char[]
279    * @param range
280    *          int[]
281    */
282   public SeqCigar(SequenceI seq, char operation[], int range[])
283   {
284     super();
285     if (seq == null)
286     {
287       throw new Error(
288               MessageManager.getString("error.implmentation_bug_seq_null"));
289     }
290     if (operation.length != range.length)
291     {
292       throw new Error(
293               MessageManager
294                       .getString("error.implementation_bug_cigar_operation_list_range_list"));
295     }
296
297     if (operation != null)
298     {
299       this.operation = new char[operation.length + _inc_length];
300       this.range = new int[operation.length + _inc_length];
301
302       if (_setSeq(seq, false, 0, 0))
303       {
304         throw new Error(
305                 MessageManager
306                         .getString("error.not_yet_implemented_cigar_object_from_cigar_string"));
307       }
308       for (int i = this.length, j = 0; j < operation.length; i++, j++)
309       {
310         char op = operation[j];
311         if (op != M && op != I && op != D)
312         {
313           throw new Error(MessageManager.formatMessage(
314                   "error.implementation_bug_cigar_operation", new String[] {
315                       Integer.valueOf(j).toString(),
316                       Integer.valueOf(op).toString(),
317                       Integer.valueOf(M).toString(),
318                       Integer.valueOf(I).toString(),
319                       Integer.valueOf(D).toString() }));
320         }
321         this.operation[i] = op;
322         this.range[i] = range[j];
323       }
324       this.length += operation.length;
325     }
326     else
327     {
328       this.operation = null;
329       this.range = null;
330       this.length = 0;
331       if (_setSeq(seq, false, 0, 0))
332       {
333         throw new Error(
334                 MessageManager
335                         .getString("error.not_yet_implemented_cigar_object_from_cigar_string"));
336       }
337     }
338   }
339
340   /**
341    * add range matched residues to cigar string
342    * 
343    * @param range
344    *          int
345    */
346   public void addMatch(int range)
347   {
348     this.addOperation(M, range);
349   }
350
351   /**
352    * Adds insertion and match operations based on seq to the cigar up to the
353    * endpos column of seq.
354    * 
355    * @param cigar
356    *          CigarBase
357    * @param seq
358    *          SequenceI
359    * @param startpos
360    *          int
361    * @param endpos
362    *          int
363    * @param initialDeletions
364    *          if true then initial deletions will be added from start of seq to
365    *          startpos
366    */
367   protected static void addSequenceOps(CigarBase cigar, SequenceI seq,
368           int startpos, int endpos, boolean initialDeletions)
369   {
370     char op = '\0';
371     int range = 0;
372     int p = 0, res = seq.getLength();
373
374     if (!initialDeletions)
375     {
376       p = startpos;
377     }
378
379     while (p <= endpos)
380     {
381       boolean isGap = (p < res) ? jalview.util.Comparison.isGap(seq
382               .getCharAt(p)) : true;
383       if ((startpos <= p) && (p <= endpos))
384       {
385         if (isGap)
386         {
387           if (range > 0 && op != I)
388           {
389             cigar.addOperation(op, range);
390             range = 0;
391           }
392           op = I;
393           range++;
394         }
395         else
396         {
397           if (range > 0 && op != M)
398           {
399             cigar.addOperation(op, range);
400             range = 0;
401           }
402           op = M;
403           range++;
404         }
405       }
406       else
407       {
408         if (!isGap)
409         {
410           if (range > 0 && op != D)
411           {
412             cigar.addOperation(op, range);
413             range = 0;
414           }
415           op = D;
416           range++;
417         }
418         else
419         {
420           // do nothing - insertions are not made in flanking regions
421         }
422       }
423       p++;
424     }
425     if (range > 0)
426     {
427       cigar.addOperation(op, range);
428     }
429   }
430
431   /**
432    * create a cigar string for given sequence
433    * 
434    * @param seq
435    *          SequenceI
436    */
437   public SeqCigar(SequenceI seq)
438   {
439     super();
440     if (seq == null)
441     {
442       throw new Error(
443               MessageManager
444                       .getString("error.implementation_error_for_new_cigar"));
445     }
446     _setSeq(seq, false, 0, 0);
447     // there is still work to do
448     addSequenceOps(this, seq, 0, seq.getLength() - 1, false);
449   }
450
451   /**
452    * Create Cigar from a range of gaps and residues on a sequence object
453    * 
454    * @param seq
455    *          SequenceI
456    * @param start
457    *          int - first column in range
458    * @param end
459    *          int - last column in range
460    */
461   public SeqCigar(SequenceI seq, int start, int end)
462   {
463     super();
464     if (seq == null)
465     {
466       throw new Error(
467               MessageManager
468                       .getString("error.implementation_error_for_new_cigar"));
469     }
470     _setSeq(seq, false, start, end + 1);
471     // there is still work to do
472     addSequenceOps(this, seq, start, end, false);
473   }
474
475   /**
476    * Create a cigar object from a cigar string like '[<I|D|M><range>]+' Will
477    * fail if the given seq already contains gaps (JBPNote: future implementation
478    * will fix)
479    * 
480    * @param seq
481    *          SequenceI object resolvable to a dataset sequence
482    * @param cigarString
483    *          String
484    * @return Cigar
485    */
486   public static SeqCigar parseCigar(SequenceI seq, String cigarString)
487           throws Exception
488   {
489     Object[] opsandrange = parseCigarString(cigarString);
490     return new SeqCigar(seq, (char[]) opsandrange[0],
491             (int[]) opsandrange[1]);
492   }
493
494   /**
495    * create an alignment from the given array of cigar sequences and gap
496    * character, and marking the given segments as visible in the given
497    * columselection.
498    * 
499    * @param alseqs
500    * @param gapCharacter
501    * @param colsel
502    *          - columnSelection where hidden regions are marked
503    * @param segments
504    *          - visible regions of alignment
505    * @return SequenceI[]
506    */
507   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
508           char gapCharacter, ColumnSelection colsel, int[] segments)
509   {
510     SequenceI[] seqs = new SequenceI[alseqs.length];
511     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
512     String[] alseqs_string = new String[alseqs.length];
513     Object[] gs_regions = new Object[alseqs.length];
514     for (int i = 0; i < alseqs.length; i++)
515     {
516       alseqs_string[i] = alseqs[i].getRefSeq().getSequenceAsString(
517               alseqs[i].start, alseqs[i].end);
518       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
519               gapCharacter); // gapped sequence, {start, start col, end.
520       // endcol}, hidden regions {{start, end, col}})
521       if (gs_regions[i] == null)
522       {
523         throw new Error(MessageManager.formatMessage(
524                 "error.implementation_error_cigar_seq_no_operations",
525                 new String[] { Integer.valueOf(i).toString() }));
526       }
527       g_seqs[i] = new StringBuffer((String) ((Object[]) gs_regions[i])[0]); // the
528       // visible
529       // gapped
530       // sequence
531     }
532     // Now account for insertions. (well - deletions)
533     // this is complicated because we must keep track of shifted positions in
534     // each sequence
535     ShiftList shifts = new ShiftList();
536     for (int i = 0; i < alseqs.length; i++)
537     {
538       Object[] gs_region = ((Object[]) ((Object[]) gs_regions[i])[2]);
539       if (gs_region != null)
540
541       {
542         for (int hr = 0; hr < gs_region.length; hr++)
543         {
544           int[] region = (int[]) gs_region[hr];
545           char[] insert = new char[region[1] - region[0] + 1];
546           for (int s = 0; s < insert.length; s++)
547           {
548             insert[s] = gapCharacter;
549           }
550           int inspos = shifts.shift(region[2]); // resolve insertion position in
551           // current alignment frame of
552           // reference
553           for (int s = 0; s < alseqs.length; s++)
554           {
555             if (s != i)
556             {
557               if (g_seqs[s].length() <= inspos)
558               {
559                 // prefix insertion with more gaps.
560                 for (int l = inspos - g_seqs[s].length(); l > 0; l--)
561                 {
562                   g_seqs[s].append(gapCharacter); // to debug - use a diffferent
563                   // gap character here
564                 }
565               }
566               g_seqs[s].insert(inspos, insert);
567             }
568             else
569             {
570               g_seqs[s].insert(inspos,
571                       alseqs_string[i].substring(region[0], region[1] + 1));
572             }
573           }
574           shifts.addShift(region[2], insert.length); // update shift in
575           // alignment frame of
576           // reference
577           if (segments == null)
578           {
579             // add a hidden column for this deletion
580             colsel.hideColumns(inspos, inspos + insert.length - 1);
581           }
582         }
583       }
584     }
585     for (int i = 0; i < alseqs.length; i++)
586     {
587       int[] bounds = ((int[]) ((Object[]) gs_regions[i])[1]);
588       SequenceI ref = alseqs[i].getRefSeq();
589       seqs[i] = new Sequence(ref.getName(), g_seqs[i].toString(),
590               ref.getStart() + alseqs[i].start + bounds[0], ref.getStart()
591                       + alseqs[i].start + (bounds[2] == 0 ? -1 : bounds[2]));
592       seqs[i].setDatasetSequence(ref);
593       seqs[i].setDescription(ref.getDescription());
594     }
595     if (segments != null)
596     {
597       for (int i = 0; i < segments.length; i += 3)
598       {
599         // int start=shifts.shift(segments[i]-1)+1;
600         // int end=shifts.shift(segments[i]+segments[i+1]-1)-1;
601         colsel.hideColumns(segments[i + 1], segments[i + 1]
602                 + segments[i + 2] - 1);
603       }
604     }
605     return seqs;
606   }
607
608   /**
609    * references to entities that this sequence cigar is associated with.
610    */
611   private Hashtable selGroups = null;
612
613   public void setGroupMembership(Object group)
614   {
615     if (selGroups == null)
616     {
617       selGroups = new Hashtable();
618     }
619     selGroups.put(group, new int[0]);
620   }
621
622   /**
623    * Test for and if present remove association to group.
624    * 
625    * @param group
626    * @return true if group was associated and it was removed
627    */
628   public boolean removeGroupMembership(Object group)
629   {
630     if (selGroups != null && selGroups.containsKey(group))
631     {
632       selGroups.remove(group);
633       return true;
634     }
635     return false;
636   }
637
638   /**
639    * forget all associations for this sequence.
640    */
641   public void clearMemberships()
642   {
643     if (selGroups != null)
644     {
645       selGroups.clear();
646     }
647     selGroups = null;
648   }
649
650   /**
651    * 
652    * @return null or array of all associated entities
653    */
654   public Object[] getAllMemberships()
655   {
656     if (selGroups == null)
657     {
658       return null;
659     }
660     Object[] mmbs = new Object[selGroups.size()];
661     Enumeration en = selGroups.keys();
662     for (int i = 0; en.hasMoreElements(); i++)
663     {
664       mmbs[i] = en.nextElement();
665     }
666     return mmbs;
667   }
668
669   /**
670    * Test for group membership
671    * 
672    * @param sgr
673    *          - a selection group or some other object that may be associated
674    *          with seqCigar
675    * @return true if sgr is associated with this seqCigar
676    */
677   public boolean isMemberOf(Object sgr)
678   {
679     return (selGroups != null) && selGroups.get(sgr) != null;
680   }
681 }