b552c21413d761a427ed227c8876ffad801e5b7e
[jalview.git] / MappingUtils.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ 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.util;
22
23 import jalview.analysis.AlignmentSorter;
24 import jalview.api.AlignViewportI;
25 import jalview.commands.CommandI;
26 import jalview.commands.EditCommand;
27 import jalview.commands.EditCommand.Action;
28 import jalview.commands.EditCommand.Edit;
29 import jalview.commands.OrderCommand;
30 import jalview.datamodel.AlignedCodonFrame;
31 import jalview.datamodel.AlignmentI;
32 import jalview.datamodel.AlignmentOrder;
33 import jalview.datamodel.ColumnSelection;
34 import jalview.datamodel.HiddenColumns;
35 import jalview.datamodel.SearchResultMatchI;
36 import jalview.datamodel.SearchResults;
37 import jalview.datamodel.SearchResultsI;
38 import jalview.datamodel.Sequence;
39 import jalview.datamodel.SequenceGroup;
40 import jalview.datamodel.SequenceI;
41
42 import java.util.ArrayList;
43 import java.util.Arrays;
44 import java.util.HashMap;
45 import java.util.Iterator;
46 import java.util.List;
47 import java.util.Map;
48
49 /**
50  * Helper methods for manipulations involving sequence mappings.
51  * 
52  * @author gmcarstairs
53  *
54  */
55 public final class MappingUtils
56 {
57
58   /**
59    * Helper method to map a CUT or PASTE command.
60    * 
61    * @param edit
62    *          the original command
63    * @param undo
64    *          if true, the command is to be undone
65    * @param targetSeqs
66    *          the mapped sequences to apply the mapped command to
67    * @param result
68    *          the mapped EditCommand to add to
69    * @param mappings
70    */
71   protected static void mapCutOrPaste(Edit edit, boolean undo,
72           List<SequenceI> targetSeqs, EditCommand result,
73           List<AlignedCodonFrame> mappings)
74   {
75     Action action = edit.getAction();
76     if (undo)
77     {
78       action = action.getUndoAction();
79     }
80     // TODO write this
81     System.err.println("MappingUtils.mapCutOrPaste not yet implemented");
82   }
83
84   /**
85    * Returns a new EditCommand representing the given command as mapped to the
86    * given sequences. If there is no mapping, returns null.
87    * 
88    * @param command
89    * @param undo
90    * @param mapTo
91    * @param gapChar
92    * @param mappings
93    * @return
94    */
95   public static EditCommand mapEditCommand(EditCommand command,
96           boolean undo, final AlignmentI mapTo, char gapChar,
97           List<AlignedCodonFrame> mappings)
98   {
99     /*
100      * For now, only support mapping from protein edits to cDna
101      */
102     if (!mapTo.isNucleotide())
103     {
104       return null;
105     }
106
107     /*
108      * Cache a copy of the target sequences so we can mimic successive edits on
109      * them. This lets us compute mappings for all edits in the set.
110      */
111     Map<SequenceI, SequenceI> targetCopies = new HashMap<>();
112     for (SequenceI seq : mapTo.getSequences())
113     {
114       SequenceI ds = seq.getDatasetSequence();
115       if (ds != null)
116       {
117         final SequenceI copy = new Sequence(seq);
118         copy.setDatasetSequence(ds);
119         targetCopies.put(ds, copy);
120       }
121     }
122
123     /*
124      * Compute 'source' sequences as they were before applying edits:
125      */
126     Map<SequenceI, SequenceI> originalSequences = command.priorState(undo);
127
128     EditCommand result = new EditCommand();
129     Iterator<Edit> edits = command.getEditIterator(!undo);
130     while (edits.hasNext())
131     {
132       Edit edit = edits.next();
133       if (edit.getAction() == Action.CUT
134               || edit.getAction() == Action.PASTE)
135       {
136         mapCutOrPaste(edit, undo, mapTo.getSequences(), result, mappings);
137       }
138       else if (edit.getAction() == Action.INSERT_GAP
139               || edit.getAction() == Action.DELETE_GAP)
140       {
141         mapInsertOrDelete(edit, undo, originalSequences,
142                 mapTo.getSequences(), targetCopies, gapChar, result,
143                 mappings);
144       }
145     }
146     return result.getSize() > 0 ? result : null;
147   }
148
149   /**
150    * Helper method to map an edit command to insert or delete gaps.
151    * 
152    * @param edit
153    *          the original command
154    * @param undo
155    *          if true, the action is to undo the command
156    * @param originalSequences
157    *          the sequences the command acted on
158    * @param targetSeqs
159    * @param targetCopies
160    * @param gapChar
161    * @param result
162    *          the new EditCommand to add mapped commands to
163    * @param mappings
164    */
165   protected static void mapInsertOrDelete(Edit edit, boolean undo,
166           Map<SequenceI, SequenceI> originalSequences,
167           final List<SequenceI> targetSeqs,
168           Map<SequenceI, SequenceI> targetCopies, char gapChar,
169           EditCommand result, List<AlignedCodonFrame> mappings)
170   {
171     Action action = edit.getAction();
172
173     /*
174      * Invert sense of action if an Undo.
175      */
176     if (undo)
177     {
178       action = action.getUndoAction();
179     }
180     final int count = edit.getNumber();
181     final int editPos = edit.getPosition();
182     for (SequenceI seq : edit.getSequences())
183     {
184       /*
185        * Get residue position at (or to right of) edit location. Note we use our
186        * 'copy' of the sequence before editing for this.
187        */
188       SequenceI ds = seq.getDatasetSequence();
189       if (ds == null)
190       {
191         continue;
192       }
193       final SequenceI actedOn = originalSequences.get(ds);
194       final int seqpos = actedOn.findPosition(editPos);
195
196       /*
197        * Determine all mappings from this position to mapped sequences.
198        */
199       SearchResultsI sr = buildSearchResults(seq, seqpos, mappings);
200
201       if (!sr.isEmpty())
202       {
203         for (SequenceI targetSeq : targetSeqs)
204         {
205           ds = targetSeq.getDatasetSequence();
206           if (ds == null)
207           {
208             continue;
209           }
210           SequenceI copyTarget = targetCopies.get(ds);
211           final int[] match = sr.getResults(copyTarget, 0,
212                   copyTarget.getLength());
213           if (match != null)
214           {
215             final int ratio = 3; // TODO: compute this - how?
216             final int mappedCount = count * ratio;
217
218             /*
219              * Shift Delete start position left, as it acts on positions to its
220              * right.
221              */
222             int mappedEditPos = action == Action.DELETE_GAP
223                     ? match[0] - mappedCount
224                     : match[0];
225             Edit e = result.new Edit(action, new SequenceI[] { targetSeq },
226                     mappedEditPos, mappedCount, gapChar);
227             result.addEdit(e);
228
229             /*
230              * and 'apply' the edit to our copy of its target sequence
231              */
232             if (action == Action.INSERT_GAP)
233             {
234               copyTarget.setSequence(new String(
235                       StringUtils.insertCharAt(copyTarget.getSequence(),
236                               mappedEditPos, mappedCount, gapChar)));
237             }
238             else if (action == Action.DELETE_GAP)
239             {
240               copyTarget.setSequence(new String(
241                       StringUtils.deleteChars(copyTarget.getSequence(),
242                               mappedEditPos, mappedEditPos + mappedCount)));
243             }
244           }
245         }
246       }
247       /*
248        * and 'apply' the edit to our copy of its source sequence
249        */
250       if (action == Action.INSERT_GAP)
251       {
252         actedOn.setSequence(new String(StringUtils.insertCharAt(
253                 actedOn.getSequence(), editPos, count, gapChar)));
254       }
255       else if (action == Action.DELETE_GAP)
256       {
257         actedOn.setSequence(new String(StringUtils.deleteChars(
258                 actedOn.getSequence(), editPos, editPos + count)));
259       }
260     }
261   }
262
263   /**
264    * Returns a SearchResults object describing the mapped region corresponding
265    * to the specified sequence position.
266    * 
267    * @param seq
268    * @param index
269    * @param seqmappings
270    * @return
271    */
272   public static SearchResultsI buildSearchResults(SequenceI seq, int index,
273           List<AlignedCodonFrame> seqmappings)
274   {
275     SearchResultsI results = new SearchResults();
276     addSearchResults(results, seq, index, seqmappings);
277     return results;
278   }
279
280   /**
281    * Adds entries to a SearchResults object describing the mapped region
282    * corresponding to the specified sequence position.
283    * 
284    * @param results
285    * @param seq
286    * @param index
287    * @param seqmappings
288    */
289   public static void addSearchResults(SearchResultsI results, SequenceI seq,
290           int index, List<AlignedCodonFrame> seqmappings)
291   {
292     if (index >= seq.getStart() && index <= seq.getEnd())
293     {
294       for (AlignedCodonFrame acf : seqmappings)
295       {
296         acf.markMappedRegion(seq, index, results);
297       }
298     }
299   }
300
301   /**
302    * Returns a (possibly empty) SequenceGroup containing any sequences in the
303    * mapped viewport corresponding to the given group in the source viewport.
304    * 
305    * @param sg
306    * @param mapFrom
307    * @param mapTo
308    * @return
309    */
310   public static SequenceGroup mapSequenceGroup(final SequenceGroup sg,
311           final AlignViewportI mapFrom, final AlignViewportI mapTo)
312   {
313     /*
314      * Note the SequenceGroup holds aligned sequences, the mappings hold dataset
315      * sequences.
316      */
317     boolean targetIsNucleotide = mapTo.isNucleotide();
318     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
319     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
320             .getCodonFrames();
321     /*
322      * Copy group name, colours etc, but not sequences or sequence colour scheme
323      */
324     SequenceGroup mappedGroup = new SequenceGroup(sg);
325     mappedGroup.setColourScheme(mapTo.getGlobalColourScheme());
326     mappedGroup.clear();
327
328     int minStartCol = -1;
329     int maxEndCol = -1;
330     final int selectionStartRes = sg.getStartRes();
331     final int selectionEndRes = sg.getEndRes();
332     for (SequenceI selected : sg.getSequences())
333     {
334       /*
335        * Find the widest range of non-gapped positions in the selection range
336        */
337       int firstUngappedPos = selectionStartRes;
338       while (firstUngappedPos <= selectionEndRes
339               && Comparison.isGap(selected.getCharAt(firstUngappedPos)))
340       {
341         firstUngappedPos++;
342       }
343
344       /*
345        * If this sequence is only gaps in the selected range, skip it
346        */
347       if (firstUngappedPos > selectionEndRes)
348       {
349         continue;
350       }
351
352       int lastUngappedPos = selectionEndRes;
353       while (lastUngappedPos >= selectionStartRes
354               && Comparison.isGap(selected.getCharAt(lastUngappedPos)))
355       {
356         lastUngappedPos--;
357       }
358
359       /*
360        * Find the selected start/end residue positions in sequence
361        */
362       int startResiduePos = selected.findPosition(firstUngappedPos);
363       int endResiduePos = selected.findPosition(lastUngappedPos);
364
365       for (AlignedCodonFrame acf : codonFrames)
366       {
367         SequenceI mappedSequence = targetIsNucleotide
368                 ? acf.getDnaForAaSeq(selected)
369                 : acf.getAaForDnaSeq(selected);
370         if (mappedSequence != null)
371         {
372           for (SequenceI seq : mapTo.getAlignment().getSequences())
373           {
374             int mappedStartResidue = 0;
375             int mappedEndResidue = 0;
376             if (seq.getDatasetSequence() == mappedSequence)
377             {
378               /*
379                * Found a sequence mapping. Locate the start/end mapped residues.
380                */
381               List<AlignedCodonFrame> mapping = Arrays
382                       .asList(new AlignedCodonFrame[]
383                       { acf });
384               SearchResultsI sr = buildSearchResults(selected,
385                       startResiduePos, mapping);
386               for (SearchResultMatchI m : sr.getResults())
387               {
388                 mappedStartResidue = m.getStart();
389                 mappedEndResidue = m.getEnd();
390               }
391               sr = buildSearchResults(selected, endResiduePos, mapping);
392               for (SearchResultMatchI m : sr.getResults())
393               {
394                 mappedStartResidue = Math.min(mappedStartResidue,
395                         m.getStart());
396                 mappedEndResidue = Math.max(mappedEndResidue, m.getEnd());
397               }
398
399               /*
400                * Find the mapped aligned columns, save the range. Note findIndex
401                * returns a base 1 position, SequenceGroup uses base 0
402                */
403               int mappedStartCol = seq.findIndex(mappedStartResidue) - 1;
404               minStartCol = minStartCol == -1 ? mappedStartCol
405                       : Math.min(minStartCol, mappedStartCol);
406               int mappedEndCol = seq.findIndex(mappedEndResidue) - 1;
407               maxEndCol = maxEndCol == -1 ? mappedEndCol
408                       : Math.max(maxEndCol, mappedEndCol);
409               mappedGroup.addSequence(seq, false);
410               break;
411             }
412           }
413         }
414       }
415     }
416     mappedGroup.setStartRes(minStartCol < 0 ? 0 : minStartCol);
417     mappedGroup.setEndRes(maxEndCol < 0 ? 0 : maxEndCol);
418     return mappedGroup;
419   }
420
421   /**
422    * Returns an OrderCommand equivalent to the given one, but acting on mapped
423    * sequences as described by the mappings, or null if no mapping can be made.
424    * 
425    * @param command
426    *          the original order command
427    * @param undo
428    *          if true, the action is to undo the sort
429    * @param mapTo
430    *          the alignment we are mapping to
431    * @param mappings
432    *          the mappings available
433    * @return
434    */
435   public static CommandI mapOrderCommand(OrderCommand command, boolean undo,
436           AlignmentI mapTo, List<AlignedCodonFrame> mappings)
437   {
438     SequenceI[] sortOrder = command.getSequenceOrder(undo);
439     List<SequenceI> mappedOrder = new ArrayList<>();
440     int j = 0;
441
442     /*
443      * Assumption: we are only interested in a cDNA/protein mapping; refactor in
444      * future if we want to support sorting (c)dna as (c)dna or protein as
445      * protein
446      */
447     boolean mappingToNucleotide = mapTo.isNucleotide();
448     for (SequenceI seq : sortOrder)
449     {
450       for (AlignedCodonFrame acf : mappings)
451       {
452         SequenceI mappedSeq = mappingToNucleotide ? acf.getDnaForAaSeq(seq)
453                 : acf.getAaForDnaSeq(seq);
454         if (mappedSeq != null)
455         {
456           for (SequenceI seq2 : mapTo.getSequences())
457           {
458             if (seq2.getDatasetSequence() == mappedSeq)
459             {
460               mappedOrder.add(seq2);
461               j++;
462               break;
463             }
464           }
465         }
466       }
467     }
468
469     /*
470      * Return null if no mappings made.
471      */
472     if (j == 0)
473     {
474       return null;
475     }
476
477     /*
478      * Add any unmapped sequences on the end of the sort in their original
479      * ordering.
480      */
481     if (j < mapTo.getHeight())
482     {
483       for (SequenceI seq : mapTo.getSequences())
484       {
485         if (!mappedOrder.contains(seq))
486         {
487           mappedOrder.add(seq);
488         }
489       }
490     }
491
492     /*
493      * Have to sort the sequences before constructing the OrderCommand - which
494      * then resorts them?!?
495      */
496     final SequenceI[] mappedOrderArray = mappedOrder
497             .toArray(new SequenceI[mappedOrder.size()]);
498     SequenceI[] oldOrder = mapTo.getSequencesArray();
499     AlignmentSorter.sortBy(mapTo, new AlignmentOrder(mappedOrderArray));
500     final OrderCommand result = new OrderCommand(command.getDescription(),
501             oldOrder, mapTo);
502     return result;
503   }
504
505   /**
506    * Returns a ColumnSelection in the 'mapTo' view which corresponds to the
507    * given selection in the 'mapFrom' view. We assume one is nucleotide, the
508    * other is protein (and holds the mappings from codons to protein residues).
509    * 
510    * @param colsel
511    * @param mapFrom
512    * @param mapTo
513    * @return
514    */
515   public static void mapColumnSelection(ColumnSelection colsel,
516           HiddenColumns hiddencols, AlignViewportI mapFrom,
517           AlignViewportI mapTo, ColumnSelection newColSel,
518           HiddenColumns newHidden)
519   {
520     boolean targetIsNucleotide = mapTo.isNucleotide();
521     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
522     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
523             .getCodonFrames();
524
525     if (colsel == null)
526     {
527       return; // mappedColumns;
528     }
529
530     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
531
532     /*
533      * For each mapped column, find the range of columns that residues in that
534      * column map to.
535      */
536     List<SequenceI> fromSequences = mapFrom.getAlignment().getSequences();
537     List<SequenceI> toSequences = mapTo.getAlignment().getSequences();
538
539     for (Integer sel : colsel.getSelected())
540     {
541       mapColumn(sel.intValue(), codonFrames, newColSel, fromSequences,
542               toSequences, fromGapChar);
543     }
544
545     Iterator<int[]> regions = hiddencols.iterator();
546     while (regions.hasNext())
547     {
548       mapHiddenColumns(regions.next(), codonFrames, newHidden,
549               fromSequences,
550               toSequences, fromGapChar);
551     }
552     return; // mappedColumns;
553   }
554
555   /**
556    * Helper method that maps a [start, end] hidden column range to its mapped
557    * equivalent
558    * 
559    * @param hidden
560    * @param mappings
561    * @param mappedColumns
562    * @param fromSequences
563    * @param toSequences
564    * @param fromGapChar
565    */
566   protected static void mapHiddenColumns(int[] hidden,
567           List<AlignedCodonFrame> mappings, HiddenColumns mappedColumns,
568           List<SequenceI> fromSequences, List<SequenceI> toSequences,
569           char fromGapChar)
570   {
571     for (int col = hidden[0]; col <= hidden[1]; col++)
572     {
573       int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
574               toSequences, fromGapChar);
575
576       /*
577        * Add the range of hidden columns to the mapped selection (converting
578        * base 1 to base 0).
579        */
580       if (mappedTo != null)
581       {
582         mappedColumns.hideColumns(mappedTo[0] - 1, mappedTo[1] - 1);
583       }
584     }
585   }
586
587   /**
588    * Helper method to map one column selection
589    * 
590    * @param col
591    *          the column number (base 0)
592    * @param mappings
593    *          the sequence mappings
594    * @param mappedColumns
595    *          the mapped column selections to add to
596    * @param fromSequences
597    * @param toSequences
598    * @param fromGapChar
599    */
600   protected static void mapColumn(int col, List<AlignedCodonFrame> mappings,
601           ColumnSelection mappedColumns, List<SequenceI> fromSequences,
602           List<SequenceI> toSequences, char fromGapChar)
603   {
604     int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
605             toSequences, fromGapChar);
606
607     /*
608      * Add the range of mapped columns to the mapped selection (converting
609      * base 1 to base 0). Note that this may include intron-only regions which
610      * lie between the start and end ranges of the selection.
611      */
612     if (mappedTo != null)
613     {
614       for (int i = mappedTo[0]; i <= mappedTo[1]; i++)
615       {
616         mappedColumns.addElement(i - 1);
617       }
618     }
619   }
620
621   /**
622    * Helper method to find the range of columns mapped to from one column.
623    * Returns the maximal range of columns mapped to from all sequences in the
624    * source column, or null if no mappings were found.
625    * 
626    * @param col
627    * @param mappings
628    * @param fromSequences
629    * @param toSequences
630    * @param fromGapChar
631    * @return
632    */
633   protected static int[] findMappedColumns(int col,
634           List<AlignedCodonFrame> mappings, List<SequenceI> fromSequences,
635           List<SequenceI> toSequences, char fromGapChar)
636   {
637     int[] mappedTo = new int[] { Integer.MAX_VALUE, Integer.MIN_VALUE };
638     boolean found = false;
639
640     /*
641      * For each sequence in the 'from' alignment
642      */
643     for (SequenceI fromSeq : fromSequences)
644     {
645       /*
646        * Ignore gaps (unmapped anyway)
647        */
648       if (fromSeq.getCharAt(col) == fromGapChar)
649       {
650         continue;
651       }
652
653       /*
654        * Get the residue position and find the mapped position.
655        */
656       int residuePos = fromSeq.findPosition(col);
657       SearchResultsI sr = buildSearchResults(fromSeq, residuePos, mappings);
658       for (SearchResultMatchI m : sr.getResults())
659       {
660         int mappedStartResidue = m.getStart();
661         int mappedEndResidue = m.getEnd();
662         SequenceI mappedSeq = m.getSequence();
663
664         /*
665          * Locate the aligned sequence whose dataset is mappedSeq. TODO a
666          * datamodel that can do this efficiently.
667          */
668         for (SequenceI toSeq : toSequences)
669         {
670           if (toSeq.getDatasetSequence() == mappedSeq)
671           {
672             int mappedStartCol = toSeq.findIndex(mappedStartResidue);
673             int mappedEndCol = toSeq.findIndex(mappedEndResidue);
674             mappedTo[0] = Math.min(mappedTo[0], mappedStartCol);
675             mappedTo[1] = Math.max(mappedTo[1], mappedEndCol);
676             found = true;
677             break;
678             // note: remove break if we ever want to map one to many sequences
679           }
680         }
681       }
682     }
683     return found ? mappedTo : null;
684   }
685
686   /**
687    * Returns the mapped codon or codons for a given aligned sequence column
688    * position (base 0).
689    * 
690    * @param seq
691    *          an aligned peptide sequence
692    * @param col
693    *          an aligned column position (base 0)
694    * @param mappings
695    *          a set of codon mappings
696    * @return the bases of the mapped codon(s) in the cDNA dataset sequence(s),
697    *         or an empty list if none found
698    */
699   public static List<char[]> findCodonsFor(SequenceI seq, int col,
700           List<AlignedCodonFrame> mappings)
701   {
702     List<char[]> result = new ArrayList<>();
703     int dsPos = seq.findPosition(col);
704     for (AlignedCodonFrame mapping : mappings)
705     {
706       if (mapping.involvesSequence(seq))
707       {
708         List<char[]> codons = mapping
709                 .getMappedCodons(seq.getDatasetSequence(), dsPos);
710         if (codons != null)
711         {
712           result.addAll(codons);
713         }
714       }
715     }
716     return result;
717   }
718
719   /**
720    * Converts a series of [start, end] range pairs into an array of individual
721    * positions. This also caters for 'reverse strand' (start > end) cases.
722    * 
723    * @param ranges
724    * @return
725    */
726   public static int[] flattenRanges(int[] ranges)
727   {
728     /*
729      * Count how many positions altogether
730      */
731     int count = 0;
732     for (int i = 0; i < ranges.length - 1; i += 2)
733     {
734       count += Math.abs(ranges[i + 1] - ranges[i]) + 1;
735     }
736
737     int[] result = new int[count];
738     int k = 0;
739     for (int i = 0; i < ranges.length - 1; i += 2)
740     {
741       int from = ranges[i];
742       final int to = ranges[i + 1];
743       int step = from <= to ? 1 : -1;
744       do
745       {
746         result[k++] = from;
747         from += step;
748       } while (from != to + step);
749     }
750     return result;
751   }
752
753   /**
754    * Returns a list of any mappings that are from or to the given (aligned or
755    * dataset) sequence.
756    * 
757    * @param sequence
758    * @param mappings
759    * @return
760    */
761   public static List<AlignedCodonFrame> findMappingsForSequence(
762           SequenceI sequence, List<AlignedCodonFrame> mappings)
763   {
764     return findMappingsForSequenceAndOthers(sequence, mappings, null);
765   }
766
767   /**
768    * Returns a list of any mappings that are from or to the given (aligned or
769    * dataset) sequence, optionally limited to mappings involving one of a given
770    * list of sequences.
771    * 
772    * @param sequence
773    * @param mappings
774    * @param filterList
775    * @return
776    */
777   public static List<AlignedCodonFrame> findMappingsForSequenceAndOthers(
778           SequenceI sequence, List<AlignedCodonFrame> mappings,
779           List<SequenceI> filterList)
780   {
781     List<AlignedCodonFrame> result = new ArrayList<>();
782     if (sequence == null || mappings == null)
783     {
784       return result;
785     }
786     for (AlignedCodonFrame mapping : mappings)
787     {
788       if (mapping.involvesSequence(sequence))
789       {
790         if (filterList != null)
791         {
792           for (SequenceI otherseq : filterList)
793           {
794             SequenceI otherDataset = otherseq.getDatasetSequence();
795             if (otherseq == sequence
796                     || otherseq == sequence.getDatasetSequence()
797                     || (otherDataset != null && (otherDataset == sequence
798                             || otherDataset == sequence
799                                     .getDatasetSequence())))
800             {
801               // skip sequences in subset which directly relate to sequence
802               continue;
803             }
804             if (mapping.involvesSequence(otherseq))
805             {
806               // selected a mapping contained in subselect alignment
807               result.add(mapping);
808               break;
809             }
810           }
811         }
812         else
813         {
814           result.add(mapping);
815         }
816       }
817     }
818     return result;
819   }
820
821   /**
822    * Returns the total length of the supplied ranges, which may be as single
823    * [start, end] or multiple [start, end, start, end ...]
824    * 
825    * @param ranges
826    * @return
827    */
828   public static int getLength(List<int[]> ranges)
829   {
830     if (ranges == null)
831     {
832       return 0;
833     }
834     int length = 0;
835     for (int[] range : ranges)
836     {
837       if (range.length % 2 != 0)
838       {
839         System.err.println(
840                 "Error unbalance start/end ranges: " + ranges.toString());
841         return 0;
842       }
843       for (int i = 0; i < range.length - 1; i += 2)
844       {
845         length += Math.abs(range[i + 1] - range[i]) + 1;
846       }
847     }
848     return length;
849   }
850
851   /**
852    * Answers true if any range includes the given value
853    * 
854    * @param ranges
855    * @param value
856    * @return
857    */
858   public static boolean contains(List<int[]> ranges, int value)
859   {
860     if (ranges == null)
861     {
862       return false;
863     }
864     for (int[] range : ranges)
865     {
866       if (range[1] >= range[0] && value >= range[0] && value <= range[1])
867       {
868         /*
869          * value within ascending range
870          */
871         return true;
872       }
873       if (range[1] < range[0] && value <= range[0] && value >= range[1])
874       {
875         /*
876          * value within descending range
877          */
878         return true;
879       }
880     }
881     return false;
882   }
883
884   /**
885    * Removes a specified number of positions from the start of a ranges list.
886    * For example, could be used to adjust cds ranges to allow for an incomplete
887    * start codon. Subranges are removed completely, or their start positions
888    * adjusted, until the required number of positions has been removed from the
889    * range. Reverse strand ranges are supported. The input array is not
890    * modified.
891    * 
892    * @param removeCount
893    * @param ranges
894    *          an array of [start, end, start, end...] positions
895    * @return a new array with the first removeCount positions removed
896    */
897   public static int[] removeStartPositions(int removeCount,
898           final int[] ranges)
899   {
900     if (removeCount <= 0)
901     {
902       return ranges;
903     }
904
905     int[] copy = Arrays.copyOf(ranges, ranges.length);
906     int sxpos = -1;
907     int cdspos = 0;
908     for (int x = 0; x < copy.length && sxpos == -1; x += 2)
909     {
910       cdspos += Math.abs(copy[x + 1] - copy[x]) + 1;
911       if (removeCount < cdspos)
912       {
913         /*
914          * we have removed enough, time to finish
915          */
916         sxpos = x;
917
918         /*
919          * increment start of first exon, or decrement if reverse strand
920          */
921         if (copy[x] <= copy[x + 1])
922         {
923           copy[x] = copy[x + 1] - cdspos + removeCount + 1;
924         }
925         else
926         {
927           copy[x] = copy[x + 1] + cdspos - removeCount - 1;
928         }
929         break;
930       }
931     }
932
933     if (sxpos > 0)
934     {
935       /*
936        * we dropped at least one entire sub-range - compact the array
937        */
938       int[] nxon = new int[copy.length - sxpos];
939       System.arraycopy(copy, sxpos, nxon, 0, copy.length - sxpos);
940       return nxon;
941     }
942     return copy;
943   }
944
945   /**
946    * Answers true if range's start-end positions include those of queryRange,
947    * where either range might be in reverse direction, else false
948    * 
949    * @param range
950    *          a start-end range
951    * @param queryRange
952    *          a candidate subrange of range (start2-end2)
953    * @return
954    */
955   public static boolean rangeContains(int[] range, int[] queryRange)
956   {
957     if (range == null || queryRange == null || range.length != 2
958             || queryRange.length != 2)
959     {
960       /*
961        * invalid arguments
962        */
963       return false;
964     }
965
966     int min = Math.min(range[0], range[1]);
967     int max = Math.max(range[0], range[1]);
968   
969     return (min <= queryRange[0] && max >= queryRange[0]
970             && min <= queryRange[1] && max >= queryRange[1]);
971   }
972
973   /**
974    * Removes the specified number of positions from the given ranges. Provided
975    * to allow a stop codon to be stripped from a CDS sequence so that it matches
976    * the peptide translation length.
977    * 
978    * @param positions
979    * @param ranges
980    *          a list of (single) [start, end] ranges
981    * @return
982    */
983   public static void removeEndPositions(int positions,
984           List<int[]> ranges)
985   {
986     int toRemove = positions;
987     Iterator<int[]> it = new ReverseListIterator<>(ranges);
988     while (toRemove > 0)
989     {
990       int[] endRange = it.next();
991       if (endRange.length != 2)
992       {
993         /*
994          * not coded for [start1, end1, start2, end2, ...]
995          */
996         System.err
997                 .println("MappingUtils.removeEndPositions doesn't handle multiple  ranges");
998         return;
999       }
1000
1001       int length = endRange[1] - endRange[0] + 1;
1002       if (length <= 0)
1003       {
1004         /*
1005          * not coded for a reverse strand range (end < start)
1006          */
1007         System.err
1008                 .println("MappingUtils.removeEndPositions doesn't handle reverse strand");
1009         return;
1010       }
1011       if (length > toRemove)
1012       {
1013         endRange[1] -= toRemove;
1014         toRemove = 0;
1015       }
1016       else
1017       {
1018         toRemove -= length;
1019         it.remove();
1020       }
1021     }
1022   }
1023 }