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