Send gapChar to print original Data
[jalview.git] / src / jalview / datamodel / SeqCigar.java
1 package jalview.datamodel;
2
3 import jalview.analysis.AlignSeq;
4
5 public class SeqCigar
6     extends CigarSimple
7 {
8
9   private SequenceI refseq=null;
10   /**
11    * Reference dataset sequence for the cigar string
12    * @return SequenceI
13    */
14   public SequenceI getRefSeq() {
15     return refseq;
16   }
17   /**
18    * Returns sequence as a string with cigar operations applied to it
19    * @return String
20    */
21   public String getSequenceString(char GapChar)
22   {
23     return (length==0) ? "" : (String) getSequenceAndDeletions(refseq.getSequence(), GapChar)[0];
24   }
25
26   /**
27    * recreates a gapped and edited version of RefSeq or null for an empty cigar string
28    * @return SequenceI
29    */
30   public SequenceI getSeq(char GapChar) {
31     Sequence seq;
32     if (refseq==null || length==0)
33       return null;
34     Object[] edit_result=getSequenceAndDeletions(refseq.getSequence(), GapChar);
35     if (edit_result==null)
36       throw new Error("Implementation Error - unexpected null from getSequenceAndDeletions");
37
38     seq = new Sequence(refseq.getName(), (String) edit_result[0], refseq.getStart()+((int[]) edit_result[1])[0], refseq.getStart()+((int[]) edit_result[1])[2]);
39     seq.setDatasetSequence(refseq);
40     return seq;
41   }
42   /*
43   We don't allow this - refseq is given at construction time only
44    public void setSeq(SequenceI seq) {
45     this.seq = seq;
46   }
47   */
48   /**
49    * internal constructor - sets seq to a gapless sequence derived from seq
50    * and prepends any 'D' operations needed to get to the first residue of seq.
51    * @param seq SequenceI
52    * @return true if gaps are present in seq
53    */
54   private boolean _setSeq(SequenceI seq) {
55     boolean hasgaps=false;
56
57     if (seq==null)
58       throw new Error("Implementation Error - _setSeq(null)");
59
60     // Find correct sequence to reference and add initial hidden offset
61     SequenceI ds = seq.getDatasetSequence();
62     if (ds==null) {
63       ds = new Sequence(seq.getName(),
64                         AlignSeq.extractGaps(jalview.util.Comparison.GapChars, new String(seq.getSequence())),
65                         seq.getStart(),
66                         seq.getEnd());
67     }
68     // check that we haven't just duplicated an ungapped sequence.
69       if (ds.getLength()==seq.getLength()) {
70         ds = seq;
71       } else {
72         hasgaps = true;
73       }
74     this.refseq = ds;
75     // Adjust offset
76     if (ds.getStart()<seq.getStart()) {
77       addDeleted(seq.getStart()-ds.getStart());
78     }
79     return hasgaps;
80   }
81   /**
82    * directly initialise a cigar object with a sequence of range, operation pairs and a sequence to apply it to.
83    * operation and range should be relative to the seq.getStart()'th residue of the dataset seq resolved from seq.
84    * @param seq SequenceI
85    * @param operation char[]
86    * @param range int[]
87    */
88   public SeqCigar(SequenceI seq, char operation[], int range[]) {
89     super();
90     if (seq==null)
91       throw new Error("Implementation Bug. Null seq !");
92     if (operation.length!=range.length) {
93       throw new Error("Implementation Bug. Cigar Operation list!= range list");
94     }
95
96     if (operation!=null) {
97       this.operation = new char[operation.length+_inc_length];
98       this.range = new int[operation.length+_inc_length];
99
100       if (_setSeq(seq)) {
101         throw new Error("NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
102       }
103       for (int i = this.length, j=0; j < operation.length; i++,j++)
104       {
105         char op = operation[j];
106         if (op != M && op != I && op != D)
107         {
108           throw new Error(
109               "Implementation Bug. Cigar Operation '"+j+"' '"+op+"' not one of '"+M+"', '"+I+"', or '"+D+"'.");
110         }
111         this.operation[i] = op;
112         this.range[i] = range[j];
113       }
114       this.length+=operation.length;
115     } else {
116       this.operation = null;
117       this.range = null;
118       this.length=0;
119       if (_setSeq(seq)) {
120         throw new Error("NOT YET Implemented: Constructing a Cigar object from a cigar string and a gapped sequence.");
121       }
122     }
123   }
124   /**
125    * add range matched residues to cigar string
126    * @param range int
127    */
128   public void addMatch(int range) {
129     this.addOperation(M, range);
130   }
131   /**
132    * Deleted regions mean that there will be discontinuous sequence numbering in the
133    * sequence returned by getSeq(char).
134    * @return true if there are non-terminal deletions
135    */
136   public boolean hasDeletedRegions() {
137     for (int i=1, l=length-1; i<l; i++) {
138       if (operation[i]==D)
139         return true;
140     }
141     return false;
142   }
143   protected static void addSequenceOps(CigarBase cigar, SequenceI seq, int startpos, int endpos) {
144   char op = '\0';
145   int range = 0;
146   int p = 0, res = seq.getLength();
147   startpos++;
148   endpos++;
149   while (p<res)
150   {
151     boolean isGap = jalview.util.Comparison.isGap(seq.getCharAt(p++));
152     if ( (startpos <= p) && (p < endpos))
153     {
154       if (isGap)
155       {
156         if (range > 0 && op != I)
157         {
158           cigar.addOperation(op, range);
159           range = 0;
160         }
161         op = I;
162         range++;
163       }
164       else
165       {
166         if (range > 0 && op != M)
167         {
168           cigar.addOperation(op, range);
169           range = 0;
170         }
171         op = M;
172         range++;
173       }
174     }
175     else
176     {
177       if (!isGap)
178       {
179         if (range > 0 && op != D)
180         {
181           cigar.addOperation(op, range);
182           range = 0;
183         }
184         op = D;
185         range++;
186       }
187       else
188       {
189         // do nothing - insertions are not recorded in flanking regions.
190       }
191     }
192   }
193   if (range > 0)
194   {
195     cigar.addOperation(op, range);
196   }
197   }
198   /**
199    * create a cigar string for given sequence
200    * @param seq SequenceI
201    */
202   public SeqCigar(SequenceI seq) {
203     super();
204     if (seq == null)
205       throw new Error("Implementation error for new Cigar(SequenceI)");
206     if (_setSeq(seq))
207     {
208         // there is still work to do
209       addSequenceOps(this, seq, 0, seq.getLength());
210     }
211   }
212   public SeqCigar(SequenceI seq, int start, int end) {
213     super();
214     if (seq == null)
215       throw new Error("Implementation error for new Cigar(SequenceI)");
216     if (_setSeq(seq))
217     {
218       // there is still work to do
219       addSequenceOps(this, seq, start, end);
220     }
221   }
222
223   /**
224    * Create a cigar object from a cigar string like '[<I|D|M><range>]+'
225    * Will fail if the given seq already contains gaps (JBPNote: future implementation will fix)
226    * @param seq SequenceI object resolvable to a dataset sequence
227    * @param cigarString String
228    * @return Cigar
229    */
230   public static SeqCigar parseCigar(SequenceI seq, String cigarString)
231       throws Exception
232   {
233     Object[] opsandrange = parseCigarString(cigarString);
234     return new SeqCigar(seq, (char[]) opsandrange[0], (int[]) opsandrange[1]);
235   }
236   /**
237    * non rigorous testing
238    */
239   /**
240    *
241    * @param seq Sequence
242    * @param ex_cs_gapped String
243    * @return String
244    */
245   public static String testCigar_string(Sequence seq, String ex_cs_gapped) {
246     SeqCigar c_sgapped = new SeqCigar(seq);
247     String cs_gapped = c_sgapped.getCigarstring();
248     if (!cs_gapped.equals(ex_cs_gapped))
249       System.err.println("Failed getCigarstring: incorect string '"+cs_gapped+"' != "+ex_cs_gapped);
250     return cs_gapped;
251   }
252   public static boolean testSeqRecovery(SeqCigar gen_sgapped, SequenceI s_gapped) {
253     SequenceI gen_sgapped_s = gen_sgapped.getSeq('-');
254     if (!gen_sgapped_s.getSequence().equals(s_gapped.getSequence())) {
255       System.err.println("Couldn't reconstruct sequence.\n" +
256                          gen_sgapped_s.getSequence() + "\n" +
257                          s_gapped.getSequence());
258       return false;
259     }
260     return true;
261   }
262   public static void main(String argv[]) throws Exception {
263     Sequence s=new Sequence("MySeq", "asdfktryasdtqwrtsaslldddptyipqqwaslchvhttt",39,80);
264     String orig_gapped;
265     Sequence s_gapped=new Sequence("MySeq", orig_gapped="----asdf------ktryas---dtqwrtsasll----dddptyipqqwa----slchvhttt", 39,80);
266     String ex_cs_gapped="4I4M6I6M3I11M4I12M4I9M";
267     s_gapped.setDatasetSequence(s);
268     String sub_gapped_s;
269     Sequence s_subsequence_gapped=new Sequence("MySeq", sub_gapped_s="------ktryas---dtqwrtsasll----dddptyipqqwa----slchvh", 43,77);
270
271     s_subsequence_gapped.setDatasetSequence(s);
272     SeqCigar c_null = new SeqCigar(s);
273     String cs_null = c_null.getCigarstring();
274     if (cs_null.length()>0)
275       System.err.println("Failed getCigarstring: Unexpected cigar operations:"+cs_null);
276     testCigar_string(s_gapped, ex_cs_gapped);
277     SeqCigar gen_sgapped = SeqCigar.parseCigar(s, ex_cs_gapped);
278     if (!gen_sgapped.getCigarstring().equals(ex_cs_gapped))
279       System.err.println("Failed parseCigar("+ex_cs_gapped+")->getCigarString()->'"+gen_sgapped.getCigarstring()+"'");
280     testSeqRecovery(gen_sgapped, s_gapped);
281     // Test dataset resolution
282     SeqCigar sub_gapped = new SeqCigar(s_subsequence_gapped);
283     if (!testSeqRecovery(sub_gapped, s_subsequence_gapped))
284         System.err.println("Failed recovery for subsequence of dataset sequence");
285     // width functions
286     if (sub_gapped.getWidth()!=sub_gapped_s.length())
287       System.err.println("Failed getWidth()");
288
289     sub_gapped.getFullWidth();
290     if (sub_gapped.hasDeletedRegions())
291       System.err.println("hasDeletedRegions is incorrect.");
292     // Test start-end region SeqCigar
293     SeqCigar sub_se_gp= new SeqCigar(s_subsequence_gapped, 8, 48);
294     if (sub_se_gp.getWidth()!=40)
295       System.err.println("SeqCigar(seq, start, end) not properly clipped alignsequence.");
296     System.out.println("Original sequence align:\n"+sub_gapped_s+"\nReconstructed window from 8 to 48\n"+"XXXXXXXX"+sub_se_gp.getSequenceString('-')+"...."+"\nCigar String:"+sub_se_gp.getCigarstring()+"");
297   }
298
299 }