src/jalview/datamodel/SeqCigar.java

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