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) ? ""
 117             : (String) getSequenceAndDeletions(
 118                     refseq.getSequenceAsString(start, end), GapChar)[0];
 119   }
 120
 121   /**
 122    * recreates a gapped and edited version of RefSeq or null for an empty cigar
 123    * string
 124    *
 125    * @return SequenceI
 126    */
 127   public SequenceI getSeq(char GapChar)
 128   {
 129     Sequence seq;
 130     if (refseq == null || length == 0)
 131     {
 132       return null;
 133     }
 134     Object[] edit_result = getSequenceAndDeletions(
 135             refseq.getSequenceAsString(start, end), GapChar);
 136     if (edit_result == null)
 137     {
 138       throw new Error(MessageManager.getString(
 139               "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() + start
 144                     + ((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(MessageManager
 188               .getString("error.implementation_error_set_seq_null"));
 189     }
 190     if (_s < 0)
 191     {
 192       throw new Error(MessageManager
 193               .formatMessage("error.implementation_error_s", new String[]
 194               { Integer.valueOf(_s).toString() }));
 195     }
 196     String seq_string = seq.getSequenceAsString();
 197     if (_e == 0 || _e < _s || _e > seq_string.length())
 198     {
 199       _e = seq_string.length();
 200     }
 201     // resolve start and end positions relative to ungapped reference sequence
 202     start = seq.findPosition(_s) - seq.getStart();
 203     end = seq.findPosition(_e) - seq.getStart();
 204     int l_ungapped = end - start;
 205     // Find correct sequence to reference and correct start and end - if
 206     // necessary
 207     SequenceI ds = seq.getDatasetSequence();
 208     if (ds == null)
 209     {
 210       // make a new dataset sequence
 211       String ungapped = AlignSeq.extractGaps(
 212               jalview.util.Comparison.GapChars, new String(seq_string));
 213       l_ungapped = ungapped.length();
 214       // check that we haven't just duplicated an ungapped sequence.
 215       if (l_ungapped == seq.getLength())
 216       {
 217         ds = seq;
 218       }
 219       else
 220       {
 221         ds = new Sequence(seq.getName(), ungapped, seq.getStart(),
 222                 seq.getStart() + ungapped.length() - 1);
 223         // JBPNote: this would be consistent but may not be useful
 224         // seq.setDatasetSequence(ds);
 225       }
 226     }
 227     // add in offset between seq and the dataset sequence
 228     if (ds.getStart() < seq.getStart())
 229     {
 230       int offset = seq.getStart() - ds.getStart();
 231       if (initialDeletion)
 232       {
 233         // absolute cigar string
 234         addDeleted(_s + offset);
 235         start = 0;
 236         end += offset;
 237       }
 238       else
 239       {
 240         // normal behaviour - just mark start and end subsequence
 241         start += offset;
 242         end += offset;
 243
 244       }
 245
 246     }
 247
 248     // any gaps to process ?
 249     if (l_ungapped != (_e - _s))
 250     {
 251       hasgaps = true;
 252     }
 253
 254     refseq = ds;
 255     // copy over local properties for the sequence instance of the refseq
 256     seqProps = SeqsetUtils.SeqCharacterHash(seq);
 257     // Check offsets
 258     if (end > ds.getLength())
 259     {
 260       throw new Error(MessageManager
 261               .getString("error.implementation_error_seqcigar_possible"));
 262       // end = ds.getLength();
 263     }
 264
 265     return hasgaps;
 266   }
 267
 268   /**
 269    * directly initialise a cigar object with a sequence of range, operation
 270    * pairs and a sequence to apply it to. operation and range should be relative
 271    * to the seq.getStart()'th residue of the dataset seq resolved from seq.
 272    *
 273    * @param seq
 274    *          SequenceI
 275    * @param operation
 276    *          char[]
 277    * @param range
 278    *          int[]
 279    */
 280   public SeqCigar(SequenceI seq, char operation[], int range[])
 281   {
 282     super();
 283     if (seq == null)
 284     {
 285       throw new Error(
 286               MessageManager.getString("error.implmentation_bug_seq_null"));
 287     }
 288     if (operation.length != range.length)
 289     {
 290       throw new Error(MessageManager.getString(
 291               "error.implementation_bug_cigar_operation_list_range_list"));
 292     }
 293
 294     if (operation != null)
 295     {
 296       this.operation = new char[operation.length + _inc_length];
 297       this.range = new int[operation.length + _inc_length];
 298
 299       if (_setSeq(seq, false, 0, 0))
 300       {
 301         throw new Error(MessageManager.getString(
 302                 "error.not_yet_implemented_cigar_object_from_cigar_string"));
 303       }
 304       for (int i = this.length, j = 0; j < operation.length; i++, j++)
 305       {
 306         char op = operation[j];
 307         if (op != M && op != I && op != D)
 308         {
 309           throw new Error(MessageManager.formatMessage(
 310                   "error.implementation_bug_cigar_operation", new String[]
 311                   { Integer.valueOf(j).toString(),
 312                       Integer.valueOf(op).toString(),
 313                       Integer.valueOf(M).toString(),
 314                       Integer.valueOf(I).toString(),
 315                       Integer.valueOf(D).toString() }));
 316         }
 317         this.operation[i] = op;
 318         this.range[i] = range[j];
 319       }
 320       this.length += operation.length;
 321     }
 322     else
 323     {
 324       this.operation = null;
 325       this.range = null;
 326       this.length = 0;
 327       if (_setSeq(seq, false, 0, 0))
 328       {
 329         throw new Error(MessageManager.getString(
 330                 "error.not_yet_implemented_cigar_object_from_cigar_string"));
 331       }
 332     }
 333   }
 334
 335   /**
 336    * add range matched residues to cigar string
 337    *
 338    * @param range
 339    *          int
 340    */
 341   public void addMatch(int range)
 342   {
 343     this.addOperation(M, range);
 344   }
 345
 346   /**
 347    * Adds insertion and match operations based on seq to the cigar up to the
 348    * endpos column of seq.
 349    *
 350    * @param cigar
 351    *          CigarBase
 352    * @param seq
 353    *          SequenceI
 354    * @param startpos
 355    *          int
 356    * @param endpos
 357    *          int
 358    * @param initialDeletions
 359    *          if true then initial deletions will be added from start of seq to
 360    *          startpos
 361    */
 362   protected static void addSequenceOps(CigarBase cigar, SequenceI seq,
 363           int startpos, int endpos, boolean initialDeletions)
 364   {
 365     char op = '\0';
 366     int range = 0;
 367     int p = 0, res = seq.getLength();
 368
 369     if (!initialDeletions)
 370     {
 371       p = startpos;
 372     }
 373
 374     while (p <= endpos)
 375     {
 376       boolean isGap = (p < res)
 377               ? jalview.util.Comparison.isGap(seq.getCharAt(p))
 378               : true;
 379       if ((startpos <= p) && (p <= endpos))
 380       {
 381         if (isGap)
 382         {
 383           if (range > 0 && op != I)
 384           {
 385             cigar.addOperation(op, range);
 386             range = 0;
 387           }
 388           op = I;
 389           range++;
 390         }
 391         else
 392         {
 393           if (range > 0 && op != M)
 394           {
 395             cigar.addOperation(op, range);
 396             range = 0;
 397           }
 398           op = M;
 399           range++;
 400         }
 401       }
 402       else
 403       {
 404         if (!isGap)
 405         {
 406           if (range > 0 && op != D)
 407           {
 408             cigar.addOperation(op, range);
 409             range = 0;
 410           }
 411           op = D;
 412           range++;
 413         }
 414         else
 415         {
 416           // do nothing - insertions are not made in flanking regions
 417         }
 418       }
 419       p++;
 420     }
 421     if (range > 0)
 422     {
 423       cigar.addOperation(op, range);
 424     }
 425   }
 426
 427   /**
 428    * create a cigar string for given sequence
 429    *
 430    * @param seq
 431    *          SequenceI
 432    */
 433   public SeqCigar(SequenceI seq)
 434   {
 435     super();
 436     if (seq == null)
 437     {
 438       throw new Error(MessageManager
 439               .getString("error.implementation_error_for_new_cigar"));
 440     }
 441     _setSeq(seq, false, 0, 0);
 442     // there is still work to do
 443     addSequenceOps(this, seq, 0, seq.getLength() - 1, false);
 444   }
 445
 446   /**
 447    * Create Cigar from a range of gaps and residues on a sequence object
 448    *
 449    * @param seq
 450    *          SequenceI
 451    * @param start
 452    *          int - first column in range
 453    * @param end
 454    *          int - last column in range
 455    */
 456   public SeqCigar(SequenceI seq, int start, int end)
 457   {
 458     super();
 459     if (seq == null)
 460     {
 461       throw new Error(MessageManager
 462               .getString("error.implementation_error_for_new_cigar"));
 463     }
 464     _setSeq(seq, false, start, end + 1);
 465     // there is still work to do
 466     addSequenceOps(this, seq, start, end, false);
 467   }
 468
 469   /**
 470    * create an alignment from the given array of cigar sequences and gap
 471    * character, and marking the given segments as visible in the given
 472    * hiddenColumns.
 473    *
 474    * @param alseqs
 475    * @param gapCharacter
 476    * @param hidden
 477    *          - hiddenColumns where hidden regions are marked
 478    * @param segments
 479    *          - visible regions of alignment
 480    * @return SequenceI[]
 481    */
 482   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
 483           char gapCharacter, HiddenColumns hidden, int[] segments)
 484   {
 485     SequenceI[] seqs = new SequenceI[alseqs.length];
 486     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
 487     String[] alseqs_string = new String[alseqs.length];
 488     Object[] gs_regions = new Object[alseqs.length];
 489     for (int i = 0; i < alseqs.length; i++)
 490     {
 491       alseqs_string[i] = alseqs[i].getRefSeq()
 492               .getSequenceAsString(alseqs[i].start, alseqs[i].end);
 493       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
 494               gapCharacter); // gapped sequence, {start, start col, end.
 495       // endcol}, hidden regions {{start, end, col}})
 496       if (gs_regions[i] == null)
 497       {
 498         throw new Error(MessageManager.formatMessage(
 499                 "error.implementation_error_cigar_seq_no_operations",
 500                 new String[]
 501                 { Integer.valueOf(i).toString() }));
 502       }
 503       g_seqs[i] = new StringBuffer((String) ((Object[]) gs_regions[i])[0]); // the
 504       // visible
 505       // gapped
 506       // sequence
 507     }
 508     // Now account for insertions. (well - deletions)
 509     // this is complicated because we must keep track of shifted positions in
 510     // each sequence
 511     ShiftList shifts = new ShiftList();
 512     for (int i = 0; i < alseqs.length; i++)
 513     {
 514       Object[] gs_region = ((Object[]) ((Object[]) gs_regions[i])[2]);
 515       if (gs_region != null)
 516
 517       {
 518         for (int hr = 0; hr < gs_region.length; hr++)
 519         {
 520           int[] region = (int[]) gs_region[hr];
 521           char[] insert = new char[region[1] - region[0] + 1];
 522           for (int s = 0; s < insert.length; s++)
 523           {
 524             insert[s] = gapCharacter;
 525           }
 526           int inspos = shifts.shift(region[2]); // resolve insertion position in
 527           // current alignment frame of
 528           // reference
 529           for (int s = 0; s < alseqs.length; s++)
 530           {
 531             if (s != i)
 532             {
 533               if (g_seqs[s].length() <= inspos)
 534               {
 535                 // prefix insertion with more gaps.
 536                 for (int l = inspos - g_seqs[s].length(); l > 0; l--)
 537                 {
 538                   g_seqs[s].append(gapCharacter); // to debug - use a diffferent
 539                   // gap character here
 540                 }
 541               }
 542               g_seqs[s].insert(inspos, insert);
 543             }
 544             else
 545             {
 546               g_seqs[s].insert(inspos,
 547                       alseqs_string[i].substring(region[0], region[1] + 1));
 548             }
 549           }
 550           shifts.addShift(region[2], insert.length); // update shift in
 551           // alignment frame of
 552           // reference
 553           if (segments == null)
 554           {
 555             // add a hidden column for this deletion
 556             hidden.hideColumns(inspos, inspos + insert.length - 1);
 557           }
 558         }
 559       }
 560     }
 561     for (int i = 0; i < alseqs.length; i++)
 562     {
 563       int[] bounds = ((int[]) ((Object[]) gs_regions[i])[1]);
 564       SequenceI ref = alseqs[i].getRefSeq();
 565       seqs[i] = new Sequence(ref.getName(), g_seqs[i].toString(),
 566               ref.getStart() + alseqs[i].start + bounds[0],
 567               ref.getStart() + alseqs[i].start
 568                       + (bounds[2] == 0 ? -1 : bounds[2]));
 569       seqs[i].setDatasetSequence(ref);
 570       seqs[i].setDescription(ref.getDescription());
 571     }
 572     if (segments != null)
 573     {
 574       for (int i = 0; i < segments.length; i += 3)
 575       {
 576         // int start=shifts.shift(segments[i]-1)+1;
 577         // int end=shifts.shift(segments[i]+segments[i+1]-1)-1;
 578         hidden.hideColumns(segments[i + 1],
 579                 segments[i + 1] + segments[i + 2] - 1);
 580       }
 581     }
 582     return seqs;
 583   }
 584
 585   /**
 586    * references to entities that this sequence cigar is associated with.
 587    */
 588   private Hashtable selGroups = null;
 589
 590   public void setGroupMembership(Object group)
 591   {
 592     if (selGroups == null)
 593     {
 594       selGroups = new Hashtable();
 595     }
 596     selGroups.put(group, new int[0]);
 597   }
 598
 599   /**
 600    * Test for and if present remove association to group.
 601    *
 602    * @param group
 603    * @return true if group was associated and it was removed
 604    */
 605   public boolean removeGroupMembership(Object group)
 606   {
 607     if (selGroups != null && selGroups.containsKey(group))
 608     {
 609       selGroups.remove(group);
 610       return true;
 611     }
 612     return false;
 613   }
 614
 615   /**
 616    * forget all associations for this sequence.
 617    */
 618   public void clearMemberships()
 619   {
 620     if (selGroups != null)
 621     {
 622       selGroups.clear();
 623     }
 624     selGroups = null;
 625   }
 626
 627   /**
 628    *
 629    * @return null or array of all associated entities
 630    */
 631   public Object[] getAllMemberships()
 632   {
 633     if (selGroups == null)
 634     {
 635       return null;
 636     }
 637     Object[] mmbs = new Object[selGroups.size()];
 638     Enumeration en = selGroups.keys();
 639     for (int i = 0; en.hasMoreElements(); i++)
 640     {
 641       mmbs[i] = en.nextElement();
 642     }
 643     return mmbs;
 644   }
 645
 646   /**
 647    * Test for group membership
 648    *
 649    * @param sgr
 650    *          - a selection group or some other object that may be associated
 651    *          with seqCigar
 652    * @return true if sgr is associated with this seqCigar
 653    */
 654   public boolean isMemberOf(Object sgr)
 655   {
 656     return (selGroups != null) && selGroups.get(sgr) != null;
 657   }
 658 }