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