f65f068d55d0ff18843e34702560ef3910392f86
[jalview.git] / src / jalview / datamodel / AlignedCodonFrame.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2)
3  * Copyright (C) 2014 The Jalview Authors
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.datamodel;
22
23 import jalview.util.MapList;
24
25 import java.util.ArrayList;
26 import java.util.List;
27
28 /**
29  * Stores mapping between the columns of a protein alignment and a DNA alignment
30  * and a list of individual codon to amino acid mappings between sequences.
31  */
32 public class AlignedCodonFrame
33 {
34
35   /**
36    * array of nucleotide positions for aligned codons at column of aligned
37    * proteins.
38    */
39   public int[][] codons = null;
40
41   /**
42    * width of protein sequence alignement implicit assertion that codons.length
43    * >= aaWidth
44    */
45   public int aaWidth = 0;
46
47   /*
48    * TODO: not an ideal solution - we reference the aligned amino acid sequences
49    * in order to make insertions on them Better would be dnaAlignment and
50    * aaAlignment reference....
51    */
52   private List<SequenceI> a_aaSeqs = new ArrayList<SequenceI>();
53
54   /*
55    * tied array of na Sequence objects.
56    */
57   private SequenceI[] dnaSeqs = null;
58
59   /*
60    * tied array of Mappings to protein sequence Objects and SequenceI[]
61    * aaSeqs=null; MapLists where eac maps from the corresponding dnaSeqs element
62    * to corresponding aaSeqs element
63    */
64   private Mapping[] dnaToProt = null;
65
66   /**
67    * initialise codon frame with a nominal alignment width
68    * 
69    * @param aWidth
70    */
71   public AlignedCodonFrame(int aWidth)
72   {
73     if (aWidth <= 0)
74     {
75       codons = null;
76       return;
77     }
78     codons = new int[aWidth][];
79     for (int res = 0; res < aWidth; res++)
80     {
81       codons[res] = null;
82     }
83   }
84
85   /**
86    * Construct a 'near copy' of the given AlignedCodonFrame, that references the
87    * same dataset sequences, but the given protein aligned sequences.
88    * 
89    * @param acf
90    * @param alignment
91    * @throws IllegalStateException
92    *           if the copied mapping references any dataset not in the alignment
93    */
94   public AlignedCodonFrame(AlignedCodonFrame acf, SequenceI[] alignment)
95   {
96     this.codons = acf.codons;
97     this.dnaSeqs = acf.dnaSeqs;
98     this.dnaToProt = acf.dnaToProt;
99
100     for (SequenceI seq : acf.a_aaSeqs)
101     {
102       boolean found = false;
103       for (SequenceI newseq : alignment)
104       {
105         if (seq.getDatasetSequence() == newseq.getDatasetSequence())
106         {
107           this.a_aaSeqs.add(newseq);
108           found = true;
109           break;
110         }
111       }
112       if (!found)
113       {
114         throw new IllegalStateException("Copying codon mapping for"
115                 + seq.getSequenceAsString());
116       }
117     }
118   }
119
120   /**
121    * ensure that codons array is at least as wide as aslen residues
122    * 
123    * @param aslen
124    * @return (possibly newly expanded) codon array
125    */
126   public int[][] checkCodonFrameWidth(int aslen)
127   {
128     if (codons.length <= aslen + 1)
129     {
130       // probably never have to do this ?
131       int[][] c = new int[codons.length + 10][];
132       for (int i = 0; i < codons.length; i++)
133       {
134         c[i] = codons[i];
135         codons[i] = null;
136       }
137       codons = c;
138     }
139     return codons;
140   }
141
142   /**
143    * @return width of aligned translated amino acid residues
144    */
145   public int getaaWidth()
146   {
147     return aaWidth;
148   }
149
150   /**
151    * increase aaWidth by one and insert a new aligned codon position space at
152    * aspos.
153    * 
154    * @param aspos
155    */
156   public void insertAAGap(int aspos, char gapCharacter)
157   {
158     // this aa appears before the aligned codons at aspos - so shift them in
159     // each pair of mapped sequences
160     aaWidth++;
161     // we actually have to modify the aligned sequences here, so use the
162     // a_aaSeqs vector
163     for (SequenceI seq : a_aaSeqs)
164     {
165       seq.insertCharAt(aspos, gapCharacter);
166     }
167
168     checkCodonFrameWidth(aspos);
169     if (aspos < aaWidth)
170     {
171       aaWidth++;
172       System.arraycopy(codons, aspos, codons, aspos + 1, codons.length
173               - aspos - 1);
174       codons[aspos] = null; // clear so new codon position can be marked.
175     }
176   }
177
178   public void setAaWidth(int aapos)
179   {
180     aaWidth = aapos;
181   }
182
183   /**
184    * add a mapping between the dataset sequences for the associated dna and
185    * protein sequence objects
186    * 
187    * @param dnaseq
188    * @param aaseq
189    * @param map
190    */
191   public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map)
192   {
193     int nlen = 1;
194     if (dnaSeqs != null)
195     {
196       nlen = dnaSeqs.length + 1;
197     }
198     SequenceI[] ndna = new SequenceI[nlen];
199     Mapping[] ndtp = new Mapping[nlen];
200     if (dnaSeqs != null)
201     {
202       System.arraycopy(dnaSeqs, 0, ndna, 0, dnaSeqs.length);
203       System.arraycopy(dnaToProt, 0, ndtp, 0, dnaSeqs.length);
204     }
205     dnaSeqs = ndna;
206     dnaToProt = ndtp;
207     nlen--;
208     dnaSeqs[nlen] = (dnaseq.getDatasetSequence() == null) ? dnaseq : dnaseq
209             .getDatasetSequence();
210     Mapping mp = new Mapping(map);
211     // JBPNote DEBUG! THIS !
212     // dnaseq.transferAnnotation(aaseq, mp);
213     // aaseq.transferAnnotation(dnaseq, new Mapping(map.getInverse()));
214     mp.to = (aaseq.getDatasetSequence() == null) ? aaseq : aaseq
215             .getDatasetSequence();
216     a_aaSeqs.add(aaseq);
217     dnaToProt[nlen] = mp;
218   }
219
220   public SequenceI[] getdnaSeqs()
221   {
222     return dnaSeqs;
223   }
224
225   public SequenceI[] getAaSeqs()
226   {
227     if (dnaToProt == null)
228     {
229       return null;
230     }
231     SequenceI[] sqs = new SequenceI[dnaToProt.length];
232     for (int sz = 0; sz < dnaToProt.length; sz++)
233     {
234       sqs[sz] = dnaToProt[sz].to;
235     }
236     return sqs;
237   }
238
239   public MapList[] getdnaToProt()
240   {
241     if (dnaToProt == null)
242     {
243       return null;
244     }
245     MapList[] sqs = new MapList[dnaToProt.length];
246     for (int sz = 0; sz < dnaToProt.length; sz++)
247     {
248       sqs[sz] = dnaToProt[sz].map;
249     }
250     return sqs;
251   }
252
253   public Mapping[] getProtMappings()
254   {
255     return dnaToProt;
256   }
257
258   /**
259    * 
260    * @param sequenceRef
261    * @return null or corresponding aaSeq dataset sequence for dnaSeq entry
262    */
263   public SequenceI getAaForDnaSeq(SequenceI dnaSeqRef)
264   {
265     return getAaForDnaSeq(dnaSeqRef, true);
266   }
267
268   /**
269    * Return the corresponding aligned or dataset aa sequence for given dna
270    * sequence, null if not found.
271    * 
272    * @param sequenceRef
273    * @param returnDataset
274    *          if true, return the aa dataset, else the aligned sequence
275    * @return
276    */
277   public SequenceI getAaForDnaSeq(SequenceI dnaSeqRef, boolean returnDataset)
278   {
279     if (dnaSeqs == null)
280     {
281       return null;
282     }
283     SequenceI dnads = dnaSeqRef.getDatasetSequence();
284     for (int ds = 0; ds < dnaSeqs.length; ds++)
285     {
286       if (dnaSeqs[ds] == dnaSeqRef || dnaSeqs[ds] == dnads)
287       {
288         if (returnDataset)
289         {
290           return dnaToProt[ds].to;
291         }
292         else
293         {
294           // TODO very fragile - depends on dnaSeqs, dnaToProt, a_aaSeqs moving
295           // in parallel; revise data model to guarantee this
296           return a_aaSeqs.get(ds);
297         }
298       }
299     }
300     return null;
301   }
302
303   /**
304    * 
305    * @param sequenceRef
306    * @return null or corresponding aaSeq entry for dnaSeq entry
307    */
308   public SequenceI getDnaForAaSeq(SequenceI aaSeqRef)
309   {
310     if (dnaToProt == null)
311     {
312       return null;
313     }
314     SequenceI aads = aaSeqRef.getDatasetSequence();
315     for (int as = 0; as < dnaToProt.length; as++)
316     {
317       if (dnaToProt[as].to == aaSeqRef || dnaToProt[as].to == aads)
318       {
319         return dnaSeqs[as];
320       }
321     }
322     return null;
323   }
324
325   /**
326    * test to see if codon frame involves seq in any way
327    * 
328    * @param seq
329    *          a nucleotide or protein sequence
330    * @return true if a mapping exists to or from this sequence to any translated
331    *         sequence
332    */
333   public boolean involvesSequence(SequenceI seq)
334   {
335     return getAaForDnaSeq(seq) != null || getDnaForAaSeq(seq) != null;
336   }
337
338   /**
339    * Add search results for regions in other sequences that translate or are
340    * translated from a particular position in seq
341    * 
342    * @param seq
343    * @param index
344    *          position in seq
345    * @param results
346    *          where highlighted regions go
347    */
348   public void markMappedRegion(SequenceI seq, int index,
349           SearchResults results)
350   {
351     if (dnaToProt == null)
352     {
353       return;
354     }
355     int[] codon;
356     SequenceI ds = seq.getDatasetSequence();
357     for (int mi = 0; mi < dnaToProt.length; mi++)
358     {
359       if (dnaSeqs[mi] == seq || dnaSeqs[mi] == ds)
360       {
361         // DEBUG System.err.println("dna pos "+index);
362         codon = dnaToProt[mi].map.locateInTo(index, index);
363         if (codon != null)
364         {
365           for (int i = 0; i < codon.length; i += 2)
366           {
367             results.addResult(dnaToProt[mi].to, codon[i], codon[i + 1]);
368           }
369         }
370       }
371       else if (dnaToProt[mi].to == seq || dnaToProt[mi].to == ds)
372       {
373         // DEBUG System.err.println("aa pos "+index);
374         {
375           codon = dnaToProt[mi].map.locateInFrom(index, index);
376           if (codon != null)
377           {
378             for (int i = 0; i < codon.length; i += 2)
379             {
380               results.addResult(dnaSeqs[mi], codon[i], codon[i + 1]);
381             }
382           }
383         }
384       }
385     }
386   }
387
388   /**
389    * Returns the DNA codon positions (base 1) for the given position (base 1) in
390    * a mapped protein sequence, or null if no mapping is found.
391    * 
392    * Intended for use in aligning cDNA to match aligned protein. Only the first
393    * mapping found is returned, so not suitable for use if multiple protein
394    * sequences are mapped to the same cDNA (but aligning cDNA as protein is
395    * ill-defined for this case anyway).
396    * 
397    * @param seq
398    *          the DNA dataset sequence
399    * @param aaPos
400    *          residue position (base 1) in a protein sequence
401    * @return
402    */
403   public int[] getDnaPosition(SequenceI seq, int aaPos)
404   {
405     /*
406      * Adapted from markMappedRegion().
407      */
408     MapList ml = null;
409     for (int i = 0; i < dnaToProt.length; i++)
410     {
411       if (dnaSeqs[i] == seq)
412       {
413         ml = getdnaToProt()[i];
414         break;
415       }
416     }
417     return ml == null ? null : ml.locateInFrom(aaPos, aaPos);
418   }
419 }