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