JAL-845 implement alignment of protein to match cDNA alignment
[jalview.git] / src / jalview / util / MappingUtils.java
1 package jalview.util;
2
3 import jalview.analysis.AlignmentSorter;
4 import jalview.api.AlignViewportI;
5 import jalview.commands.CommandI;
6 import jalview.commands.EditCommand;
7 import jalview.commands.EditCommand.Action;
8 import jalview.commands.EditCommand.Edit;
9 import jalview.commands.OrderCommand;
10 import jalview.datamodel.AlignedCodonFrame;
11 import jalview.datamodel.AlignmentI;
12 import jalview.datamodel.AlignmentOrder;
13 import jalview.datamodel.ColumnSelection;
14 import jalview.datamodel.SearchResults;
15 import jalview.datamodel.SearchResults.Match;
16 import jalview.datamodel.Sequence;
17 import jalview.datamodel.SequenceGroup;
18 import jalview.datamodel.SequenceI;
19
20 import java.util.ArrayList;
21 import java.util.HashMap;
22 import java.util.Iterator;
23 import java.util.List;
24 import java.util.Map;
25 import java.util.Set;
26
27 /**
28  * Helper methods for manipulations involving sequence mappings.
29  * 
30  * @author gmcarstairs
31  *
32  */
33 public final class MappingUtils
34 {
35
36   /**
37    * Helper method to map a CUT or PASTE command.
38    * 
39    * @param edit
40    *          the original command
41    * @param undo
42    *          if true, the command is to be undone
43    * @param targetSeqs
44    *          the mapped sequences to apply the mapped command to
45    * @param result
46    *          the mapped EditCommand to add to
47    * @param mappings
48    */
49   protected static void mapCutOrPaste(Edit edit, boolean undo,
50           List<SequenceI> targetSeqs, EditCommand result,
51           Set<AlignedCodonFrame> mappings)
52   {
53     Action action = edit.getAction();
54     if (undo)
55     {
56       action = action.getUndoAction();
57     }
58     // TODO write this
59     System.err.println("MappingUtils.mapCutOrPaste not yet implemented");
60   }
61
62   /**
63    * Returns a new EditCommand representing the given command as mapped to the
64    * given sequences. If there is no mapping, returns null.
65    * 
66    * @param command
67    * @param undo
68    * @param mapTo
69    * @param gapChar
70    * @param mappings
71    * @return
72    */
73   public static EditCommand mapEditCommand(EditCommand command,
74           boolean undo, final AlignmentI mapTo, char gapChar,
75           Set<AlignedCodonFrame> mappings)
76   {
77     /*
78      * For now, only support mapping from protein edits to cDna
79      */
80     if (!mapTo.isNucleotide())
81     {
82       return null;
83     }
84
85     /*
86      * Cache a copy of the target sequences so we can mimic successive edits on
87      * them. This lets us compute mappings for all edits in the set.
88      */
89     Map<SequenceI, SequenceI> targetCopies = new HashMap<SequenceI, SequenceI>();
90     for (SequenceI seq : mapTo.getSequences())
91     {
92       SequenceI ds = seq.getDatasetSequence();
93       if (ds != null)
94       {
95         final SequenceI copy = new Sequence("", new String(
96                 seq.getSequence()));
97         copy.setDatasetSequence(ds);
98         targetCopies.put(ds, copy);
99       }
100     }
101
102     /*
103      * Compute 'source' sequences as they were before applying edits:
104      */
105     Map<SequenceI, SequenceI> originalSequences = command.priorState(undo);
106
107     EditCommand result = new EditCommand();
108     Iterator<Edit> edits = command.getEditIterator(!undo);
109     while (edits.hasNext())
110     {
111       Edit edit = edits.next();
112       if (edit.getAction() == Action.CUT
113               || edit.getAction() == Action.PASTE)
114       {
115         mapCutOrPaste(edit, undo, mapTo.getSequences(), result, mappings);
116       }
117       else if (edit.getAction() == Action.INSERT_GAP
118               || edit.getAction() == Action.DELETE_GAP)
119       {
120         mapInsertOrDelete(edit, undo, originalSequences,
121                 mapTo.getSequences(), targetCopies, gapChar, result,
122                 mappings);
123       }
124     }
125     return result.getSize() > 0 ? result : null;
126   }
127
128   /**
129    * Helper method to map an edit command to insert or delete gaps.
130    * 
131    * @param edit
132    *          the original command
133    * @param undo
134    *          if true, the action is to undo the command
135    * @param originalSequences
136    *          the sequences the command acted on
137    * @param targetSeqs
138    * @param targetCopies
139    * @param gapChar
140    * @param result
141    *          the new EditCommand to add mapped commands to
142    * @param mappings
143    */
144   protected static void mapInsertOrDelete(Edit edit, boolean undo,
145           Map<SequenceI, SequenceI> originalSequences,
146           final List<SequenceI> targetSeqs,
147           Map<SequenceI, SequenceI> targetCopies, char gapChar,
148           EditCommand result, Set<AlignedCodonFrame> mappings)
149   {
150     Action action = edit.getAction();
151
152     /*
153      * Invert sense of action if an Undo.
154      */
155     if (undo)
156     {
157       action = action.getUndoAction();
158     }
159     final int count = edit.getNumber();
160     final int editPos = edit.getPosition();
161     for (SequenceI seq : edit.getSequences())
162     {
163       /*
164        * Get residue position at (or to right of) edit location. Note we use our
165        * 'copy' of the sequence before editing for this.
166        */
167       SequenceI ds = seq.getDatasetSequence();
168       if (ds == null)
169       {
170         continue;
171       }
172       final SequenceI actedOn = originalSequences.get(ds);
173       final int seqpos = actedOn.findPosition(editPos);
174
175       /*
176        * Determine all mappings from this position to mapped sequences.
177        */
178       SearchResults sr = buildSearchResults(seq, seqpos, mappings);
179
180       if (!sr.isEmpty())
181       {
182         for (SequenceI targetSeq : targetSeqs)
183         {
184           ds = targetSeq.getDatasetSequence();
185           if (ds == null)
186           {
187             continue;
188           }
189           SequenceI copyTarget = targetCopies.get(ds);
190           final int[] match = sr.getResults(copyTarget, 0,
191                   copyTarget.getLength());
192           if (match != null)
193           {
194             final int ratio = 3; // TODO: compute this - how?
195             final int mappedCount = count * ratio;
196
197             /*
198              * Shift Delete start position left, as it acts on positions to its
199              * right.
200              */
201             int mappedEditPos = action == Action.DELETE_GAP ? match[0]
202                     - mappedCount : match[0];
203             Edit e = result.new Edit(action, new SequenceI[]
204             { targetSeq }, mappedEditPos, mappedCount, gapChar);
205             result.addEdit(e);
206
207             /*
208              * and 'apply' the edit to our copy of its target sequence
209              */
210             if (action == Action.INSERT_GAP)
211             {
212               copyTarget.setSequence(new String(StringUtils.insertCharAt(
213                       copyTarget.getSequence(), mappedEditPos, mappedCount,
214                       gapChar)));
215             }
216             else if (action == Action.DELETE_GAP)
217             {
218               copyTarget.setSequence(new String(StringUtils.deleteChars(
219                       copyTarget.getSequence(), mappedEditPos,
220                       mappedEditPos + mappedCount)));
221             }
222           }
223         }
224       }
225       /*
226        * and 'apply' the edit to our copy of its source sequence
227        */
228       if (action == Action.INSERT_GAP)
229       {
230         actedOn.setSequence(new String(StringUtils.insertCharAt(
231                 actedOn.getSequence(), editPos, count, gapChar)));
232       }
233       else if (action == Action.DELETE_GAP)
234       {
235         actedOn.setSequence(new String(StringUtils.deleteChars(
236                 actedOn.getSequence(), editPos, editPos + count)));
237       }
238     }
239   }
240
241   /**
242    * Returns a SearchResults object describing the mapped region corresponding
243    * to the specified sequence position.
244    * 
245    * @param seq
246    * @param index
247    * @param seqmappings
248    * @return
249    */
250   public static SearchResults buildSearchResults(SequenceI seq, int index,
251           Set<AlignedCodonFrame> seqmappings)
252   {
253     SearchResults results;
254     results = new SearchResults();
255     if (index >= seq.getStart() && index <= seq.getEnd())
256     {
257       for (AlignedCodonFrame acf : seqmappings)
258       {
259         acf.markMappedRegion(seq, index, results);
260       }
261     }
262     return results;
263   }
264
265   /**
266    * Returns a (possibly empty) SequenceGroup containing any sequences in the
267    * mapped viewport corresponding to the given group in the source viewport.
268    * 
269    * @param sg
270    * @param mapFrom
271    * @param mapTo
272    * @return
273    */
274   public static SequenceGroup mapSequenceGroup(SequenceGroup sg,
275           AlignViewportI mapFrom, AlignViewportI mapTo)
276   {
277     /*
278      * Note the SequenceGroup holds aligned sequences, the mappings hold dataset
279      * sequences.
280      */
281     boolean targetIsNucleotide = mapTo.isNucleotide();
282     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
283     Set<AlignedCodonFrame> codonFrames = protein.getAlignment()
284             .getCodonFrames();
285
286     /*
287      * Copy group name, name colours, but not sequences or sequence colour
288      * scheme
289      */
290     SequenceGroup mappedGroup = new SequenceGroup(sg);
291     sg.cs = mapTo.getGlobalColourScheme();
292     mappedGroup.clear();
293     // TODO set width of mapped group
294
295     for (SequenceI selected : sg.getSequences())
296     {
297       for (AlignedCodonFrame acf : codonFrames)
298       {
299         SequenceI mappedSequence = targetIsNucleotide ? acf
300                 .getDnaForAaSeq(selected) : acf.getAaForDnaSeq(selected);
301         if (mappedSequence != null)
302         {
303           for (SequenceI seq : mapTo.getAlignment().getSequences())
304           {
305             if (seq.getDatasetSequence() == mappedSequence)
306             {
307               mappedGroup.addSequence(seq, false);
308               break;
309             }
310           }
311         }
312       }
313     }
314     return mappedGroup;
315   }
316
317   /**
318    * Returns an OrderCommand equivalent to the given one, but acting on mapped
319    * sequences as described by the mappings, or null if no mapping can be made.
320    * 
321    * @param command
322    *          the original order command
323    * @param undo
324    *          if true, the action is to undo the sort
325    * @param mapTo
326    *          the alignment we are mapping to
327    * @param mappings
328    *          the mappings available
329    * @return
330    */
331   public static CommandI mapOrderCommand(OrderCommand command,
332           boolean undo, AlignmentI mapTo, Set<AlignedCodonFrame> mappings)
333   {
334     SequenceI[] sortOrder = command.getSequenceOrder(undo);
335     List<SequenceI> mappedOrder = new ArrayList<SequenceI>();
336     int j = 0;
337     for (SequenceI seq : sortOrder)
338     {
339       for (AlignedCodonFrame acf : mappings)
340       {
341         /*
342          * Try protein-to-Dna, failing that try dna-to-protein
343          */
344         SequenceI mappedSeq = acf.getDnaForAaSeq(seq);
345         if (mappedSeq == null)
346         {
347           mappedSeq = acf.getAaForDnaSeq(seq);
348         }
349         if (mappedSeq != null)
350         {
351           for (SequenceI seq2 : mapTo.getSequences())
352           {
353             if (seq2.getDatasetSequence() == mappedSeq)
354             {
355               mappedOrder.add(seq2);
356               j++;
357               break;
358             }
359           }
360         }
361       }
362     }
363
364     /*
365      * Return null if no mappings made.
366      */
367     if (j == 0)
368     {
369       return null;
370     }
371
372     /*
373      * Add any unmapped sequences on the end of the sort in their original
374      * ordering.
375      */
376     if (j < mapTo.getHeight())
377     {
378       for (SequenceI seq : mapTo.getSequences())
379       {
380         if (!mappedOrder.contains(seq))
381         {
382           mappedOrder.add(seq);
383         }
384       }
385     }
386
387     /*
388      * Have to sort the sequences before constructing the OrderCommand - which
389      * then resorts them?!?
390      */
391     final SequenceI[] mappedOrderArray = mappedOrder
392             .toArray(new SequenceI[mappedOrder.size()]);
393     SequenceI[] oldOrder = mapTo.getSequencesArray();
394     AlignmentSorter.sortBy(mapTo, new AlignmentOrder(mappedOrderArray));
395     final OrderCommand result = new OrderCommand(command.getDescription(),
396             oldOrder, mapTo);
397     return result;
398   }
399
400   /**
401    * Returns a ColumnSelection in the 'mapTo' view which corresponds to the
402    * given selection in the 'mapFrom' view. We assume one is nucleotide, the
403    * other is protein (and holds the mappings from codons to protein residues).
404    * 
405    * @param colsel
406    * @param mapFrom
407    * @param mapTo
408    * @return
409    */
410   public static ColumnSelection mapColumnSelection(ColumnSelection colsel,
411           AlignViewportI mapFrom, AlignViewportI mapTo)
412   {
413     boolean targetIsNucleotide = mapTo.isNucleotide();
414     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
415     Set<AlignedCodonFrame> codonFrames = protein.getAlignment()
416             .getCodonFrames();
417     ColumnSelection mappedColumns = new ColumnSelection();
418     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
419
420     // FIXME allow for hidden columns
421
422     /*
423      * For each mapped column, find the range of columns that residues in that
424      * column map to.
425      */
426     for (Object obj : colsel.getSelected())
427     {
428       int col = ((Integer) obj).intValue();
429       int mappedToMin = Integer.MAX_VALUE;
430       int mappedToMax = Integer.MIN_VALUE;
431
432       /*
433        * For each sequence in the 'from' alignment
434        */
435       for (SequenceI fromSeq : mapFrom.getAlignment().getSequences())
436       {
437         /*
438          * Ignore gaps (unmapped anyway)
439          */
440         if (fromSeq.getCharAt(col) == fromGapChar)
441         {
442           continue;
443         }
444
445         /*
446          * Get the residue position and find the mapped position.
447          */
448         int residuePos = fromSeq.findPosition(col);
449         SearchResults sr = buildSearchResults(fromSeq, residuePos,
450                 codonFrames);
451         for (Match m : sr.getResults())
452         {
453           int mappedStartResidue = m.getStart();
454           int mappedEndResidue = m.getEnd();
455           SequenceI mappedSeq = m.getSequence();
456
457           /*
458            * Locate the aligned sequence whose dataset is mappedSeq. TODO a
459            * datamodel that can do this efficiently.
460            */
461           for (SequenceI toSeq : mapTo.getAlignment().getSequences())
462           {
463             if (toSeq.getDatasetSequence() == mappedSeq)
464             {
465               int mappedStartCol = toSeq.findIndex(mappedStartResidue);
466               int mappedEndCol = toSeq.findIndex(mappedEndResidue);
467               mappedToMin = Math.min(mappedToMin, mappedStartCol);
468               mappedToMax = Math.max(mappedToMax, mappedEndCol);
469               // System.out.println(fromSeq.getName() + " mapped to cols "
470               // + mappedStartCol + ":" + mappedEndCol);
471               break;
472               // note: remove break if we ever want to map one to many sequences
473             }
474           }
475         }
476       }
477       /*
478        * Add the range of mapped columns to the mapped selection (converting
479        * base 1 to base 0). Note that this may include intron-only regions which
480        * lie between the start and end ranges of the selection.
481        */
482       for (int i = mappedToMin; i <= mappedToMax; i++)
483       {
484         mappedColumns.addElement(i - 1);
485       }
486     }
487     return mappedColumns;
488   }
489
490 }