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 }