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