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