src/jalview/datamodel/SeqCigar.java

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