src/jalview/datamodel/SeqCigar.java

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