src/jalview/datamodel/SeqCigar.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.4.0.b2)
   3  * Copyright (C) 2009 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
   4  *
   5  * This program is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU General Public License
   7  * as published by the Free Software Foundation; either version 2
   8  * of the License, or (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  18  */
  19 package jalview.datamodel;
  20
  21 import java.util.Hashtable;
  22
  23 import jalview.analysis.*;
  24 import jalview.util.*;
  25
  26 public class SeqCigar extends CigarSimple
  27 {
  28   /**
  29    * start(inclusive) and end(exclusive) of subsequence on refseq
  30    */
  31   private int start, end;
  32
  33   private SequenceI refseq = null;
  34
  35   private Hashtable seqProps;
  36
  37   /**
  38    * Reference dataset sequence for the cigar string
  39    *
  40    * @return SequenceI
  41    */
  42   public SequenceI getRefSeq()
  43   {
  44     return refseq;
  45   }
  46
  47   /**
  48    *
  49    * @return int start index of cigar ops on refSeq
  50    */
  51   public int getStart()
  52   {
  53     return start;
  54   }
  55
  56   /**
  57    *
  58    * @return int end index (exclusive) of cigar ops on refSeq
  59    */
  60   public int getEnd()
  61   {
  62     return end;
  63   }
  64
  65   /**
  66    * Returns sequence as a string with cigar operations applied to it
  67    *
  68    * @return String
  69    */
  70   public String getSequenceString(char GapChar)
  71   {
  72     return (length == 0) ? "" : (String) getSequenceAndDeletions(refseq
  73             .getSequenceAsString(start, end), GapChar)[0];
  74   }
  75
  76   /**
  77    * recreates a gapped and edited version of RefSeq or null for an empty cigar
  78    * string
  79    *
  80    * @return SequenceI
  81    */
  82   public SequenceI getSeq(char GapChar)
  83   {
  84     Sequence seq;
  85     if (refseq == null || length == 0)
  86     {
  87       return null;
  88     }
  89     Object[] edit_result = getSequenceAndDeletions(refseq
  90             .getSequenceAsString(start, end), GapChar);
  91     if (edit_result == null)
  92     {
  93       throw new Error(
  94               "Implementation Error - unexpected null from getSequenceAndDeletions");
  95     }
  96     int bounds[] = (int[]) edit_result[1];
  97     seq = new Sequence(refseq.getName(), (String) edit_result[0], refseq
  98             .getStart()
  99             + start + bounds[0], refseq.getStart() + start
 100             + ((bounds[2] == 0) ? -1 : bounds[2]));
 101     seq.setDescription(refseq.getDescription());
 102     int sstart = seq.getStart(), send = seq.getEnd();
 103     // seq.checkValidRange(); probably not needed
 104     // recover local properties if present
 105     if (seqProps != null)
 106     {
 107       // this recovers dataset sequence reference as well as local features,
 108       // names, start/end settings.
 109       SeqsetUtils.SeqCharacterUnhash(seq, seqProps);
 110     }
 111     // ensure dataset sequence is up to date from local reference
 112     seq.setDatasetSequence(refseq);
 113     seq.setStart(sstart);
 114     seq.setEnd(send);
 115     return seq;
 116   }
 117
 118   /*
 119    * We don't allow this - refseq is given at construction time only public void
 120    * setSeq(SequenceI seq) { this.seq = seq; }
 121    */
 122   /**
 123    * internal constructor - sets seq to a gapless sequence derived from seq and
 124    * prepends any 'D' operations needed to get to the first residue of seq.
 125    *
 126    * @param seq
 127    *                SequenceI
 128    * @param initialDeletion
 129    *                true to mark initial dataset sequence residues as deleted in
 130    *                subsequence
 131    * @param _s
 132    *                index of first position in seq
 133    * @param _e
 134    *                index after last position in (possibly gapped) seq
 135    * @return true if gaps are present in seq
 136    */
 137   private boolean _setSeq(SequenceI seq, boolean initialDeletion, int _s,
 138           int _e)
 139   {
 140     boolean hasgaps = false;
 141     if (seq == null)
 142     {
 143       throw new Error("Implementation Error - _setSeq(null,...)");
 144     }
 145     if (_s < 0)
 146     {
 147       throw new Error("Implementation Error: _s=" + _s);
 148     }
 149     String seq_string = seq.getSequenceAsString();
 150     if (_e == 0 || _e < _s || _e > seq_string.length())
 151     {
 152       _e = seq_string.length();
 153     }
 154     // resolve start and end positions relative to ungapped reference sequence
 155     start = seq.findPosition(_s) - seq.getStart();
 156     end = seq.findPosition(_e) - seq.getStart();
 157     int l_ungapped = end - start;
 158     // Find correct sequence to reference and correct start and end - if
 159     // necessary
 160     SequenceI ds = seq.getDatasetSequence();
 161     if (ds == null)
 162     {
 163       // make a new dataset sequence
 164       String ungapped = AlignSeq.extractGaps(
 165               jalview.util.Comparison.GapChars, new String(seq_string));
 166       l_ungapped = ungapped.length();
 167       // check that we haven't just duplicated an ungapped sequence.
 168       if (l_ungapped == seq.getLength())
 169       {
 170         ds = seq;
 171       }
 172       else
 173       {
 174         ds = new Sequence(seq.getName(), ungapped, seq.getStart(), seq
 175                 .getStart()
 176                 + ungapped.length() - 1);
 177         // JBPNote: this would be consistent but may not be useful
 178         // seq.setDatasetSequence(ds);
 179       }
 180     }
 181     // add in offset between seq and the dataset sequence
 182     if (ds.getStart() < seq.getStart())
 183     {
 184       int offset = seq.getStart() - ds.getStart();
 185       if (initialDeletion)
 186       {
 187         // absolute cigar string
 188         addDeleted(_s + offset);
 189         start = 0;
 190         end += offset;
 191       }
 192       else
 193       {
 194         // normal behaviour - just mark start and end subsequence
 195         start += offset;
 196         end += offset;
 197
 198       }
 199
 200     }
 201
 202     // any gaps to process ?
 203     if (l_ungapped != (_e - _s))
 204     {
 205       hasgaps = true;
 206     }
 207
 208     refseq = ds;
 209     // copy over local properties for the sequence instance of the refseq
 210     seqProps = SeqsetUtils.SeqCharacterHash(seq);
 211     // Check offsets
 212     if (end > ds.getLength())
 213     {
 214       throw new Error(
 215               "SeqCigar: Possible implementation error: sequence is longer than dataset sequence");
 216       // end = ds.getLength();
 217     }
 218
 219     return hasgaps;
 220   }
 221
 222   /**
 223    * directly initialise a cigar object with a sequence of range, operation
 224    * pairs and a sequence to apply it to. operation and range should be relative
 225    * to the seq.getStart()'th residue of the dataset seq resolved from seq.
 226    *
 227    * @param seq
 228    *                SequenceI
 229    * @param operation
 230    *                char[]
 231    * @param range
 232    *                int[]
 233    */
 234   public SeqCigar(SequenceI seq, char operation[], int range[])
 235   {
 236     super();
 237     if (seq == null)
 238     {
 239       throw new Error("Implementation Bug. Null seq !");
 240     }
 241     if (operation.length != range.length)
 242     {
 243       throw new Error(
 244               "Implementation Bug. Cigar Operation list!= range list");
 245     }
 246
 247     if (operation != null)
 248     {
 249       this.operation = new char[operation.length + _inc_length];
 250       this.range = new int[operation.length + _inc_length];
 251
 252       if (_setSeq(seq, false, 0, 0))
 253       {
 254         throw new Error(
 255                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 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("Implementation Bug. Cigar Operation '" + j
 263                   + "' '" + op + "' not one of '" + M + "', '" + I
 264                   + "', or '" + D + "'.");
 265         }
 266         this.operation[i] = op;
 267         this.range[i] = range[j];
 268       }
 269       this.length += operation.length;
 270     }
 271     else
 272     {
 273       this.operation = null;
 274       this.range = null;
 275       this.length = 0;
 276       if (_setSeq(seq, false, 0, 0))
 277       {
 278         throw new Error(
 279                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 280       }
 281     }
 282   }
 283
 284   /**
 285    * add range matched residues to cigar string
 286    *
 287    * @param range
 288    *                int
 289    */
 290   public void addMatch(int range)
 291   {
 292     this.addOperation(M, range);
 293   }
 294
 295   /**
 296    * Adds insertion and match operations based on seq to the cigar up to the
 297    * endpos column of seq.
 298    *
 299    * @param cigar
 300    *                CigarBase
 301    * @param seq
 302    *                SequenceI
 303    * @param startpos
 304    *                int
 305    * @param endpos
 306    *                int
 307    * @param initialDeletions
 308    *                if true then initial deletions will be added from start of
 309    *                seq to startpos
 310    */
 311   protected static void addSequenceOps(CigarBase cigar, SequenceI seq,
 312           int startpos, int endpos, boolean initialDeletions)
 313   {
 314     char op = '\0';
 315     int range = 0;
 316     int p = 0, res = seq.getLength();
 317
 318     if (!initialDeletions)
 319     {
 320       p = startpos;
 321     }
 322
 323     while (p <= endpos)
 324     {
 325       boolean isGap = (p < res) ? jalview.util.Comparison.isGap(seq
 326               .getCharAt(p)) : true;
 327       if ((startpos <= p) && (p <= endpos))
 328       {
 329         if (isGap)
 330         {
 331           if (range > 0 && op != I)
 332           {
 333             cigar.addOperation(op, range);
 334             range = 0;
 335           }
 336           op = I;
 337           range++;
 338         }
 339         else
 340         {
 341           if (range > 0 && op != M)
 342           {
 343             cigar.addOperation(op, range);
 344             range = 0;
 345           }
 346           op = M;
 347           range++;
 348         }
 349       }
 350       else
 351       {
 352         if (!isGap)
 353         {
 354           if (range > 0 && op != D)
 355           {
 356             cigar.addOperation(op, range);
 357             range = 0;
 358           }
 359           op = D;
 360           range++;
 361         }
 362         else
 363         {
 364           // do nothing - insertions are not made in flanking regions
 365         }
 366       }
 367       p++;
 368     }
 369     if (range > 0)
 370     {
 371       cigar.addOperation(op, range);
 372     }
 373   }
 374
 375   /**
 376    * create a cigar string for given sequence
 377    *
 378    * @param seq
 379    *                SequenceI
 380    */
 381   public SeqCigar(SequenceI seq)
 382   {
 383     super();
 384     if (seq == null)
 385     {
 386       throw new Error("Implementation error for new Cigar(SequenceI)");
 387     }
 388     _setSeq(seq, false, 0, 0);
 389     // there is still work to do
 390     addSequenceOps(this, seq, 0, seq.getLength() - 1, false);
 391   }
 392
 393   /**
 394    * Create Cigar from a range of gaps and residues on a sequence object
 395    *
 396    * @param seq
 397    *                SequenceI
 398    * @param start
 399    *                int - first column in range
 400    * @param end
 401    *                int - last column in range
 402    */
 403   public SeqCigar(SequenceI seq, int start, int end)
 404   {
 405     super();
 406     if (seq == null)
 407     {
 408       throw new Error("Implementation error for new Cigar(SequenceI)");
 409     }
 410     _setSeq(seq, false, start, end + 1);
 411     // there is still work to do
 412     addSequenceOps(this, seq, start, end, false);
 413   }
 414
 415   /**
 416    * Create a cigar object from a cigar string like '[<I|D|M><range>]+' Will
 417    * fail if the given seq already contains gaps (JBPNote: future implementation
 418    * will fix)
 419    *
 420    * @param seq
 421    *                SequenceI object resolvable to a dataset sequence
 422    * @param cigarString
 423    *                String
 424    * @return Cigar
 425    */
 426   public static SeqCigar parseCigar(SequenceI seq, String cigarString)
 427           throws Exception
 428   {
 429     Object[] opsandrange = parseCigarString(cigarString);
 430     return new SeqCigar(seq, (char[]) opsandrange[0],
 431             (int[]) opsandrange[1]);
 432   }
 433
 434   /**
 435    * createAlignment
 436    *
 437    * @param alseqs
 438    *                SeqCigar[]
 439    * @param gapCharacter
 440    *                char
 441    * @return SequenceI[]
 442    */
 443   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
 444           char gapCharacter, ColumnSelection colsel, int[] segments)
 445   {
 446     SequenceI[] seqs = new SequenceI[alseqs.length];
 447     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
 448     String[] alseqs_string = new String[alseqs.length];
 449     Object[] gs_regions = new Object[alseqs.length];
 450     for (int i = 0; i < alseqs.length; i++)
 451     {
 452       alseqs_string[i] = alseqs[i].getRefSeq().getSequenceAsString(
 453               alseqs[i].start, alseqs[i].end);
 454       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
 455               gapCharacter); // gapped sequence, {start, start col, end.
 456                               // endcol}, hidden regions {{start, end, col}})
 457       if (gs_regions[i] == null)
 458       {
 459         throw new Error("Implementation error: " + i
 460                 + "'th sequence Cigar has no operations.");
 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, alseqs_string[i].substring(
 506                       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(), ref
 525               .getStart()
 526               + alseqs[i].start + bounds[0], ref.getStart()
 527               + alseqs[i].start + (bounds[2] == 0 ? -1 : bounds[2]));
 528       seqs[i].setDatasetSequence(ref);
 529       seqs[i].setDescription(ref.getDescription());
 530     }
 531     if (segments != null)
 532     {
 533       for (int i = 0; i < segments.length; i += 3)
 534       {
 535         // int start=shifts.shift(segments[i]-1)+1;
 536         // int end=shifts.shift(segments[i]+segments[i+1]-1)-1;
 537         colsel.hideColumns(segments[i + 1], segments[i + 1]
 538                 + segments[i + 2] - 1);
 539       }
 540     }
 541     return seqs;
 542   }
 543
 544   /**
 545    * non rigorous testing
 546    */
 547   /**
 548    *
 549    * @param seq
 550    *                Sequence
 551    * @param ex_cs_gapped
 552    *                String
 553    * @return String
 554    */
 555   public static String testCigar_string(Sequence seq, String ex_cs_gapped)
 556   {
 557     SeqCigar c_sgapped = new SeqCigar(seq);
 558     String cs_gapped = c_sgapped.getCigarstring();
 559     if (!cs_gapped.equals(ex_cs_gapped))
 560     {
 561       System.err.println("Failed getCigarstring: incorect string '"
 562               + cs_gapped + "' != " + ex_cs_gapped);
 563     }
 564     return cs_gapped;
 565   }
 566
 567   public static boolean testSeqRecovery(SeqCigar gen_sgapped,
 568           SequenceI s_gapped)
 569   {
 570     // this is non-rigorous - start and end recovery is not tested.
 571     SequenceI gen_sgapped_s = gen_sgapped.getSeq('-');
 572     if (!gen_sgapped_s.getSequence().equals(s_gapped.getSequence()))
 573     {
 574       System.err.println("Couldn't reconstruct sequence.\n"
 575               + gen_sgapped_s.getSequenceAsString() + "\n"
 576               + s_gapped.getSequenceAsString());
 577       return false;
 578     }
 579     return true;
 580   }
 581
 582   public static void main(String argv[]) throws Exception
 583   {
 584     String o_seq;
 585     Sequence s = new Sequence("MySeq",
 586             o_seq = "asdfktryasdtqwrtsaslldddptyipqqwaslchvhttt", 39, 80);
 587     String orig_gapped;
 588     Sequence s_gapped = new Sequence(
 589             "MySeq",
 590             orig_gapped = "----asdf------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhttt",
 591             39, 80);
 592     String ex_cs_gapped = "4I4M6I6M3I11M4I12M4I9M";
 593     s_gapped.setDatasetSequence(s);
 594     String sub_gapped_s;
 595     Sequence s_subsequence_gapped = new Sequence(
 596             "MySeq",
 597             sub_gapped_s = "------ktryas---dtqwrtsasll----dddptyipqqwa----slchvh",
 598             43, 77);
 599
 600     s_subsequence_gapped.setDatasetSequence(s);
 601     SeqCigar c_null = new SeqCigar(s);
 602     String cs_null = c_null.getCigarstring();
 603     if (!cs_null.equals("42M"))
 604     {
 605       System.err
 606               .println("Failed to recover ungapped sequence cigar operations:"
 607                       + ((cs_null == "") ? "empty string" : cs_null));
 608     }
 609     testCigar_string(s_gapped, ex_cs_gapped);
 610     SeqCigar gen_sgapped = SeqCigar.parseCigar(s, ex_cs_gapped);
 611     if (!gen_sgapped.getCigarstring().equals(ex_cs_gapped))
 612     {
 613       System.err.println("Failed parseCigar(" + ex_cs_gapped
 614               + ")->getCigarString()->'" + gen_sgapped.getCigarstring()
 615               + "'");
 616     }
 617     testSeqRecovery(gen_sgapped, s_gapped);
 618     // Test dataset resolution
 619     SeqCigar sub_gapped = new SeqCigar(s_subsequence_gapped);
 620     if (!testSeqRecovery(sub_gapped, s_subsequence_gapped))
 621     {
 622       System.err
 623               .println("Failed recovery for subsequence of dataset sequence");
 624     }
 625     // width functions
 626     if (sub_gapped.getWidth() != sub_gapped_s.length())
 627     {
 628       System.err.println("Failed getWidth()");
 629     }
 630
 631     sub_gapped.getFullWidth();
 632     if (sub_gapped.hasDeletedRegions())
 633     {
 634       System.err.println("hasDeletedRegions is incorrect.");
 635     }
 636     // Test start-end region SeqCigar
 637     SeqCigar sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 638     if (sub_se_gp.getWidth() != 41)
 639     {
 640       System.err
 641               .println("SeqCigar(seq, start, end) not properly clipped alignsequence.");
 642     }
 643     System.out.println("Original sequence align:\n" + sub_gapped_s
 644             + "\nReconstructed window from 8 to 48\n" + "XXXXXXXX"
 645             + sub_se_gp.getSequenceString('-') + "..." + "\nCigar String:"
 646             + sub_se_gp.getCigarstring() + "\n");
 647     SequenceI ssgp = sub_se_gp.getSeq('-');
 648     System.out.println("\t " + ssgp.getSequenceAsString());
 649     for (int r = 0; r < 10; r++)
 650     {
 651       sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 652       int sl = sub_se_gp.getWidth();
 653       int st = sl - 1 - r;
 654       for (int rs = 0; rs < 10; rs++)
 655       {
 656         int e = st + rs;
 657         sub_se_gp.deleteRange(st, e);
 658         String ssgapedseq = sub_se_gp.getSeq('-').getSequenceAsString();
 659         System.out.println(st + "," + e + "\t:" + ssgapedseq);
 660         st -= 3;
 661       }
 662     }
 663     {
 664       SeqCigar[] set = new SeqCigar[]
 665       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 666           new SeqCigar(s_gapped) };
 667       Alignment al = new Alignment(set);
 668       for (int i = 0; i < al.getHeight(); i++)
 669       {
 670         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 671                 + al.getSequenceAt(i).getStart() + "\t"
 672                 + al.getSequenceAt(i).getEnd() + "\t"
 673                 + al.getSequenceAt(i).getSequenceAsString());
 674       }
 675     }
 676     {
 677       System.out.println("Gapped.");
 678       SeqCigar[] set = new SeqCigar[]
 679       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 680           new SeqCigar(s_gapped) };
 681       set[0].deleteRange(20, 25);
 682       Alignment al = new Alignment(set);
 683       for (int i = 0; i < al.getHeight(); i++)
 684       {
 685         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 686                 + al.getSequenceAt(i).getStart() + "\t"
 687                 + al.getSequenceAt(i).getEnd() + "\t"
 688                 + al.getSequenceAt(i).getSequenceAsString());
 689       }
 690     }
 691     // if (!ssgapedseq.equals("ryas---dtqqwa----slchvh"))
 692     // System.err.println("Subseqgaped\n------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhryas---dtqwrtsasll--qwa----slchvh\n"+ssgapedseq+"\n"+sub_se_gp.getCigarstring());
 693   }
 694
 695 }