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.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.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       boolean allGapped = (firstUngappedPos > selectionEndRes);
 347
 348       int lastUngappedPos = selectionEndRes;
 349       if (!allGapped)
 350       {
 351         while (lastUngappedPos >= selectionStartRes
 352                 && Comparison.isGap(selected.getCharAt(lastUngappedPos)))
 353         {
 354           lastUngappedPos--;
 355         }
 356       }
 357
 358       /*
 359        * Find the selected start/end residue positions in sequence
 360        */
 361       int startResiduePos = allGapped ? 0 : selected.findPosition(firstUngappedPos);
 362       int endResiduePos = allGapped ? 0 : selected.findPosition(lastUngappedPos);
 363
 364       for (AlignedCodonFrame acf : codonFrames)
 365       {
 366         for (SequenceI seq : mapTo.getAlignment().getSequences())
 367         {
 368           SequenceI peptide = targetIsNucleotide ? selected : seq;
 369           SequenceI cds = targetIsNucleotide ? seq : selected;
 370           SequenceToSequenceMapping s2s = acf.getCoveringMapping(cds,
 371                   peptide);
 372           if (s2s == null)
 373           {
 374             continue;
 375           }
 376           mappedGroup.addSequence(seq, false);
 377           if (allGapped)
 378           {
 379             /*
 380              * sequence is mapped but includes no mapped residues
 381              */
 382             continue;
 383           }
 384           int mappedStartResidue = 0;
 385           int mappedEndResidue = 0;
 386           List<AlignedCodonFrame> mapping = Arrays.asList(acf);
 387           SearchResultsI sr = buildSearchResults(selected, startResiduePos,
 388                   mapping);
 389           for (SearchResultMatchI m : sr.getResults())
 390           {
 391             mappedStartResidue = m.getStart();
 392             mappedEndResidue = m.getEnd();
 393           }
 394           sr = buildSearchResults(selected, endResiduePos, mapping);
 395           for (SearchResultMatchI m : sr.getResults())
 396           {
 397             mappedStartResidue = Math.min(mappedStartResidue, m.getStart());
 398             mappedEndResidue = Math.max(mappedEndResidue, m.getEnd());
 399           }
 400
 401           /*
 402            * Find the mapped aligned columns, save the range. Note findIndex
 403            * returns a base 1 position, SequenceGroup uses base 0
 404            */
 405           int mappedStartCol = seq.findIndex(mappedStartResidue) - 1;
 406           minStartCol = minStartCol == -1 ? mappedStartCol
 407                   : Math.min(minStartCol, mappedStartCol);
 408           int mappedEndCol = seq.findIndex(mappedEndResidue) - 1;
 409           maxEndCol = maxEndCol == -1 ? mappedEndCol
 410                   : Math.max(maxEndCol, mappedEndCol);
 411           break;
 412         }
 413       }
 414     }
 415     mappedGroup.setStartRes(minStartCol < 0 ? 0 : minStartCol);
 416     mappedGroup.setEndRes(maxEndCol < 0 ? 0 : maxEndCol);
 417     return mappedGroup;
 418   }
 419
 420   /**
 421    * Returns an OrderCommand equivalent to the given one, but acting on mapped
 422    * sequences as described by the mappings, or null if no mapping can be made.
 423    *
 424    * @param command
 425    *          the original order command
 426    * @param undo
 427    *          if true, the action is to undo the sort
 428    * @param mapTo
 429    *          the alignment we are mapping to
 430    * @param mappings
 431    *          the mappings available
 432    * @return
 433    */
 434   public static CommandI mapOrderCommand(OrderCommand command, boolean undo,
 435           AlignmentI mapTo, List<AlignedCodonFrame> mappings)
 436   {
 437     SequenceI[] sortOrder = command.getSequenceOrder(undo);
 438     List<SequenceI> mappedOrder = new ArrayList<>();
 439     int j = 0;
 440
 441     /*
 442      * Assumption: we are only interested in a cDNA/protein mapping; refactor in
 443      * future if we want to support sorting (c)dna as (c)dna or protein as
 444      * protein
 445      */
 446     boolean mappingToNucleotide = mapTo.isNucleotide();
 447     for (SequenceI seq : sortOrder)
 448     {
 449       for (AlignedCodonFrame acf : mappings)
 450       {
 451           for (SequenceI seq2 : mapTo.getSequences())
 452           {
 453             /*
 454              * the corresponding peptide / CDS is the one for which there is
 455              * a complete ('covering') mapping to 'seq'
 456              */
 457             SequenceI peptide = mappingToNucleotide ? seq2 : seq;
 458             SequenceI cds = mappingToNucleotide ? seq : seq2;
 459             SequenceToSequenceMapping s2s = acf.getCoveringMapping(cds,
 460                     peptide);
 461             if (s2s != null)
 462             {
 463               mappedOrder.add(seq2);
 464               j++;
 465               break;
 466             }
 467           }
 468       }
 469     }
 470
 471     /*
 472      * Return null if no mappings made.
 473      */
 474     if (j == 0)
 475     {
 476       return null;
 477     }
 478
 479     /*
 480      * Add any unmapped sequences on the end of the sort in their original
 481      * ordering.
 482      */
 483     if (j < mapTo.getHeight())
 484     {
 485       for (SequenceI seq : mapTo.getSequences())
 486       {
 487         if (!mappedOrder.contains(seq))
 488         {
 489           mappedOrder.add(seq);
 490         }
 491       }
 492     }
 493
 494     /*
 495      * Have to sort the sequences before constructing the OrderCommand - which
 496      * then resorts them?!?
 497      */
 498     final SequenceI[] mappedOrderArray = mappedOrder
 499             .toArray(new SequenceI[mappedOrder.size()]);
 500     SequenceI[] oldOrder = mapTo.getSequencesArray();
 501     AlignmentSorter.sortBy(mapTo, new AlignmentOrder(mappedOrderArray));
 502     final OrderCommand result = new OrderCommand(command.getDescription(),
 503             oldOrder, mapTo);
 504     return result;
 505   }
 506
 507   /**
 508    * Returns a ColumnSelection in the 'mapTo' view which corresponds to the
 509    * given selection in the 'mapFrom' view. We assume one is nucleotide, the
 510    * other is protein (and holds the mappings from codons to protein residues).
 511    *
 512    * @param colsel
 513    * @param mapFrom
 514    * @param mapTo
 515    * @return
 516    */
 517   public static void mapColumnSelection(ColumnSelection colsel,
 518           HiddenColumns hiddencols, AlignViewportI mapFrom,
 519           AlignViewportI mapTo, ColumnSelection newColSel,
 520           HiddenColumns newHidden)
 521   {
 522     boolean targetIsNucleotide = mapTo.isNucleotide();
 523     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 524     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
 525             .getCodonFrames();
 526
 527     if (colsel == null)
 528     {
 529       return;
 530     }
 531
 532     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
 533
 534     /*
 535      * For each mapped column, find the range of columns that residues in that
 536      * column map to.
 537      */
 538     List<SequenceI> fromSequences = mapFrom.getAlignment().getSequences();
 539     List<SequenceI> toSequences = mapTo.getAlignment().getSequences();
 540
 541     for (Integer sel : colsel.getSelected())
 542     {
 543       mapColumn(sel.intValue(), codonFrames, newColSel, fromSequences,
 544               toSequences, fromGapChar);
 545     }
 546
 547     Iterator<int[]> regions = hiddencols.iterator();
 548     while (regions.hasNext())
 549     {
 550       mapHiddenColumns(regions.next(), codonFrames, newHidden,
 551               fromSequences, toSequences, fromGapChar);
 552     }
 553     return;
 554   }
 555
 556   /**
 557    * Helper method that maps a [start, end] hidden column range to its mapped
 558    * equivalent
 559    *
 560    * @param hidden
 561    * @param mappings
 562    * @param mappedColumns
 563    * @param fromSequences
 564    * @param toSequences
 565    * @param fromGapChar
 566    */
 567   protected static void mapHiddenColumns(int[] hidden,
 568           List<AlignedCodonFrame> mappings, HiddenColumns mappedColumns,
 569           List<SequenceI> fromSequences, List<SequenceI> toSequences,
 570           char fromGapChar)
 571   {
 572     for (int col = hidden[0]; col <= hidden[1]; col++)
 573     {
 574       int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
 575               toSequences, fromGapChar);
 576
 577       /*
 578        * Add the range of hidden columns to the mapped selection (converting
 579        * base 1 to base 0).
 580        */
 581       if (mappedTo != null)
 582       {
 583         mappedColumns.hideColumns(mappedTo[0] - 1, mappedTo[1] - 1);
 584       }
 585     }
 586   }
 587
 588   /**
 589    * Helper method to map one column selection
 590    *
 591    * @param col
 592    *          the column number (base 0)
 593    * @param mappings
 594    *          the sequence mappings
 595    * @param mappedColumns
 596    *          the mapped column selections to add to
 597    * @param fromSequences
 598    * @param toSequences
 599    * @param fromGapChar
 600    */
 601   protected static void mapColumn(int col, List<AlignedCodonFrame> mappings,
 602           ColumnSelection mappedColumns, List<SequenceI> fromSequences,
 603           List<SequenceI> toSequences, char fromGapChar)
 604   {
 605     int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
 606             toSequences, fromGapChar);
 607
 608     /*
 609      * Add the range of mapped columns to the mapped selection (converting
 610      * base 1 to base 0). Note that this may include intron-only regions which
 611      * lie between the start and end ranges of the selection.
 612      */
 613     if (mappedTo != null)
 614     {
 615       for (int i = mappedTo[0]; i <= mappedTo[1]; i++)
 616       {
 617         mappedColumns.addElement(i - 1);
 618       }
 619     }
 620   }
 621
 622   /**
 623    * Helper method to find the range of columns mapped to from one column.
 624    * Returns the maximal range of columns mapped to from all sequences in the
 625    * source column, or null if no mappings were found.
 626    *
 627    * @param col
 628    * @param mappings
 629    * @param fromSequences
 630    * @param toSequences
 631    * @param fromGapChar
 632    * @return
 633    */
 634   protected static int[] findMappedColumns(int col,
 635           List<AlignedCodonFrame> mappings, List<SequenceI> fromSequences,
 636           List<SequenceI> toSequences, char fromGapChar)
 637   {
 638     int[] mappedTo = new int[] { Integer.MAX_VALUE, Integer.MIN_VALUE };
 639     boolean found = false;
 640
 641     /*
 642      * For each sequence in the 'from' alignment
 643      */
 644     for (SequenceI fromSeq : fromSequences)
 645     {
 646       /*
 647        * Ignore gaps (unmapped anyway)
 648        */
 649       if (fromSeq.getCharAt(col) == fromGapChar)
 650       {
 651         continue;
 652       }
 653
 654       /*
 655        * Get the residue position and find the mapped position.
 656        */
 657       int residuePos = fromSeq.findPosition(col);
 658       SearchResultsI sr = buildSearchResults(fromSeq, residuePos, mappings);
 659       for (SearchResultMatchI m : sr.getResults())
 660       {
 661         int mappedStartResidue = m.getStart();
 662         int mappedEndResidue = m.getEnd();
 663         SequenceI mappedSeq = m.getSequence();
 664
 665         /*
 666          * Locate the aligned sequence whose dataset is mappedSeq. TODO a
 667          * datamodel that can do this efficiently.
 668          */
 669         for (SequenceI toSeq : toSequences)
 670         {
 671           if (toSeq.getDatasetSequence() == mappedSeq
 672                   && mappedStartResidue >= toSeq.getStart()
 673                   && mappedEndResidue <= toSeq.getEnd())
 674           {
 675             int mappedStartCol = toSeq.findIndex(mappedStartResidue);
 676             int mappedEndCol = toSeq.findIndex(mappedEndResidue);
 677             mappedTo[0] = Math.min(mappedTo[0], mappedStartCol);
 678             mappedTo[1] = Math.max(mappedTo[1], mappedEndCol);
 679             found = true;
 680             break;
 681             // note: remove break if we ever want to map one to many sequences
 682           }
 683         }
 684       }
 685     }
 686     return found ? mappedTo : null;
 687   }
 688
 689   /**
 690    * Returns the mapped codon or codons for a given aligned sequence column
 691    * position (base 0).
 692    *
 693    * @param seq
 694    *          an aligned peptide sequence
 695    * @param col
 696    *          an aligned column position (base 0)
 697    * @param mappings
 698    *          a set of codon mappings
 699    * @return the bases of the mapped codon(s) in the cDNA dataset sequence(s),
 700    *         or an empty list if none found
 701    */
 702   public static List<char[]> findCodonsFor(SequenceI seq, int col,
 703           List<AlignedCodonFrame> mappings)
 704   {
 705     List<char[]> result = new ArrayList<>();
 706     int dsPos = seq.findPosition(col);
 707     for (AlignedCodonFrame mapping : mappings)
 708     {
 709       if (mapping.involvesSequence(seq))
 710       {
 711         List<char[]> codons = mapping
 712                 .getMappedCodons(seq.getDatasetSequence(), dsPos);
 713         if (codons != null)
 714         {
 715           result.addAll(codons);
 716         }
 717       }
 718     }
 719     return result;
 720   }
 721
 722   /**
 723    * Converts a series of [start, end] range pairs into an array of individual
 724    * positions. This also caters for 'reverse strand' (start > end) cases.
 725    *
 726    * @param ranges
 727    * @return
 728    */
 729   public static int[] flattenRanges(int[] ranges)
 730   {
 731     /*
 732      * Count how many positions altogether
 733      */
 734     int count = 0;
 735     for (int i = 0; i < ranges.length - 1; i += 2)
 736     {
 737       count += Math.abs(ranges[i + 1] - ranges[i]) + 1;
 738     }
 739
 740     int[] result = new int[count];
 741     int k = 0;
 742     for (int i = 0; i < ranges.length - 1; i += 2)
 743     {
 744       int from = ranges[i];
 745       final int to = ranges[i + 1];
 746       int step = from <= to ? 1 : -1;
 747       do
 748       {
 749         result[k++] = from;
 750         from += step;
 751       } while (from != to + step);
 752     }
 753     return result;
 754   }
 755
 756   /**
 757    * Returns a list of any mappings that are from or to the given (aligned or
 758    * dataset) sequence.
 759    *
 760    * @param sequence
 761    * @param mappings
 762    * @return
 763    */
 764   public static List<AlignedCodonFrame> findMappingsForSequence(
 765           SequenceI sequence, List<AlignedCodonFrame> mappings)
 766   {
 767     return findMappingsForSequenceAndOthers(sequence, mappings, null);
 768   }
 769
 770   /**
 771    * Returns a list of any mappings that are from or to the given (aligned or
 772    * dataset) sequence, optionally limited to mappings involving one of a given
 773    * list of sequences.
 774    *
 775    * @param sequence
 776    * @param mappings
 777    * @param filterList
 778    * @return
 779    */
 780   public static List<AlignedCodonFrame> findMappingsForSequenceAndOthers(
 781           SequenceI sequence, List<AlignedCodonFrame> mappings,
 782           List<SequenceI> filterList)
 783   {
 784     List<AlignedCodonFrame> result = new ArrayList<>();
 785     if (sequence == null || mappings == null)
 786     {
 787       return result;
 788     }
 789     for (AlignedCodonFrame mapping : mappings)
 790     {
 791       if (mapping.involvesSequence(sequence))
 792       {
 793         if (filterList != null)
 794         {
 795           for (SequenceI otherseq : filterList)
 796           {
 797             SequenceI otherDataset = otherseq.getDatasetSequence();
 798             if (otherseq == sequence
 799                     || otherseq == sequence.getDatasetSequence()
 800                     || (otherDataset != null && (otherDataset == sequence
 801                             || otherDataset == sequence
 802                                     .getDatasetSequence())))
 803             {
 804               // skip sequences in subset which directly relate to sequence
 805               continue;
 806             }
 807             if (mapping.involvesSequence(otherseq))
 808             {
 809               // selected a mapping contained in subselect alignment
 810               result.add(mapping);
 811               break;
 812             }
 813           }
 814         }
 815         else
 816         {
 817           result.add(mapping);
 818         }
 819       }
 820     }
 821     return result;
 822   }
 823
 824   /**
 825    * Returns the total length of the supplied ranges, which may be as single
 826    * [start, end] or multiple [start, end, start, end ...]
 827    *
 828    * @param ranges
 829    * @return
 830    */
 831   public static int getLength(List<int[]> ranges)
 832   {
 833     if (ranges == null)
 834     {
 835       return 0;
 836     }
 837     int length = 0;
 838     for (int[] range : ranges)
 839     {
 840       if (range.length % 2 != 0)
 841       {
 842         Cache.log.error(
 843                 "Error unbalance start/end ranges: " + ranges.toString());
 844         return 0;
 845       }
 846       for (int i = 0; i < range.length - 1; i += 2)
 847       {
 848         length += Math.abs(range[i + 1] - range[i]) + 1;
 849       }
 850     }
 851     return length;
 852   }
 853
 854   /**
 855    * Answers true if any range includes the given value
 856    *
 857    * @param ranges
 858    * @param value
 859    * @return
 860    */
 861   public static boolean contains(List<int[]> ranges, int value)
 862   {
 863     if (ranges == null)
 864     {
 865       return false;
 866     }
 867     for (int[] range : ranges)
 868     {
 869       if (range[1] >= range[0] && value >= range[0] && value <= range[1])
 870       {
 871         /*
 872          * value within ascending range
 873          */
 874         return true;
 875       }
 876       if (range[1] < range[0] && value <= range[0] && value >= range[1])
 877       {
 878         /*
 879          * value within descending range
 880          */
 881         return true;
 882       }
 883     }
 884     return false;
 885   }
 886
 887   /**
 888    * Removes a specified number of positions from the start of a ranges list.
 889    * For example, could be used to adjust cds ranges to allow for an incomplete
 890    * start codon. Subranges are removed completely, or their start positions
 891    * adjusted, until the required number of positions has been removed from the
 892    * range. Reverse strand ranges are supported. The input array is not
 893    * modified.
 894    *
 895    * @param removeCount
 896    * @param ranges
 897    *          an array of [start, end, start, end...] positions
 898    * @return a new array with the first removeCount positions removed
 899    */
 900   public static int[] removeStartPositions(int removeCount,
 901           final int[] ranges)
 902   {
 903     if (removeCount <= 0)
 904     {
 905       return ranges;
 906     }
 907
 908     int[] copy = Arrays.copyOf(ranges, ranges.length);
 909     int sxpos = -1;
 910     int cdspos = 0;
 911     for (int x = 0; x < copy.length && sxpos == -1; x += 2)
 912     {
 913       cdspos += Math.abs(copy[x + 1] - copy[x]) + 1;
 914       if (removeCount < cdspos)
 915       {
 916         /*
 917          * we have removed enough, time to finish
 918          */
 919         sxpos = x;
 920
 921         /*
 922          * increment start of first exon, or decrement if reverse strand
 923          */
 924         if (copy[x] <= copy[x + 1])
 925         {
 926           copy[x] = copy[x + 1] - cdspos + removeCount + 1;
 927         }
 928         else
 929         {
 930           copy[x] = copy[x + 1] + cdspos - removeCount - 1;
 931         }
 932         break;
 933       }
 934     }
 935
 936     if (sxpos > 0)
 937     {
 938       /*
 939        * we dropped at least one entire sub-range - compact the array
 940        */
 941       int[] nxon = new int[copy.length - sxpos];
 942       System.arraycopy(copy, sxpos, nxon, 0, copy.length - sxpos);
 943       return nxon;
 944     }
 945     return copy;
 946   }
 947
 948   /**
 949    * Answers true if range's start-end positions include those of queryRange,
 950    * where either range might be in reverse direction, else false
 951    *
 952    * @param range
 953    *          a start-end range
 954    * @param queryRange
 955    *          a candidate subrange of range (start2-end2)
 956    * @return
 957    */
 958   public static boolean rangeContains(int[] range, int[] queryRange)
 959   {
 960     if (range == null || queryRange == null || range.length != 2
 961             || queryRange.length != 2)
 962     {
 963       /*
 964        * invalid arguments
 965        */
 966       return false;
 967     }
 968
 969     int min = Math.min(range[0], range[1]);
 970     int max = Math.max(range[0], range[1]);
 971
 972     return (min <= queryRange[0] && max >= queryRange[0]
 973             && min <= queryRange[1] && max >= queryRange[1]);
 974   }
 975
 976   /**
 977    * Removes the specified number of positions from the given ranges. Provided
 978    * to allow a stop codon to be stripped from a CDS sequence so that it matches
 979    * the peptide translation length.
 980    *
 981    * @param positions
 982    * @param ranges
 983    *          a list of (single) [start, end] ranges
 984    * @return
 985    */
 986   public static void removeEndPositions(int positions, List<int[]> ranges)
 987   {
 988     int toRemove = positions;
 989     Iterator<int[]> it = new ReverseListIterator<>(ranges);
 990     while (toRemove > 0)
 991     {
 992       int[] endRange = it.next();
 993       if (endRange.length != 2)
 994       {
 995         /*
 996          * not coded for [start1, end1, start2, end2, ...]
 997          */
 998         Cache.log.error(
 999                 "MappingUtils.removeEndPositions doesn't handle multiple  ranges");
1000         return;
1001       }
1002
1003       int length = endRange[1] - endRange[0] + 1;
1004       if (length <= 0)
1005       {
1006         /*
1007          * not coded for a reverse strand range (end < start)
1008          */
1009         Cache.log.error(
1010                 "MappingUtils.removeEndPositions doesn't handle reverse strand");
1011         return;
1012       }
1013       if (length > toRemove)
1014       {
1015         endRange[1] -= toRemove;
1016         toRemove = 0;
1017       }
1018       else
1019       {
1020         toRemove -= length;
1021         it.remove();
1022       }
1023     }
1024   }
1025
1026   /**
1027    * Converts a list of [start, end] ranges to a single array of [start, end,
1028    * start, end ...]
1029    *
1030    * @param ranges
1031    * @return
1032    */
1033   public static int[] listToArray(List<int[]> ranges)
1034   {
1035     int[] result = new int[ranges.size() * 2];
1036     int i = 0;
1037     for (int[] range : ranges)
1038     {
1039       result[i++] = range[0];
1040       result[i++] = range[1];
1041     }
1042     return result;
1043   }
1044
1045   /**
1046    * Returns the maximal start-end positions in the given (ordered) list of
1047    * ranges which is overlapped by the given begin-end range, or null if there
1048    * is no overlap.
1049    *
1050    * <pre>
1051    * Examples:
1052    *   if ranges is {[4, 8], [10, 12], [16, 19]}
1053    * then
1054    *   findOverlap(ranges, 1, 20) == [4, 19]
1055    *   findOverlap(ranges, 6, 11) == [6, 11]
1056    *   findOverlap(ranges, 9, 15) == [10, 12]
1057    *   findOverlap(ranges, 13, 15) == null
1058    * </pre>
1059    *
1060    * @param ranges
1061    * @param begin
1062    * @param end
1063    * @return
1064    */
1065   protected static int[] findOverlap(List<int[]> ranges, final int begin,
1066           final int end)
1067   {
1068     boolean foundStart = false;
1069     int from = 0;
1070     int to = 0;
1071
1072     /*
1073      * traverse the ranges to find the first position (if any) >= begin,
1074      * and the last position (if any) <= end
1075      */
1076     for (int[] range : ranges)
1077     {
1078       if (!foundStart)
1079       {
1080         if (range[0] >= begin)
1081         {
1082           /*
1083            * first range that starts with, or follows, begin
1084            */
1085           foundStart = true;
1086           from = Math.max(range[0], begin);
1087         }
1088         else if (range[1] >= begin)
1089         {
1090           /*
1091            * first range that contains begin
1092            */
1093           foundStart = true;
1094           from = begin;
1095         }
1096       }
1097
1098       if (range[0] <= end)
1099       {
1100         to = Math.min(end, range[1]);
1101       }
1102     }
1103
1104     return foundStart && to >= from ? new int[] { from, to } : null;
1105   }
1106 }