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