src/jalview/datamodel/SeqCigar.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)
   3  * Copyright (C) 2010 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle
   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 of the License, or (at your option) any later version.
  10  *
  11  * Jalview is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14  * PURPOSE.  See the GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  17  */
  18 package jalview.datamodel;
  19
  20 import java.util.Hashtable;
  21
  22 import jalview.analysis.*;
  23 import jalview.util.*;
  24
  25 public class SeqCigar extends CigarSimple
  26 {
  27   /**
  28    * start(inclusive) and end(exclusive) of subsequence on refseq
  29    */
  30   private int start, end;
  31
  32   private SequenceI refseq = null;
  33
  34   private Hashtable seqProps;
  35
  36   /**
  37    * Reference dataset sequence for the cigar string
  38    *
  39    * @return SequenceI
  40    */
  41   public SequenceI getRefSeq()
  42   {
  43     return refseq;
  44   }
  45
  46   /**
  47    *
  48    * @return int start index of cigar ops on refSeq
  49    */
  50   public int getStart()
  51   {
  52     return start;
  53   }
  54
  55   /**
  56    *
  57    * @return int end index (exclusive) of cigar ops on refSeq
  58    */
  59   public int getEnd()
  60   {
  61     return end;
  62   }
  63
  64   /**
  65    * Returns sequence as a string with cigar operations applied to it
  66    *
  67    * @return String
  68    */
  69   public String getSequenceString(char GapChar)
  70   {
  71     return (length == 0) ? "" : (String) getSequenceAndDeletions(
  72             refseq.getSequenceAsString(start, end), GapChar)[0];
  73   }
  74
  75   /**
  76    * recreates a gapped and edited version of RefSeq or null for an empty cigar
  77    * string
  78    *
  79    * @return SequenceI
  80    */
  81   public SequenceI getSeq(char GapChar)
  82   {
  83     Sequence seq;
  84     if (refseq == null || length == 0)
  85     {
  86       return null;
  87     }
  88     Object[] edit_result = getSequenceAndDeletions(
  89             refseq.getSequenceAsString(start, end), GapChar);
  90     if (edit_result == null)
  91     {
  92       throw new Error(
  93               "Implementation Error - unexpected null from getSequenceAndDeletions");
  94     }
  95     int bounds[] = (int[]) edit_result[1];
  96     seq = new Sequence(refseq.getName(), (String) edit_result[0],
  97             refseq.getStart() + start + bounds[0], refseq.getStart()
  98                     + start + ((bounds[2] == 0) ? -1 : bounds[2]));
  99     seq.setDescription(refseq.getDescription());
 100     int sstart = seq.getStart(), send = seq.getEnd();
 101     // seq.checkValidRange(); probably not needed
 102     // recover local properties if present
 103     if (seqProps != null)
 104     {
 105       // this recovers dataset sequence reference as well as local features,
 106       // names, start/end settings.
 107       SeqsetUtils.SeqCharacterUnhash(seq, seqProps);
 108     }
 109     // ensure dataset sequence is up to date from local reference
 110     seq.setDatasetSequence(refseq);
 111     seq.setStart(sstart);
 112     seq.setEnd(send);
 113     return seq;
 114   }
 115
 116   /*
 117    * We don't allow this - refseq is given at construction time only public void
 118    * setSeq(SequenceI seq) { this.seq = seq; }
 119    */
 120   /**
 121    * internal constructor - sets seq to a gapless sequence derived from seq and
 122    * prepends any 'D' operations needed to get to the first residue of seq.
 123    *
 124    * @param seq
 125    *          SequenceI
 126    * @param initialDeletion
 127    *          true to mark initial dataset sequence residues as deleted in
 128    *          subsequence
 129    * @param _s
 130    *          index of first position in seq
 131    * @param _e
 132    *          index after last position in (possibly gapped) seq
 133    * @return true if gaps are present in seq
 134    */
 135   private boolean _setSeq(SequenceI seq, boolean initialDeletion, int _s,
 136           int _e)
 137   {
 138     boolean hasgaps = false;
 139     if (seq == null)
 140     {
 141       throw new Error("Implementation Error - _setSeq(null,...)");
 142     }
 143     if (_s < 0)
 144     {
 145       throw new Error("Implementation Error: _s=" + _s);
 146     }
 147     String seq_string = seq.getSequenceAsString();
 148     if (_e == 0 || _e < _s || _e > seq_string.length())
 149     {
 150       _e = seq_string.length();
 151     }
 152     // resolve start and end positions relative to ungapped reference sequence
 153     start = seq.findPosition(_s) - seq.getStart();
 154     end = seq.findPosition(_e) - seq.getStart();
 155     int l_ungapped = end - start;
 156     // Find correct sequence to reference and correct start and end - if
 157     // necessary
 158     SequenceI ds = seq.getDatasetSequence();
 159     if (ds == null)
 160     {
 161       // make a new dataset sequence
 162       String ungapped = AlignSeq.extractGaps(
 163               jalview.util.Comparison.GapChars, new String(seq_string));
 164       l_ungapped = ungapped.length();
 165       // check that we haven't just duplicated an ungapped sequence.
 166       if (l_ungapped == seq.getLength())
 167       {
 168         ds = seq;
 169       }
 170       else
 171       {
 172         ds = new Sequence(seq.getName(), ungapped, seq.getStart(),
 173                 seq.getStart() + ungapped.length() - 1);
 174         // JBPNote: this would be consistent but may not be useful
 175         // seq.setDatasetSequence(ds);
 176       }
 177     }
 178     // add in offset between seq and the dataset sequence
 179     if (ds.getStart() < seq.getStart())
 180     {
 181       int offset = seq.getStart() - ds.getStart();
 182       if (initialDeletion)
 183       {
 184         // absolute cigar string
 185         addDeleted(_s + offset);
 186         start = 0;
 187         end += offset;
 188       }
 189       else
 190       {
 191         // normal behaviour - just mark start and end subsequence
 192         start += offset;
 193         end += offset;
 194
 195       }
 196
 197     }
 198
 199     // any gaps to process ?
 200     if (l_ungapped != (_e - _s))
 201     {
 202       hasgaps = true;
 203     }
 204
 205     refseq = ds;
 206     // copy over local properties for the sequence instance of the refseq
 207     seqProps = SeqsetUtils.SeqCharacterHash(seq);
 208     // Check offsets
 209     if (end > ds.getLength())
 210     {
 211       throw new Error(
 212               "SeqCigar: Possible implementation error: sequence is longer than dataset sequence");
 213       // end = ds.getLength();
 214     }
 215
 216     return hasgaps;
 217   }
 218
 219   /**
 220    * directly initialise a cigar object with a sequence of range, operation
 221    * pairs and a sequence to apply it to. operation and range should be relative
 222    * to the seq.getStart()'th residue of the dataset seq resolved from seq.
 223    *
 224    * @param seq
 225    *          SequenceI
 226    * @param operation
 227    *          char[]
 228    * @param range
 229    *          int[]
 230    */
 231   public SeqCigar(SequenceI seq, char operation[], int range[])
 232   {
 233     super();
 234     if (seq == null)
 235     {
 236       throw new Error("Implementation Bug. Null seq !");
 237     }
 238     if (operation.length != range.length)
 239     {
 240       throw new Error(
 241               "Implementation Bug. Cigar Operation list!= range list");
 242     }
 243
 244     if (operation != null)
 245     {
 246       this.operation = new char[operation.length + _inc_length];
 247       this.range = new int[operation.length + _inc_length];
 248
 249       if (_setSeq(seq, false, 0, 0))
 250       {
 251         throw new Error(
 252                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 253       }
 254       for (int i = this.length, j = 0; j < operation.length; i++, j++)
 255       {
 256         char op = operation[j];
 257         if (op != M && op != I && op != D)
 258         {
 259           throw new Error("Implementation Bug. Cigar Operation '" + j
 260                   + "' '" + op + "' not one of '" + M + "', '" + I
 261                   + "', or '" + D + "'.");
 262         }
 263         this.operation[i] = op;
 264         this.range[i] = range[j];
 265       }
 266       this.length += operation.length;
 267     }
 268     else
 269     {
 270       this.operation = null;
 271       this.range = null;
 272       this.length = 0;
 273       if (_setSeq(seq, false, 0, 0))
 274       {
 275         throw new Error(
 276                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 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) ? jalview.util.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("Implementation error for new Cigar(SequenceI)");
 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("Implementation error for new Cigar(SequenceI)");
 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    * createAlignment
 433    *
 434    * @param alseqs
 435    *          SeqCigar[]
 436    * @param gapCharacter
 437    *          char
 438    * @return SequenceI[]
 439    */
 440   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
 441           char gapCharacter, ColumnSelection colsel, int[] segments)
 442   {
 443     SequenceI[] seqs = new SequenceI[alseqs.length];
 444     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
 445     String[] alseqs_string = new String[alseqs.length];
 446     Object[] gs_regions = new Object[alseqs.length];
 447     for (int i = 0; i < alseqs.length; i++)
 448     {
 449       alseqs_string[i] = alseqs[i].getRefSeq().getSequenceAsString(
 450               alseqs[i].start, alseqs[i].end);
 451       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
 452               gapCharacter); // gapped sequence, {start, start col, end.
 453       // endcol}, hidden regions {{start, end, col}})
 454       if (gs_regions[i] == null)
 455       {
 456         throw new Error("Implementation error: " + i
 457                 + "'th sequence Cigar has no operations.");
 458       }
 459       g_seqs[i] = new StringBuffer((String) ((Object[]) gs_regions[i])[0]); // the
 460       // visible
 461       // gapped
 462       // sequence
 463     }
 464     // Now account for insertions. (well - deletions)
 465     // this is complicated because we must keep track of shifted positions in
 466     // each sequence
 467     ShiftList shifts = new ShiftList();
 468     for (int i = 0; i < alseqs.length; i++)
 469     {
 470       Object[] gs_region = ((Object[]) ((Object[]) gs_regions[i])[2]);
 471       if (gs_region != null)
 472
 473       {
 474         for (int hr = 0; hr < gs_region.length; hr++)
 475         {
 476           int[] region = (int[]) gs_region[hr];
 477           char[] insert = new char[region[1] - region[0] + 1];
 478           for (int s = 0; s < insert.length; s++)
 479           {
 480             insert[s] = gapCharacter;
 481           }
 482           int inspos = shifts.shift(region[2]); // resolve insertion position in
 483           // current alignment frame of
 484           // reference
 485           for (int s = 0; s < alseqs.length; s++)
 486           {
 487             if (s != i)
 488             {
 489               if (g_seqs[s].length() <= inspos)
 490               {
 491                 // prefix insertion with more gaps.
 492                 for (int l = inspos - g_seqs[s].length(); l > 0; l--)
 493                 {
 494                   g_seqs[s].append(gapCharacter); // to debug - use a diffferent
 495                   // gap character here
 496                 }
 497               }
 498               g_seqs[s].insert(inspos, insert);
 499             }
 500             else
 501             {
 502               g_seqs[s].insert(inspos,
 503                       alseqs_string[i].substring(region[0], region[1] + 1));
 504             }
 505           }
 506           shifts.addShift(region[2], insert.length); // update shift in
 507           // alignment frame of
 508           // reference
 509           if (segments == null)
 510           {
 511             // add a hidden column for this deletion
 512             colsel.hideColumns(inspos, inspos + insert.length - 1);
 513           }
 514         }
 515       }
 516     }
 517     for (int i = 0; i < alseqs.length; i++)
 518     {
 519       int[] bounds = ((int[]) ((Object[]) gs_regions[i])[1]);
 520       SequenceI ref = alseqs[i].getRefSeq();
 521       seqs[i] = new Sequence(ref.getName(), g_seqs[i].toString(),
 522               ref.getStart() + alseqs[i].start + bounds[0], ref.getStart()
 523                       + alseqs[i].start + (bounds[2] == 0 ? -1 : bounds[2]));
 524       seqs[i].setDatasetSequence(ref);
 525       seqs[i].setDescription(ref.getDescription());
 526     }
 527     if (segments != null)
 528     {
 529       for (int i = 0; i < segments.length; i += 3)
 530       {
 531         // int start=shifts.shift(segments[i]-1)+1;
 532         // int end=shifts.shift(segments[i]+segments[i+1]-1)-1;
 533         colsel.hideColumns(segments[i + 1], segments[i + 1]
 534                 + segments[i + 2] - 1);
 535       }
 536     }
 537     return seqs;
 538   }
 539
 540   /**
 541    * non rigorous testing
 542    */
 543   /**
 544    *
 545    * @param seq
 546    *          Sequence
 547    * @param ex_cs_gapped
 548    *          String
 549    * @return String
 550    */
 551   public static String testCigar_string(Sequence seq, String ex_cs_gapped)
 552   {
 553     SeqCigar c_sgapped = new SeqCigar(seq);
 554     String cs_gapped = c_sgapped.getCigarstring();
 555     if (!cs_gapped.equals(ex_cs_gapped))
 556     {
 557       System.err.println("Failed getCigarstring: incorect string '"
 558               + cs_gapped + "' != " + ex_cs_gapped);
 559     }
 560     return cs_gapped;
 561   }
 562
 563   public static boolean testSeqRecovery(SeqCigar gen_sgapped,
 564           SequenceI s_gapped)
 565   {
 566     // this is non-rigorous - start and end recovery is not tested.
 567     SequenceI gen_sgapped_s = gen_sgapped.getSeq('-');
 568     if (!gen_sgapped_s.getSequence().equals(s_gapped.getSequence()))
 569     {
 570       System.err.println("Couldn't reconstruct sequence.\n"
 571               + gen_sgapped_s.getSequenceAsString() + "\n"
 572               + s_gapped.getSequenceAsString());
 573       return false;
 574     }
 575     return true;
 576   }
 577
 578   public static void main(String argv[]) throws Exception
 579   {
 580     String o_seq;
 581     Sequence s = new Sequence("MySeq",
 582             o_seq = "asdfktryasdtqwrtsaslldddptyipqqwaslchvhttt", 39, 80);
 583     String orig_gapped;
 584     Sequence s_gapped = new Sequence(
 585             "MySeq",
 586             orig_gapped = "----asdf------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhttt",
 587             39, 80);
 588     String ex_cs_gapped = "4I4M6I6M3I11M4I12M4I9M";
 589     s_gapped.setDatasetSequence(s);
 590     String sub_gapped_s;
 591     Sequence s_subsequence_gapped = new Sequence(
 592             "MySeq",
 593             sub_gapped_s = "------ktryas---dtqwrtsasll----dddptyipqqwa----slchvh",
 594             43, 77);
 595
 596     s_subsequence_gapped.setDatasetSequence(s);
 597     SeqCigar c_null = new SeqCigar(s);
 598     String cs_null = c_null.getCigarstring();
 599     if (!cs_null.equals("42M"))
 600     {
 601       System.err
 602               .println("Failed to recover ungapped sequence cigar operations:"
 603                       + ((cs_null == "") ? "empty string" : cs_null));
 604     }
 605     testCigar_string(s_gapped, ex_cs_gapped);
 606     SeqCigar gen_sgapped = SeqCigar.parseCigar(s, ex_cs_gapped);
 607     if (!gen_sgapped.getCigarstring().equals(ex_cs_gapped))
 608     {
 609       System.err.println("Failed parseCigar(" + ex_cs_gapped
 610               + ")->getCigarString()->'" + gen_sgapped.getCigarstring()
 611               + "'");
 612     }
 613     testSeqRecovery(gen_sgapped, s_gapped);
 614     // Test dataset resolution
 615     SeqCigar sub_gapped = new SeqCigar(s_subsequence_gapped);
 616     if (!testSeqRecovery(sub_gapped, s_subsequence_gapped))
 617     {
 618       System.err
 619               .println("Failed recovery for subsequence of dataset sequence");
 620     }
 621     // width functions
 622     if (sub_gapped.getWidth() != sub_gapped_s.length())
 623     {
 624       System.err.println("Failed getWidth()");
 625     }
 626
 627     sub_gapped.getFullWidth();
 628     if (sub_gapped.hasDeletedRegions())
 629     {
 630       System.err.println("hasDeletedRegions is incorrect.");
 631     }
 632     // Test start-end region SeqCigar
 633     SeqCigar sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 634     if (sub_se_gp.getWidth() != 41)
 635     {
 636       System.err
 637               .println("SeqCigar(seq, start, end) not properly clipped alignsequence.");
 638     }
 639     System.out.println("Original sequence align:\n" + sub_gapped_s
 640             + "\nReconstructed window from 8 to 48\n" + "XXXXXXXX"
 641             + sub_se_gp.getSequenceString('-') + "..." + "\nCigar String:"
 642             + sub_se_gp.getCigarstring() + "\n");
 643     SequenceI ssgp = sub_se_gp.getSeq('-');
 644     System.out.println("\t " + ssgp.getSequenceAsString());
 645     for (int r = 0; r < 10; r++)
 646     {
 647       sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 648       int sl = sub_se_gp.getWidth();
 649       int st = sl - 1 - r;
 650       for (int rs = 0; rs < 10; rs++)
 651       {
 652         int e = st + rs;
 653         sub_se_gp.deleteRange(st, e);
 654         String ssgapedseq = sub_se_gp.getSeq('-').getSequenceAsString();
 655         System.out.println(st + "," + e + "\t:" + ssgapedseq);
 656         st -= 3;
 657       }
 658     }
 659     {
 660       SeqCigar[] set = new SeqCigar[]
 661       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 662           new SeqCigar(s_gapped) };
 663       Alignment al = new Alignment(set);
 664       for (int i = 0; i < al.getHeight(); i++)
 665       {
 666         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 667                 + al.getSequenceAt(i).getStart() + "\t"
 668                 + al.getSequenceAt(i).getEnd() + "\t"
 669                 + al.getSequenceAt(i).getSequenceAsString());
 670       }
 671     }
 672     {
 673       System.out.println("Gapped.");
 674       SeqCigar[] set = new SeqCigar[]
 675       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 676           new SeqCigar(s_gapped) };
 677       set[0].deleteRange(20, 25);
 678       Alignment al = new Alignment(set);
 679       for (int i = 0; i < al.getHeight(); i++)
 680       {
 681         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 682                 + al.getSequenceAt(i).getStart() + "\t"
 683                 + al.getSequenceAt(i).getEnd() + "\t"
 684                 + al.getSequenceAt(i).getSequenceAsString());
 685       }
 686     }
 687     // if (!ssgapedseq.equals("ryas---dtqqwa----slchvh"))
 688     // System.err.println("Subseqgaped\n------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhryas---dtqwrtsasll--qwa----slchvh\n"+ssgapedseq+"\n"+sub_se_gp.getCigarstring());
 689   }
 690
 691 }