src/jalview/datamodel/SeqCigar.java

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