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 jalview.analysis.AlignmentSorter;
  24 import jalview.api.AlignViewportI;
  25 import jalview.commands.CommandI;
  26 import jalview.commands.EditCommand;
  27 import jalview.commands.EditCommand.Action;
  28 import jalview.commands.EditCommand.Edit;
  29 import jalview.commands.OrderCommand;
  30 import jalview.datamodel.AlignedCodonFrame;
  31 import jalview.datamodel.AlignmentI;
  32 import jalview.datamodel.AlignmentOrder;
  33 import jalview.datamodel.ColumnSelection;
  34 import jalview.datamodel.HiddenColumns;
  35 import jalview.datamodel.SearchResultMatchI;
  36 import jalview.datamodel.SearchResults;
  37 import jalview.datamodel.SearchResultsI;
  38 import jalview.datamodel.Sequence;
  39 import jalview.datamodel.SequenceGroup;
  40 import jalview.datamodel.SequenceI;
  41
  42 import java.util.ArrayList;
  43 import java.util.Arrays;
  44 import java.util.HashMap;
  45 import java.util.Iterator;
  46 import java.util.List;
  47 import java.util.Map;
  48
  49 /**
  50  * Helper methods for manipulations involving sequence mappings.
  51  *
  52  * @author gmcarstairs
  53  *
  54  */
  55 public final class MappingUtils
  56 {
  57
  58   /**
  59    * Helper method to map a CUT or PASTE command.
  60    *
  61    * @param edit
  62    *          the original command
  63    * @param undo
  64    *          if true, the command is to be undone
  65    * @param targetSeqs
  66    *          the mapped sequences to apply the mapped command to
  67    * @param result
  68    *          the mapped EditCommand to add to
  69    * @param mappings
  70    */
  71   protected static void mapCutOrPaste(Edit edit, boolean undo,
  72           List<SequenceI> targetSeqs, EditCommand result,
  73           List<AlignedCodonFrame> mappings)
  74   {
  75     Action action = edit.getAction();
  76     if (undo)
  77     {
  78       action = action.getUndoAction();
  79     }
  80     // TODO write this
  81     System.err.println("MappingUtils.mapCutOrPaste not yet implemented");
  82   }
  83
  84   /**
  85    * Returns a new EditCommand representing the given command as mapped to the
  86    * given sequences. If there is no mapping, returns null.
  87    *
  88    * @param command
  89    * @param undo
  90    * @param mapTo
  91    * @param gapChar
  92    * @param mappings
  93    * @return
  94    */
  95   public static EditCommand mapEditCommand(EditCommand command,
  96           boolean undo, final AlignmentI mapTo, char gapChar,
  97           List<AlignedCodonFrame> mappings)
  98   {
  99     /*
 100      * For now, only support mapping from protein edits to cDna
 101      */
 102     if (!mapTo.isNucleotide())
 103     {
 104       return null;
 105     }
 106
 107     /*
 108      * Cache a copy of the target sequences so we can mimic successive edits on
 109      * them. This lets us compute mappings for all edits in the set.
 110      */
 111     Map<SequenceI, SequenceI> targetCopies = new HashMap<>();
 112     for (SequenceI seq : mapTo.getSequences())
 113     {
 114       SequenceI ds = seq.getDatasetSequence();
 115       if (ds != null)
 116       {
 117         final SequenceI copy = new Sequence(seq);
 118         copy.setDatasetSequence(ds);
 119         targetCopies.put(ds, copy);
 120       }
 121     }
 122
 123     /*
 124      * Compute 'source' sequences as they were before applying edits:
 125      */
 126     Map<SequenceI, SequenceI> originalSequences = command.priorState(undo);
 127
 128     EditCommand result = new EditCommand();
 129     Iterator<Edit> edits = command.getEditIterator(!undo);
 130     while (edits.hasNext())
 131     {
 132       Edit edit = edits.next();
 133       if (edit.getAction() == Action.CUT
 134               || edit.getAction() == Action.PASTE)
 135       {
 136         mapCutOrPaste(edit, undo, mapTo.getSequences(), result, mappings);
 137       }
 138       else if (edit.getAction() == Action.INSERT_GAP
 139               || edit.getAction() == Action.DELETE_GAP)
 140       {
 141         mapInsertOrDelete(edit, undo, originalSequences,
 142                 mapTo.getSequences(), targetCopies, gapChar, result,
 143                 mappings);
 144       }
 145     }
 146     return result.getSize() > 0 ? result : null;
 147   }
 148
 149   /**
 150    * Helper method to map an edit command to insert or delete gaps.
 151    *
 152    * @param edit
 153    *          the original command
 154    * @param undo
 155    *          if true, the action is to undo the command
 156    * @param originalSequences
 157    *          the sequences the command acted on
 158    * @param targetSeqs
 159    * @param targetCopies
 160    * @param gapChar
 161    * @param result
 162    *          the new EditCommand to add mapped commands to
 163    * @param mappings
 164    */
 165   protected static void mapInsertOrDelete(Edit edit, boolean undo,
 166           Map<SequenceI, SequenceI> originalSequences,
 167           final List<SequenceI> targetSeqs,
 168           Map<SequenceI, SequenceI> targetCopies, char gapChar,
 169           EditCommand result, List<AlignedCodonFrame> mappings)
 170   {
 171     Action action = edit.getAction();
 172
 173     /*
 174      * Invert sense of action if an Undo.
 175      */
 176     if (undo)
 177     {
 178       action = action.getUndoAction();
 179     }
 180     final int count = edit.getNumber();
 181     final int editPos = edit.getPosition();
 182     for (SequenceI seq : edit.getSequences())
 183     {
 184       /*
 185        * Get residue position at (or to right of) edit location. Note we use our
 186        * 'copy' of the sequence before editing for this.
 187        */
 188       SequenceI ds = seq.getDatasetSequence();
 189       if (ds == null)
 190       {
 191         continue;
 192       }
 193       final SequenceI actedOn = originalSequences.get(ds);
 194       final int seqpos = actedOn.findPosition(editPos);
 195
 196       /*
 197        * Determine all mappings from this position to mapped sequences.
 198        */
 199       SearchResultsI sr = buildSearchResults(seq, seqpos, mappings);
 200
 201       if (!sr.isEmpty())
 202       {
 203         for (SequenceI targetSeq : targetSeqs)
 204         {
 205           ds = targetSeq.getDatasetSequence();
 206           if (ds == null)
 207           {
 208             continue;
 209           }
 210           SequenceI copyTarget = targetCopies.get(ds);
 211           final int[] match = sr.getResults(copyTarget, 0,
 212                   copyTarget.getLength());
 213           if (match != null)
 214           {
 215             final int ratio = 3; // TODO: compute this - how?
 216             final int mappedCount = count * ratio;
 217
 218             /*
 219              * Shift Delete start position left, as it acts on positions to its
 220              * right.
 221              */
 222             int mappedEditPos = action == Action.DELETE_GAP ? match[0]
 223                     - mappedCount : match[0];
 224             Edit e = result.new Edit(action, new SequenceI[] { targetSeq },
 225                     mappedEditPos, mappedCount, gapChar);
 226             result.addEdit(e);
 227
 228             /*
 229              * and 'apply' the edit to our copy of its target sequence
 230              */
 231             if (action == Action.INSERT_GAP)
 232             {
 233               copyTarget.setSequence(new String(StringUtils.insertCharAt(
 234                       copyTarget.getSequence(), mappedEditPos, mappedCount,
 235                       gapChar)));
 236             }
 237             else if (action == Action.DELETE_GAP)
 238             {
 239               copyTarget.setSequence(new String(StringUtils.deleteChars(
 240                       copyTarget.getSequence(), mappedEditPos,
 241                       mappedEditPos + mappedCount)));
 242             }
 243           }
 244         }
 245       }
 246       /*
 247        * and 'apply' the edit to our copy of its source sequence
 248        */
 249       if (action == Action.INSERT_GAP)
 250       {
 251         actedOn.setSequence(new String(StringUtils.insertCharAt(
 252                 actedOn.getSequence(), editPos, count, gapChar)));
 253       }
 254       else if (action == Action.DELETE_GAP)
 255       {
 256         actedOn.setSequence(new String(StringUtils.deleteChars(
 257                 actedOn.getSequence(), editPos, editPos + count)));
 258       }
 259     }
 260   }
 261
 262   /**
 263    * Returns a SearchResults object describing the mapped region corresponding
 264    * to the specified sequence position.
 265    *
 266    * @param seq
 267    * @param index
 268    * @param seqmappings
 269    * @return
 270    */
 271   public static SearchResultsI buildSearchResults(SequenceI seq, int index,
 272           List<AlignedCodonFrame> seqmappings)
 273   {
 274     SearchResultsI results = new SearchResults();
 275     addSearchResults(results, seq, index, seqmappings);
 276     return results;
 277   }
 278
 279   /**
 280    * Adds entries to a SearchResults object describing the mapped region
 281    * corresponding to the specified sequence position.
 282    *
 283    * @param results
 284    * @param seq
 285    * @param index
 286    * @param seqmappings
 287    */
 288   public static void addSearchResults(SearchResultsI results, SequenceI seq,
 289           int index, List<AlignedCodonFrame> seqmappings)
 290   {
 291     if (index >= seq.getStart() && index <= seq.getEnd())
 292     {
 293       for (AlignedCodonFrame acf : seqmappings)
 294       {
 295         acf.markMappedRegion(seq, index, results);
 296       }
 297     }
 298   }
 299
 300   /**
 301    * Returns a (possibly empty) SequenceGroup containing any sequences in the
 302    * mapped viewport corresponding to the given group in the source viewport.
 303    *
 304    * @param sg
 305    * @param mapFrom
 306    * @param mapTo
 307    * @return
 308    */
 309   public static SequenceGroup mapSequenceGroup(final SequenceGroup sg,
 310           final AlignViewportI mapFrom, final AlignViewportI mapTo)
 311   {
 312     /*
 313      * Note the SequenceGroup holds aligned sequences, the mappings hold dataset
 314      * sequences.
 315      */
 316     boolean targetIsNucleotide = mapTo.isNucleotide();
 317     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 318     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
 319             .getCodonFrames();
 320     /*
 321      * Copy group name, colours etc, but not sequences or sequence colour scheme
 322      */
 323     SequenceGroup mappedGroup = new SequenceGroup(sg);
 324     mappedGroup.setColourScheme(mapTo.getGlobalColourScheme());
 325     mappedGroup.clear();
 326
 327     int minStartCol = -1;
 328     int maxEndCol = -1;
 329     final int selectionStartRes = sg.getStartRes();
 330     final int selectionEndRes = sg.getEndRes();
 331     for (SequenceI selected : sg.getSequences())
 332     {
 333       /*
 334        * Find the widest range of non-gapped positions in the selection range
 335        */
 336       int firstUngappedPos = selectionStartRes;
 337       while (firstUngappedPos <= selectionEndRes
 338               && Comparison.isGap(selected.getCharAt(firstUngappedPos)))
 339       {
 340         firstUngappedPos++;
 341       }
 342
 343       /*
 344        * If this sequence is only gaps in the selected range, skip it
 345        */
 346       if (firstUngappedPos > selectionEndRes)
 347       {
 348         continue;
 349       }
 350
 351       int lastUngappedPos = selectionEndRes;
 352       while (lastUngappedPos >= selectionStartRes
 353               && Comparison.isGap(selected.getCharAt(lastUngappedPos)))
 354       {
 355         lastUngappedPos--;
 356       }
 357
 358       /*
 359        * Find the selected start/end residue positions in sequence
 360        */
 361       int startResiduePos = selected.findPosition(firstUngappedPos);
 362       int endResiduePos = selected.findPosition(lastUngappedPos);
 363
 364       for (AlignedCodonFrame acf : codonFrames)
 365       {
 366         SequenceI mappedSequence = targetIsNucleotide ? acf
 367                 .getDnaForAaSeq(selected) : acf.getAaForDnaSeq(selected);
 368         if (mappedSequence != null)
 369         {
 370           for (SequenceI seq : mapTo.getAlignment().getSequences())
 371           {
 372             int mappedStartResidue = 0;
 373             int mappedEndResidue = 0;
 374             if (seq.getDatasetSequence() == mappedSequence)
 375             {
 376               /*
 377                * Found a sequence mapping. Locate the start/end mapped residues.
 378                */
 379               List<AlignedCodonFrame> mapping = Arrays
 380                       .asList(new AlignedCodonFrame[] { acf });
 381               SearchResultsI sr = buildSearchResults(selected,
 382                       startResiduePos, mapping);
 383               for (SearchResultMatchI m : sr.getResults())
 384               {
 385                 mappedStartResidue = m.getStart();
 386                 mappedEndResidue = m.getEnd();
 387               }
 388               sr = buildSearchResults(selected, endResiduePos, mapping);
 389               for (SearchResultMatchI m : sr.getResults())
 390               {
 391                 mappedStartResidue = Math.min(mappedStartResidue,
 392                         m.getStart());
 393                 mappedEndResidue = Math.max(mappedEndResidue, m.getEnd());
 394               }
 395
 396               /*
 397                * Find the mapped aligned columns, save the range. Note findIndex
 398                * returns a base 1 position, SequenceGroup uses base 0
 399                */
 400               int mappedStartCol = seq.findIndex(mappedStartResidue) - 1;
 401               minStartCol = minStartCol == -1 ? mappedStartCol : Math.min(
 402                       minStartCol, mappedStartCol);
 403               int mappedEndCol = seq.findIndex(mappedEndResidue) - 1;
 404               maxEndCol = maxEndCol == -1 ? mappedEndCol : Math.max(
 405                       maxEndCol, mappedEndCol);
 406               mappedGroup.addSequence(seq, false);
 407               break;
 408             }
 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,
 433           boolean undo, 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     // ColumnSelection mappedColumns = new ColumnSelection();
 522
 523     if (colsel == null)
 524     {
 525       return; // mappedColumns;
 526     }
 527
 528     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
 529
 530     /*
 531      * For each mapped column, find the range of columns that residues in that
 532      * column map to.
 533      */
 534     List<SequenceI> fromSequences = mapFrom.getAlignment().getSequences();
 535     List<SequenceI> toSequences = mapTo.getAlignment().getSequences();
 536
 537     for (Integer sel : colsel.getSelected())
 538     {
 539       mapColumn(sel.intValue(), codonFrames, newColSel, fromSequences,
 540               toSequences, fromGapChar);
 541     }
 542
 543     for (int[] hidden : hiddencols)
 544     {
 545       mapHiddenColumns(hidden, codonFrames, newHidden, fromSequences,
 546               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,
 597           List<AlignedCodonFrame> mappings, ColumnSelection mappedColumns,
 598           List<SequenceI> fromSequences, List<SequenceI> toSequences,
 599           char fromGapChar)
 600   {
 601     int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
 602             toSequences, fromGapChar);
 603
 604     /*
 605      * Add the range of mapped columns to the mapped selection (converting
 606      * base 1 to base 0). Note that this may include intron-only regions which
 607      * lie between the start and end ranges of the selection.
 608      */
 609     if (mappedTo != null)
 610     {
 611       for (int i = mappedTo[0]; i <= mappedTo[1]; i++)
 612       {
 613         mappedColumns.addElement(i - 1);
 614       }
 615     }
 616   }
 617
 618   /**
 619    * Helper method to find the range of columns mapped to from one column.
 620    * Returns the maximal range of columns mapped to from all sequences in the
 621    * source column, or null if no mappings were found.
 622    *
 623    * @param col
 624    * @param mappings
 625    * @param fromSequences
 626    * @param toSequences
 627    * @param fromGapChar
 628    * @return
 629    */
 630   protected static int[] findMappedColumns(int col,
 631           List<AlignedCodonFrame> mappings, List<SequenceI> fromSequences,
 632           List<SequenceI> toSequences, char fromGapChar)
 633   {
 634     int[] mappedTo = new int[] { Integer.MAX_VALUE, Integer.MIN_VALUE };
 635     boolean found = false;
 636
 637     /*
 638      * For each sequence in the 'from' alignment
 639      */
 640     for (SequenceI fromSeq : fromSequences)
 641     {
 642       /*
 643        * Ignore gaps (unmapped anyway)
 644        */
 645       if (fromSeq.getCharAt(col) == fromGapChar)
 646       {
 647         continue;
 648       }
 649
 650       /*
 651        * Get the residue position and find the mapped position.
 652        */
 653       int residuePos = fromSeq.findPosition(col);
 654       SearchResultsI sr = buildSearchResults(fromSeq, residuePos, mappings);
 655       for (SearchResultMatchI m : sr.getResults())
 656       {
 657         int mappedStartResidue = m.getStart();
 658         int mappedEndResidue = m.getEnd();
 659         SequenceI mappedSeq = m.getSequence();
 660
 661         /*
 662          * Locate the aligned sequence whose dataset is mappedSeq. TODO a
 663          * datamodel that can do this efficiently.
 664          */
 665         for (SequenceI toSeq : toSequences)
 666         {
 667           if (toSeq.getDatasetSequence() == mappedSeq)
 668           {
 669             int mappedStartCol = toSeq.findIndex(mappedStartResidue);
 670             int mappedEndCol = toSeq.findIndex(mappedEndResidue);
 671             mappedTo[0] = Math.min(mappedTo[0], mappedStartCol);
 672             mappedTo[1] = Math.max(mappedTo[1], mappedEndCol);
 673             found = true;
 674             break;
 675             // note: remove break if we ever want to map one to many sequences
 676           }
 677         }
 678       }
 679     }
 680     return found ? mappedTo : null;
 681   }
 682
 683   /**
 684    * Returns the mapped codon or codons for a given aligned sequence column
 685    * position (base 0).
 686    *
 687    * @param seq
 688    *          an aligned peptide sequence
 689    * @param col
 690    *          an aligned column position (base 0)
 691    * @param mappings
 692    *          a set of codon mappings
 693    * @return the bases of the mapped codon(s) in the cDNA dataset sequence(s),
 694    *         or an empty list if none found
 695    */
 696   public static List<char[]> findCodonsFor(SequenceI seq, int col,
 697           List<AlignedCodonFrame> mappings)
 698   {
 699     List<char[]> result = new ArrayList<>();
 700     int dsPos = seq.findPosition(col);
 701     for (AlignedCodonFrame mapping : mappings)
 702     {
 703       if (mapping.involvesSequence(seq))
 704       {
 705         List<char[]> codons = mapping.getMappedCodons(
 706                 seq.getDatasetSequence(), dsPos);
 707         if (codons != null)
 708         {
 709           result.addAll(codons);
 710         }
 711       }
 712     }
 713     return result;
 714   }
 715
 716   /**
 717    * Converts a series of [start, end] range pairs into an array of individual
 718    * positions. This also caters for 'reverse strand' (start > end) cases.
 719    *
 720    * @param ranges
 721    * @return
 722    */
 723   public static int[] flattenRanges(int[] ranges)
 724   {
 725     /*
 726      * Count how many positions altogether
 727      */
 728     int count = 0;
 729     for (int i = 0; i < ranges.length - 1; i += 2)
 730     {
 731       count += Math.abs(ranges[i + 1] - ranges[i]) + 1;
 732     }
 733
 734     int[] result = new int[count];
 735     int k = 0;
 736     for (int i = 0; i < ranges.length - 1; i += 2)
 737     {
 738       int from = ranges[i];
 739       final int to = ranges[i + 1];
 740       int step = from <= to ? 1 : -1;
 741       do
 742       {
 743         result[k++] = from;
 744         from += step;
 745       } while (from != to + step);
 746     }
 747     return result;
 748   }
 749
 750   /**
 751    * Returns a list of any mappings that are from or to the given (aligned or
 752    * dataset) sequence.
 753    *
 754    * @param sequence
 755    * @param mappings
 756    * @return
 757    */
 758   public static List<AlignedCodonFrame> findMappingsForSequence(
 759           SequenceI sequence, List<AlignedCodonFrame> mappings)
 760   {
 761     return findMappingsForSequenceAndOthers(sequence, mappings, null);
 762   }
 763
 764   /**
 765    * Returns a list of any mappings that are from or to the given (aligned or
 766    * dataset) sequence, optionally limited to mappings involving one of a given
 767    * list of sequences.
 768    *
 769    * @param sequence
 770    * @param mappings
 771    * @param filterList
 772    * @return
 773    */
 774   public static List<AlignedCodonFrame> findMappingsForSequenceAndOthers(
 775           SequenceI sequence, List<AlignedCodonFrame> mappings,
 776           List<SequenceI> filterList)
 777   {
 778     List<AlignedCodonFrame> result = new ArrayList<>();
 779     if (sequence == null || mappings == null)
 780     {
 781       return result;
 782     }
 783     for (AlignedCodonFrame mapping : mappings)
 784     {
 785       if (mapping.involvesSequence(sequence))
 786       {
 787         if (filterList != null)
 788         {
 789           for (SequenceI otherseq : filterList)
 790           {
 791             SequenceI otherDataset = otherseq.getDatasetSequence();
 792             if (otherseq == sequence
 793                     || otherseq == sequence.getDatasetSequence()
 794                     || (otherDataset != null && (otherDataset == sequence || 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         System.err.println("Error unbalance start/end ranges: "
 836                 + 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 }