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(refseq
  72             .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(refseq
  89             .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], refseq
  97             .getStart()
  98             + start + bounds[0], refseq.getStart() + start
  99             + ((bounds[2] == 0) ? -1 : bounds[2]));
 100     seq.setDescription(refseq.getDescription());
 101     int sstart = seq.getStart(), send = seq.getEnd();
 102     // seq.checkValidRange(); probably not needed
 103     // recover local properties if present
 104     if (seqProps != null)
 105     {
 106       // this recovers dataset sequence reference as well as local features,
 107       // names, start/end settings.
 108       SeqsetUtils.SeqCharacterUnhash(seq, seqProps);
 109     }
 110     // ensure dataset sequence is up to date from local reference
 111     seq.setDatasetSequence(refseq);
 112     seq.setStart(sstart);
 113     seq.setEnd(send);
 114     return seq;
 115   }
 116
 117   /*
 118    * We don't allow this - refseq is given at construction time only public void
 119    * setSeq(SequenceI seq) { this.seq = seq; }
 120    */
 121   /**
 122    * internal constructor - sets seq to a gapless sequence derived from seq and
 123    * prepends any 'D' operations needed to get to the first residue of seq.
 124    *
 125    * @param seq
 126    *          SequenceI
 127    * @param initialDeletion
 128    *          true to mark initial dataset sequence residues as deleted in
 129    *          subsequence
 130    * @param _s
 131    *          index of first position in seq
 132    * @param _e
 133    *          index after last position in (possibly gapped) seq
 134    * @return true if gaps are present in seq
 135    */
 136   private boolean _setSeq(SequenceI seq, boolean initialDeletion, int _s,
 137           int _e)
 138   {
 139     boolean hasgaps = false;
 140     if (seq == null)
 141     {
 142       throw new Error("Implementation Error - _setSeq(null,...)");
 143     }
 144     if (_s < 0)
 145     {
 146       throw new Error("Implementation Error: _s=" + _s);
 147     }
 148     String seq_string = seq.getSequenceAsString();
 149     if (_e == 0 || _e < _s || _e > seq_string.length())
 150     {
 151       _e = seq_string.length();
 152     }
 153     // resolve start and end positions relative to ungapped reference sequence
 154     start = seq.findPosition(_s) - seq.getStart();
 155     end = seq.findPosition(_e) - seq.getStart();
 156     int l_ungapped = end - start;
 157     // Find correct sequence to reference and correct start and end - if
 158     // necessary
 159     SequenceI ds = seq.getDatasetSequence();
 160     if (ds == null)
 161     {
 162       // make a new dataset sequence
 163       String ungapped = AlignSeq.extractGaps(
 164               jalview.util.Comparison.GapChars, new String(seq_string));
 165       l_ungapped = ungapped.length();
 166       // check that we haven't just duplicated an ungapped sequence.
 167       if (l_ungapped == seq.getLength())
 168       {
 169         ds = seq;
 170       }
 171       else
 172       {
 173         ds = new Sequence(seq.getName(), ungapped, seq.getStart(), seq
 174                 .getStart()
 175                 + ungapped.length() - 1);
 176         // JBPNote: this would be consistent but may not be useful
 177         // seq.setDatasetSequence(ds);
 178       }
 179     }
 180     // add in offset between seq and the dataset sequence
 181     if (ds.getStart() < seq.getStart())
 182     {
 183       int offset = seq.getStart() - ds.getStart();
 184       if (initialDeletion)
 185       {
 186         // absolute cigar string
 187         addDeleted(_s + offset);
 188         start = 0;
 189         end += offset;
 190       }
 191       else
 192       {
 193         // normal behaviour - just mark start and end subsequence
 194         start += offset;
 195         end += offset;
 196
 197       }
 198
 199     }
 200
 201     // any gaps to process ?
 202     if (l_ungapped != (_e - _s))
 203     {
 204       hasgaps = true;
 205     }
 206
 207     refseq = ds;
 208     // copy over local properties for the sequence instance of the refseq
 209     seqProps = SeqsetUtils.SeqCharacterHash(seq);
 210     // Check offsets
 211     if (end > ds.getLength())
 212     {
 213       throw new Error(
 214               "SeqCigar: Possible implementation error: sequence is longer than dataset sequence");
 215       // end = ds.getLength();
 216     }
 217
 218     return hasgaps;
 219   }
 220
 221   /**
 222    * directly initialise a cigar object with a sequence of range, operation
 223    * pairs and a sequence to apply it to. operation and range should be relative
 224    * to the seq.getStart()'th residue of the dataset seq resolved from seq.
 225    *
 226    * @param seq
 227    *          SequenceI
 228    * @param operation
 229    *          char[]
 230    * @param range
 231    *          int[]
 232    */
 233   public SeqCigar(SequenceI seq, char operation[], int range[])
 234   {
 235     super();
 236     if (seq == null)
 237     {
 238       throw new Error("Implementation Bug. Null seq !");
 239     }
 240     if (operation.length != range.length)
 241     {
 242       throw new Error(
 243               "Implementation Bug. Cigar Operation list!= range list");
 244     }
 245
 246     if (operation != null)
 247     {
 248       this.operation = new char[operation.length + _inc_length];
 249       this.range = new int[operation.length + _inc_length];
 250
 251       if (_setSeq(seq, false, 0, 0))
 252       {
 253         throw new Error(
 254                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 255       }
 256       for (int i = this.length, j = 0; j < operation.length; i++, j++)
 257       {
 258         char op = operation[j];
 259         if (op != M && op != I && op != D)
 260         {
 261           throw new Error("Implementation Bug. Cigar Operation '" + j
 262                   + "' '" + op + "' not one of '" + M + "', '" + I
 263                   + "', or '" + D + "'.");
 264         }
 265         this.operation[i] = op;
 266         this.range[i] = range[j];
 267       }
 268       this.length += operation.length;
 269     }
 270     else
 271     {
 272       this.operation = null;
 273       this.range = null;
 274       this.length = 0;
 275       if (_setSeq(seq, false, 0, 0))
 276       {
 277         throw new Error(
 278                 "NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
 279       }
 280     }
 281   }
 282
 283   /**
 284    * add range matched residues to cigar string
 285    *
 286    * @param range
 287    *          int
 288    */
 289   public void addMatch(int range)
 290   {
 291     this.addOperation(M, range);
 292   }
 293
 294   /**
 295    * Adds insertion and match operations based on seq to the cigar up to the
 296    * endpos column of seq.
 297    *
 298    * @param cigar
 299    *          CigarBase
 300    * @param seq
 301    *          SequenceI
 302    * @param startpos
 303    *          int
 304    * @param endpos
 305    *          int
 306    * @param initialDeletions
 307    *          if true then initial deletions will be added from start of seq to
 308    *          startpos
 309    */
 310   protected static void addSequenceOps(CigarBase cigar, SequenceI seq,
 311           int startpos, int endpos, boolean initialDeletions)
 312   {
 313     char op = '\0';
 314     int range = 0;
 315     int p = 0, res = seq.getLength();
 316
 317     if (!initialDeletions)
 318     {
 319       p = startpos;
 320     }
 321
 322     while (p <= endpos)
 323     {
 324       boolean isGap = (p < res) ? jalview.util.Comparison.isGap(seq
 325               .getCharAt(p)) : true;
 326       if ((startpos <= p) && (p <= endpos))
 327       {
 328         if (isGap)
 329         {
 330           if (range > 0 && op != I)
 331           {
 332             cigar.addOperation(op, range);
 333             range = 0;
 334           }
 335           op = I;
 336           range++;
 337         }
 338         else
 339         {
 340           if (range > 0 && op != M)
 341           {
 342             cigar.addOperation(op, range);
 343             range = 0;
 344           }
 345           op = M;
 346           range++;
 347         }
 348       }
 349       else
 350       {
 351         if (!isGap)
 352         {
 353           if (range > 0 && op != D)
 354           {
 355             cigar.addOperation(op, range);
 356             range = 0;
 357           }
 358           op = D;
 359           range++;
 360         }
 361         else
 362         {
 363           // do nothing - insertions are not made in flanking regions
 364         }
 365       }
 366       p++;
 367     }
 368     if (range > 0)
 369     {
 370       cigar.addOperation(op, range);
 371     }
 372   }
 373
 374   /**
 375    * create a cigar string for given sequence
 376    *
 377    * @param seq
 378    *          SequenceI
 379    */
 380   public SeqCigar(SequenceI seq)
 381   {
 382     super();
 383     if (seq == null)
 384     {
 385       throw new Error("Implementation error for new Cigar(SequenceI)");
 386     }
 387     _setSeq(seq, false, 0, 0);
 388     // there is still work to do
 389     addSequenceOps(this, seq, 0, seq.getLength() - 1, false);
 390   }
 391
 392   /**
 393    * Create Cigar from a range of gaps and residues on a sequence object
 394    *
 395    * @param seq
 396    *          SequenceI
 397    * @param start
 398    *          int - first column in range
 399    * @param end
 400    *          int - last column in range
 401    */
 402   public SeqCigar(SequenceI seq, int start, int end)
 403   {
 404     super();
 405     if (seq == null)
 406     {
 407       throw new Error("Implementation error for new Cigar(SequenceI)");
 408     }
 409     _setSeq(seq, false, start, end + 1);
 410     // there is still work to do
 411     addSequenceOps(this, seq, start, end, false);
 412   }
 413
 414   /**
 415    * Create a cigar object from a cigar string like '[<I|D|M><range>]+' Will
 416    * fail if the given seq already contains gaps (JBPNote: future implementation
 417    * will fix)
 418    *
 419    * @param seq
 420    *          SequenceI object resolvable to a dataset sequence
 421    * @param cigarString
 422    *          String
 423    * @return Cigar
 424    */
 425   public static SeqCigar parseCigar(SequenceI seq, String cigarString)
 426           throws Exception
 427   {
 428     Object[] opsandrange = parseCigarString(cigarString);
 429     return new SeqCigar(seq, (char[]) opsandrange[0],
 430             (int[]) opsandrange[1]);
 431   }
 432
 433   /**
 434    * createAlignment
 435    *
 436    * @param alseqs
 437    *          SeqCigar[]
 438    * @param gapCharacter
 439    *          char
 440    * @return SequenceI[]
 441    */
 442   public static SequenceI[] createAlignmentSequences(SeqCigar[] alseqs,
 443           char gapCharacter, ColumnSelection colsel, int[] segments)
 444   {
 445     SequenceI[] seqs = new SequenceI[alseqs.length];
 446     StringBuffer[] g_seqs = new StringBuffer[alseqs.length];
 447     String[] alseqs_string = new String[alseqs.length];
 448     Object[] gs_regions = new Object[alseqs.length];
 449     for (int i = 0; i < alseqs.length; i++)
 450     {
 451       alseqs_string[i] = alseqs[i].getRefSeq().getSequenceAsString(
 452               alseqs[i].start, alseqs[i].end);
 453       gs_regions[i] = alseqs[i].getSequenceAndDeletions(alseqs_string[i],
 454               gapCharacter); // gapped sequence, {start, start col, end.
 455       // endcol}, hidden regions {{start, end, col}})
 456       if (gs_regions[i] == null)
 457       {
 458         throw new Error("Implementation error: " + i
 459                 + "'th sequence Cigar has no operations.");
 460       }
 461       g_seqs[i] = new StringBuffer((String) ((Object[]) gs_regions[i])[0]); // the
 462       // visible
 463       // gapped
 464       // sequence
 465     }
 466     // Now account for insertions. (well - deletions)
 467     // this is complicated because we must keep track of shifted positions in
 468     // each sequence
 469     ShiftList shifts = new ShiftList();
 470     for (int i = 0; i < alseqs.length; i++)
 471     {
 472       Object[] gs_region = ((Object[]) ((Object[]) gs_regions[i])[2]);
 473       if (gs_region != null)
 474
 475       {
 476         for (int hr = 0; hr < gs_region.length; hr++)
 477         {
 478           int[] region = (int[]) gs_region[hr];
 479           char[] insert = new char[region[1] - region[0] + 1];
 480           for (int s = 0; s < insert.length; s++)
 481           {
 482             insert[s] = gapCharacter;
 483           }
 484           int inspos = shifts.shift(region[2]); // resolve insertion position in
 485           // current alignment frame of
 486           // reference
 487           for (int s = 0; s < alseqs.length; s++)
 488           {
 489             if (s != i)
 490             {
 491               if (g_seqs[s].length() <= inspos)
 492               {
 493                 // prefix insertion with more gaps.
 494                 for (int l = inspos - g_seqs[s].length(); l > 0; l--)
 495                 {
 496                   g_seqs[s].append(gapCharacter); // to debug - use a diffferent
 497                   // gap character here
 498                 }
 499               }
 500               g_seqs[s].insert(inspos, insert);
 501             }
 502             else
 503             {
 504               g_seqs[s].insert(inspos, alseqs_string[i].substring(
 505                       region[0], region[1] + 1));
 506             }
 507           }
 508           shifts.addShift(region[2], insert.length); // update shift in
 509           // alignment frame of
 510           // reference
 511           if (segments == null)
 512           {
 513             // add a hidden column for this deletion
 514             colsel.hideColumns(inspos, inspos + insert.length - 1);
 515           }
 516         }
 517       }
 518     }
 519     for (int i = 0; i < alseqs.length; i++)
 520     {
 521       int[] bounds = ((int[]) ((Object[]) gs_regions[i])[1]);
 522       SequenceI ref = alseqs[i].getRefSeq();
 523       seqs[i] = new Sequence(ref.getName(), g_seqs[i].toString(), ref
 524               .getStart()
 525               + 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    * non rigorous testing
 545    */
 546   /**
 547    *
 548    * @param seq
 549    *          Sequence
 550    * @param ex_cs_gapped
 551    *          String
 552    * @return String
 553    */
 554   public static String testCigar_string(Sequence seq, String ex_cs_gapped)
 555   {
 556     SeqCigar c_sgapped = new SeqCigar(seq);
 557     String cs_gapped = c_sgapped.getCigarstring();
 558     if (!cs_gapped.equals(ex_cs_gapped))
 559     {
 560       System.err.println("Failed getCigarstring: incorect string '"
 561               + cs_gapped + "' != " + ex_cs_gapped);
 562     }
 563     return cs_gapped;
 564   }
 565
 566   public static boolean testSeqRecovery(SeqCigar gen_sgapped,
 567           SequenceI s_gapped)
 568   {
 569     // this is non-rigorous - start and end recovery is not tested.
 570     SequenceI gen_sgapped_s = gen_sgapped.getSeq('-');
 571     if (!gen_sgapped_s.getSequence().equals(s_gapped.getSequence()))
 572     {
 573       System.err.println("Couldn't reconstruct sequence.\n"
 574               + gen_sgapped_s.getSequenceAsString() + "\n"
 575               + s_gapped.getSequenceAsString());
 576       return false;
 577     }
 578     return true;
 579   }
 580
 581   public static void main(String argv[]) throws Exception
 582   {
 583     String o_seq;
 584     Sequence s = new Sequence("MySeq",
 585             o_seq = "asdfktryasdtqwrtsaslldddptyipqqwaslchvhttt", 39, 80);
 586     String orig_gapped;
 587     Sequence s_gapped = new Sequence(
 588             "MySeq",
 589             orig_gapped = "----asdf------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhttt",
 590             39, 80);
 591     String ex_cs_gapped = "4I4M6I6M3I11M4I12M4I9M";
 592     s_gapped.setDatasetSequence(s);
 593     String sub_gapped_s;
 594     Sequence s_subsequence_gapped = new Sequence(
 595             "MySeq",
 596             sub_gapped_s = "------ktryas---dtqwrtsasll----dddptyipqqwa----slchvh",
 597             43, 77);
 598
 599     s_subsequence_gapped.setDatasetSequence(s);
 600     SeqCigar c_null = new SeqCigar(s);
 601     String cs_null = c_null.getCigarstring();
 602     if (!cs_null.equals("42M"))
 603     {
 604       System.err
 605               .println("Failed to recover ungapped sequence cigar operations:"
 606                       + ((cs_null == "") ? "empty string" : cs_null));
 607     }
 608     testCigar_string(s_gapped, ex_cs_gapped);
 609     SeqCigar gen_sgapped = SeqCigar.parseCigar(s, ex_cs_gapped);
 610     if (!gen_sgapped.getCigarstring().equals(ex_cs_gapped))
 611     {
 612       System.err.println("Failed parseCigar(" + ex_cs_gapped
 613               + ")->getCigarString()->'" + gen_sgapped.getCigarstring()
 614               + "'");
 615     }
 616     testSeqRecovery(gen_sgapped, s_gapped);
 617     // Test dataset resolution
 618     SeqCigar sub_gapped = new SeqCigar(s_subsequence_gapped);
 619     if (!testSeqRecovery(sub_gapped, s_subsequence_gapped))
 620     {
 621       System.err
 622               .println("Failed recovery for subsequence of dataset sequence");
 623     }
 624     // width functions
 625     if (sub_gapped.getWidth() != sub_gapped_s.length())
 626     {
 627       System.err.println("Failed getWidth()");
 628     }
 629
 630     sub_gapped.getFullWidth();
 631     if (sub_gapped.hasDeletedRegions())
 632     {
 633       System.err.println("hasDeletedRegions is incorrect.");
 634     }
 635     // Test start-end region SeqCigar
 636     SeqCigar sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 637     if (sub_se_gp.getWidth() != 41)
 638     {
 639       System.err
 640               .println("SeqCigar(seq, start, end) not properly clipped alignsequence.");
 641     }
 642     System.out.println("Original sequence align:\n" + sub_gapped_s
 643             + "\nReconstructed window from 8 to 48\n" + "XXXXXXXX"
 644             + sub_se_gp.getSequenceString('-') + "..." + "\nCigar String:"
 645             + sub_se_gp.getCigarstring() + "\n");
 646     SequenceI ssgp = sub_se_gp.getSeq('-');
 647     System.out.println("\t " + ssgp.getSequenceAsString());
 648     for (int r = 0; r < 10; r++)
 649     {
 650       sub_se_gp = new SeqCigar(s_subsequence_gapped, 8, 48);
 651       int sl = sub_se_gp.getWidth();
 652       int st = sl - 1 - r;
 653       for (int rs = 0; rs < 10; rs++)
 654       {
 655         int e = st + rs;
 656         sub_se_gp.deleteRange(st, e);
 657         String ssgapedseq = sub_se_gp.getSeq('-').getSequenceAsString();
 658         System.out.println(st + "," + e + "\t:" + ssgapedseq);
 659         st -= 3;
 660       }
 661     }
 662     {
 663       SeqCigar[] set = new SeqCigar[]
 664       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 665           new SeqCigar(s_gapped) };
 666       Alignment al = new Alignment(set);
 667       for (int i = 0; i < al.getHeight(); i++)
 668       {
 669         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 670                 + al.getSequenceAt(i).getStart() + "\t"
 671                 + al.getSequenceAt(i).getEnd() + "\t"
 672                 + al.getSequenceAt(i).getSequenceAsString());
 673       }
 674     }
 675     {
 676       System.out.println("Gapped.");
 677       SeqCigar[] set = new SeqCigar[]
 678       { new SeqCigar(s), new SeqCigar(s_subsequence_gapped, 8, 48),
 679           new SeqCigar(s_gapped) };
 680       set[0].deleteRange(20, 25);
 681       Alignment al = new Alignment(set);
 682       for (int i = 0; i < al.getHeight(); i++)
 683       {
 684         System.out.println("" + al.getSequenceAt(i).getName() + "\t"
 685                 + al.getSequenceAt(i).getStart() + "\t"
 686                 + al.getSequenceAt(i).getEnd() + "\t"
 687                 + al.getSequenceAt(i).getSequenceAsString());
 688       }
 689     }
 690     // if (!ssgapedseq.equals("ryas---dtqqwa----slchvh"))
 691     // System.err.println("Subseqgaped\n------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhryas---dtqwrtsasll--qwa----slchvh\n"+ssgapedseq+"\n"+sub_se_gp.getCigarstring());
 692   }
 693
 694 }